「国立科学博物館で学ぶ物理学」メニュー

ジャイロホイール

回折

HR図

国立科学博物館で学ぶ物理学　＜数値計算法＞

[概説]

[適用例]

[Ⅰ]微分方程式の数値解法の概要

物理方程式では、一般的に２階微分方程式を扱うので、簡単のため以下では２階微分方程式の数値解法について説明する
常微分方程式の数値解析に対し最も多く用いられる一般的方法は離散変数法と呼ばれるもので、①オイラー法　②ルンゲクッタ法　などがあります。その違いは関数をテーラー展開したときの微少量の次数をどこまで取るかに拠ります。②は４次の微少量まで取るので①より精度が高く、多用されています。
一方、偏微分方程式式に対し最も多く用いられる一般的方法は、③差分法と呼ばれるものすが、常微分方程式に使用することも出来ます。

３つの方法の比較は下図のようになります。

オイラー法（赤）は誤差が短時間に現れるが、ルンゲクッタ法（赤紫）・差分化法（青）はいずれも精度よく解が得られている。

[Ⅱ] オイラー法

(a)従属変数が１個の１階常微分方程式
　　　　　 $dx/dt=f(t,x)$
の場合、現在値 x_iに対し、⊿t 後の値 x_i+1は、単純に積分した値を加えるだけです。
　　　　　 $x_{i+1}=x_i +f(t_i,x_i)\Delta t$

(b)従属変数が１個の２階微分方程式
　　　　　 $dx^2/dt^2=f(t,x,dx/dt)$
の場合は、これを以下のような１階の連立常微分方程式
　　　　　 $dx/dt=v=f_1(t,x,v)$
　　　　　 $dv/dt=f_2(t,x,v)$
に置き換えて、夫々の式に対し１階常微分方程式と同じ操作を行います。
　　　　　 $v_{i+1}=v_i +f_1(t_i,x_i,v_i)\Delta t$
　　　　　 $x_{i+1}=x_i +f_2(t_i,x_i,v_i)\Delta t$
これを単純に繰り返すだけの方法です。極短時間の計算では、この方法で十分でしょう。

[Ⅲ] ルンゲクッタ法

(a)従属変数が１個の１階常微分方程式
　　　　　 $dx/dt=f(t,x)$
の場合、
　　　　　 $k_1=\Delta t\cdot f(t_i ,\hspace{5mm} x_i)$
　　　　　 $k_2=\Delta t\cdot f(t_i+\Delta t/2\hspace{5mm} ,\hspace{5mm} x_i+k_1/2)$
　　　　　 $k_3=\Delta t\cdot f(t_i+\Delta t/2\hspace{5mm} ,\hspace{5mm} x_i+k_2/2)$
　　　　　 $k_4=\Delta t\cdot f(t_i+\Delta t\hspace{5mm} ,\hspace{5mm} x_i+k_3)$
として、⊿t 後の値x_i+1を次式で求めます。
　　　　　 $x_{i+1}=x_i+(k_1+2k_2+2k_3+k_4)/6$

(b)従属変数が１個の２階微分方程式
　　　　　 $d^2x/dt^2=f(t,x,dx/dt)$
の場合は、[Ⅰ](b)と同様に１階の連立常微分方程式
　　　　　 $dx/dt=v=f_1(t,x,v)$
　　　　　 $dv/dt=f_2(t,x,v)$
に分け、上の１階微分方程式の解法を適用します。すなわち、
　　　　　 $k_{x1}=\Delta t\cdot f_1( v_i)\hspace{90mm}k_{v1}=\Delta t\cdot f_2(t_i ,\hspace{5mm} x_i,\hspace{5mm} v_i)$
　　　　　 $k_{x2}=\Delta t\cdot f_1(v_i+k_{v1}/2)\hspace{20mm}k_{v2}=\Delta t\cdot f_2(t_i+\Delta t/2\hspace{5mm} ,\hspace{5mm} x_i+k_{x1}/2,\hspace{5mm} v_i+k_{v1}/2)$
　　　　　 $k_{x3}=\Delta t\cdot f_1(v_i+k_{v2}/2)\hspace{20mm}k_{v3}=\Delta t\cdot f_2(t_i+\Delta t/2\hspace{5mm} ,\hspace{5mm} x_i+k_{x2}/2,\hspace{5mm} v_i+k_{v2}/2)$
　　　　　 $k_{x4}=\Delta t\cdot f_1(v_i+k_{v3})\hspace{40mm}k_{v4}=\Delta t\cdot f_2(t_i+\Delta t\hspace{5mm} ,\hspace{5mm} x_i+k_{x3},\hspace{5mm} v_i+k_{v3})$
として、⊿t 後の値v_i+1、x_i+1を次式で求めます。
　　　　　 $v_{i+1}=v_i+(k_{v1}+2k_{v2}+2k_{v3}+k_{v4})/6\hspace{100mm}$
　　　　　 $x_{i+1}=x_i+(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})/6$

