HOME

特集:MATLABを活かした技術開発

数値解析を見直してみよう!  その2:微分代数方程式

1 はじめに

ダイナミックシステムの問題を考えるとき、しばしば“微分代数方程式(Differential Algebraic Equation:以下DAE)”と呼ばれる問題に遭遇します。DAEは機構解析など、一部の分野を除いてはあまり意識されませんでしたが、ダイナミクスを表す微分方程式と、動きに拘束を与える代数方程式を同時に満足する方程式のことで、図1に示すように、機構系以外にも、流体系、電気系など多くのシステムに自然に存在します。


図1 実世界のDAEの例

図1スライダークランク機構では、各リンクは力を受けて運動(ダイナミクス)しますが、各リンクはそれぞれ独立に動くことはできず、ピン結合により互いの節点を共有しているので、位置が等しい(拘束)必要があります。サーボ弁では、シンダーの圧力差によりピストンが動き(ダイナミクス)ますが、シリンダに流入する流体と流出する流体の量(質量)は等しくなる(拘束)必要があります。電気回路でも同様なことが言えます。

このように多くのシステムは、某かの拘束の下に動いており、多くの場合はDAE問題となります。しかしこのDAE問題は、数値シミュレーションにおいて注意すべき側面を持っています。ここでは、普段はシミュレーションツールが自動的に計算しているDAE問題についてもう少し掘り下げるとともに、数値シミュレーション上の留意点を簡単なベンチマークで示し、理解を深めます。

2 DAEとは?

良く知られた多次元の常微分方程式(Ordinary Differential Equation: 以下ODE)は、(1)式のような陽的多次元形式で表現できます。

これは、Runge-Kuttaなどの陽的解法(現在の値から次のステップの値を逐次的に求める)で容易に答えを得ることができます。これに対しDAE問題は、一般的に(2)式のような陰的多次元形式で表現されます。

これは特別な場合、(3)式のように書き直すことができます。

これは、半陽的多次元形式または、制約のあるODEと呼ばれます

ここで、INDEX(指数)という概念を導入します。INDEXは、DAEとODEの間のメトリック(距離の概念:どれくらい離れているか)で、以下のように定義されます。

一般の多次元DAE(2)に対してINDEXは、y をyとtについて解くときに(すなわち、y に対する常微分方程式を定義るすために)必要な、最小の微分回数pのことである。

ここで、簡単な例によりINDEXの概念を具体的に示します。(5)式のDAE(常微分方程式と代数方程式の連立)を考えます。

y2´に対する陽的形式を得るために、(6)式のようにq(t)を2回微分する必要があるので、INDEXは2となります。

従って、通常のODEのINDEXは0となります。

3 INDEXリダクション


図2 振り子

INDEXが2以上のDAEを高INDEX DAEと呼びます。数値シミュレーションを用いる場合は、INDEXが0もしくは1まで低くすることが望まれます。これはODEソルバで高速に解いたり、あるいはDASSLに代表される多くのDAEソルバはINDEX 1に対応しており、これを利用して正確にシミュレーションできるからです。これについては次の節で紹介するとして、ここでは高INDEX DAEを低INDEX DAEに置き換える、「INDEX リダクション」というテクニックを紹介します。
図2に示す振り子を用いて、INDEXリダクションのしくみを説明します。


Lagrange乗数λを使って直交座標系でNewtonの方程式を立てると(7)式となります。

(7b)式を微分し、(7a)式を代入すると以下の式を得ます。

さらに微分(2回目)し、(7a)式を代入すると以下の式を得ます。

さらに微分(3回目)し、(7a)式を代入すると以下の式を得ます。

3回の微分後、拘束方程式(7b)を陽的微分方程式に変更できました。よってこのシステムはINDEX 3となります。これより(7)式は(8)式の陽的多次元ODEとなります。


図3 Maple/MapleSimによるDAE問題のシミュレーション手順

この例に示すように、INDEXリダクションは数式処理を伴います。INDEXリダクションは人手で行ってもできなくはありませんが、非常に煩雑な数式処理が伴うため、ツールによる自動化が望まれます。その代表がMaple/MapleSimで、これにより、数式処理/数値解析を融合したシミュレーションを簡単に実施できます。図3にMaple/MapleSimを用いたDAE問題のシミュレーション手順を示します。MapleSimでは高INDEX DAE問題を自動的にINDEX 1までINDEXリダクションし、高精度にシミュレーションします。なお、MapleSimの詳細についてはここでは触れませんので、以下のURLを参照下さい。
http://www.cybernet.co.jp/ad/maplesim/


4 ベンチマーク:DAE vs ODE

多くの場合、図2に示す振り子は(9)式のような極座標系で回転の運動方程式を立てて、ODE問題として解かれます。

しかしこの場合、リンクの長さによる拘束を陽に解いていない(回転系の積分演算だけで解く)ため、シミュレーション時間が長くなると積分誤差が累積し、リンクの長さが変わったように計算されるなど、物理的に矛盾のある結果を導く可能性があります。


図4 ODEとDAEのシミュレーション結果

(9)式をSimulinkでODE問題としてモデル化した場合と、図2をMapleSimでモデル化し、INDEX 1 DAE問題としてモデル化した場合のシミュレーション結果を図4に示します。どちらも振り子の初期角度を水平位置とした場合の自由振動です。a),b)は0〜5秒までのx位置とy位置を示しており、ODEでもDAEでも大きな差は認められません。しかし、c), d)は0〜105秒と長時間シミュレーションした場合のy位置ですが、ODEでは振り子の角度が小さくなっていくのが確認されます。

このように、ODEソルバでは拘束が陽に計算されないので、累積積分誤差に注意を払う必要があります。


関連製品・サービス

Mapleシリーズ

数式処理・数式モデル設計環境

お問い合わせは

アドバンスドソリューション統括部 モデルベース開発推進室

E-mail:infomaple@cybernet.co.jp

2/2

(5/14)