ホイン（Heun）法

ホイン（Heun）法

 

§１　オイラー法

 

微分方程式

の解は

である。

（２）式の右辺の積分を

で近似すると、（２）式は

と近似することにより、x₀を計算の起点とし、x₁、x₂、・・・と築時的に計算をすることができる。

とあらわすことにすると、

になる。

とくに、点が等間隔hに並んでいるとき、

である。

この方法をオイラー法という。

 

問題　オイラー法を用いて、次の微分方程式のx=0.1、0.2、0.3におけるyの値を求めよ。

【解】

x₀=0、y₀=y(x₀)=1、h=0.1とおくと、

x₁=0.1、y₁=1.05だから

x₂=0.2、y₂=1.163だから

したがって、

x=0.1のとき、y=1.05

x=0.2のとき、y=1.1106

x=0.3のとき、y=1.1846

（解答終）

 

表計算ソフトを用いた計算結果を以下に示す。

 

 

 

オイラー法だと、計算を進めると、誤差が蓄積して、その結果、厳密解

との食い違いが大きくなっていくことが理解できると思う。

 

 

§２　ホイン法

まず、オイラー法

を用いてを計算し、次に台形公式

と修正する。

以下同様に、所定の精度が得られるまで

と反復計算をする方法をホイン法という。

 

ホイン法を用いて、問題の微分方程式をh=0.1と同一の条件で解いた結果は次の通り。

 

 

 

 

この計算に用いたプログラムは以下の通り。

 

#include <stdio.h>
#include <math.h>

double f(double x, double y) {
    return 0.5*(1+x)*y*y;
}

void Heun(double *x , double *y , double h , int n) {
    double eps = 1.0e-6, err;
    double yn = 0.;
    int i,j, limit = 100;

    for(i=0; i< n; i++) {
        x[i+1]=x[i]+h;
        // Euler法で予備計算
        y[i+1] = y[i] + h*f(x[i],y[i]);
        for (j=1; j <= limit; j++) { // 所定の精度が得られるまで反復
            // 台形公式で修正
            yn = y[i]+0.5*(f(x[i],y[i])+f(x[i+1],y[i+1]))*h;
            err =fabs(y[i+1]-yn);
            if (err <= eps) break;
            // y[i+1]の更新
            y[i+1]=yn;
        }

    }
}

double Exact(double x) {
    return 4./(4 - 2*x -x*x);
}

main() {
    double x[11], y[11];
    double h;
    int i, n;

    h=0.1;
    n=10;

    for (i = 0; i<=n; i++) {
        x[i] = 0; y[i]=0;
    }
   
    x[0]=0; y[0]=1.; // xとyの初期値
   
    Heun(x,y,h, n); // Heun法を呼び出す

    printf("   x      y(Heun)  y(exact) error(％) \n");
    for (i=0; i <= n; i++) {
        printf("%f %f %f %f\n", x[i],y[i],Exact(x[i]),fabs(1.0-y[i]/Exact(x[i]))*100);
    }
}

 

オイラー法より、Heun（予測修正子法ともいう）のほうが精度よく計算できていることがわかるだろう。

タグ：数値解析