(c)従属変数が2個の２階連立微分方程式
　　　　　 $d^2x/dt^2=f_x(t,x,y,dx/dt,dy/dt)$
　　　　　 $d^2y/dt^2=f_y(t,x,y,dx/dt,dy/dt)$
の場合、これを更に分けて１階の連立常微分方程式

　　　　　 $dx/dt=v_x=f_{x1}(t,x,y,v_x,v_y)$ 　　　　　　 $dv_x/dt=f_{x2}(t,x,y,v_x,v_y)$
　　　　　 $dy/dt=v_y=f_{y1}(t,x,y,v_x,v_y)$ 　　　　　　 $dv_y/dt=f_{y2}(t,x,y,v_x,v_y)$
に分け、上の１階微分方程式の解法を適用します。すなわち、
　　　　　 $k_{x1}=\Delta t\cdot f_{x1}(t_i,x_i,y_i, v_{xi},v_{yi})\hspace{90mm}k_{vx1}=\Delta t\cdot f_{x2}(t_i ,\hspace{5mm} x_i,\hspace{5mm} y_i ,\hspace{5mm} v_{xi},\hspace{5mm} v_{yi})$
　　　　　 $k_{y1}=\Delta t\cdot f_{y1}(t_i,x_i,y_i, v_{xi},v_{yi})\hspace{90mm}k_{vy1}=\Delta t\cdot f_{y2}(t_i ,\hspace{5mm} x_i,\hspace{5mm} y_i ,\hspace{5mm} v_{xi},\hspace{5mm} v_{yi})$

　　　　　 $k_{x2}=\Delta t\cdot f_{x1}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x1}/2\hspace{5mm},\hspace{5mm}y_i+k_{y1}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx1}/2\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy1}/2)$
　　　　　 $k_{y2}=\Delta t\cdot f_{y1}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x1}/2\hspace{5mm},\hspace{5mm}y_i+k_{y1}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx1}/2\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy1}/2)$
　　　　　 $k_{vx2}=\Delta t\cdot f_{x2}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x1}/2\hspace{5mm},\hspace{5mm}y_i+k_{y1}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx1}/2,v_{yi}+k_{vy1}/2)$
　　　　　 $k_{vy2}=\Delta t\cdot f_{y2}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x1}/2\hspace{5mm},\hspace{5mm}y_i+k_{y1}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx1}/2,v_{yi}+k_{vy1}/2)$

　　　　　 $k_{x3}=\Delta t\cdot f_{x1}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x2}/2\hspace{5mm},\hspace{5mm}y_i+k_{y2}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx2}/2\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy2}/2)$
　　　　　 $k_{y3}=\Delta t\cdot f_{y1}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x2}/2\hspace{5mm},\hspace{5mm}y_i+k_{y2}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx2}/2\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy2}/2)$
　　　　　 $k_{vx3}=\Delta t\cdot f_{x2}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x2}/2\hspace{5mm},\hspace{5mm}y_i+k_{y2}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx2}/2,v_{yi}+k_{vy2}/2)$
　　　　　 $k_{vy3}=\Delta t\cdot f_{y2}(t_i+\Delta t/2\hspace{5mm},\hspace{5mm}x_i+k_{x2}/2\hspace{5mm},\hspace{5mm}y_i+k_{y2}/2\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx2}/2,v_{yi}+k_{vy2}/2)$

　　　　　 $k_{x4}=\Delta t\cdot f_{x1}(t_i+\Delta t\hspace{5mm},\hspace{5mm}x_i+k_{x3}\hspace{5mm},\hspace{5mm}y_i+k_{y3}\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx3}\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy3})$
　　　　　 $k_{y4}=\Delta t\cdot f_{y1}(t_i+\Delta t\hspace{5mm},\hspace{5mm}x_i+k_{x3}\hspace{5mm},\hspace{5mm}y_i+k_{y3}\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx3}\hspace{5mm},\hspace{5mm}v_{yi}+k_{vy3})$
　　　　　 $k_{vx4}=\Delta t\cdot f_{x2}(t_i+\Delta t\hspace{5mm},\hspace{5mm}x_i+k_{x3}\hspace{5mm},\hspace{5mm}y_i+k_{y3}\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx3},v_{yi}+k_{vy3})$
　　　　　 $k_{vy4}=\Delta t\cdot f_{y2}(t_i+\Delta t\hspace{5mm},\hspace{5mm}x_i+k_{x3}\hspace{5mm},\hspace{5mm}y_i+k_{y3}\hspace{5mm}, \hspace{5mm}v_{xi}+k_{vx3},v_{yi}+k_{vy3})$

