GNU Octave を用いた,離散時間の外力に対する常微分方程式の数値解析

はじめに

GNU Octave の常微分方程式のソルバーとしては,標準で組み込まれている lsodeode45 があります. これらのソルバーでは,変数 xvt を引数として導関数 xv′ を返すという,解くべき常微分方程式を表す関数を自分で定義して,各ソルバーの引数として与える必要があります.

なお,解くべき方程式が運動方程式の場合には変位の2階微分である加速度の微分方程式となりますが,常微分方程式のソルバーは1階微分方程式しか解けないので,解くべき変数に速度 v を,そして解くべき微分方程式に dx/dt = v を加えて,2本の1階微分方程式に変換する必要があります. また,外力が時間 t の関数として与えられる場合には,常微分方程式を表す関数の中に外力の関数を記述することができますが,配列データとして外力が与えられた場合には,この外力もパラメータとして常微分方程式を表す関数に渡す必要があります.

一般的なプログラムではグローバル変数は敬遠されることが多いですが,このような常微分方程式を表す関数にパラメータを渡す場合にはグローバル変数が一番簡便な方法であると思います. 以上より,グローバル変数を用いて任意外力を常微分方程式の関数に渡す方法を記します.

前提条件

考えている任意外力はサンプリング周期 Δt で離散化されたデータであるはずなので,任意外力だけでなくサンプリング周期 Δt もグローバル変数として渡す必要があります. また,外力と応答間で何らかの計算を行う場合には,応答 xv も同じサンプリング周期 Δt である必要があります. これらより,以下の条件で解くものとします.

なお, global によるグローバル変数の宣言は,常微分方程式のソルバーを呼び出すプログラム内で対象の変数を初めて用いる時に行う必要があり,常微分方程式を表す関数内でも宣言する必要があります. global 宣言では名前で変数の照合を行うため,これら2つの global 宣言で変数の順序が一致していなくても構わないのですが,あえて変える必要もないので,特に理由がない限り同じ順序にしておく方が良いでしょう.

ソルバーに lsode を用いる場合

GNU Octave に標準で組み込まれている lsode を用いた場合の手順は,解くべき常微分方程式を与える関数を lsode_target_equation とした場合,以下となります. (...の部分には,左辺の変数に与える値が入ります.)

      global dt Excitation
      dt = ...
      Excitation = ...
      N = ...
      xv0 = ...
      t = transpose(dt * [0 : N - 1]);
      xv = lsode(@lsode_target_equation, xv0, t);
    

解くべき方程式が無次元化された減衰のない1自由度振動系の場合,lsode_target_equation.m は以下となります. なお,区分求積法のように dt の間で Excitation を一定値として取り扱うと後述する crkrkf で結果に差が出ることがありますが,線形加速度法のように区間 dt の間で線形補間すると crkrkf の差が非常に小さくなります. そのため,以下の関数では dt の間で Excitation を一次関数で補間しています. 以降のソルバーでも同様です.

      function xvdot = lsode_target_equation(xv, t)
        global dt Excitation
        idx = floor(t / dt) + 1;
        if (idx < numel(Excitation))
          F = (idx - t / dt) * Excitation(idx) + (1.0 + t / dt - idx) * Excitation(idx + 1);
        else
          F = Excitation(idx);
        endif
        xvdot(1) = xv(2);
        xvdot(2) = - xv(1) + F;
      endfunction
    

ただし,この lsode における計算精度の自動調整機能は,配列データとして与えられる離散時間の外力と相性が悪いようで,他のソルバーに比べると計算時間が非常に長くなります.

ソルバーに ode45 を用いる場合

GNU Octave に標準で組み込まれてるもう一つのソルバーの ode45 を用いた場合の手順は,解くべき常微分方程式を与える関数を target_equation とした場合,以下となります. (...の部分には,左辺の変数に与える値が入ります.)

      global dt Excitation
      dt = ...
      Excitation = ...
      N = ...
      xv0 = ...
      t = transpose(dt * [0 : N - 1]);
      [t1, xv] = ode45(@target_equation, t, xv0);
    

