微分代数方程式の解析と操作

この例では、Symbolic Math Toolbox™ を使用して微分指数の高い微分代数方程式 (DAE) を求解する方法を説明します。

エンジニアは多くの場合、物体 (機械システム、電子機器など) の振る舞いを微分方程式と代数方程式の組み合わせで指定します。MATLAB® では、このような DAE を積分することのできる ode15i や ode15s などの特別な数値ソルバーが提供されています (それらの '微分指数' が 1 を上回らない場合に限る)。

この例では、代数制約をもつ微分方程式系としてモデルをセットアップするところから、数値シミュレーションを行うまでのワークフローを説明します。以下の Symbolic Math Toolbox 関数を使用します。

daeFunction
findDecoupledBlocks
incidenceMatrix
isOfLowDAEIndex
reduceDifferentialOrder
massMatrixForm
reduceDAEIndex
reduceDAEToODE
reduceRedundancies
sym/decic

モデルのパラメーターの定義

2 次元の実体振子について考えます。この振子は一定の長さ r の紐で原点に取り付けられた重り m で構成されています。重りに作用するのは重力加速度 g = 9.81 m/s^2 のみです。このモデルは、重りの位置 (x(t), y(t)) と重りを円上に留めている紐の力 F(t) から成る 2 階微分方程式で構成されます。力は紐に沿って掛かっています。

syms x(t) y(t) F(t) m g r
eqs = [m*diff(x(t), t, t) == F(t)/r*x(t);
       m*diff(y(t), t, t) == F(t)/r*y(t) - m*g;
       x(t)^2 + y(t)^2 == r^2]

eqs = 
 $(\begin{array}{c} m \frac{\partial^{2}}{\partial t^{2}} x (t) = \frac{F (t) x (t)}{r} \\ m \frac{\partial^{2}}{\partial t^{2}} y (t) = \frac{F (t) y (t)}{r} - g m \\ {x (t)}^{2} + {y (t)}^{2} = r^{2} \end{array})$

vars = [x(t), y(t), F(t)]

vars =  $(\begin{array}{c} x (t) & y (t) & F (t) \end{array})$

この DAE 系を 1 階微分代数方程式系に書き換えます。

[eqs, vars, newVars] = reduceDifferentialOrder(eqs, vars)

eqs = 
 $(\begin{array}{c} m \frac{\partial}{\partial t} Dxt (t) - \frac{F (t) x (t)}{r} \\ g m + m \frac{\partial}{\partial t} Dyt (t) - \frac{F (t) y (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ Dxt (t) - \frac{\partial}{\partial t} x (t) \\ Dyt (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \end{array})$

newVars = 
 $(\begin{array}{c} Dxt (t) & \frac{\partial}{\partial t} x (t) \\ Dyt (t) & \frac{\partial}{\partial t} y (t) \end{array})$

高指数 DAE 系を解く

ode15i などの MATLAB の数値ソルバーを使用可能にするため、以下の手順を行う必要があります。

DAE 系を MATLAB 関数ハンドルに変換する。
系のシンボリックパラメーターの数値を選択する。
矛盾のない初期条件を設定する。

DAE 系を MATLAB 関数ハンドルに変換するには、daeFunction を使用します。

F = daeFunction(eqs, vars, [m, g, r])

F = function_handle with value:
    @(t,in2,in3,in4)[in3(4,:).*in4(:,1)-(in2(3,:).*in2(1,:))./in4(:,3);in3(5,:).*in4(:,1)+in4(:,1).*in4(:,2)-(in2(3,:).*in2(2,:))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:)-in3(1,:);in2(5,:)-in3(2,:)]

系のシンボリックパラメーターに次の数値 m = 1kg、g = 9.18m/s^2 および r = 1m を代入します。

f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])

f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

関数ハンドル f は数値ソルバー ode15i に適した入力です。次の手順では矛盾のない初期条件を計算します。odeset を使用して数値の許容誤差を設定します。それから、MATLAB 関数 decic を使用して、時間 t0 = 0 における位置と導関数に整合する初期条件 y0, yp0 を計算します。

opt = odeset('RelTol', 10.0^(-4), 'AbsTol' , 10.0^(-4));
t0 = 0;
[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(3,1)], [], zeros(5,1), [], opt)

初期条件をテストします。

f(t0, y0, yp0)

これで、ode15i を使用して系の求解が実施できます。ode15i を呼び出すと、積分は即座に停止し、次の警告を発します。

Warning: Matrix is singular, close to singular or badly scaled.

Results may be inaccurate. RCOND = NaN.

Warning: Failure at t=0.000000e+00.

Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (0.000000e+00) at time t.

