2010. 7.27.
石立 喬
Visual C++ 2010 Express の易しい使い方（１５）
―――― 4次Runge-Kutta法による微分方程式の求解とLorenz Attractor図形の描画 ――――
概要
微分方程式の数値解法は、現代のコンピュータ出現以前から、数多く研究され、実用化されてきた。初期のコンピュータの科学技術計算への応用は、ほとんどが微分方程式を解くことであったと言っても過言ではない。ここでは、微分方程式の解法で最も良く使われているRunge-Kutta法を用いて、カオス図形で著名なLorenz
Attractorの図形を、三元連立微分方程式の解として描画する。
微分方程式
微分方程式は、微分係数が含まれる方程式で、これを解くと、数値ではなく、微分係数を含まない数式が得られる。微分方程式は、常微分方程式と偏微分方程式に分けられる。単独方程式、連立方程式があり、一階微分方程式、二階微分方程式と高階微分方程式、線形微分方程式と非線形微分方程式などがある。
微分方程式の数値解法
すべての微分方程式が解析的に解けるわけではなく、自然界を支配する多くの微分方程式は、数値的に近似的にしか解けない。微分方程式の数値解法には、Euler法、改良Euler法、Heun法、Runge-Kutta法などがある。ここでは、精度が良く広く実用されている４次Runge-Kutta法を紹介する。
４次Runge-Kutta法
微分方程式dx/dt=f(t,x)をdx=h*f(t,x)のように変形して計算するのは、Eular法などと同じであるが、下記のような漸化式を用いる。刻み幅hを小さくしなくても精度が上がるので、時間が掛からず誤差が少ない。
４次Runge-Kutta法の連立微分方程式への応用
解くべき連立微分方程式を一般化して
dx/dt=fx(t,x,y,z)
dy/dt=fy(t,x,y,z)
dz/dt=fz(t,x,y,z)
とすると、次の漸化式によって解が求まる。
k1=h*fx(t,x,y,z)
l1=h*fy(t,x,y,z)
m1=h*fz(t,x,y,z)
k2=h*fx(t+h/2,x+k1/2,y+l1/2,z+m1/2)
l2=h*fy(t+h/2,x+k1/2,y+l1/2,z+m1/2)
m2=h*fz(t+h/2,x+k1/2,y+l1/2,z+m1/2)
k3=h*fx(t+h/2,x+k2/2,y+l2/2,z+m2/2)
l3=h*fy(t+h/2,x+k2/2,y+l2/2,z+m2/2)
m3=h*fz(t+h/2,x+k2/2,y+l2/2,z+m2/2)
k4=h*fx(t+h,x+k3,y+l3,z+m3)
l4=h*fy(t+h,x+k3,y+l3,z+m3)
m4=h*fz(t+h,x+k3,y+l3,z+m3)
x+=(k1+2*k2+2*k3+k4)/6
y+=(l1+2*l2+2*l3+l4)/6
z+=(m1+2*m2+2*m3+m4)/6
Lorenz Attractor（代表的なカオス図形）
微分方程式を解くのであれば、無味乾燥な図形を描くよりは、多少芸術的？な描画をしてみたくなる。Lorenzの方程式と呼ばれる著名な連立微分方程式があるので、これを試みる。Edward
N. Lorenzはカオス現象の発見者の一人として有名で、多くの文献に紹介されている。Attractorとは「引き込むもの」の意味で、図形が一定値に集まって行くところから命名されている。
Lorenz自身およびその他の研究者により、面白いとされた連立微分方程式は、下記の通りである。
初期条件についても、x=-8、y=8、z=27が推奨されているので、ここではこれに倣った。
Lorenz Attractorのプログラムの概要
図形の美しさに重点を置いたため、ここでは邪魔になる座標軸や目盛は省略した。
１） 連立微分方程式の解は三次元の関数として得られるので、３通りの方向から見た図形を描くため、見る角度（thetaとphi）を３通り定義し、有名なButterfly、Owl（ふくろう）、筆者が比較的綺麗だと思った図形を、コンボボックスで切り替えて表示できるようにした。
コンボボックスを単純に使用すると、コンボボックスを表示するよりも先にAttractor図形の描画が始まってしまい、見苦しいので、これを避けるために、コンボボックスの初期設定をtype=0（コンボボックス上のコンテンツは、「---タイプを選択してください---」としておく）とし、このままでは描画しないようにした。
２） 三次元表示に必要な係数（sint,sinp,cost,cosp）をあらかじめ計算しておく。
３） 連立微分方程式の解の初期値を与える。
４） Runge-Kuttaの漸化式によって解を求める。ただし、Lorenz の式には、右辺に t が含まれていないので、関数function_x、function_y、function_zはtを含まない。また、function_xには変数zを必要としない。
５） 時間tに伴って表示の色を変えながらパソコン画面に解を三次元表示する。時間tによって色を変えるには、「Visual
C++ 2010 Expressの易しい使い方（１２）」で用いたメソッドを流用した。
６）描画に時間がかかるので、描画中には「描画中」、終了したら「描画終了」と表示した。
Lorenz Attractorのプログラム
int type;
private: System::Void Form1_Load(System::Object^ sender, System::EventArgs^ e) {
type=0;
comboBox1->SelectedIndex=0;
}
private: System::Void Form1_Paint(System::Object^ sender, System::Windows::Forms::PaintEventArgs^
e) {
Graphics^ g=e->Graphics;
g->FillRectangle(Brushes::Black,20,50,400,320); //背景を黒にする
if(type==0) return; //タイプが選択されていないときは描画をしない
int X0=220,Y0=340; //微分方程式の解を描画する原点
//4次Runge-Kutta法に使用
double t,x,y,z;
double k1,k2,k3,k4;
double l1,l2,l3,l4;
double m1,m2,m3,m4;
double h=0.00001; //精度的には、もっと粗くても良いが、意図的に描画時間をかけるため
//色を付けての三次元表示に使用
double sint,sinp,cost,cosp;
int xx,yy,xx_old,yy_old;
Pen^ pen1;
array<double>^ theta={-90.0,0.0,-10.0};
array<double>^ phi={90.0,90.0,150.0};
sint=Math::Sin(theta[type-1]*Math::PI/180);
sinp=Math::Sin(phi[type-1]*Math::PI/180);
cost=Math::Cos(theta[type-1]*Math::PI/180);
cosp=Math::Cos(phi[type-1]*Math::PI/180);
//初期値の設定
x=-8.0;
y=8.0;
z=27.0;
//描画中と表示
g->DrawString("描画中",Font,Brushes::White,350,340);
//4次Runge-Kutta法
for(t=0.0;t<40.0;t+=h){
//Runge-Kutta法による連立微分方程式の求解
k1=h*function_x(x,y);
l1=h*function_y(x,y,z);
m1=h*function_z(x,y,z);
k2=h*function_x(x+k1/2,y+l1/2);
l2=h*function_y(x+k1/2,y+l1/2,z+m1/2);
m2=h*function_z(x+k1/2,y+l1/2,z+m1/2);
k3=h*function_x(x+k2/2,y+l2/2);
l3=h*function_y(x+k2/2,y+l2/2,z+m2/2);
m3=h*function_z(x+k2/2,y+l2/2,z+m2/2);
k4=h*function_x(x+k3,y+l3);
l4=h*function_y(x+k3,y+l3,z+m3);
m4=h*function_z(x+k3,y+l3,z+m3);
x+=(k1+2*k2+2*k3+k4)/6;
y+=(l1+2*l2+2*l3+l4)/6;
z+=(m1+2*m2+2*m3+m4)/6;
//三次元表示への変換
xx=(int)((-sint*x+cost*y)*6+0.5);
yy=(int)((-cost*cosp*x-sint*cosp*y+sinp*z)*6+0.5);
//表示色の設定
pen1=gcnew Pen(countToColor((int)t,40));
//解の描画
if(t<=0.0){
xx_old=xx;
yy_old=yy;
}
else{
g->DrawLine(pen1,X0+xx_old,Y0-yy_old,X0+xx,Y0-yy);
xx_old=xx;
yy_old=yy;
}
}
//描画終了と表示
g->FillRectangle(Brushes::Black,345,335,60,20);
g->DrawString("描画終了",Font,Brushes::White,350,340);
}
//求解のための数式
private: double function_x(double x,double y){
return -10.0*x+10.0*y;
}
private: double function_y(double x,double y,double z){
return -x*z+28.0*x-y;
}
private: double function_z(double x,double y,double z){
return x*y-2.666666667*z;
}
private: Color countToColor(int n,int base){
int d=n % base;
d*=(256/base);
return NumberToRGBColor(d);
}
private: Color NumberToRGBColor(int n){
int h=n/43;
int r,g,b;
switch(h){
case 0:r=255; g=n*6; b=0; break; //n= 0～ 42
case 1:r=(85-n)*6; g=255; b=0; break; //n= 43～ 85
case 2:r=0; g=255; b=(n-86)*6; break; //n= 86～128
case 3:r=0; g=(171-n)*6; b=255; break; //n=129～171
case 4:r=(n-172)*6; g=0; b=255; break; //n=172～214
case 5:r=255; g=0; b=(255-n)*6; break; //n=215～255
};
return Color::FromArgb(r,g,b);
}
private: System::Void comboBox1_SelectedIndexChanged(System::Object^ sender, System::EventArgs^
e) {
type=comboBox1->SelectedIndex;
Invalidate();
}
実行結果
連立微分方程式の解は三次元図形として得られるので、見る方向によって感じが変わる。Butterfly、Owlの名称で知られているもの、綺麗だと思った角度から見たものを示す。
図1 Butterfly（蝶の羽根のように見える）と呼ばれるLorenz図形
図2 Owl（フクロウの眼のように見える）と呼ばれるLorenz図形
図3 筆者がきれいと思った角度から見たLorenz図形
「Visual C++ の勉強部屋」（目次）へ