ここでは常微分方程式の近似解法としてよく知られている オイラー法と 4 次のルンゲ=クッタ法について簡単に説明し、 そのプログラムを書く。 y=y(x) とし、常微分方程式
y'=f(x,y)
y(x0)=y0
の近似解を求める。 オイラー法、ルンゲ=クッタ法のいずれも、適当な刻み幅 h を定めて xn=x0+nh とし yn=y(xn) を順次求める方法である。
オイラー法は簡単で
yn+1=yn +hf(xn,yn)
とする。 すなわち x=xn における微分係数を 区間 [xn, xn+1] で用いる。 したがって近似の精度はそれほど良くない。
オイラー法の改良は色々と考えられるが、 ここでは比較的精度の良い 4 次のルンゲ=クッタ法を説明する。
f1=hf (xn,yn)
f2=hf (xn+h/2,yn+f1/2)
f3=hf (xn+h/2,yn+f2/2)
f4=hf (xn+h,yn+f3)
として
yn+1=yn +(f1+2f2+2f3+f4)/6
とするのである。 理論については適当な文献を参照して欲しい。
実際のプログラムを見てみよう。
#include <stdio.h>#define STEP 16float myfunct(float x, float y){ return 1 - y * y;}float euler(float x, float y, float targ, int step, float (*funct)(float, float)){ int i; float h; h = (targ - x) / step; for (i = 0; i < step; i++) { y = y + h * funct(x, y); x += h; } return y;}float rungekutta(float x, float y, float targ, int step, float (*funct)(float, float)){ int i; float h, f1, f2, f3, f4; h = (targ - x) / step; for (i = 0; i < step; i++) { f1 = h * funct(x, y); f2 = h * funct(x + h/2, y + f1/2); f3 = h * funct(x + h/2, y + f2/2); f4 = h * funct(x + h, y + f3); y = y + (f1 + 2*f2 + 2*f3 + f4)/6; x += h; } return y;}int main(int argc, char *argv[]){ float x, y, t; printf("Input initial values : "); scanf("%f %f", &x, &y); printf("Input a target value : "); scanf("%f", &t); printf("Euler : %f\n", euler(x, y, t, STEP, myfunct)); printf("Runge-Kutta : %f\n", rungekutta(x, y, t, STEP, myfunct)); return 0;}オイラー法を行う関数 euler とルンゲ=クッタ法を行う関数 rungekutta の引数は同じで、
euler( "x の初期値", "y の初期値", "関数の値を求めたい x の値", "区間を区切る数", "関数名" )
である。 初めの myfunct で微分方程式を定義している。 また、この例では STEP で区切りの数を固定して計算している。
課題 10.01
厳密解を求めることが出来る常微分方程式で上のプログラムを実行し、 近似解の誤差を調べよ。
課題 10.02
関数 rungekutta を書き換え、 一定間隔の近似解を出力するようにせよ。 また、出力されたデータから GNUPLOT, maxima, 表計算ソフトなどを用いてグラフを作成せよ。