JP5692739B2

JP5692739B2 - 常微分方程式を解くための方法、プログラム及びシステム

Info

Publication number: JP5692739B2
Application number: JP2010054251A
Authority: JP
Inventors: アルキメデス・カネド; 吉澤　武朗; 武朗吉澤; 小松　秀昭; 秀昭小松
Original assignee: International Business Machines Corp
Current assignee: International Business Machines Corp
Priority date: 2010-03-11
Filing date: 2010-03-11
Publication date: 2015-04-01
Anticipated expiration: 2030-03-11
Also published as: US20110225225A1; JP2011186991A

Description

この発明は、コンピュータを利用したシミュレーション・システム等において使用される常微分方程式を解くための方法、プログラム及びシステムに関し、特に、常微分方程式を解くための計算を高速化し、あるいは計算量を減らすための技法に関するものである。

従来より、科学技術計算、シミュレーションなどの分野で、コンピュータが利用されている。

最近になって特に盛んに開発されるようになってきたシミュレーションの分野として、ロボット、自動車、飛行機などのメトカトロニクスのプラントのシミュレーション用ソフトウェアがある。電子部品とソフトウェア技術に発展の恩恵により、ロボット、自動車、飛行機などでは、神経のように張り巡らされたワイヤ結線や無線ＬＡＮなどを利用して、大部分の制御が電子的に行われる。

それらは、本来的には機械的装置であるのに、大量の制御ソフトウェアをも内蔵している。そのため、製品の開発に当たっては、制御プログラムの開発とそのテストに、長い時間と、膨大な費用と、多数の人員を費やす必要が出てきた。

このようなテストにために従来行われている技法として、ＨＩＬＳ(Hardware In the Loop Simulation)がある。特に、自動車全体の電子制御ユニット（ＥＣＵ）をテストする環境は、フルビークルＨＩＬＳと呼ばれる。フルビークルＨＩＬＳにおいては、実験室内で、本物のＥＣＵが、エンジン、トランスミッション機構などをエミュレーションする専用のハードウェア装置に接続され、所定のシナリオに従って、テストが行われる。ＥＣＵからの出力は、監視用のコンピュータに入力され、さらにはディスプレイに表示されて、テスト担当者がディスプレイを眺めながら、異常動作がないかどうか、チェックする。

そこで近年、高価なエミュレーション用ハードウェア装置を使うことなく、ソフトウェアで構成する手法が提案されている。この手法は、ＳＩＬＳ(Software In the Loop Simulation)と呼ばれ、ＥＣＵに搭載されるマイクロコンピュータ、入出力回路、制御のシナリオ、エンジンやトランスミッションなどのプラントを全て、ソフトウェア・シミュレータで構成する技法である。これによれば、ＥＣＵのハードウェアが存在しなくても、テストを実行可能である。

このようなＳＩＬＳの構築を支援するシステムとして例えば、MathWorks社から入手可能なシミュレーション・モデリング・システムである、MATLAB(R)/Simulink(R)がある。MATLAB(R)/Simulink(R)を使用すると、図１に示すように、画面上にグラフィカル・インターフェースによって、矩形で示す機能ブロックを配置し、矢印のようにその処理の流れを指定することによって、シミュレーション・プログラムを作成することができる。これらのブロック線図は、シミュレーションの１タイムステップ分の処理を表しており、これが所定回繰り返されることにより、シミュレーション対象となるシステムの時系列における振る舞いを得ることができる。こうして、MATLAB(R)/Simulink(R)上で、機能ブロックなどのブロック線図が作成されると、Real-Time Workshop(R)の機能により、等価な機能のＣ言語など既知のコンピュータ言語のソース・コードに変換することができる。このＣ言語のソース・コードをコンパイルすることにより、別のコンピュータ・システムで、ＳＩＬＳとして、シミュレーションを実行することができる。

そこで、図２に示すように、機能ブロックを、Ａ、Ｂ、Ｃ及びＤのように、機能ブロックの複数の集合に分け、それらをＣＰＵによって計算する技法が従来より実施されている。

このようにして、機能ブロックの個々の集合毎のコードをコンパイルして実行可能コードを生成し、プロセッサ上で実行させることができる。

