タイトル

計算物理学

 常備分方程式編

タイトル

1,常備分方程式 Euler法 Runge−Kutta法 レポート
2,偏微分方程式 (13KB) Euler法 Crank-Nicholson法 レポート
3,確率論 (13KB) モンテカルロ法 レポート

葉

1,常微分方程式

物理の法則は、微分方程式の形で表されることが多いです。たとえば、


F=m 2
    dt2

というのは、2階微分方程式になっています。

微分方程式の種類は物理数学2(17KB)でも説明したのですが
常微分方程式と偏微分方程式の2種類があります。
偏微分方程式は「∂」という記号をつかい、解くのが複雑になります。
(例としてSch.eg.や波動方程式)

一般に、常微分方程式は次のようになります。


dx(t) = f(x,t)
dt

これは、xの次元が1の時なので、n次元の常微分方程式は
今の式がn個つくる必要があります。

ボーダー

Euler法

ここでは簡単のため、1次元の微分方程式について考えます。


dx(t) = f(x,t)
dt

という関係から


x(t+凾) = x(t) + (凾)f(x,t)

これは、xをTalor展開して、1次の項まで取り入れたものです。
ただ、1次の項までしか取り入れていないので誤差は大きくなります。


問題
x+ωx=0
で表される運動の、位相空間での動きを表しなさい。
ただし、質量はm=1
p=dx とする。
  dt

解答
今回は、xとpの2次元の常微分方程式を解くことになります。

dx=p    dp=−ω2x
dt      dt

以上の関係からEuler法では次のようになります。

x(t+凾) = x(t) +  p(t)凾
p(t+凾) = p(t) −ω2x(t)凾


プログラムは以下の通りです。
(SGLを使っています。色の宣言などは省略しているのでこのままでは動きません。)

main()
{
      char c;
      float   omega;
      float   x, p, xnew, pnew;
      float   delt, etot;
      int     itime,ntime;
      x = 1.0;
      p = 0.0;
      ntime = 5000;
      printf("Enter omega & delt \n");
      scanf("%f %f",&omega, &delt); 

      for(itime = 1; itime <= ntime; itime++) {
         xnew  = x + delt*p;
         pnew  = p - delt*omega*omega*x;

         x    = xnew;
         p    = pnew;
 
         sgl_color(SGL_YELLOW);
         sgl_disc(x,p,0.05);
         sgl_display2();
         if(itime - (itime/10)*10 == 0)
           {  
            etot = p*p/2.0 + omega*omega*x*x/2.0;
            printf("itime = %d  etot = %f \n",itime,etot);
           }
        }
      printf("Hit return to terminate \n");
      scanf("%c",&c);
}

これを実行すると、円にならず広がっていってしまいます。
これは、大きな誤差が存在するためです。

実行結果

ボーダー

Runge−Kutta法

Runge−Kutta法は、Euler法と比べると式は複雑なのですが
精度は1000倍になります。(有効数値が3桁アップしたので)
式は次のようになります。


x(t+凾)=x(t)+(a+2b+2c+d)凾煤^6

   a=f(x,t)
   b=f(x+a凾煤^2,t+凾煤^2)
   c=f(x+b凾煤^2,t+凾煤^2)
   d=f(x+c凾煤Ct+凾)

例 (1変数の場合)
f(x,t)=2x の場合

a=2x
b=2(x+a凾煤^2)
c=2(x+b凾煤^2)
d=2(x+c凾)

x(t+凾) = x(t) + (a+2b+2c+d)凾煤^6

また、多変数の場合は次のようになります。


xi(t+凾)=xi(t)+(ai+2bi+2ci+di)凾煤^6

   ai=fi(x1,・・,xn,t)
   bi=fi(x1+a1凾煤^2,・・,xn+an凾煤^2,t+凾煤^2)
   ci=fi(x1+b1凾煤^2,・・,xn+bn凾煤^2,t+凾煤^2)
   di=fi(x1+c1凾煤@ ,・・,xn+cn凾煤@ ,t+凾)

問題 (多変数の場合)
「Euler法」の問題をRunge−Kutta法で解きなさい。

解答
プログラムは先ほどの黄色の部分を次のように書き換えます。
ここで注意するところは、a1,a2,b1,・・・d1,d2の順に書くことです。