この例では、ode15i はこれらの警告を複数回発します。見やすくするため、ode15i を呼び出す前に warning('off','all') を使用して警告を無効にし、呼び出した後に再び有効にします。

tfinal = 0.5;
s = warning('off','all');
ode15i(f, [t0, tfinal], y0, yp0, opt);

Figure contains an axes object. The axes object contains 10 objects of type line.

warning(s)

DAE 系の解析と調整

DAE 系の微分指数を確認します。

isLowIndexDAE(eqs, vars)

ans = logical
   0

この結果は ode15i でこの系が解けない理由を示しています。この関数では入力の DAE 系の微分指数が 0 または 1 である必要があります。いくつかの隠れた代数制約を含む同等のより大きな DAE 系にモデルを拡張することにより微分指数を簡約します。

[eqs, vars, newVars, index] = reduceDAEIndex(eqs, vars)

eqs = 
 $(\begin{array}{c} m Dxtt (t) - \frac{F (t) x (t)}{r} \\ g m + m Dytt (t) - \frac{F (t) y (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ Dxt (t) - {Dxt}_{1} (t) \\ Dyt (t) - {Dyt}_{1} (t) \\ 2 {Dxt}_{1} (t) x (t) + 2 {Dyt}_{1} (t) y (t) \\ 2 y (t) \frac{\partial}{\partial t} {Dyt}_{1} (t) + 2 {{Dxt}_{1} (t)}^{2} + 2 {{Dyt}_{1} (t)}^{2} + 2 Dxt1t (t) x (t) \\ Dxtt (t) - Dxt1t (t) \\ Dytt (t) - \frac{\partial}{\partial t} {Dyt}_{1} (t) \\ {Dyt}_{1} (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \\ {Dxt}_{1} (t) \\ {Dyt}_{1} (t) \\ Dxt1t (t) \end{array})$

newVars = 
 $(\begin{array}{c} Dytt (t) & \frac{\partial}{\partial t} Dyt (t) \\ Dxtt (t) & \frac{\partial}{\partial t} Dxt (t) \\ {Dxt}_{1} (t) & \frac{\partial}{\partial t} x (t) \\ {Dyt}_{1} (t) & \frac{\partial}{\partial t} y (t) \\ Dxt1t (t) & \frac{\partial^{2}}{\partial t^{2}} x (t) \end{array})$

index = 3

4 番目の出力からは元のモデルの微分指数が 3 であることが分かります。新しい系を単純化します。

[eqs, vars, S] = reduceRedundancies(eqs, vars)

eqs = 
 $(\begin{array}{c} - \frac{F (t) x (t) - m r Dxtt (t)}{r} \\ \frac{g m r - F (t) y (t) + m r Dytt (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ 2 Dxt (t) x (t) + 2 Dyt (t) y (t) \\ 2 {Dxt (t)}^{2} + 2 {Dyt (t)}^{2} + 2 Dxtt (t) x (t) + 2 Dytt (t) y (t) \\ Dytt (t) - \frac{\partial}{\partial t} Dyt (t) \\ Dyt (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \end{array})$

S = struct with fields:
      solvedEquations: [3x1 sym]
    constantVariables: [0x2 sym]
    replacedVariables: [3x2 sym]
       otherEquations: [0x1 sym]

新しい系の微分指数が低い (0 または 1) ことを確認します。

isLowIndexDAE(eqs, vars)

ans = logical
   1

低指数 DAE 系を解く

シンボリックパラメーターを数値で置き換える MATLAB 関数ハンドルを生成します。

F = daeFunction(eqs, vars, [m, g, r])

F = function_handle with value:
    @(t,in2,in3,in4)[-(in2(3,:).*in2(1,:)-in2(7,:).*in4(:,1).*in4(:,3))./in4(:,3);(-in2(3,:).*in2(2,:)+in2(6,:).*in4(:,1).*in4(:,3)+in4(:,1).*in4(:,2).*in4(:,3))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:).*in2(1,:).*2.0+in2(5,:).*in2(2,:).*2.0;in2(7,:).*in2(1,:).*2.0+in2(6,:).*in2(2,:).*2.0+in2(4,:).^2.*2.0+in2(5,:).^2.*2.0;in2(6,:)-in3(5,:);in2(5,:)-in3(2,:)]

f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])

f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

MATLAB 関数 decic によって簡約された指数に整合する初期条件を計算します。ここで、opt は数値の許容誤差を設定するオプション構造体です。これは odeset を使用して計算済みです。

[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(5,1)], [], zeros(7,1), [], opt)

系を解いて解をプロットします。

ode15i(f, [t0, tfinal], y0, yp0, opt)

Figure contains an axes object. The axes object contains 14 objects of type line.