そのような実行コードを、コンピュータ・システムの上で走らせて、シミュレーションの結果を得るということは、多くの場合、コンピュータ・システムに、微分方程式を解かせることに帰着される。特に、自動車やロボットなどのメトカトロニクスと呼ばれるシステムでは、解くべき微分方程式は、電気回路の応答システム、電子回路のフィードバック制御システム、機械的な力学方程式など、常微分方程式で記述され、従って、シミュレーションの結果を得るためには、コンピュータ・システムによって、連立常微分方程式を解く必要がある。そのような連立常微分方程式は、次のように定式化される。
y'(t) = f(t,y(t))

但しここで、tは時間、y'(t)は時間tに関する一階微分であり、
yとfは一般的にはn次元ベクトルであり、
y(t) ≡ (y^[1](t), y^[2](t), ..., y^[N](t))^T
f(t,y(t)) ≡ (f^[1](t,y(t)), f^[2](t,y(t)), ...,f^[N](t,y(t)))^T
これは、y, f ∈ Ｒ^Nとも書かれる。

すなわち、並べて書くと、
y'^[1](t) = f^[1](t,y(t))
y'^[2](t) = f^[2](t,y(t))
y'^[3](t) = f^[3](t,y(t))
...
y'^[N](t) = f^[N](t,y(t))
という連立常微分方程式になる。
典型的には、上記の１つの機能ブロックの集合に対するコードがそれぞれ、y'^[i](t) = f^[i](t,y(t)) (但し、i = 1,...,N)という１つの式の右辺に対応する。この明細書では、その計算処理の単位をストランドと呼ぶことにする。

より高速なシミュレーションを達成するために、使用するコンピュータとして高性能のものを用いるだけではなく、より合理的な計算量の計算アルゴリズムを使用するという要望が強い。

これに関して、コンピュータ上の数値計算により近似的に常微分方程式を解くためのアルゴリズムとして、ルンゲ・クッタ(Runge-Kutta)法が知られている。これは、１９００年に考案されたものであり、計算のステップは固定であり、応用例によっては、計算精度の点で、十分でない場合があった。

その後、計算精度と計算速度の両立のため、１９６０年代に、適応的に可変ステップを用いるルンゲ・クッタ・フェールベルク(Runge-Kutta-Fehlberg)法が考案された。ここでは、説明の便宜上、y'(t) = f(t,y(t))を１変数の場合に見立てて、ルンゲ・クッタ・フェールベルク法を説明する。すなわち、ここでは、y(t)も、f(t,y(t))もスカラーである。

ルンゲ・クッタ・フェールベルク法のうち、特に、例として、ＯＤＥ４５として知られているスキームについて説明する。ここで、ＯＤＥとは、ordinary differential equation = 常微分方程式の略である。

