Edited at

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

More than 3 years have passed since last update.


サンプル

微分方程式

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

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


sample.c

#include <stdio.h>

#include <math.h>

int main(void) {

double f, f1, f2, g1, g2, dx, x;

f1 = 2;
g1 = 0;
dx = 0.01;

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

f = 2*cos(x);

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

f1 = f2;
g1 = g2;
}
return 0;
}



実行結果


問題

微分方程式

\frac{d^2y}{dx^2} - 6\frac{dy}{dx} + 9y = 0 (初期条件:x = 0で y = 0, \frac{dy}{dx} = 1)

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


problem.c

#include <stdio.h>

#include <math.h>

int main(void) {

double f, f1, f2, g1, g2, dx, x;

f1 = 0;
g1 = 1;
dx = 0.01;

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

f2 = f1 + dx*g1;
g2 = g1 + dx*(6*g1 - 9*f1);

f = x*exp(3*x);

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

f1 = f2;
g1 = g2;
}
return 0;
}



実行結果


計算のメモ

\frac{d^2y}{dx^2} + a\frac{dy}{dx} + by = R(x)

において

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

とすると

\\

\begin{equation}

\begin{cases}

\frac{dg(x)}{dx} + a\frac{dg(x)}{dx} + bf(x) = R(x) \\
\frac{df(x)}{dx} = g(x)

\end{cases}

\end{equation}

となる。

今、

\frac{dg(x)}{dx} - 6g(x) + 9f(x) = 0

なので、

\frac{dg(x)}{dx} = 6g(x) - 9f(x)

となり、これがコード中のg2


ちなみに

もっと細かく計算すると


problem.c

#include <stdio.h>

#include <math.h>

int main(void) {

double f, f1, f2, g1, g2, dx, x;

f1 = 0;
g1 = 1;
dx = 0.001;

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

f2 = f1 + dx*g1;
g2 = g1 + dx*(6*g1 - 9*f1);

f = x*exp(3*x);

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

f1 = f2;
g1 = g2;
}
return 0;
}


こんな感じ