1,常備分方程式 | Euler法 Runge−Kutta法 レポート |
2,偏微分方程式 (13KB) | Euler法 Crank-Nicholson法 レポート |
3,確率論 (13KB) | モンテカルロ法 レポート |
物理の法則は、微分方程式の形で表されることが多いです。たとえば、
F=m d2x dt2
というのは、2階微分方程式になっています。
微分方程式の種類は物理数学2(17KB)でも説明したのですが
常微分方程式と偏微分方程式の2種類があります。
偏微分方程式は「∂」という記号をつかい、解くのが複雑になります。
(例としてSch.eg.や波動方程式)
一般に、常微分方程式は次のようになります。
dx(t) = f(x,t) dt
これは、xの次元が1の時なので、n次元の常微分方程式は
今の式がn個つくる必要があります。
ここでは簡単のため、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法は、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法 |
ここでは、「計算物理学」の講義で実際に出題されたレポートの問題の解放を解説します。
気象学者のローレンツは大気を簡単なモデルとして
dx=A(y−x) dt dy=−xz+Rx−y dt dz=xy−Bz dt
という微分方程式を考えた。(A,B,Rはパラメータ)
この系をシミュレーションせよ。
ただし、A=4 B=8/3 R=80という値を代入するといい。
※実際の課題は、
の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(); }
シミュレーション結果は以上のようになりました。
なお、プログラムのバグや、初期値が不正なこともあるので
これが正解かは不確かです。(ミスを発見したらメールをください。)
今回レポート用に制作したアプリケーションがあります。
・制作日 1997/11/26(Wed) ・Version 0.01 ・著作権 小林成徳
このSoftはリリース版なので単独で動作します。