Edited at

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

More than 3 years have passed since last update.


サンプル

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

の解をオイラー法で求める。


オイラー法の概要

\frac{dy}{dx} = xy, y = f(x)\ とする。 

\\
\frac{dy}{dx} = \frac{f(x+Δx)-f(x)}{Δx} …(☆)
\\
= xf(x)
\\
⇔f(x+Δx) = f(x)(1+xΔx)\ …①
\\
初期条件 x=0, y=1より
\\
f(Δx) = 1\ …②
\\
x = Δx のとき①に代入すると
\\
f(2Δx) = f(Δx)\{1+(Δx)^2\}
\\
= {1+(Δx)^2}\ (∵②)

という具合に、1点が分かれば次の点が分かる。


正確な数値を求める

オイラー法で計算した値と比較するために正確な数値を求める(手計算)。

\frac{dy}{dx} = xy

\\
⇔\frac{1}{y}dy = xdx
\\
⇔\int \frac{1}{y}dy = \int xdx
\\
⇔logy = \frac{1}{2}x^2 + C
\\
⇔y = exp(\frac{1}{2}x^2 + C)
\\
初期条件より
\\
1 = exp(C)
\\
⇔1 = e^C ⇔ C = 0
\\
ゆえに
\\
y = exp(\frac{1}{2}x^2)


プログラムで使うf2を求める

式①のこと


プログラム


sample.c


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

int main(void) {
double f, f1, f2, dx, x;

f1 = 1;
dx = 0.01;

for (int i=0; i<=100; ++i) {
x = i*dx;
f2 = f1 + dx*x*f1;

f = exp(x*x/2);

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

f1 = f2;
}
}


実行結果


問題

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

の数値解をオイラー法で求め、正確な値と比較する。


正確な値を求める

\int \frac{1}{y^2}dy = \int -3xdx

\\
⇔-\frac{1}{y} = -\frac{3}{2}x^2 + C
\\
初期条件より
\\
-\frac{1}{2} = 0 + C
\\
よって
\\
-\frac{1}{y} = -\frac{3}{2}x^2 - \frac{1}{2}
\\
⇔ \frac{1}{y} = \frac{3x^2+1}{2}
\\
⇔ y = \frac{2}{3x^2+1}


f2を求める

\frac{dy}{dx} = \frac{f(x+Δx)-f(x)}{Δx} …(☆)

\\
= -3x{f(x)}^2
\\
⇔ f(x+Δx) = Δx(-3x{f(x)}^2)+f(x)
\\
= f(x) -3xΔx{f(x)}^2


プログラム


problem.c

#include <stdio.h>

#include <math.h>

int main(void) {
double f, f1, f2, dx, x;

f1 = 2;
dx = 0.01;

for (int i=0; i<=100; ++i) {
x = i*dx;

f2 = f1 - 3*x*dx*pow(f1, 2);

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

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

f1 = f2;
}
}


実行結果