x_nを時間、
k_i(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
a_i, b_ij, c_i,c^* _i (ここで、i,jは添字)を所定の定数とすると、
k₁ = hf(x_n,y_n)
k₂ = hf(x_n+a₂h,y_n+b₂₁k₁)
k₃ = hf(x_n+a₃h,y_n+b₃₁k₁+b₃₂k₂)
k₄ = hf(x_n+a₄h,y_n+b₄₁k₁+b₄₂k₂+b₄₃k₃)
k₅ = hf(x_n+a₅h,y_n+b₅₁k₁+b₅₂k₂+b₅₃k₃+b₅₄k₄)
k₆ = hf(x_n+a₆h,y_n+b₆₁k₁+b₆₂k₂+b₆₃k₃+b₆₄k₄+b₆₅k₅)

ルンゲ・クッタ・フェールベルク法においては、計算精度を所望以内に保つために、上記のようにして計算された誤差Δが、所定の閾値Δ₀の範囲内かどうかが判断され、もしそうなら、上記の式により、ステップh₀が計算され、それが次のループのステップとなる。

もし、上記のようにして計算された誤差Δが、所定の閾値Δ₀の範囲を超えると、ルンゲ・クッタ・フェールベルク法においては、エラーと見なして、もともと使用されたステップhよりも小さいステップに変えて、計算がやり直される。

ところが、y'(t) = f(t,y(t))、y, f ∈ Ｒⁿという連立常微分方程式においては、y(t)は実は、y^[1](t), y^[2](t), ..., y^[n](t)を含むため、例えば、y'^[i](t) = f^[i](t,y(t)) (iは、1からnまでの間のある数)でのみ上記のエラーが生じたとしても、y'(t) = f(t,y(t))、y, f ∈ Ｒ^N全体の計算のやり直しが必要となる。このようなy'(t) = f(t,y(t))全体の計算のやり直しは、やり直しが多く発生する場合や、全体の計算量が大きい場合には、とても大きいオーバーヘッドとなる。

特開平５−３３４３３７号公報は、常微分方程式の並列処理方式において、マルチプロセッサ環境下で、異なる初期刻み幅を複数プロセッサの各々に割付け、数値積分部で近似値を求め、局所的誤差計算部で局所的誤差を、許容誤差計算部で許容誤差を求め、この二つの誤差の関係から、刻み幅調節部で各プロセッサの刻み幅を調節すると共に、評価装置では複数プロセッサで処理された刻み幅の最大値を選択して各プロセッサに再設定することを開示する。

W. H. Enright, K. R. Jackson, S. P. Norsett, P. G. Thomsen, "Interpolants for Runge-Kutta formulas" ACM Transactions on Mathematical Software, Volume 12, Issue 3 (September 1986) は、ルンゲ・クッタ法及び、ルンゲ・クッタ・フェールベルク法における補間式を構成する方法について記述する。

特開平５−３３４３３７号公報

W. H. Enright, K. R. Jackson, S. P. Norsett, P. G. Thomsen, "Interpolants for Runge-Kutta formulas" ACM Transactions on Mathematical Software, Volume 12, Issue 3 (September 1986)

上記従来技術は、コンピュータによる常微分方程式の計算における、刻み量あるいは補間量の計算技法については記述するものの、ルンゲ・クッタ法または、ルンゲ・クッタ・フェールベルク法等に従い、連立常微分方程式を解く際の計算量の低減技法については、示唆するものではない。

従って、本発明の目的は、シミュレーション・システムなどで使用される連立常微分方程式をコンピュータで解く計算処理において、計算量を低減する技法を提供することにある。

本発明は、上記課題を解決するためになされたものであり、以下に記述する処理によって、コンピュータによる連立常微分方程式の計算量を低減させる。

本発明は、上記背景技術で述べた、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法を用いた、コンピュータによる常微分方程式の数値解法を前提とする。すなわち、連立常微分方程式の個々の常微分方程式を解く際に、次数Ｎ（Ｎは所定の整数）の計算項と、次数Ｎ＋１の計算項の差Δを計算して、その差が所定の閾値Δ₀よりも小さいかどうか判断し、もしΔ =< Δ₀であるなら、Δ₀/Δで決まる所定の計算式に従い、ステップ・サイズを決定して、次の計算に進む。ここまでは、従来のルンゲ・クッタ・フェールベルク法による連立常微分方程式の計算方法と同じである。

ところが、次数Ｎの計算項と、次数Ｎ＋１の計算項の差Δが閾値Δ₀よりも大きい場合に、本発明の計算処理は異なる。すなわち、従来技術では、連立常微分方程式を構成する常微分方程式のうちの１つでもエラーになると、ステップを戻って、ステップ・サイズを減らし、連立常微分方程式全体の計算をやり直す必要があった。

しかし、本発明においては、連立常微分方程式を構成する常微分方程式のうち、Δ > Δ₀となったエラーを生じた常微分方程式を計算するストランドに再計算が命じられる。再計算のため、エラーを生じたストランドには、Δ₀/Δで決まるステップサイズがセットされ、その分だけ時間ステップが先に進む。そこから更に、再計算を行っていない他のストランドのステップに追いつくまで時間ステップを進めるためには、再計算を行っているストランドに関連する他のストランドの計算結果の、当該ステップ開始時点での値が必要となる。この値は、再計算を行っていない他のストランドの計算を再実行するのではなく、他のストランドの前回ステップにおける値と、最終ステップにおける値との補間値として計算される。

そのような常微分方程式の補間値の好適な計算アルゴリズムは、上述のW. H. Enright, K. R. Jackson, S. P. Norsett, P. G. Thomsen, "Interpolants for Runge-Kutta formulas" ACM Transactions on Mathematical Software, Volume 12, Issue 3 (September 1986)に記述されている。しかし、本発明で利用可能な補間値の計算アルゴリズムは、これには限定されず、所望の精度を満たす任意の補間値の計算アルゴリズムを使用することができる。

補間値で以って再計算することにより、エラーを生じた常微分方程式を計算するストランドにおいて、エラーを所定の閾値Δ₀よりも小さく保ちながら、エラーが十分少なかった他のストランドと足並みが揃うように、時間ステップを進める。その後、連立常微分方程式全体の計算が進められる。

この発明によれば、コンピュータによる、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法に従う連立常微分方程式を計算する際に、閾値よりも大きい誤差を生じたストランドのみ再計算すればよくなるので、計算量が減り、処理を高速化することが可能となる。

機能ブロックの例を示す図である。機能ブロックのストランド作成の例を示す図である。ハードウェア構成を示すブロック図である。論理構成の機能ブロック図である。連立常微分方程式を並列処理で計算するステップを示す概要フローチャートである。ＯＤＥ４５のスキームに本発明を適用した処理で常微分方程式を解く処理を示すフローチャートである。連立常微分方程式を並列処理で計算する際の、補間値のタイミングを示す図である。

以下、図面を参照して、本発明の一実施例の構成及び処理を説明する。以下の記述では、特に断わらない限り、図面に亘って、同一の要素は同一の符号で参照されるものとする。なお、ここで説明する構成と処理は、一実施例として説明するものであり、本発明の技術的範囲をこの実施例に限定して解釈する意図はないことを理解されたい。

図３を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・バス３０２には、ＣＰＵ３０４と、主記憶（ＲＡＭ）３０６と、ハードディスク・ドライブ（ＨＤＤ）３０８と、キーボード３１０と、マウス３１２と、ディスプレイ３１４が接続されている。ＣＰＵ３０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）４、インテル社のＣｏｒｅ（商標）２ＤＵＯ、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶３０６は、好適には、１ＧＢ以上の容量、より好ましくは、２ＧＢ以上の容量をもつものである。

ハードディスク・ドライブ３０８には、オペレーティング・システムが、格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ７、ＷｉｎｄｏｗｓＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃＯＳ（商標）などの、ＣＰＵ３０４に適合する任意のものでよい。

ハードティスク・ドライブ３０８にはさらに、MATLAB(R)/Simulink(R)、Ｃコンパイラまたは、Ｃ＋＋コンパイラ、後述する本発明に係る解析、ストランド作成のためのモジュール、ＣＰＵ割り当て用コード生成モジュールなどが格納されており、オペレータのキーボードやマウス操作に応答して、メイン・メモリ３０６にロードされて実行される。

尚、使用可能なシミュレーション・モデリング・ツールは、MATLAB(R)/Simulink(R)に限定されず、オープンソースのScilab/Scicosなど任意のシミュレーション・モデリング・ツールを使用することが可能である。

あるいは、場合によっては、シミュレーション・モデリング・ツールを使わず、直接、Ｃ、Ｃ＋＋などでシミュレーション・システムのソース・コードを書くことも可能であり、その場合にも、個々の機能の少なくとも一部が、常微分方程式の計算として記述できるなら、本発明は適用可能である。

キーボード１１０及びマウス１１２は、オペレーティング・システムが提供するグラフィック・ユーザ・インターフェースに従い、シミュレーションを開始したり、所定の式のパラメータを提供するために使用される。

ディスプレイ３１４は、連立常微分方程式の解としてのシミュレーションの挙動を、適宜表示するために使用される。

図４は、本発明の実施例に係る一実施例の機能ブロック図である。各々のブロックは、基本的に、ハードティスク・ドライブ３０８に格納されている。

図４において、シミュレーション・モデリング・ツール４０２は、MATLAB(R)/Simulink(R)、Scilab/Scicosなどの既存の任意のモデリング・ツールを使用することができる。シミュレーション・モデリング・ツール４０２は、基本的には、オペレータが、ディスプレイ３１４上でＧＵＩ的に機能ブロックを配置し、数式など必要な属性を記述し、必要に応じて、機能ブロック間を関連付けてブロック線図を記述することを可能ならしめるような機能をもつ。シミュレーション・モデリング・ツール４０２はさらに、記述されたブロック線図に等価な機能を記述するＣのソースコードを出力する機能をもつ。Ｃ以外にも、Ｃ＋＋、ＦＯＲＴＲＡＮなどを使用することができる。特に、ＭＤＬファイルは、Simulink(R)独自のフォーマットであり、機能ブロック間の依存関係を記述するために、ＭＤＬファイルを生成することができる。

なお、シミュレーション・モデリング・ツールは、別のパーソナル・コンピュータに導入して、そこで生成されたソース・コードを、ネットワークなどを経由して、ハードティスク・ドライブ３１６にダウンロードするようにすることもできる。

こうして出力されたソース・コード４０４は、ハードティスク・ドライブ３０８に保存される。なお、ソース・コード４０４以外に、機能ブロック間の依存関係を記述するためのＭＤＬファイルを保存してもよい。

解析モジュール４０６は、ソースコード４０４を入力して構文解析し、ブロックのつながりを、グラフ表現に変換する。グラフ表現のデータは、好適には、ハードディスク・ドライブ３１６に格納される。なお、コンピュータ上のグラフ表現のデータ構造は周知であるので、ここでは説明を省略する。

ストランド作成モジュール４０８は、解析モジュール４０６によって作成されたグラフ表現を読み取って、これには限定しないが、次のような方法で、各積分ブロックごとにストランドを構成する。すなわち、ブロック線図中の各積分ブロックごとに、その積分ブロックからフローの順方向に辿って、再び（自身および他の）積分ブロックが辿られるまでの全てのブロックの集合を見つけ、そのブロックの集合を一つのストランドとする。この操作は、ブロック線図から各常微分方程式を取り出すことに相当する。

コード生成モジュール４１０は、ストランド作成モジュール４０８が生成したストランドの情報に基づき、コンパイラ４１２がコンパイルするためのソースコードを生成する。コンパイラ４１２が想定するプログラミング言語としては、好適にはＣ、Ｃ＋＋、Ｃ＃、Ｊａｖａ（商標）などの任意のプログラミング言語を使用することができ、コード生成モジュール４１０はそれに対応して、ストランド毎に、ソースコードを生成することになる。

コンパイラ４１２が生成した実行可能バイナリ・コード（図示しない）は、ストランド毎に、オペレーティング・システムの作用により、実行環境４１４で実行される。

本発明においては、シミュレーションの処理全体は、下記のような、Ｎ個の常微分方程式からなる、連立常微分方程式として記述される。

y'(t) = f(t,y(t)) y, f ∈ Ｒ^N

但しここで、tは時間、y'(t)は時間tに関する一階微分であり、
yとfは一般的にはn次元ベクトルであり、
y(t) ≡ (y^[1](t), y^[2](t), ..., y^[N](t))^T
f(t,y(t)) ≡ (f^[1](t,y(t)), f^[2](t,y(t)), ...,f^[N](t,y(t)))^T

並べて書くと、
y'^[1](t) = f^[1](t,y(t))
y'^[2](t) = f^[2](t,y(t))
y'^[3](t) = f^[3](t,y(t))
...
y'^[N](t) = f^[N](t,y(t))
という連立常微分方程式になる。
典型的には、上記の１つのストランドが、y'^[i](t) = f^[i](t,y(t)) (但し、i = 1,...,N)という１つの式の右辺に対応する。

すなわち、個々のストランドが、連立常微分方程式のうちの１つの常微分方程式である、y'^[i](t) = f^[i](t,y(t))に関する右辺の数値計算を実行する処理を行う。

なお、シミュレーションの基礎となる制御理論によれば、制御系の特性方程式や応答関数は、ラプラス変換のパラメータsで記述され、これらは多くの場合常微分方程式に帰着されるので、シミュレーション・モデルを連立常微分方程式で記述するというこの実施例の前提は、十分に一般性を有すると言える。

図５は、実行環境４１４で、コンパイラ４１２によって生成された実行可能モジュール（図示しない）が、連立常微分方程式を解くための数値計算を実行する処理を示す概要フローチャートである。

図５の処理を実行するために、この実施例によれば、ストランドに対応する個々の常微分方程式を解くための個別のコード以外に、全体の処理を統合する実行コード（図示しない）が生成される。

図５のステップ５０２において、管理用の実行コードは、y'(t) = f(t,y(t)) y, f ∈ Ｒ^Nという連立常微分方程式における、t=t₀での初期値y(t₀)を与える。これは実際は、実行させるシミュレーションに応じて、オペレータが事前にセットしておくものである。

ステップ５０４では、管理用の実行コードは、個々のストランドの常微分方程式を解くストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎを一斉にスタートさせ、並列に動作させる。

起動されると、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは各々、本発明の特徴を含むＯＤＥ４５のスキームにより、常微分方程式を解く計算を行う。各々のストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの処理の詳細は、図６のフローチャートで詳しく説明する。

管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎが全て終了するのを待つ。というのは、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは、計算量が同等とは限らず、また、本発明の特徴によれば、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎのどのストランドも、エラーが検出されたことに応答して再計算を行うことがあるので、その分、遅延する可能性があるからである。

管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの１ステップ分の計算ステップが終了したのを確認すると、ステップ５１０で、シミュレーションの終了かどうかの判断を行う。

シミュレーションの終了は、オペレータの操作、予定したシナリオの終了、または予定時間経過などによって決定される。管理用の実行コードは、シミュレーションの終了でないと判断すると、ステップ５０４に戻り、次のステップの計算に進む。もしシミュレーションの終了であるなら、管理用の実行コードは、シミュレーションを終了させる。

次に、図６のフローチャートを参照して、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎで個々に実行される計算処理について説明する。

背景技術でも説明したが、この実施例では、ルンゲ・クッタ・フェールベルク法において、ＯＤＥ４５のスキームで、連立常微分方程式をコンピュータによる数値計算により解くものとする。なお、本発明は、ルンゲ・クッタ・フェールベルク法の特定の次数に限定されず、ＯＤＥ３４、ＯＤＥ５６などのその他のスキームも使用可能である。なお、ルンゲ・クッタ・フェールベルク法のより詳しい説明は、Ward Cheney, David Kincaid, "Numerical Mathematics and Computing", Brooks/Cole Pub Co; 6版 (2007/8/3) などを参照されたい。

ステップ６０２では、図５のステップ５０２で設定された初期値が、この特定のストランドに対応する常微分方程式の処理条件としてセットされる。

ステップ６０４では、シミュレーションの終了かどうかが判断され、もしそうなら、シミュレーションは終了される。ステップ６０４は実質的に、図５のステップ５１０と同一であるかまたは、連動する。

シミュレーションの終了でないと判断されると、ステップ６０６に進み、ＯＤＥ４５のスキームに従い、計算が行われる。その具体的計算は次のとおりである。但し、ここでは、j番目のストランドに関連する、常微分方程式y'^[j](t) = f^[j](t,y(t))に着目するものとする。

ここでは、x₀ = x_nとして、処理を開始する。
x_nを時間、
k_i(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
a_i, b_ij, c_i,c^* _i (ここで、i,jは添字)を所定の定数とすると、
k₁ = hf^[j](x_n,y_n)
k₂ = hf^[j](x_n+a₂h,y_n+b₂₁k₁)
k₃ = hf^[j](x_n+a₃h,y_n+b₃₁k₁+b₃₂k₂)
k₄ = hf^[j](x_n+a₄h,y_n+b₄₁k₁+b₄₂k₂+b₄₃k₃)
k₅ = hf^[j](x_n+a₅h,y_n+b₅₁k₁+b₅₂k₂+b₅₃k₃+b₅₄k₄)
k₆ = hf^[j](x_n+a₆h,y_n+b₆₁k₁+b₆₂k₂+b₆₃k₃+b₆₄k₄+b₆₅k₅)

ステップ６０８では、次の式により誤差Δが計算される。

ステップ６１０では、計算した誤差Δが、予め定めた閾値Δ₀より大きいかどうか、判断される。

ステップ６１０で、Δ > Δ₀なら、ステップ６１２で、次の式により、h_pを計算する。ただし、計算されたh_pとx₀を加えた値がx_n+1を超える場合には、x₀ + h_p = x_n+1とする。

次に、ステップ６１４では、次の式に従い、i = 1, .., Nの５次エルミート補間値U^[i]が計算される。

但し、

ここで、f^[i] _lは、次の式であらわされる。

さらに、

ところで、k₁〜k₆の右辺にあらわれるy_nは実際は、
y_n ≡ y_n(y^[1],y^[2],...,y^[N])のように、
y^[1],y^[2],...,y^[N]の関数である。

そこで、数５で決定される補間式Ｕ^[i](x_n+h_p)を、y^[i]に置き換えて(但し、i = 1, ,,,, Nの各々に亘る)、ステップ６１６では、x₀+h_pから出発して、ステップ６０６と同様のＯＤＥ４５の計算が行われる。なお、補間式Ｕ^[i](x_n+h_p)のより詳しい説明については、前述のW. H. Enright, K. R. Jackson, S. P. Norsett, P. G. Thomsen, "Interpolants for Runge-Kutta formulas" ACM Transactions on Mathematical Software, Volume 12, Issue 3 (September 1986)などを参照されたい。

エラーを生じたストランドｊは、対応するストランドi(但し、i = 1, ,,,, Nで、jは除く)にx_n+h_pを渡し、ストランドiが補間式Ｕ^[i](x_n+h_p)を計算して、その値を、ストランドｊに送り返してくるという処理により、ストランドｊは、Ｕ^[i](x_n+h_p)(但し、i = 1, ,,,, N)の値を集める。この処理は、図７に関連して、後で説明する。

ステップ６１８では、ステップ６０８と同様の処理で誤差が計算され、ステップ６２０では、ステップ６１０と同様の処理で、誤差と閾値が比較され、誤差が閾値より大きいと判断されると、ステップ６１２に戻る。

ステップ６２０で、誤差が閾値より小さいか閾値と等しいと判断されると、ステップ６２２に進む。

ステップ６２２では、x₀+h_p < x_n+1かどうかが判断され、そうであれれば、ステップ６１２に戻って、h_pが再計算され、その再計算されたh_pに基づき、ステップ６１４で補間値Ｕが計算され、ステップ６１６でＯＤＥ４５の計算が行われる。

ルンゲ・クッタ・フェールベルク法の誤差式の性質により、ステップ６１０からステップ６１２に分岐した直後は必ず、ステップ６２２での判断はx₀+h_p < x_n+1であり、その後、ステップ６１２、ステップ６１４、ステップ６１６、ステップ６１８、ステップ６２０、ステップ６２２のループが廻る毎にh_pの値は次第に増加する。

そうして、x₀+h_pがx_n+1に等しくなると、ステップ６２４に進み、下記の式で次ステップhの計算をする。ここで、式右辺のhは、h_pで読み替えるものとする。

ステップ６２６では、x₀ = x_n+1によって時間が進められ、ステップ６２８では、
状態y₀が、y_n+1で更新される。そして処理は、判断ステップ６０４に戻る。

判断ステップ６１０に戻って、誤差Δが閾値Δ₀に等しいか、閾値Δ₀より小さい場合は、既に説明したステップ６２４、６２６、及び６２８を経て、判断ステップ６０４に戻る。

図７は、ストランド１〜Ｎの動作を示す概要タイミング図である。ストランド１〜Ｎに対応する実行コードによる計算の結果、ストランドｊだけが、エラー、すなわち、ステップ６１０で誤差Δ > 閾値Δ₀となったとする。すると、ストランド１〜ｊ−１と、ストランドｊ＋１〜Ｎが、ステップ６２４、ステップ処理６２６、及びステップ６２８を経てx_n+1に到達している一方で、ストランドｊは、ステップ６１２から再計算の処理に入る。ストランドｊは、ステップ６１４で補間値を計算するが、その際、自己充足的に計算できるのは自ストランドの補間値だけなので、その常微分方程式の右辺に現れるy^[i]に対応するそれぞれのストランドに、x₀+h_pの値を渡して、それらのストランドから、計算された補間値を受け取る。それぞれのストランドでは、ステップ６１４で説明した補間の計算式が用いられる。

図７では、そのようにして、対応するそれぞれのストランドで計算された補間値は、それぞれ、y^*[1] _n,y^*[2] _n,...,y^*[N] _nとして示されている。

再計算のため、対応するそれぞれのストランドは、ストランドｊに補間値を通信するが、それ以外は、ストランドｊだけが再計算のループを廻すので、それぞれのストランドにはほとんど計算の負荷がかからない。

こうして、ストランドｊだけが再計算を終了すると、図５のステップ５０８で、ストランドｊが他のストランドとステップを合わせることになり、次のステップを計算する準備ができたことになる。尚、ストランドｊのステップを揃えるための計算は、実質的に、ステップ６１２で実行される。

以上、特定の実施例に基づき本発明を説明してきたが、これには限定されず、例えば、ルンゲ・クッタ・フェールベルク法は、ＯＤＥ４５以外に、ＯＤＥ３４、ＯＤＥ５６など、任意の次数のスキームを使用することができる。その際の補間の計算式は、例えば、前述のW. H. Enright, K. R. Jackson, S. P. Norsett, P. G. Thomsen, "Interpolants for Runge-Kutta formulas" ACM Transactions on Mathematical Software, Volume 12, Issue 3 (September 1986)に記載されている式から計算することができる。

また、本発明は、ルンゲ・クッタ・フェールベルク法に限定されず、それを一般化した埋め込み型ルンゲ・クッタ法、例えば、Dormand prince法にも適用可能である。

さらに、補間の計算式は、実施例に示されている５次のエルミート補間式に限定されず、必要な精度を満たすものであるなら、線形補間など、ニュートンの補間式など、任意の補間式を用いることができる。

また、上記実施例は、シングル・プロセッサを例として説明されたが、これには限定されず、本発明の技法は、マルチプロセッサ、マルチコアなどの複数のＣＰＵをもつシステムにでも同様に適用することができる。

３１６ハードディスク・ドライブ
４０２シミュレーション・モデリング・ツール
４０４ソース・コード
４０６解析モジュール
４０４ソースコード
４０８ストランド作成モジュール
４１０コード生成モジュール
４１２コンパイラ
４１４実行環境

Claims

複数の常微分方程式からなる連立常微分方程式を埋め込み型ルンゲ・クッタ法を用いて解いて、シミュレーションの結果を得るための方法であって、コンピュータが、
前記連立常微分方程式の個々の常微分方程式を解く際に、
当該個々の常微分方程式それぞれについて、次数Ｎ（ここで、前記次数Ｎは前記埋め込み型ルンゲ・クッタ法の次数によって定まる整数の定数である）の状態値と次数Ｎ＋１の状態値との差を計算するステップと、
前記差が所定の閾値を超えたことに応答して、前記閾値を前記差で除して決まるステップ・サイズを計算するステップと、
前記差が所定の閾値を超えた常微分方程式それぞれについて、
前記差が所定の閾値を超えた当該常微分方程式において計算されたステップ・サイズに対応し且つ当該ステップ・サイズが計算された常微分方程式の補間値を計算するステップと、
前記差が所定の閾値を超えた当該常微分方程式において計算されたステップ・サイズと前記計算された補間値を用いて、前記差が所定の閾値を超えた前記常微分方程式を解くことを再実行するステップと、
前記差が所定の閾値を超えた常微分方程式の再実行された値を含む前記連立常微分方程式の解としてのシミュレーションの挙動をディスプレイ上に表示するステップと
を実行することを含む、前記方法。
前記コンピュータが、前記補間値を計算するステップと前記再実行するステップとの間において、
前記差が所定の閾値を超えたことに応答して、前記計算された補間値を前記連立常微分方程式のうちの前記差が所定の閾値を超えている連立常微分方程式を構成する他の計算式において受け取るステップ
を実行することをさらに含む、請求項１に記載の方法。
前記埋め込み型ルンゲ・クッタ法において、常微分方程式（ＯＤＥ；ordinary differential equation）４５、ＯＤＥ３４、又は、ＯＤＥ５６が使用される、請求項１又は２に記載の方法。
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項１又は２に記載の方法。
前記補間値が、５次エルミート補間値である、請求項１〜４のいずれか一項に記載の方法。
複数の常微分方程式からなる連立常微分方程式を埋め込み型ルンゲ・クッタ法を用いて解いて、シミュレーションの結果を得るためのシステムあって、
前記連立常微分方程式の個々の常微分方程式を解く際に、
当該個々の常微分方程式それぞれについて、次数Ｎ（ここで、前記次数Ｎは前記埋め込み型ルンゲ・クッタ法の次数によって定まる整数の定数である）の状態値と次数Ｎ＋１の状態値との差を計算する手段と、
前記差が所定の閾値を超えたことに応答して、前記閾値を前記差で除して決まるステップ・サイズを計算する手段と、
前記差が所定の閾値を超えた常微分方程式それぞれについて、
前記差が所定の閾値を超えた当該常微分方程式において計算されたステップ・サイズに対応し且つ当該ステップ・サイズが計算された常微分方程式の補間値を計算する手段と、
前記差が所定の閾値を超えた当該常微分方程式において計算されたステップ・サイズと前記計算された補間値を用いて、前記差が所定の閾値を超えた前記常微分方程式を再計算するステップを実行させる手段と、
前記差が所定の閾値を超えた常微分方程式の再計算された値を含む前記連立常微分方程式の解としてのシミュレーションの挙動をディスプレイ上に表示する手段と
を備えている、前記システム。
前記差が所定の閾値を超えたことに応答して、前記計算された補間値を前記連立常微分方程式のうちの前記差が所定の閾値を超えている連立常微分方程式を構成する他の計算式において受け取る手段
をさらに備えており、当該受け取りが前記補間値を計算することと前記常微分方程式を解くことを再実行することとの処理の間において行われる、請求項６に記載のシステム。
前記埋め込み型ルンゲ・クッタ法において、常微分方程式（ＯＤＥ；ordinary differential equation）４５、ＯＤＥ３４、又は、ＯＤＥ５６が使用される、請求項６又は７に記載のシステム。
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項６又は７に記載のシステム。
前記補間値が、５次エルミート補間値である、請求項６〜９のいずれか一項に記載のシステム。
複数の常微分方程式からなる連立常微分方程式を埋め込み型ルンゲ・クッタ法を用いて解いて、シミュレーションの結果を得るためのコンピュータ・プログラムあって、コンピュータに、請求項１〜５のいずれか一項に記載の方法の各ステップを実行させる、前記コンピュータ・プログラム。