main()
{
      char c;
      float   omega;
      float   x, p, xnew, pnew;
      float   delt, etot;
      int     itime,ntime;
      float a1,a2,b1,b2,c1,c2,d1,d2;

      x = 1.0;
      p = 0.0;
      ntime = 5000;
      printf("Enter omega & delt \n");
      scanf("%f %f",&omega, &delt); 

      for(itime = 1; itime <= ntime; itime++) {
		a1=p;
		a2=-omega*omega*x;
		b1=p+a2*delt/2.0;
		b2=-omega*omega*(x+a1*delt/2.0);
		c1=p+b2*delt/2.0;
		c2=-omega*omega*(x+b1*delt/2.0);
		d1=p+c2*delt;
		d2=-omega*omega*(x+c1*delt);

		xnew = x + (a1+2*b1+2*c1+d1)*delt/6.0;
		pnew = p + (a2+2*b2+2*c2+d2)*delt/6.0;

         x    = xnew;
         p    = pnew;
 
         sgl_color(SGL_YELLOW);
         sgl_disc(x,p,0.05);
         sgl_display2();
         if(itime - (itime/10)*10 == 0)
           {  
            etot = p*p/2.0 + omega*omega*x*x/2.0;
            printf("itime = %d  etot = %f \n",itime,etot);
           }
        }
      printf("Hit return to terminate \n");
      scanf("%c",&c);
}

これを実行すると以下の左のような円になります。
これは、先ほどのEuler法と 倍率、計算回数、初期条件は同じです。

Runge−Kutta法 Euler法
Runge-Kutta Euler

上上へ

ボーダー

レポート

ここでは、「計算物理学」の講義で実際に出題されたレポートの問題の解放を解説します。

●課題

気象学者のローレンツは大気を簡単なモデルとして

dx=A(y−x)
dt

dy=−xz+Rx−y
dt

dz=xy−Bz
dt

という微分方程式を考えた。(A,B,Rはパラメータ)
この系をシミュレーションせよ。
ただし、A=4 B=8/3 R=80という値を代入するといい。

※実際の課題は、

  1. x-y,y-z,z-x平面に投影下図の作成
  2. パラメータに値によりCAOS的になることを確かめよ
  3. どれか一つの変数と時間の関係をグラフにせよ
  4. 初期値を微少にずらしたとき、元と大きく違うことを確かめよ

の4題です。

●解答

ここでは、計算部分のプログラムを示しています。
描画はSGLを使っていますが設定部分は省略しています。

double A=4.0,B=8.0/3.0,R=80.0;
double a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3,x,y,z,xnew,ynew,znew;
double delt=0.01;
x=1.0;	y=0.0;	z=0.0;

for (long i=0; i<5000; i++){
    a1=A*(y-x);
    a2=-x*z+R*x-y;
    a3=x*y-B*z;
    b1=A*(y+a2*delt/2-(x+a1*delt/2));
    b2=-(x+a1*delt/2)*(z+a3*delt/2)+R*(x+a1*delt/2)-(y+a2*delt/2);
    b3=(x+a1*delt/2)*(y+a2*delt/2)-B*(z+a3*delt/2);
    c1=A*(y+b2*delt/2-(x+b1*delt/2));
    c2=-(x+b1*delt/2)*(z+b3*delt/2)+R*(x+b1*delt/2)-(y+b2*delt/2);
    c3=(x+b1*delt/2)*(y+b2*delt/2)-B*(z+b3*delt/2);
    d1=A*(y+c2*delt-(x+c1*delt));
    d2=-(x+c1*delt)*(z+c3*delt)+R*(x+c1*delt)-(y+c2*delt);
    d3=(x+c1*delt)*(y+c2*delt)-B*(z+c3*delt);
    xnew=x + (a1+2*b1+2*c1+d1)*delt/6.0;
    ynew=y + (a2+2*b2+2*c2+d2)*delt/6.0;
    znew=z + (a3+2*b3+2*c3+d3)*delt/6.0;
    x=xnew;  y=ynew;  z=znew;
    sgl_color(SGL_YELLOW);
    sgl_disc(x,y,1);
    sgl_display2();
}

●シミュレーション結果

ローレンツの微分方程式のシミュレーション結果

シミュレーション結果は以上のようになりました。
なお、プログラムのバグや、初期値が不正なこともあるので
これが正解かは不確かです。(ミスを発見したらメールをください。)

●Soft

今回レポート用に制作したアプリケーションがあります。

KeisanPhy.exe(22KB)


・制作日 1997/11/26(Wed)
・Version 0.01
・著作権 小林成徳

このSoftはリリース版なので単独で動作します。

葉

上上へ

戻るTopに戻る