Compose-4010：常微分方程式を解く

Tutorial Level: Intermediate

このチュートリアルでは、ComposeのOML言語で常微分方程式（ODE）を解くことを検討します。わかりやすいように、広く使用されていることから2次微分方程式を使用します。この方程式で、機械的なばね質量減衰系や直列RLC回路などを記述できます。どちらにも、振幅をAとする任意の正弦波入力を印加するものとします。

状態変数をΦとして、ばね質量減衰系では各変数を次のように定義します。
  • Φ：位置
  • α：質量で除算した減衰係数
  • β：質量で除算したばね剛性
  • 励振は力です。
直列RLC回路系では各変数を次のように定義します。
  • Φ：電流
  • α：インダクタンスで除算した抵抗値
  • β：静電容量とインダクタンスの積の逆数
  • 励振は、A*sin(ωt)を導関数とする電圧です。

1次方程式への変換

ODEソルバーでは1次の陽的方程式のみを解くことができるので、ここで説明する手順は必須です。この変換を実現するには、元の状態変数Фの導関数に等しい状態変数ξを新たに導入します。次に示す2つの連立1次微分方程式を考えることができます。

新しい状態変数の関係方程式：

新しい状態変数の側から記述した元の方程式：

これらの方程式は陽的な形態で表現されています（左辺には導関数のみが記述されています）。ここでは、これら2つの連立方程式を解いて、2つの状態変数を時間の関数として求める必要があります。

わかりやすくするために、すべてのパラメーター（α、β、A、ω）は1に等しいものとします。

2つの方程式の実行

  1. Composeを開始します。
  2. ファイルメニューで、開くを選択して、<installation_dir>/tutorials/フォルダでファイルode45Demo.omlを探します。

    ここで取り上げている系の関数は次のようになります。このコードでは、ベクトル[φ, ξ]をxと記述しています。

    この系の関数では、上記で得られた2つの方程式を実行します。1番目の引数は現在の時間、2番目の引数は、ベクトルで記述した状態変数値です。

ODEソルバーの使用

次は、OMLの組み込みODEソルバーを使用します。このチュートリアルでは、次のように入力を設定してode45関数を使用します。

1番目の引数で、ODEによる系の関数を指定します。2番目の引数は、この方程式で得る解の時刻暦を記述したベクトルです。3番目の引数は、状態変数の初期条件ベクトルです。わかりやすいようにoptions引数でデフォルトの閾値を明示的に指定します。

timeベクトルは、1番目の出力引数として再生成されます。このチュートリアルでは入力と同じです。2番目の出力引数であるxには、ベクトルtに記述された時間ごとに、φとξの値が列単位で記述されます。

結果のプロット

次のコードでは、この結果を抽出してプロットします。

ばね質量減衰系の場合、このプロットは位置と速度を表しています。RLC回路の場合、このプロットは電流とその変化速度を表しています。