×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

物理の小さい余白 物理の大事な骨の部分だけを独自にまとめましたサイトです
ご意見ご感想はこちらまで(「」内に@gmail.comを付けてください) 「ysdswallow」
top page 微分・積分のまとめ ニュートン力学のまとめ
電磁気のまとめ 解析力学のまとめ 数値計算(微分方程式)のまとめ
おすすめ参考書等紹介

Chapter 3
高次の Runge-Kutta 法

2次より高次の RK 法は同じように

y    = y + hF (x  ,y ,h)
 n+1 m  n       n n
F = ∑   αiki
     i=1
           ∑m         m∑
ki = f(xn +    Pij,yn +    Pijkj)
           j=1        j=1
というふうにf(x,y)  を適当に組み合わせることで得られます。ここで以下
(                     )
   P11  P12  ⋅⋅⋅  P1m
||  P21  P22  ⋅⋅⋅  P2m ||
||    ..    ..  ...   ..  ||
|(    .    .        .  |)
   Pm1   ⋅⋅⋅  ⋅⋅⋅ Pmm
    α1   ⋅⋅⋅  ⋅⋅⋅  αm

といふうにm+1行m列の行列の形(Stetter の行列表現)で高次の RK 法を紹介していきます。

3.1 2次の RK 法

先ほどやったやつですがこの行列で表すと

(          )
   0    0
(  1    0  )
  1∕2  1∕2

また、さきほどβ = 1  の場合は

(  0   0 )
( 1∕2  0 )
   0   1

となります。

3.1.1 3次の RK 法

(   0    0    0 )
|  1∕2   0    0 |
|(  - 1   2    0 |)
   1∕6  4∕6  1∕6

が3次の RK 法です。もちろん係数は代表的なものであって、これが唯一ではありません。ただ、3次はあまり使われませ ん。それは手間や精度などを考慮すると4次のほうがお得だからです。

3.2 4次の RK 法

4次では例えば

(                    )
    0   0    0    0
||  1∕2  0    0    0  ||
||   0  1∕2   0    0  ||
(   0   0    1    0  )
   1∕6 2∕6  2∕6  1∕6

のようなものが有名です。4次の RK 法はかなりよく使われるものです。ためしにこれであの

 d2x    x
 dt2 = - r3
 d2y    y
 --2 = --3
 dt     r       dx(0)     dy(0)
x (0) = 1,y(0) = 0,---- = 0,-----= 1
    ∘ -------     dt       dt
r =   x2 + y2
を解いてみましょう。具体的な計算手順は
k  = f(x ,y )
 1      n  n
k2 = f(xn + h∕2,yn + k1h∕2)
k3 = f(xn + h∕2,yn + k2h∕2)
k4 = f(xn + h,yn + k3h)
y    = y + (k + 2k + 2k + k )h∕6
 n+1    n    1    2    3   4
となります。このソースは
/*四次RK法*/  
 
#include <stdio.h>  
#include <math.h>  
 
int main(void){  
 double x =1;  
 double kx[5] ={0};  
 double u =0;  
 double ku[5] ={0};  
 double y =0;  
 double ky[5] ={0};  
 double v =1;  
 double kv[5] ={0};  
 double h;  
 int k,t;  
 
 for(k=0;k < 100; k++){  
  for(t=0;t <=3;t++){  
   if(t<=2) h=0.1;  
   else h=0.2;  
   ku[t+1] =- (x+h *kx[t]) / sqrt((y+h *ky[t])*(y+h*ky[t])+(x+h *kx[t])*(x+h *kx[t]));  
   kx[t+1] = u+h *ku[t];  
   kv[t+1] = -(y+h *ky[t]) / sqrt((y+h *ky[t])*(y+h*ky[t])+(x+h *kx[t])*(x+h *kx[t]));  
   ky[t+1] =  v+h *kv[t];  
  }  
 
  u =u+(ku[1]+2*ku[2]+2*ku[3]+ku[4])*0.2/6;  
  x =x+(kx[1]+2*kx[2]+2*kx[3]+kx[4])*0.2/6;  
  v =v+(kv[1]+2*kv[2]+2*kv[3]+kv[4])*0.2/6;  
  y =y+(ky[1]+2*ky[2]+2*ky[3]+ky[4])*0.2/6;  
 
  if(k%2==0) printf("%f %f\n",x,y);  
 }  
 return 0;  
}

です。この結果が


PIC

Figure 3.1: 4次の RK 法


です。点の位置はほぼ理論値と一致していますね。確かにこの4次の RK 法がよく使われる理由がわかっていただける かと思います。