HOME

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

数値解析を見直してみよう!  その1:スティフ

はじめに

私たちはシミュレーションにおいて、様々な場面で数値解析の恩恵を受けています。長年の研究成果と技術競争により、商用化されたツールは普段はその計算の仕組みを意識しなくても、“かなり” 精度の高い答えを出してくれます。また数値計算の利点として、入力さえ間違えなければ「基本的には何がしかの答え」を出してくれます。しかしながら、それ故「本当に合っているのか?」という視点を持つことが重要なのも事実です。数値計算の中身を一つ一つ追うのは非常にタフな作業で非現実的ですが、その性質を理解するだけでも数値解析をより身近に感じ、これまで以上にうまく付き合うことができると思います。

ここではそのような観点から中身が陽に分かっている問題に対して、数値解析の性質を2回に分けて見直してみようと思います。話を具体的にするため、ユーザが自由に問題設定できる数値解析ツールの代表としてMATLABを数式処理(厳密解)ツールとしてMapleを用いますが、両者のツールの優劣を付けるのではなく数値解析の性質を再考し、より有効に使いこなすことが目的であることを、改めて付記させていただきます。

スティフな系


定義1 陽解法と陰解法


定義2 スティフな系

微分方程式を数値解析で解く場合、いわゆる“スティフな系”と呼ばれる問題に出くわすことがあります。解の安定性やシミュレーション速度に注意が必要と言われていますが、この性質を明らかにします。

スティフな系とは、微分方程式を陽解法で解く場合に重要な概念となります。定義1に陽解法の定義を、定義2にスティフな系の定義を記します。

陽解法は現在と過去の結果から未来の挙動を計算できるため、時間軸シミュレータでは都合の良い方法です。

反面、たとえ系が安定であってもシミュレーションの安定性を保証できず、精度上必要な計算刻みより遙かに小さい刻みが必要とされる場合があり、これを“スティフな系”と呼ぶというのが定義1、2の言っているところです。

スティフの概念は非線形システムにも適用できますが、 線形システムでは系の固有値で決定され、固有値の原点からの距離が大きく離れている場合に相当します。

物理的には:

  • 系の固有振動数に大きな差がある場合
  • 系の減衰に大きな差がある場合

具体的な問題としては:

  • 機械系と電気系など時定数が大きく異なる系が混在するモータのようなメカトロニクス
  • 金属とゴムなどの剛性の大きく異なる複合材料

などが挙げられます。


テストケース

1.テストモデル

図1のような固有振動数を任意に変えることができる2自由度振動系を考えます。原点からの距離の差は約660倍となってます。


図1 スティフな系のテストモデル

2.数値解析:固定刻みソルバ


図2 2自由度振動系のSimulinkモデル

図1のモデルをSimulinkでモデル化すると、図2のようになります。各質点に適当な初期変位を与え、代表的な陽解法であるRunge-Kutta系列のソルバDormand-Prince法で、刻み幅変えてシミュレーションした場合の質点1、2の変位を図3に示します。この系は本来物理的に安定であることは自明ですが、刻み幅が大きくなると解が発散することが分かります。



図3 Dormand-Prince法(陽解法)による固定刻みシミュレーション

3.数値解析:可変刻みソルバ

刻みを調整しながらシミュレーションする可変刻みソルバには、スティフな系に対応したものがあります。これは、スティフな系では一般にシミュレーションが極端に遅くなることを改善したものです。しかし、安易に用いて良いものか?について考察します。

図4のc )、d ) はスティフな系に対応したソルバです。ソルバにより大きく結果が異なることが分かります。後に説明しますが、この場合ode45が一番精度良い結果を示しています。この結果から分かるように、スティフなソルバあるいは次数の低いソルバでは、高周波の減衰が高めに見積もられた結果となります。多くの物理系、工学問題では「固有振動数が高いほど減衰しやすい」ことが多いのでスティフなソルバは有効に機能しますが、本系のように高次の減衰が小さい場合は注意を要することが分かります。


図4 可変刻みシミュレーション(スティフソルバの比較)

4.数式処理:厳密解

数式処理ツールMapleと、その上で動作する次世代物理モデラMapleSimを用いて、本問題に対する真の解を探ります。Maple/MapleSimによるシミュレーションの流れを図5に示します。


図5 Maple/MapleSimによるシミュレーション

図6にMapleSimで作成した物理モデルを、図7にそこから自動生成したシステム方程式を示します。これは線形システムであるため、Mapleの数式処理機能を用いて厳密解を求めることができます。


図6 MapleSimによる物理モデリング

図7 自動生成したシステム方程式

さらにMapleでは厳密解を数値化する際、任意の桁数で計算できます。ここではハードウェアの倍精度浮動小数点精度(10進換算で約15桁)を超える20桁で数値化した結果とSimulinkの数値解析結果を比較したものを図8に示します。


図8 厳密解と数値解の比較

Mapleの結果は厳密解を高精度に数値化しており、ほぼ真の解と考えることができるため、ここではこれを厳密解とします。これよりSimulinkのode45はわずかに高周波数で減衰するものの、他のソルバと比較(図4)し大変良好な結果を示していることが分かります。

考 察

以上の結果より以下のことが言えます。

1)スティフな系において固定刻みソルバは解の安定性を左右するため、刻み幅を慎重に選択する必要があります。

2)可変刻みソルバではトレランスにより刻みはコントロールできるので固定刻みより有効と考えられるが、ソルバにより解が異なるので注意が必要です。

3)スティフなソルバは高周波の減衰が高めに見積もられるので、減衰が少ない系に対しては適しません。

4)複雑な系に数値解析を適用するとソルバの性質がわかりにくくなってしまうので、モデルが陽に分かっているシンプルな系に対して厳密解と比較することによりソルバの性質を知ることは、数値解析をより有効に利用するために重要であると考えられます。

関連製品・サービス

Mapleシリーズ

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

お問い合わせは

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

E-mail:infomaple@cybernet.co.jp

1/2

(4/14)