前回Euler法で Euler法による常微分方程式を解きました． このとき得られる数値解において誤差が問題となりましたので もう少し正確に数値解を求めることが出来るHeun法を使ってみたいと思います． 実はこれでもあまり正確ではなく，一般的には4次近似のRunge-Kutta法が用いられます． Runge-Kutta法は次回作成してみたいと思います．
Euler法では次の図のように考えて数値解を逐次的に算出しました．
y(dt)を求める際にy'(0)のみを考慮していました． これをHeun法ではy(dt)を求めるときの傾きをy'(0)だけではなくy'(dt)も考慮しようと しています．とはいえ，y(dt)がわからなければy'(dt)もわかりませんので， y'(0)を使って出てきたy(dt)を使ってy'(dt)を求め，y'(0)とy'(dt)の平均を取って 傾きを計算し，改めてy(dt)を求めます．少しわかりにくいので次の図を見てください．
これはEuler法でy(2dt)まで求めたときのプロセスですが，このときの y'(dt)とy'(0)の平均を使ってy(dt)を求めます．つまり次のようになります．
これは式で書くと，
となりますが，k2nでfに渡している2つ目の引数を見てください． この部分は上の図でいうとy(dt)に相当します． また，k1はy'(0)に，k2はy'(dt)に相当するのがわかります． 直感的にも誤差が小さくなりそうです．
前回Euler法を作っているためそれに少々手を加えるだけで内容はそんなに変わりません．
void Heun(double (*f[])(double t, double *x),
double t0, double *x, double tn, int div, int num)
{
double h, t;
double *k1, *k2, *temp;
int i, j;
if((k1 = (double*)malloc(sizeof(double)*num)) == NULL){
fprintf(stderr, "Allocation Error\n");
exit(1);
}
if((k2 = (double*)malloc(sizeof(double)*num)) == NULL){
fprintf(stderr, "Allocation Error\n");
exit(1);
}
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++){
k1[j] = (*f[j])(t, x);
temp[j] = x[j] + h*k1[j];
}
for(j=0; j<num; j++){
k2[j] = (*f[j])(t+h, temp);
x[j] += (k1[j] + k2[j])*h/2;
}
t += h;
}
free(k1);
free(k2);
free(temp);
}
一見ずいぶん長くなってるように見えますが，mallocとそのエラーチェックが 殆どです．
変数の数や次数を可変にしたために取り得る初期値の数が 一意に決まらない為，動的に領域を確保する必要がでてきます． 領域を確保してから引数で渡してもよかったのですが，そうすると 引数の数がやたらと増えてしまいますのでこのようにしました． この方法だと，あるtにおけるy(t)の値だけを知りたいときには問題ありませんが， グラフを作成するような場合(Euler法の記事の最後のソース)にはmallocとfreeが 延々と繰り返されますので無駄が大きいです．
