Posted at

[C言語] 2階常微分方程式の解をグラフにする(ルンゲクッタ法)

More than 3 years have passed since last update.


サンプル

微分方程式

\frac{d^2y}{dx^2} = -ay (初期条件:x = 0 で y = 0, \frac{dy}{dx} = a, aは正の定数)

の解をルンゲクッタ法で求め、aの値を変えた時の解をグラフにする。


sample.c

#include <stdio.h>

#include <math.h>

int main(void) {

FILE * data;

double f, f1, f2, g1, g2, dx, x, k1, k2, k3, k4, k, h1, h2, h3, h4, h;
int a = 2;

data = fopen("sample_data.csv", "w");

f1 = 0;
g1 = a;
dx = 0.01;

for (int i=0; i<=500; i++) {

x = i*dx;

k1 = dx*g1;
h1 = dx*(-a*f1);

k2 = dx*(g1 + h1/2.);
h2 = dx*( -a*(f1 + k1/2.) );

k3 = dx*(g1 + h2/2.);
h3 = dx*( -a*(f1 + k2/2.) );

k4 = dx*(g1 + h3/2.);
h4 = dx*(-a*(f1 + k3) );

k = (k1 + 2*k2 + 2*k3 + k4)/6.;
f2 = f1 + k;

h = (h1 + 2*h2 + 2*h3 + h4)/6.;
g2 = g1 + h;

f = sqrt(a)*sin( sqrt(a)*x );

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

f1 = f2;
g1 = g2;
}

fclose(data);

return 0;
}



実行結果


問題

微分方程式

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

の解をルンゲクッタ法で求め、xとyの関係およびyと$\frac{dy}{dx}$の関係をグラフにする。aとbは正の定数でこれらの値を変えた時の変化を調べる。


problem.c

#include <stdio.h>

#include <math.h>

int main(void) {

FILE * data;

double f, f1, f2, g1, g2, dx, x, k1, k2, k3, k4, k, h1, h2, h3, h4, h;
int a, b;

data = fopen("problem_data.csv", "w");

printf("aを入力して下さい。 a = ");
scanf("%d", &a);
printf("bを入力して下さい。 b = ");
scanf("%d", &b);
printf("yを入力して下さい。 y = ");
scanf("%lf", &f1);
printf("dy/dxを入力して下さい。 dy/dx = ");
scanf("%lf", &g1);
printf("dxを入力して下さい。 dx = ");
scanf("%lf", &dx);

for (int i=0; i<=500; i++) {

x = i*dx;

k1 = dx*g1;
h1 = dx*( (-b)*g1 - a*f1);

k2 = dx*(g1 + h1/2.);
h2 = dx*( (-b)*(g1 + h1/2.) - a*(f1 + k1/2.) );

k3 = dx*(g1 + h2/2.);
h3 = dx*( (-b)*(g1 + h2/2.) - a*(f1 + k2/2.) );

k4 = dx*(g1 + h3);
h4 = dx*( (-b)*(g1 + h3) - a*(f1 + k3) );

k = (k1 + 2*k2 + 2*k3 + k4)/6.;
f2 = f1 + k;

h = (h1 + 2*h2 + 2*h3 + h4)/6.;
g2 = g1 + h;

f = 5*exp(3*x)*(1 - 3*x);

printf("x = %f, 数値解-> f1 = %f, dy/dx = %f, 正確な値-> f = %f \n", x, f1, g1, f);
fprintf(data,"%f, %f, %f\n", x, f1, g1);

f1 = f2;
g1 = g2;
}

fclose(data);

return 0;
}


手入力するようにした。


実行結果

$a = 2, b = 0, y = 5, \frac{dy}{dx} = 0, dx = 0.01$

での結果。