Edited at

[C言語] 1階常微分方程式の解法(ルンゲクッタ法)

More than 3 years have passed since last update.

前回の投稿と同じ微分方程式を使っている。

参考:[C言語] 1階常微分方程式の解法(オイラー法)


ルンゲクッタ法の概要

\frac{dy}{dx} = f(x, y)

とし、初期条件$(x_0, y_0)$とする。

$x$が$x_0$から$h$だけ増加した点 $x_1 = x_0 + h$における$y_1$を以下の計算で求める。


k_ 1 = hf(x_0, y_0)
\\
k_2 = hf(x_0 + \frac{h}{2}, y_0 + \frac{k_1}{2})
\\
k_3 = hf(x_0 + \frac{h}{2}, y_0 + \frac{k_2}{2})
\\
k_4 = hf(x_0 + h, y_0 + k_3)
\\
k = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}
\\
y_1 = y_0 + k


サンプル

\frac{dy}{dx} = xy\ (初期条件: x=0\ で\ y=1)

の解をルンゲクッタ法で求める。

#include <stdio.h>

#include <math.h>

int main(void) {

double f, f1, f2, dx, x, k1, k2, k3, k4, k;

f1 = 1;
dx = 0.01;

for (int i=0; i<=100; i++) {
x = i*dx;
k1 = dx*x*f1;
k2 = dx*(x + dx/2.)*(f1 + k1/2.);
k3 = dx*(x + dx/2.)*(f1 + k2/2.);
k4 = dx*(x + dx)*(f1 + k3);
k = (k1 + 2*k2 + 2*k3 + k4)/6.;

f2 = f1 + k;
f = exp(x*x/2.);

printf("x = %f, 計算値-> f1 = %f, 正確な値-> f = %f \n", x, f1, f);

f1 = f2;
}
return 0;
}


実行結果


問題

\frac{dy}{dx} = -3xy^2\ (初期条件: x=0\ で\ y=2)

の解をルンゲクッタ法で求める。

#include <stdio.h>

#include <math.h>

int main(void) {

double f, f1, f2, dx, x, k1, k2, k3, k4, k;

f1 = 2;
dx = 0.01;

for (int i=0; i<=100; i++) {
x = i*dx;
k1 = dx*(-3)*x*f1*f1;
k2 = dx*(-3)*(x + dx/2.)*pow( (f1 + k1/2.), 2);
k3 = dx*(-3)*(x + dx/2.)*pow( (f1 + k2/2.), 2);
k4 = dx*(-3)*(x + dx)*pow( (f1 + k3), 2);
k = (k1 + 2*k2 + 2*k3 + k4)/6.;

f2 = f1 + k;

f = 2./(3*pow(x, 2)+1);

printf("x = %f, 計算値-> f1 = %f, 正確な値-> f = %f \n", x, f1, f);

f1 = f2;
}
return 0;
}


実行結果

という感じで、オイラー法よりも精度が高い。