RIGHT:寄稿：東條遼平 * 常微分方程式とは [#u9491d2d] 常微分方程式とは一変数関数の導関数が含まれている方程式の事をいい、 例えば次の式であれば、yの導関数がtとyによって記述されています。 #ref(img100.png,nolink) 多くの物理現象は微分方程式によって記述され、 例を挙げれば、放射性物質の崩壊、空気抵抗、電気回路、減衰振動などがあります。 そのため微分方程式を解く事は大変重要なのですが、簡単な微分方程式であっても 数式で表示された解である解析解を求めるのは大変困難になります。 そこでPCを使い数値解を求めてみようと思います。 * Euler法とは [#i1b3ea9b] #ref(img101.png,nolink) このような方程式を考えてみます。y'(t)とはy(t)の傾きの事ですので、 初期値y(0)、y'(0)が判っているのであれば、 十分小さい幅dtを考える事でy(dt)がわかり、y(dt)とdtを使って y'(dt)が求められます。それを繰り返す事で任意のtでのy(t)を解こうというものです。 イメージとしては次のようになります。 #ref(euler1.png,nolink) これを式で書くと、 #ref(img102.png,nolink) であり、t=tnのとき、y(t)=yn、y'(t)=y'n とすると、yn+1は次のように表せます。 #ref(img103.png,nolink) これを使って微分方程式の数値解を求めるのがEuler法です。 * 連立微分方程式 [#dc813872] 連立微分方程式とはその名の通り、ある関数の導関数に他の関数が含まれるようなもので、 例えば次のようなものがあります。 #ref(img107.png,nolink) ここで、 #ref(img108.png,nolink) のようにおくと、 #ref(img109.png,nolink) と出来ることから、連立微分方程式も引数と変数が増えるだけで処理は同じであることが わかります。 * 高次微分方程式 [#fc86b40e] #ref(img110.png,nolink) のように高次の導関数を含む微分方程式を考えてみます。このとき、 #ref(img111.png,nolink) というふうに式変形をし、 #ref(img112.png,nolink) とすると、これは上で書いた連立微分方程式と同じ形になり、y1、y2の 初期値がわかれば同じ方法で解けることになります。 つまり、 #ref(img100.png,nolink) という一番シンプルな微分方程式さえ解ければ、連立微分方程式も、 高次微分方程式も解けるということです。 * Taylor展開における一次近似 [#i7bcdf7b] また、無限回微分出来るという事から、Taylor展開を考えると、 #ref(img104.png,nolink) #ref(img105.png,nolink) #ref(img106.png,nolink) となり、2階微分以降を無視するとEuler法の式と同じになります。 また、dt2以降を削る事になるため、誤差はdtの二乗のオーダーとなり、 精度はかなり悪いと言えます。 * ソースコード [#cd285df3] 素直にEuler法のソースコードを書くならば、tとy(t)の値を引数にとるy'(t)の関数を 定義して、その関数ポインタfとt0からtnまでの区間をdiv分割し初期値yを与えるとすると 次のようになります。 double Euler(double (*f)(double t, double y), double t0, double y, double tn, int div) { double h, t; int i; h = (tn - t0); if(div) h /= div; t = t0; for(i=0; i<div; i++){ y += h*(*f)(t, y); t += h; } return y; } これでy(tn)の値を返してくれるのですが、これでは 二階の微分方程式や、連立微分方程式を解こうとすると 新たに受け取る引数と初期値の個数が増えただけで同じ処理をする関数を 作らなければいけません。そこで関数ポインタと初期値の数を配列で渡すことで、 汎用な関数を作る事にします。 void Euler(double (*f[])(double t, double *x), double t0, double *x, double tn, int div, int num) { double h, t; double *temp; int i, j; if((temp = (double*)malloc(sizeof(double)*num)) == NULL){ fprintf(stderr, "Allocation Error\n"); exit(1); } h = (tn - t0); if(div) h /= div; t = t0; for(i=0; i<div; i++){ for(j=0; j<num; j++){ temp[j] = h*(*f[j])(t, x); } for(j=0; j<num; j++){ x[j] += temp[j]; } t += h; } free(temp); } 関数ポインタの配列と初期値を配列で渡し、その要素数をnumで渡すことで、 連立微分方程式の変数がいくつになっても、また、何次の導関数を含んでいても この関数で解を求めることができます。例えばf(t,x,y,z)であれば x[0]をxに、x[1]をyに、x[2]をzに対応させます。 double x[] = {1, 0}; double (*f[2])(double, double*); f[0] = f1; f[1] = f2; Euler(f, 0, x, 10, 1000, 2); printf("%f %f\n", x[0], x[1]); のようにすればt=10のときの値だけを求められますし、 int i; double t=0, tn=10; double x[] = {1, 0}; double (*f[2])(double, double*); double h; f[0] = f1; f[1] = f2; h=(tn-t)/1000; printf("%f %f\n", t, x[0]); for(i=0; i<1000; i++){ Euler(f, t, x, t+h, 1, 2); t += h; printf("%f %f\n", x[0], x[1]); } のようにすれば、グラフ作成用のデータも取ることができて便利です。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);