として、⊿t 後の値x_i+1を次式で求めます。
　　　　　 $x_{i+1}=x_i+(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})/6$
　　　　　 $y_{i+1}=y_i+(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})/6$
　　　　　 $vx_{i+1}=vx_i+(k_{vx1}+2k_{vx2}+2k_{vx3}+k_{vx4})/6\hspace{100mm}$
　　　　　 $vy_{i+1}=vy_i+(k_{vy1}+2k_{vy2}+2k_{vy3}+k_{vy4})/6\hspace{100mm}$

[Ⅳ] 差分方程式による解法

(a)常微分方程式　(多変数の従属関数であれ、高階微分方程式であれ計算法は同じ）
　　　　　 $d^2x/dt^2=f_x(t,x,y,z,v_x,v_y,v_z)$
　　　　　 $d^2y/dt^2=f_y(t,x,y,z,v_x,v_y,v_z)$
　　　　　 $d^2z/dt^2=f_z(t,x,y,z,v_x,v_y,v_z)$
このような場合、x について言えば、現在値をx_i、時間⊿t 前の値をx_i-1、時間⊿t 後の値をx_i+1として、dx/dt、d²x/dt²を、次のように差分化します。
　　　　　 $\frac{dx}{dt}=\frac{(x_{i+1}-x_{i-1})}{2\Delta t}$
　　　　　 $\frac{d^2x}{dt^2}=\frac{\frac{(x_{i+1}-x_{i})}{\Delta t}-\frac{(x_{i}-x_{i-1})}{\Delta t}}{\Delta t}=\frac{(x_{i+1}-2x_{i}+x_{i-1})}{(\Delta t)^2}$
y,zについても同様で、これらを２階微分方程式に代入して得られる x_i+1, y_i+1, z_i+1 に関する連立方程式を解けば、以下の形で解が得られます。
　　　　　 $x_{i+1}=F_1(x_{i},y_{i},z_{i},x_{i-1},y_{i-1},z_{i-1})$
　　　　　 $y_{i+1}=F_2(x_{i},y_{i},z_{i},x_{i-1},y_{i-1},z_{i-1})$
　　　　　 $z_{i+1}=F_3(x_{i},y_{i},z_{i},x_{i-1},y_{i-1},z_{i-1})$
この計算を繰り返せば、x,y,z の時間変化が得られます。
なお、速度の差分化には、上に述べた中央差分ではなく、次のような差分を用いて簡単化することもあります。精度との見合いで選択する必要があります。
　　　　　 $\frac{dx}{dt}=\frac{(x_{i}-x_{i-1})}{\Delta t}$

(b)偏微分方程式　例えば
　　　　　 $\partial T/\partial t=a(\partial^2T/\partial x^2+\partial T/\partial y^2)$
の場合、(x,ｙ)面を升目に分割し、
　　　　　　　 $\frac{\partial T}{\partial t}=\frac{(T_{i+1}-T_{i-1})}{2\Delta t}$
　　　　　 $\frac{\partial T}{\partial x}=\frac{(T_{i+1}-T_{i-1})}{2\Delta x}$
　　　　　 $\frac{\partial^2 T}{\partial x^2}=\frac{(T_{i+1}-2T_{i}+T_{i-1})}{(\Delta x)^2}$

[Ⅴ]数値計算例

(a)単純振り子

この方程式の厳密解は、初期値θ=θ₀とすると、
　　　　　 $\theta=\theta_0cos(\omega^{'} t)$
　　　　　 $T=(2\pi/\omega) \sum_{n=1}^\infty[(2n-1)!!/(2n)!!]^2\sin^{2n}(\theta_0/2)$
　　　　　 $\omega^{'}=2\pi /T$ ３つの方法の比較は下図のようになります。

オイラー法（赤）は誤差が短時間に現れるが、ルンゲクッタ法（赤紫）・差分化法（青）はいずれも精度よく解が得られている。
敢えて言えば、振幅の精度は極めて高いが、周期については誤差は小さいとは言え、振幅の精度ほどではない。

(b)フーコーの振り子

フーコーの振り子の近似方程式は、
　　　　　 $\ddot{x}-2\Omega{'}\dot{y}+\omega^2x =0$
　　　　　 $\ddot{y}+2\Omega{'}\dot{x}+\omega^2y =0$
この方程式の厳密解は、
　　　　　 $\xi=\mathrm{e}^{-i\Omega{'}t}\{C_1\mathrm{e}^{i\omega{t}}+C_2\mathrm{e}^{-i\omega{t}}\}$