計算したい微分方程式
m\ddot x+c\dot x+kx=F\left(t \right)
中央差分の公式
\ddot x=\frac{x^{n+1}-2x^{n}+x^{n-1}}{dt^2} \\
\dot x=\frac{x^{n+1}-x^{n}}{dt}
代入したら
m\frac{x^{n+1}-2x^{n}+x^{n-1}}{dt^2}+c\frac{x^{n+1}-x^{n}}{dt}+kx^{n+1} =F^{n+1}
右辺には$x^{n+1}$が未知数
外力Fが$ A\sin\left( 2\pi \omega t \right)$みたいな$t$さえわかれば出てくるものとかなら$x^{n+1}=$の形にすれば求まる。
外力Fが、$F\left(t,x,\dot x,\ddot x\right)$で中身が特に変な関数が入ってくると単純には解けなくなる。
教科書には、「両辺に未知数があるので解けば答えが求まる」みたいにサラっと書いてることが多いが、「そこのやり方を教えてよ」って感じになってる。
newton raphson法
ニュートン ラフソン法の一般的な式は以下 詳しくはwiki見て
x_{i+1}=x_i-\frac{f\left(x_i\right)}{f'\left(x_i\right)}
求めたいものは以下 右辺の$F^{n+1}$を左辺に持ってきて$=0$にしただけ。
f\left(x^{n+1} \right) = m\frac{x^{n+1}-2x^{n}+x^{n-1}}{dt^2}+c\frac{x^{n+1}-x^{n}}{dt}+kx^{n+1} -F^{n+1} =0
$f\left(x^{n+1} \right)$が0になる$x^{n+1}$を求めたい。
Newton Raphson法の公式をみると、${f'\left(x^{n+1}\right)}$が必要になる。関数$f$の$x^{n+1}$の微分ってこと。
$x^{n},x^{n-1}$はただの定数だから消える。
{f'\left(x^{n+1}\right)} =\frac{m}{dt^2} + \frac{c}{dt}+k -\dot F^{n+1} \left(x^{n+1} \right)
外力 $F'^{n+1}\left(t+dt,x^{n+1},\dot x^{n+1},\ddot x^{n+1} \right)$とかどうしよう・・・・
今回は確認のためだから一定力として、微分したら0になるものを使う。
外力 $F'_{i+1}\left(t,x_i,\dot x_i,\ddot x_i \right)$がマジで非線形すぎて微分出来ない時に無理やり求める方法はそのうち。
これを繰り返していって、$x_i$と$x_{i+1}$が変化しなくなるまで繰り返す。
初期値計算
あるステップでの$x^{n+1}$を求めるために$x_{i=0}^{n+1}$を仮置しなければならない。
m\frac{x^{n+1}-2x^{n}+x^{n-1}}{dt^2}+c\frac{x^{n+1}-x^{n}}{dt}+kx^{n+1} =F^{n}
右辺の$F^{n+1}$を$F^{n}$にして、陽解法的に$x^{n+1}$を求めて初期値にする。
x^{n+1}_{i=0}=\frac{F^n-m\left(\frac{-2x^n+x^{n-1}}{dt^2} \right) -c \left(- \frac{x^n}{dt} \right) }{ \left(\frac{m}{dt^2}+c+k \right)}
初期値の次の計算
ニュートン・ラフソン法の次の計算は、公式を再度みてみると以下なので
x_{i+1}=x_i-\frac{f\left(x_i\right)}{f'\left(x_i\right)}
以下のような形になって、繰り返し計算をする。
x_{2}=x_1-\frac{f\left(x_1\right)}{f'\left(x_1\right)}
$x_{i}$と$x_{i+1}$の値が変化しなくなったら$x^{n+1}$して、また連立方程式を解く。
一般的には、0.000001くらいまでやるかな。