| １，常備分方程式 | Ｅｕｌｅｒ法 Ｒｕｎｇｅ－Ｋｕｔｔａ法 レポート |
| ２，偏微分方程式 (13KB) | Ｅｕｌｅｒ法 Crank-Nicholson法 レポート |
| ３，確率論 (13KB) | モンテカルロ法 レポート |
物理の法則は、微分方程式の形で表されることが多いです。たとえば、
Ｆ＝ｍ ｄ2ｘ ｄｔ2
というのは、２階微分方程式になっています。
微分方程式の種類は物理数学２(17KB)でも説明したのですが
常微分方程式と偏微分方程式の２種類があります。
偏微分方程式は「∂」という記号をつかい、解くのが複雑になります。
(例としてSch.eg.や波動方程式)
一般に、常微分方程式は次のようになります。
ｄｘ(t) ＝ ｆ(ｘ，ｔ) ｄｔ
これは、ｘの次元が１の時なので、ｎ次元の常微分方程式は
今の式がｎ個つくる必要があります。
ここでは簡単のため、１次元の微分方程式について考えます。
ｄｘ(t) ＝ ｆ(ｘ，ｔ) ｄｔ
という関係から
ｘ(ｔ＋⊿ｔ) ＝ ｘ(ｔ) ＋ (⊿ｔ)ｆ(ｘ，ｔ)
これは、ｘをＴａｌｏｒ展開して、１次の項まで取り入れたものです。
ただ、１次の項までしか取り入れていないので誤差は大きくなります。
問題 ｘ＋ωｘ＝０ で表される運動の、位相空間での動きを表しなさい。 ただし、質量はｍ＝１ ｐ＝ｄｘ とする。 ｄｔ 解答 今回は、ｘとｐの２次元の常微分方程式を解くことになります。 ｄｘ＝ｐ ｄｐ＝－ω2ｘ ｄｔ ｄｔ 以上の関係からEuler法では次のようになります。 ｘ(ｔ＋⊿ｔ) ＝ ｘ(ｔ) ＋ ｐ(ｔ)⊿ｔ ｐ(ｔ＋⊿ｔ) ＝ ｐ(ｔ) －ω2ｘ(ｔ)⊿ｔ プログラムは以下の通りです。 (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);
}
これを実行すると、円にならず広がっていってしまいます。
これは、大きな誤差が存在するためです。
Ｒｕｎｇｅ－Ｋｕｔｔａ法は、Ｅｕｌｅｒ法と比べると式は複雑なのですが
精度は１０００倍になります。(有効数値が３桁アップしたので)
式は次のようになります。
ｘ(ｔ＋⊿ｔ)＝ｘ(ｔ)＋(ａ＋２ｂ＋２ｃ＋ｄ)⊿ｔ／６ ａ＝ｆ(ｘ，ｔ) ｂ＝ｆ(ｘ＋ａ⊿ｔ／２，ｔ＋⊿ｔ／２) ｃ＝ｆ(ｘ＋ｂ⊿ｔ／２，ｔ＋⊿ｔ／２) ｄ＝ｆ(ｘ＋ｃ⊿ｔ，ｔ＋⊿ｔ)
例 (１変数の場合) ｆ(x,t)＝２ｘ の場合 ａ＝２ｘ ｂ＝２(ｘ＋ａ⊿ｔ／２) ｃ＝２(ｘ＋ｂ⊿ｔ／２) ｄ＝２(ｘ＋ｃ⊿ｔ) ｘ(ｔ＋⊿ｔ) ＝ ｘ(ｔ) ＋ (ａ＋２ｂ＋２ｃ＋ｄ)⊿ｔ／６
また、多変数の場合は次のようになります。
ｘi(ｔ＋⊿ｔ)＝ｘi(ｔ)＋(ａi＋２ｂi＋２ｃi＋ｄi)⊿ｔ／６ ａi＝ｆi(ｘ1，・・，ｘn，ｔ) ｂi＝ｆi(ｘ1＋ａ1⊿ｔ／２，・・，ｘn＋ａn⊿ｔ／２，ｔ＋⊿ｔ／２) ｃi＝ｆi(ｘ1＋ｂ1⊿ｔ／２，・・，ｘn＋ｂn⊿ｔ／２，ｔ＋⊿ｔ／２) ｄi＝ｆi(ｘ1＋ｃ1⊿ｔ ，・・，ｘn＋ｃn⊿ｔ ，ｔ＋⊿ｔ)
問題 (多変数の場合) 「Ｅｕｌｅｒ法」の問題をＲｕｎｇｅ－Ｋｕｔｔａ法で解きなさい。 解答 プログラムは先ほどの黄色の部分を次のように書き換えます。 ここで注意するところは、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法と 倍率、計算回数、初期条件は同じです。
| Ｒｕｎｇｅ－Ｋｕｔｔａ法 | Ｅｕｌｅｒ法 |
 |
|
ここでは、「計算物理学」の講義で実際に出題されたレポートの問題の解放を解説します。
気象学者のローレンツは大気を簡単なモデルとして
ｄｘ＝Ａ(ｙ－ｘ) ｄｔ ｄｙ＝－ｘｚ＋Ｒｘ－ｙ ｄｔ ｄｚ＝ｘｙ－Ｂｚ ｄｔ
という微分方程式を考えた。(Ａ，Ｂ，Ｒはパラメータ)
この系をシミュレーションせよ。
ただし、Ａ＝４ Ｂ＝８／３ Ｒ＝８０という値を代入するといい。
※実際の課題は、
の４題です。
ここでは、計算部分のプログラムを示しています。
描画は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 ・著作権 小林成徳
このＳｏｆｔはリリース版なので単独で動作します。
KeisanPhy.exe(22KB)
Topに戻る