t1ode45 の計算に用いた時間ステップが入るベクトルであり,上記のスクリプトでは t と等しくなります. 通常の場合,ベクトル t1 は不要なので最後の行を以下とすることによって,余分な変数 t1 をなくすことができます.

      [~, xv] = ode45(@target_equation, t, xv0);
    

解くべき方程式が無次元化された減衰のない1自由度振動系の場合,target_equation.m は以下となります.

      function xvdot = target_equation(t, xv)
        global dt Excitation
        idx = floor(t / dt) + 1;
        if (idx < numel(Excitation))
          F = (idx - t / dt) * Excitation(idx) + (1.0 + t / dt - idx) * Excitation(idx + 1);
        else
          F = Excitation(idx);
        endif
        xvdot(1) = xv(2);
        xvdot(2) = - xv(1) + F;
      endfunction
    

固定時間刻みのソルバー

以上のように GNU Octave で提供される常微分方程式のソルバーを,離散化されたデータとしての外力を受ける運動方程式に適用する場合には,以下の注意点があります.

  1. lsode : 計算時間が長い
  2. ode45 : 計算に用いた時間ステップが最初の返り値となる (離散時間の外力に対する常微分方程式の場合,自分で指定したベクトルと同じものである)

これらの問題は上記の関数が時間刻みを動的に変化させるアルゴリズムを用いて実装されているためです. そのため,固定時間刻みのシンプルなソルバーとして,以下の方法の関数を作ってみました.

  1. 古典的なRunge Kutta法 (crk)
  2. Runge Kutta Fehlberg法 (rkf)

それぞれ,常微分方程式が上記のように定義された target_equation(t, xv) の順番で引数を取る関数として与えられるとき,以下のようにして使います.

      xv = crk(@target_equation, t, xv0);
      xv = rkf(@target_equation, t, xv0);
    

また,筆者が作った2つの関数 (crk, rkf),および lsodeode45 では,方程式を記述する関数の書き方やソルバーの使い方に違いがあるので,簡単に乗り換えることができません. そのため,同じ手続きでこれらのソルバーを呼び出すための関数(ode_solver_wrapper)も作ってみました. この関数を使うと,常微分方程式が上記のように定義された target_equation(t, x) の順番で引数を取る関数として与えられるとき,crk, rkf, lsode, ode45 のそれぞれのソルバーを以下の手続きで使うことができます.

      xv = ode_solver_wrapper("crk", @target_equation, t, xv0);
      xv = ode_solver_wrapper("rkf", @target_equation, t, xv0);
      xv = ode_solver_wrapper("lsode", @target_equation, t, xv0);
      xv = ode_solver_wrapper("ode45", @target_equation, t, xv0);
    

なお,この ode_solver_wrapper では他のソルバーの ode23 も同様の手続きで呼び出すことができます.

グローバル変数を使わない場合

以上ではグローバル変数を使う場合について書きましたが,グローバル変数を使わない場合は,以下のように常微分方程式を与える関数に引数として変数を渡さなければなりません.

      function xvdot = target_equation_no_global(t, xv, dt, Excitation)
        idx = floor(t / dt) + 1;
        if (idx < numel(Excitation))
          F = (idx - t / dt) * Excitation(idx) + (1.0 + t / dt - idx) * Excitation(idx + 1);
        else
          F = Excitation(idx);
        endif
        xvdot(1) = xv(2);
        xvdot(2) = - xv(1) + F;
      endfunction
    

そして,微分方程式のソルバーでは以下のように無名関数を用いて引数を変換する必要があります.

      xv = lsode(@(xv, t) target_equation_no_global(t, xv, dt, Excitation), xv0, t);
      [~, xv] = ode45(@(t, xv) target_equation_no_global(t, xv, dt, Excitation), t, xv0);
      xv = crk(@(t, xv) target_equation_no_global(t, xv, dt, Excitation), t, xv0);
      xv = rkf(@(t, xv) target_equation_no_global(t, xv, dt, Excitation), t, xv0);
    

このように大きな引数 Excitation を無名関数のデフォルト引数として渡すと,グローバル変数を使った場合に比べて計算時間が若干増えるようです.

Last update: 2024.1.25

[機械力学研究室ホームページに戻る]   [数値計算に関する情報のページに戻る]