ホイン(Heun)法 [ネコ騙し数学]
ホイン(Heun)法
§1 オイラー法
微分方程式
の解は
である。
(2)式の右辺の積分を
で近似すると、(2)式は
と近似することにより、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
(解答終)
表計算ソフトを用いた計算結果を以下に示す。
オイラー法だと、計算を進めると、誤差が蓄積して、その結果、厳密解
との食い違いが大きくなっていくことが理解できると思う。
§2 ホイン法
まず、オイラー法
を用いてを計算し、次に台形公式
と修正する。
以下同様に、所定の精度が得られるまで
と反復計算をする方法をホイン法という。
ホイン法を用いて、問題の微分方程式を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(予測修正子法ともいう)のほうが精度よく計算できていることがわかるだろう。
コメント 0