pypy.com/

Python、Unity、FX自動化などを勉強しています。あと、コーラと車も好きです。そこらへんについて、たまに記事を書きます。

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

こんにちは。

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

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

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


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

オイラー(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}

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

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;
}

これを用いて、計算を行ってみると、以下のようになります。

今回は、単振動のはずなのに、振幅が大きくなっていっていることが確認できます。

なので、もう少し刻み時間を小さくして、dt=0.01とすると、以下のようになります。

ましになったのが分かりますね。

このように、オイラー法では、刻み時間が小さいほど、精度がよくなります

今回はこれで終わります。