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

はじめに

GNU Octaveの常微分方程式のソルバとしては,標準で組み込まれているlsodeや,Octave Forgeで提供されるode45などのOdePkgがあります. これらのソルバでは,変数xtを引数として導関数x′を返すという,解くべき常微分方程式を表す関数を自分で定義して,各ソルバの引数として与える必要があります. その際,OdePkgのソルバ群では常微分方程式を表す関数に時間tと状態変数x以外のパラメータをセル配列で渡すことができますが,lsodeではパラメータを渡す方法がありません.

一方,外力が時間tの関数として与えられる場合には,常微分方程式を表す関数の中に外力の関数を記述することができますが,配列データとして外力が与えられた場合には,この外力もパラメータとして常微分方程式を表す関数に渡す必要があります.

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

前提条件

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

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

また,以下のプログラムではN>65536となる場合を考慮して,四捨五入する関数にuint32を使っています.

ソルバにlsodeを用いる場合

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

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

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

      function xdot = lsode_target_equation(x, t)
        global dt Excitation
        idx = uint32(t / dt) + 1;
        xdot(1) = x(2);
        xdot(2) = - x(1) + Excitation(idx);
      endfunction
    

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

ソルバにOdePkgを用いる場合

octave-forgeで提供されるOdePkgode45を用いた場合の手順は,解くべき常微分方程式を与える関数をtarget_equationとした場合,以下となります. (...の部分には,左辺の変数に与える値が入ります.)

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

t1ode45の計算に用いた時間ステップが入るベクトルであり,上記のスクリプトではtと等しくなります.

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

      function xdot = target_equation(t, x)
        global dt Excitation
        idx = uint32(t / dt) + 1;
        xdot(1) = x(2);
        xdot(2) = - x(1) + Excitation(idx);
      endfunction
    

なお,このOdePkgは基本的に自動時間ステップなので,上記のように時間ステップを自分で与えた場合,特定のNでは以下のようなエラーが出ることがあります. (以下は0.05の時間ステップでN=10000とした場合に出たエラーメッセージです.)

error: Solving has not been successful. The iterative integration loop exited at time t = 499.950000 before endpoint at tend = 499.950000 was reached. This may happen if the stepsize grows smaller than defined in vminstepsize. Try to reduce the value of "InitialStep" and/or "MaxStep" with the command "odeset".

このような場合はNを増減させてエラーが出ない値を探してください.

ただし,OdePkgのChangeLogを見ると以下のように書いてあるので,

      2008-10-07 21:08  treichl

        * inst/ode23.m, inst/ode45.m, inst/ode54.m, inst/ode78.m: Fixed
          stepsize bug if fixed time steps are given.
    

特定のNで計算できない上記の現象は,2008年10月以降は解決されているかもしれません. Octave 3.2.3とOdePkg 0.6.7の組み合わせで試したところ再現しませんでした.

固定時間刻みのソルバ

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

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

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

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

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

      x = crk(@target_equation, t, x0);
      x = rkf(@target_equation, t, x0);
    

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

      x = ode_solver_wrapper("crk", @target_equation, t, x0);
      x = ode_solver_wrapper("rkf", @target_equation, t, x0);
      x = ode_solver_wrapper("lsode", @target_equation, t, x0);
      x = ode_solver_wrapper("ode45", @target_equation, t, x0);
    

なお,このode_solver_wrapperではOdePkgで提供される他のソルバのode23, ode54, ode78も同様の手続きで呼び出すことができます.

Last update: 2010.7.21

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