そもそもなんでこんなことをやるのか
ただの暇つぶしです。物理をただやっていてもなんかつまんないなぁ、と思ってしまって。ならとことん攻めた物理をしようと。いまパッと思い浮かんだのでこの記事を書いています。そしてまだ勉強中の身ですので、多くの間違いがあると思いますが、その時はご指摘をよろしくお願いします。そしてめんどくさくなると色々飛ばしますが、そこはご容赦ください。
運動方程式まで
っと、いきなり運動と言われても何するかわからん。てことで三次元の直交座標系を与えよう。
三次元のベクトル空間を $ V $ とする。 $ V $ には加法・減法・スカラー積・ベクトル積・外積（楔積）・テンソル積が既に定義されているとしよう（あった方が何かと便利なので）
基底ベクトルを $ \boldsymbol{e}_x,\boldsymbol{e}_y,\boldsymbol{e}_z $ として、それぞれは直交しているとする。添え字 $i=1,2,3$ を用いて基底ベクトルを表すと
$$ \hat{e}_i = (\boldsymbol{e}_x,\boldsymbol{e}_y,\boldsymbol{e}_z)\ \ ,\ \ \hat{e}_i \cdot \hat{e}_j=0\ (i\neq j) $$
こうなりますと。でベクトルをこう表記しよう。
ベクトル$A\in V$に対して、$ A_x=A^1, A_y=A^2, A_z=A^3 $ とすると、アインシュタインの縮約記法を用いれば$ \boldsymbol{A}=A^i\hat{e}_i $ となる（めんどくさいので、縮約の時は基本的に$i=1,2,3$としよう）
ここまではベクトルを何となく定義したに過ぎない（厳密とは程遠いが、お遊びにはちょうどいいだろう）
ほいで位置ベクトル$x^i$を用いて色々定義していきたいが、単位系を決めていないではないか！世界的に有名なのがSI単位系なので、私も例に習ってSI単位系を用いることにしよう。
んじゃあ位置ベクトル$x^i[m]$を用いて速度ベクトル$v^i[m/s]$を定義する。
時間を$t[s]$とすると、
$$v^i=\cfrac{dx^i}{dt}$$
さらに加速度ベクトル$a^i[m/s^2]$もやっちまおうぜ
$$a^i=\cfrac{d^2x^i}{dt^2}=\cfrac{dv^i}{dt}$$
まあ微分なわけだ。微分の定義は、言う必要ないよな？（めんどくさいだけ）
とは言っても、運動方程式を定義してないのでさっさと定義しちゃおうよ。
力ベクトル$F^i[N]$ 質量$m[kg]$ を用いてニュートンの運動方程式を表すと
$$ma^i=F^i$$
いえーい定義できた！ついでに運動量$p^i=mv^i$なるものを定義すると、ニュートンの運動方程式は
$$F^i=\cfrac{dp^i}{dt}$$
とも表せる。
摩擦力と空気抵抗と弾性力を瞬殺したい
垂直抗力を$N[N]$とすると、最大摩擦力$F_0[N]$は静止摩擦係数$\mu$を用いて$F_0=\mu N$
動摩擦力$F'[N]$は動摩擦係数$\mu'$を用いて$F'=\mu' N$と表せる。
物体の空気抵抗力$F[N]$は空気抵抗係数$R$を用いて$F^i=-Rv^i$と表せる。
ばねを自然長から$x[m]$引っ張ったときのばねの弾性力$F[N]$はばね係数$k[N/m]$を用いて$F=kx$と表せる。
よし下準備はこれでOK。こっからめっちゃめんどくさい微分方程式立てるぞ！
- 空気抵抗を考慮した摩擦あり斜面を滑る質点の運動
重力場は鉛直下向き一様で、重力加速度をgとしよう。そして斜面は$\theta[rad] (0<\theta<\frac{\pi}{2})$で水平面と角度をなしている。ついでに言うと延々と続いている。
座標系の取り方は、今回は斜面下向きをxの正、斜面上向きをyの正とする。要は二次元で考えるって話。めんどくさいのでベクトルは添え字を使わない。
質量は$m[kg]$、空気抵抗は$R$、静止摩擦係数は$\mu$、動摩擦係数は$\mu'$を用いる。早速立式したいところだが、摩擦力が働く斜面には摩擦角という概念があるので、まずはそこから議論しなければならない。今から考えるのは最大摩擦力と重力が釣り合っている状態だ。
静止している場合には当然y方向の力は釣り合ってなければならない。y方向は重力と垂直抗力が釣り合っているため、これらについて考えればよい。重力のy成分はどのような値になっているだろう。これは三角関数を用いると$mg\cos\theta$と表せる。よって垂直抗力Nと重力の関係は$N=mg\cos\theta$で表される。
ではx方向はどうだろうか？静止しているのでこちらも釣り合っている。物体には摩擦力、重力が働いている。空気抵抗は速度が0なので働かないので、今回は考える必要はない。
さて、重力$mg$のx成分は$mg\sin\theta$、最大摩擦力は$\mu N$で表されるため、これがイコールになっていれば良い。先ほどの垂直抗力と重力の関係を用いて立式すると
$$mg\sin\theta=\mu mg\cos\theta$$
$$\theta=\arctan\mu$$
となる。ここから微分方程式を立式しよう。
$0<\theta\leq\arctan\mu$ のとき
物体は静止している。
$\arctan\mu<\theta<\frac{\pi}{2}$ のとき
$ma=F$より
$$m\cfrac{dv_x}{dt}=mg\sin\theta-\mu'mg\cos\theta-kv_x$$
こちらの式は一階線形微分方程式となっておりまっせ。
どうやって解こうかなと、思いまして。ラプラス変換か同次式からのパターンか。乱数生成を用いてどっちか決めてもらいましょう。（マジで使いました）
結果はラプラス変換でしたね。
ということで両辺ラプラス変換しまーす、の前に式を少し変形して
$\cfrac{dv_x}{dt}+\cfrac{k}{m}v_x=g\sin\theta-\mu'g\cos\theta$
そして両辺をラプラス変換変換して
$sV(s)+\cfrac{k}{m}V(s)=\cfrac{\sin\theta-\mu'\cos\theta}{s}g$
$V(s)=\cfrac{g(\sin\theta-\mu'\cos\theta)}{s(s+\cfrac{k}{m})}$
ヘヴィサイトの展開定理を用いて部分分数分解して、逆ラプラス変換をすると
$v_x(t)=\cfrac{m}{k}g(\sin\theta-\mu'\cos\theta)(1-e^{-\frac{k}{m}t})$
両辺をtで積分して、初期値が0になるように積分定数を決定すると
$x(t)=\cfrac{m}{k}g(\sin\theta-\mu'\cos\theta)t+\cfrac{m^2}{k^2}g(\sin\theta-\mu'\cos\theta)e^{-\frac{k}{m}t}-\cfrac{m^2}{k^2}g(\sin\theta-\mu'\cos\theta)$
これが解だ。このまま終わりでもいいが、もう少し分析したい。
もし仮に$\mu'=0$で$k\to0$ならば、$x=\frac{1}{2}gt^2$に近づくはずだ。今回はこいつの近似まで行きたい。
$\mu'=0$の時の式は発散しないので成立する。
$x(t)=\cfrac{m}{k}g(\sin\theta)t+\cfrac{m^2}{k^2}g(\sin\theta)e^{-\frac{k}{m}t}-\cfrac{m^2}{k^2}g(\sin\theta)$
さらに同類項をまとめる
$x(t)=\cfrac{m}{k}g\sin\theta(t+\cfrac{m}{k}e^{-\frac{k}{m}t}-\cfrac{m}{k})$
ここから試しに指数関数項をテイラー展開してみよう
$x=\cfrac{m}{k}g\sin\theta(t+\cfrac{m}{k}-t-\cfrac{m}{k}+\cfrac{k}{m}\cfrac{t^2}{2!}\cdots)$
さらに展開して
$x=g(\sin\theta)(\cfrac{t^2}{2!}-\cfrac{k}{m}\cfrac{t^3}{3!}\cdots)$
ここまでやって、ようやく$k=0$で発散しない式ができた。$k=0$を代入して
$x=\cfrac{1}{2}gt^2\sin\theta$
この式は空気抵抗無し、摩擦無しの場合の$x$の式である。
ちゃんと近似できるところも気持ちいいところだ。