c言語で二階常微分方程式を解く（オイラー法）

こんにちは。

今回は、c言語でオイラー法を使った二階常微分方程式を解く方法について、説明したいと思います。

オイラー法は以前に下記記事で説明したように、簡単に一階常微分方程式に適用できます。

c言語で微分方程式を解く（オイラー法） - pypy.com/

しかし、二階の常微分方程式を解く際には、もう少し段階を踏む必要があるので、それを説明していきます。

オイラー（Euler）法の考え方
オイラー（Euler）法のコード

オイラー（Euler）法の考え方

今回は、以下のような運動方程式を考えます。

今回は運動方程式を題材にしますが、ほかの問題でも同様にできます。

$m\frac{d^2x(t)}{dt^2}+kx(t)=0\tag{1}$

このままでは、オイラー法を適用することができないので、一階微分方程式二つに分解します。

(1)は $v(t)=\frac{dx(t)}{dt}$ とおいて、以下の二式に分解できます。

$v(t)=\frac{dx(t)}{dt}\tag{2}$
$m\frac{dv(t)}{dt}+kx(t)=0\tag{3}$

オイラー法を適用する前に、少し整理して、

$\frac{dx(t)}{dt}=v(t) \tag{4}$
$\frac{dv(t)}{dt}=-\frac{k}{m}x(t) \tag{5}$

この２式にオイラー法を適用すると、以下のようになります。

$x(t+dt)=x(t)+v(t)dt\tag{6}$
$v(t+dt)=v(t)-\frac{k}{m}x(t)dt\tag{7}$

オイラー（Euler）法のコード

それでは、コードを書いていきます。

今回は、パラメータと初期値を以下のように設定します。

$m=1.0$
$k=1.0$
$x(0)=1.0$
$v(0)=0.0$

c言語で考え方で書いたもののように、書くと以下のようになります。

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

double f_x(double x_pre, double v_pre, double dt);
double f_v(double x_pre, double v_pre, double k, double m, double dt);

int main(void){
  // パラメータ
  double m = 1.0;
  double k = 1.0;
  // 初期値
  double x = 1.0;
  double v = 0.0;
  // 解析範囲
  double t_start = 0.0;
  double t_end = 50.0;
  // 刻み時間
  double dt = 0.1;
  // 分割数
  int n = (int)((t_end-t_start) / dt) + 1;

  printf("t, x, v\n");
  for(double i=0; i<n; i++){
    double x_new = f_x(x, v, dt);
    double v_new = f_v(x, v, k, m, dt);
    x = x_new;
    v = v_new;
    printf("%5.5lf, %5.5lf, %5.5lf\n", t_start+dt*i, x, v);
  }
}

double f_x(double x_pre, double v_pre, double dt){
  return x_pre + v_pre * dt;
}

double f_v(double x_pre, double v_pre, double k, double m, double dt){
  return v_pre - k/m * x_pre * dt;
}