Heun法

寄稿:東條遼平

常微分方程式の2次近似

前回Euler法で Euler法による常微分方程式を解きました. このとき得られる数値解において誤差が問題となりましたので もう少し正確に数値解を求めることが出来るHeun法を使ってみたいと思います. 実はこれでもあまり正確ではなく,一般的には4次近似のRunge-Kutta法が用いられます. Runge-Kutta法は次回作成してみたいと思います.

Heun法とは

Euler法では次の図のように考えて数値解を逐次的に算出しました.

euler1.png

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)を求めます.少しわかりにくいので次の図を見てください.

heun1.png

これはEuler法でy(2dt)まで求めたときのプロセスですが,このときの y'(dt)とy'(0)の平均を使ってy(dt)を求めます.つまり次のようになります.

heun2.png

これは式で書くと,

img113.png
img114.png
img115.png

となりますが,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が 延々と繰り返されますので無駄が大きいです.

Valid XHTML 1.1! home > コンピュータ > プログラミング >
リロード   新規 編集 凍結 差分 添付 複製 改名   トップ 一覧 検索 最終更新 バックアップ   ヘルプ   最終更新のRSS
Modified by 物理のかぎプロジェクト PukiWiki 1.4.5_1 Copyright © 2001-2005 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-jiPowered by PHP 5.2.17HTML convert time to 0.038 sec.