常微分方程式を解くための方法、プログラム及びシステム
【課題】連立常微分方程式の計算量を低減する。
【解決手段】ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法によって、連立常微分方程式の個々の常微分方程式を解く際に、次数Ｎの計算項と、次数Ｎ＋１の計算項の差Δを計算して、その差が所定の閾値Δ0よりも小さいかどうか判断し、もしΔ=<Δ0であるなら、Δ0/Δで決まる所定の計算式に従い、ステップ・サイズを決定して次の計算に進み、Δ>Δ0となったエラーを生じた常微分方程式を計算するストランドのみに再計算が命じられる。再計算のストランドには、Δ0/Δで決まるステップ・サイズがセットされる。そこから、補間値で以って再計算することにより、エラーが所定の閾値Δ0よりも小さくなると、エラーを生じていない他の常微分方程式を計算するストランドとの足並みが揃うので、次に連立常微分方程式全体の計算が進められる。
【解決手段】ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法によって、連立常微分方程式の個々の常微分方程式を解く際に、次数Ｎの計算項と、次数Ｎ＋１の計算項の差Δを計算して、その差が所定の閾値Δ0よりも小さいかどうか判断し、もしΔ=<Δ0であるなら、Δ0/Δで決まる所定の計算式に従い、ステップ・サイズを決定して次の計算に進み、Δ>Δ0となったエラーを生じた常微分方程式を計算するストランドのみに再計算が命じられる。再計算のストランドには、Δ0/Δで決まるステップ・サイズがセットされる。そこから、補間値で以って再計算することにより、エラーが所定の閾値Δ0よりも小さくなると、エラーを生じていない他の常微分方程式を計算するストランドとの足並みが揃うので、次に連立常微分方程式全体の計算が進められる。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、コンピュータを利用したシミュレーション・システム等において使用される常微分方程式を解くための方法、プログラム及びシステムに関し、特に、常微分方程式を解くための計算を高速化し、あるいは計算量を減らすための技法に関するものである。
【背景技術】
【０００２】
従来より、科学技術計算、シミュレーションなどの分野で、コンピュータが利用されている。
【０００３】
最近になって特に盛んに開発されるようになってきたシミュレーションの分野として、ロボット、自動車、飛行機などのメトカトロニクスのプラントのシミュレーション用ソフトウェアがある。電子部品とソフトウェア技術に発展の恩恵により、ロボット、自動車、飛行機などでは、神経のように張り巡らされたワイヤ結線や無線ＬＡＮなどを利用して、大部分の制御が電子的に行われる。
【０００４】
それらは、本来的には機械的装置であるのに、大量の制御ソフトウェアをも内蔵している。そのため、製品の開発に当たっては、制御プログラムの開発とそのテストに、長い時間と、膨大な費用と、多数の人員を費やす必要が出てきた。
【０００５】
このようなテストにために従来行われている技法として、ＨＩＬＳ(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 = 常微分方程式の略である。
【００１７】
xnを時間、
ki(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
ai, bij, ci,c*i (ここで、i,jは添字)を所定の定数とすると、
k1 = hf(xn,yn)
k2 = hf(xn+a2h,yn+b21k1)
k3 = hf(xn+a3h,yn+b31k1+b32k2)
k4 = hf(xn+a4h,yn+b41k1+b42k2+b43k3)
k5 = hf(xn+a5h,yn+b51k1+b52k2+b53k3+b54k4)
k6 = hf(xn+a6h,yn+b61k1+b62k2+b63k3+b64k4+b65k5)
【数１】
【００１８】
ルンゲ・クッタ・フェールベルク法においては、計算精度を所望以内に保つために、上記のようにして計算された誤差Δが、所定の閾値Δ0の範囲内かどうかが判断され、もしそうなら、上記の式により、ステップh0が計算され、それが次のループのステップとなる。
【００１９】
もし、上記のようにして計算された誤差Δが、所定の閾値Δ0の範囲を超えると、ルンゲ・クッタ・フェールベルク法においては、エラーと見なして、もともと使用されたステップhよりも小さいステップに変えて、計算がやり直される。
【００２０】
ところが、y'(t) = f(t,y(t))、y, f ∈ Ｒnという連立常微分方程式においては、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)
【発明の概要】
【発明が解決しようとする課題】
【００２５】
上記従来技術は、コンピュータによる常微分方程式の計算における、刻み量あるいは補間量の計算技法については記述するものの、ルンゲ・クッタ法または、ルンゲ・クッタ・フェールベルク法等に従い、連立常微分方程式を解く際の計算量の低減技法については、示唆するものではない。
【００２６】
従って、本発明の目的は、シミュレーション・システムなどで使用される連立常微分方程式をコンピュータで解く計算処理において、計算量を低減する技法を提供することにある。
【課題を解決するための手段】
【００２７】
本発明は、上記課題を解決するためになされたものであり、以下に記述する処理によって、コンピュータによる連立常微分方程式の計算量を低減させる。
【００２８】
本発明は、上記背景技術で述べた、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法を用いた、コンピュータによる常微分方程式の数値解法を前提とする。すなわち、連立常微分方程式の個々の常微分方程式を解く際に、次数Ｎ（Ｎは所定の整数）の計算項と、次数Ｎ＋１の計算項の差Δを計算して、その差が所定の閾値Δ0よりも小さいかどうか判断し、もしΔ =< Δ0であるなら、Δ0/Δで決まる所定の計算式に従い、ステップ・サイズを決定して、次の計算に進む。ここまでは、従来のルンゲ・クッタ・フェールベルク法による連立常微分方程式の計算方法と同じである。
【００２９】
ところが、次数Ｎの計算項と、次数Ｎ＋１の計算項の差Δが閾値Δ0よりも大きい場合に、本発明の計算処理は異なる。すなわち、従来技術では、連立常微分方程式を構成する常微分方程式のうちの１つでもエラーになると、ステップを戻って、ステップ・サイズを減らし、連立常微分方程式全体の計算をやり直す必要があった。
【００３０】
しかし、本発明においては、連立常微分方程式を構成する常微分方程式のうち、Δ > Δ0となったエラーを生じた常微分方程式を計算するストランドに再計算が命じられる。再計算のため、エラーを生じたストランドには、Δ0/Δで決まるステップサイズがセットされ、その分だけ時間ステップが先に進む。そこから更に、再計算を行っていない他のストランドのステップに追いつくまで時間ステップを進めるためには、再計算を行っているストランドに関連する他のストランドの計算結果の、当該ステップ開始時点での値が必要となる。この値は、再計算を行っていない他のストランドの計算を再実行するのではなく、他のストランドの前回ステップにおける値と、最終ステップにおける値との補間値として計算される。
【００３１】
そのような常微分方程式の補間値の好適な計算アルゴリズムは、上述の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)に記述されている。しかし、本発明で利用可能な補間値の計算アルゴリズムは、これには限定されず、所望の精度を満たす任意の補間値の計算アルゴリズムを使用することができる。
【００３２】
補間値で以って再計算することにより、エラーを生じた常微分方程式を計算するストランドにおいて、エラーを所定の閾値Δ0よりも小さく保ちながら、エラーが十分少なかった他のストランドと足並みが揃うように、時間ステップを進める。その後、連立常微分方程式全体の計算が進められる。
【発明の効果】
【００３３】
この発明によれば、コンピュータによる、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法に従う連立常微分方程式を計算する際に、閾値よりも大きい誤差を生じたストランドのみ再計算すればよくなるので、計算量が減り、処理を高速化することが可能となる。
【図面の簡単な説明】
【００３４】
【図１】機能ブロックの例を示す図である。
【図２】機能ブロックのストランド作成の例を示す図である。
【図３】ハードウェア構成を示すブロック図である。
【図４】論理構成の機能ブロック図である。
【図５】連立常微分方程式を並列処理で計算するステップを示す概要フローチャートである。
【図６】ＯＤＥ４５のスキームに本発明を適用した処理で常微分方程式を解く処理を示すフローチャートである。
【図７】連立常微分方程式を並列処理で計算する際の、補間値のタイミングを示す図である。
【発明を実施するための形態】
【００３５】
以下、図面を参照して、本発明の一実施例の構成及び処理を説明する。以下の記述では、特に断わらない限り、図面に亘って、同一の要素は同一の符号で参照されるものとする。なお、ここで説明する構成と処理は、一実施例として説明するものであり、本発明の技術的範囲をこの実施例に限定して解釈する意図はないことを理解されたい。
【００３６】
図３を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・バス３０２には、ＣＰＵ３０４と、主記憶（ＲＡＭ）３０６と、ハードディスク・ドライブ（ＨＤＤ）３０８と、キーボード３１０と、マウス３１２と、ディスプレイ３１４が接続されている。ＣＰＵ３０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）４、インテル社のＣｏｒｅ（商標） ２ ＤＵＯ、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶３０６は、好適には、１ＧＢ以上の容量、より好ましくは、２ＧＢ以上の容量をもつものである。
【００３７】
ハードディスク・ドライブ３０８には、オペレーティング・システムが、格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ ７、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃ ＯＳ（商標）などの、ＣＰＵ３０４に適合する任意のものでよい。
【００３８】
ハードティスク・ドライブ３０８にはさらに、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=t0での初期値y(t0)を与える。これは実際は、実行させるシミュレーションに応じて、オペレータが事前にセットしておくものである。
【００６０】
ステップ５０４では、管理用の実行コードは、個々のストランドの常微分方程式を解くストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎを一斉にスタートさせ、並列に動作させる。
【００６１】
起動されると、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは各々、本発明の特徴を含むＯＤＥ４５のスキームにより、常微分方程式を解く計算を行う。各々のストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの処理の詳細は、図６のフローチャートで詳しく説明する。
【００６２】
管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎが全て終了するのを待つ。というのは、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは、計算量が同等とは限らず、また、本発明の特徴によれば、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎのどのストランドも、エラーが検出されたことに応答して再計算を行うことがあるので、その分、遅延する可能性があるからである。
【００６３】
管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの１ステップ分の計算ステップが終了したのを確認すると、ステップ５１０で、シミュレーションの終了かどうかの判断を行う。
【００６４】
シミュレーションの終了は、オペレータの操作、予定したシナリオの終了、または予定時間経過などによって決定される。管理用の実行コードは、シミュレーションの終了でないと判断すると、ステップ５０４に戻り、次のステップの計算に進む。もしシミュレーションの終了であるなら、管理用の実行コードは、シミュレーションを終了させる。
【００６５】
次に、図６のフローチャートを参照して、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎで個々に実行される計算処理について説明する。
【００６６】
背景技術でも説明したが、この実施例では、ルンゲ・クッタ・フェールベルク法において、ＯＤＥ４５のスキームで、連立常微分方程式をコンピュータによる数値計算により解くものとする。なお、本発明は、ルンゲ・クッタ・フェールベルク法の特定の次数に限定されず、ＯＤＥ３４、ＯＤＥ５６などのその他のスキームも使用可能である。なお、ルンゲ・クッタ・フェールベルク法のより詳しい説明は、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))に着目するものとする。
【００７０】
ここでは、x0 = xnとして、処理を開始する。
xnを時間、
ki(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
ai, bij, ci,c*i (ここで、i,jは添字)を所定の定数とすると、
k1 = hf[j](xn,yn)
k2 = hf[j](xn+a2h,yn+b21k1)
k3 = hf[j](xn+a3h,yn+b31k1+b32k2)
k4 = hf[j](xn+a4h,yn+b41k1+b42k2+b43k3)
k5 = hf[j](xn+a5h,yn+b51k1+b52k2+b53k3+b54k4)
k6 = hf[j](xn+a6h,yn+b61k1+b62k2+b63k3+b64k4+b65k5)
【数２】
【００７１】
ステップ６０８では、次の式により誤差Δが計算される。
【数３】
【００７２】
ステップ６１０では、計算した誤差Δが、予め定めた閾値Δ0より大きいかどうか、判断される。
【００７３】
ステップ６１０で、Δ > Δ0なら、ステップ６１２で、次の式により、hpを計算する。ただし、計算されたhpとx0を加えた値がxn+1を超える場合には、x0 + hp = xn+1とする。
【数４】
【００７４】
次に、ステップ６１４では、次の式に従い、i = 1, .., Nの５次エルミート補間値U[i]が計算される。
【数５】
但し、
【数６】
【００７５】
ここで、f[i]lは、次の式であらわされる。
【数７】
【００７６】
さらに、
【数８】
【００７７】
ところで、k1〜k6の右辺にあらわれるynは実際は、
yn ≡ yn(y[1],y[2],...,y[N])のように、
y[1],y[2],...,y[N]の関数である。
【００７８】
そこで、数５で決定される補間式Ｕ[i](xn+hp)を、y[i]に置き換えて(但し、i = 1, ,,,, Nの各々に亘る)、ステップ６１６では、x0+hpから出発して、ステップ６０６と同様のＯＤＥ４５の計算が行われる。なお、補間式Ｕ[i](xn+hp)のより詳しい説明については、前述の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は除く)にxn+hpを渡し、ストランドiが補間式Ｕ[i](xn+hp)を計算して、その値を、ストランドｊに送り返してくるという処理により、ストランドｊは、Ｕ[i](xn+hp)(但し、i = 1, ,,,, N)の値を集める。この処理は、図７に関連して、後で説明する。
【００８０】
ステップ６１８では、ステップ６０８と同様の処理で誤差が計算され、ステップ６２０では、ステップ６１０と同様の処理で、誤差と閾値が比較され、誤差が閾値より大きいと判断されると、ステップ６１２に戻る。
【００８１】
ステップ６２０で、誤差が閾値より小さいか閾値と等しいと判断されると、ステップ６２２に進む。
【００８２】
ステップ６２２では、x0+hp < xn+1かどうかが判断され、そうであれれば、ステップ６１２に戻って、hpが再計算され、その再計算されたhpに基づき、ステップ６１４で補間値Ｕが計算され、ステップ６１６でＯＤＥ４５の計算が行われる。
【００８３】
ルンゲ・クッタ・フェールベルク法の誤差式の性質により、ステップ６１０からステップ６１２に分岐した直後は必ず、ステップ６２２での判断はx0+hp < xn+1であり、その後、ステップ６１２、ステップ６１４、ステップ６１６、ステップ６１８、ステップ６２０、ステップ６２２のループが廻る毎にhpの値は次第に増加する。
【００８４】
そうして、x0+hpがxn+1に等しくなると、ステップ６２４に進み、下記の式で次ステップhの計算をする。ここで、式右辺のhは、hpで読み替えるものとする。
【数９】
【００８５】
ステップ６２６では、x0 = xn+1によって時間が進められ、ステップ６２８では、
状態y0が、yn+1で更新される。そして処理は、判断ステップ６０４に戻る。
【００８６】
判断ステップ６１０に戻って、誤差Δが閾値Δ0に等しいか、閾値Δ0より小さい場合は、既に説明したステップ６２４、６２６、及び６２８を経て、判断ステップ６０４に戻る。
【００８７】
図７は、ストランド１〜Ｎの動作を示す概要タイミング図である。ストランド１〜Ｎに対応する実行コードによる計算の結果、ストランドｊだけが、エラー、すなわち、ステップ６１０で誤差Δ > 閾値Δ0となったとする。すると、ストランド１〜ｊ−１と、ストランドｊ＋１〜Ｎが、ステップ６２４、ステップ処理６２６、及びステップ６２８を経てxn+1に到達している一方で、ストランドｊは、ステップ６１２から再計算の処理に入る。ストランドｊは、ステップ６１４で補間値を計算するが、その際、自己充足的に計算できるのは自ストランドの補間値だけなので、その常微分方程式の右辺に現れるy[i]に対応するそれぞれのストランドに、x0+hpの値を渡して、それらのストランドから、計算された補間値を受け取る。それぞれのストランドでは、ステップ６１４で説明した補間の計算式が用いられる。
【００８８】
図７では、そのようにして、対応するそれぞれのストランドで計算された補間値は、それぞれ、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法にも適用可能である。
【００９３】
さらに、補間の計算式は、実施例に示されている５次のエルミート補間式に限定されず、必要な精度を満たすものであるなら、線形補間など、ニュートンの補間式など、任意の補間式を用いることができる。
【００９４】
また、上記実施例は、シングル・プロセッサを例として説明されたが、これには限定されず、本発明の技法は、マルチプロセッサ、マルチコアなどの複数のＣＰＵをもつシステムにでも同様に適用することができる。
【符号の説明】
【００９５】
３１６ ハードディスク・ドライブ
４０２ シミュレーション・モデリング・ツール
４０４ ソース・コード
４０６ 解析モジュール
４０４ ソースコード
４０８ ストランド作成モジュール
４１０ コード生成モジュール
４１２ コンパイラ
４１４ 実行環境
【技術分野】
【０００１】
この発明は、コンピュータを利用したシミュレーション・システム等において使用される常微分方程式を解くための方法、プログラム及びシステムに関し、特に、常微分方程式を解くための計算を高速化し、あるいは計算量を減らすための技法に関するものである。
【背景技術】
【０００２】
従来より、科学技術計算、シミュレーションなどの分野で、コンピュータが利用されている。
【０００３】
最近になって特に盛んに開発されるようになってきたシミュレーションの分野として、ロボット、自動車、飛行機などのメトカトロニクスのプラントのシミュレーション用ソフトウェアがある。電子部品とソフトウェア技術に発展の恩恵により、ロボット、自動車、飛行機などでは、神経のように張り巡らされたワイヤ結線や無線ＬＡＮなどを利用して、大部分の制御が電子的に行われる。
【０００４】
それらは、本来的には機械的装置であるのに、大量の制御ソフトウェアをも内蔵している。そのため、製品の開発に当たっては、制御プログラムの開発とそのテストに、長い時間と、膨大な費用と、多数の人員を費やす必要が出てきた。
【０００５】
このようなテストにために従来行われている技法として、ＨＩＬＳ(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 = 常微分方程式の略である。
【００１７】
xnを時間、
ki(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
ai, bij, ci,c*i (ここで、i,jは添字)を所定の定数とすると、
k1 = hf(xn,yn)
k2 = hf(xn+a2h,yn+b21k1)
k3 = hf(xn+a3h,yn+b31k1+b32k2)
k4 = hf(xn+a4h,yn+b41k1+b42k2+b43k3)
k5 = hf(xn+a5h,yn+b51k1+b52k2+b53k3+b54k4)
k6 = hf(xn+a6h,yn+b61k1+b62k2+b63k3+b64k4+b65k5)
【数１】
【００１８】
ルンゲ・クッタ・フェールベルク法においては、計算精度を所望以内に保つために、上記のようにして計算された誤差Δが、所定の閾値Δ0の範囲内かどうかが判断され、もしそうなら、上記の式により、ステップh0が計算され、それが次のループのステップとなる。
【００１９】
もし、上記のようにして計算された誤差Δが、所定の閾値Δ0の範囲を超えると、ルンゲ・クッタ・フェールベルク法においては、エラーと見なして、もともと使用されたステップhよりも小さいステップに変えて、計算がやり直される。
【００２０】
ところが、y'(t) = f(t,y(t))、y, f ∈ Ｒnという連立常微分方程式においては、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)
【発明の概要】
【発明が解決しようとする課題】
【００２５】
上記従来技術は、コンピュータによる常微分方程式の計算における、刻み量あるいは補間量の計算技法については記述するものの、ルンゲ・クッタ法または、ルンゲ・クッタ・フェールベルク法等に従い、連立常微分方程式を解く際の計算量の低減技法については、示唆するものではない。
【００２６】
従って、本発明の目的は、シミュレーション・システムなどで使用される連立常微分方程式をコンピュータで解く計算処理において、計算量を低減する技法を提供することにある。
【課題を解決するための手段】
【００２７】
本発明は、上記課題を解決するためになされたものであり、以下に記述する処理によって、コンピュータによる連立常微分方程式の計算量を低減させる。
【００２８】
本発明は、上記背景技術で述べた、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法を用いた、コンピュータによる常微分方程式の数値解法を前提とする。すなわち、連立常微分方程式の個々の常微分方程式を解く際に、次数Ｎ（Ｎは所定の整数）の計算項と、次数Ｎ＋１の計算項の差Δを計算して、その差が所定の閾値Δ0よりも小さいかどうか判断し、もしΔ =< Δ0であるなら、Δ0/Δで決まる所定の計算式に従い、ステップ・サイズを決定して、次の計算に進む。ここまでは、従来のルンゲ・クッタ・フェールベルク法による連立常微分方程式の計算方法と同じである。
【００２９】
ところが、次数Ｎの計算項と、次数Ｎ＋１の計算項の差Δが閾値Δ0よりも大きい場合に、本発明の計算処理は異なる。すなわち、従来技術では、連立常微分方程式を構成する常微分方程式のうちの１つでもエラーになると、ステップを戻って、ステップ・サイズを減らし、連立常微分方程式全体の計算をやり直す必要があった。
【００３０】
しかし、本発明においては、連立常微分方程式を構成する常微分方程式のうち、Δ > Δ0となったエラーを生じた常微分方程式を計算するストランドに再計算が命じられる。再計算のため、エラーを生じたストランドには、Δ0/Δで決まるステップサイズがセットされ、その分だけ時間ステップが先に進む。そこから更に、再計算を行っていない他のストランドのステップに追いつくまで時間ステップを進めるためには、再計算を行っているストランドに関連する他のストランドの計算結果の、当該ステップ開始時点での値が必要となる。この値は、再計算を行っていない他のストランドの計算を再実行するのではなく、他のストランドの前回ステップにおける値と、最終ステップにおける値との補間値として計算される。
【００３１】
そのような常微分方程式の補間値の好適な計算アルゴリズムは、上述の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)に記述されている。しかし、本発明で利用可能な補間値の計算アルゴリズムは、これには限定されず、所望の精度を満たす任意の補間値の計算アルゴリズムを使用することができる。
【００３２】
補間値で以って再計算することにより、エラーを生じた常微分方程式を計算するストランドにおいて、エラーを所定の閾値Δ0よりも小さく保ちながら、エラーが十分少なかった他のストランドと足並みが揃うように、時間ステップを進める。その後、連立常微分方程式全体の計算が進められる。
【発明の効果】
【００３３】
この発明によれば、コンピュータによる、ルンゲ・クッタ・フェールベルク法などの埋め込み型ルンゲ・クッタ法に従う連立常微分方程式を計算する際に、閾値よりも大きい誤差を生じたストランドのみ再計算すればよくなるので、計算量が減り、処理を高速化することが可能となる。
【図面の簡単な説明】
【００３４】
【図１】機能ブロックの例を示す図である。
【図２】機能ブロックのストランド作成の例を示す図である。
【図３】ハードウェア構成を示すブロック図である。
【図４】論理構成の機能ブロック図である。
【図５】連立常微分方程式を並列処理で計算するステップを示す概要フローチャートである。
【図６】ＯＤＥ４５のスキームに本発明を適用した処理で常微分方程式を解く処理を示すフローチャートである。
【図７】連立常微分方程式を並列処理で計算する際の、補間値のタイミングを示す図である。
【発明を実施するための形態】
【００３５】
以下、図面を参照して、本発明の一実施例の構成及び処理を説明する。以下の記述では、特に断わらない限り、図面に亘って、同一の要素は同一の符号で参照されるものとする。なお、ここで説明する構成と処理は、一実施例として説明するものであり、本発明の技術的範囲をこの実施例に限定して解釈する意図はないことを理解されたい。
【００３６】
図３を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・バス３０２には、ＣＰＵ３０４と、主記憶（ＲＡＭ）３０６と、ハードディスク・ドライブ（ＨＤＤ）３０８と、キーボード３１０と、マウス３１２と、ディスプレイ３１４が接続されている。ＣＰＵ３０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のＰｅｎｔｉｕｍ（商標）４、インテル社のＣｏｒｅ（商標） ２ ＤＵＯ、ＡＭＤ社のＡｔｈｌｏｎ（商標）などを使用することができる。主記憶３０６は、好適には、１ＧＢ以上の容量、より好ましくは、２ＧＢ以上の容量をもつものである。
【００３７】
ハードディスク・ドライブ３０８には、オペレーティング・システムが、格納されている。オペレーティング・システムは、Ｌｉｎｕｘ（商標）、マイクロソフト社のＷｉｎｄｏｗｓ ７、Ｗｉｎｄｏｗｓ ＸＰ（商標）、Ｗｉｎｄｏｗｓ（商標）２０００、アップルコンピュータのＭａｃ ＯＳ（商標）などの、ＣＰＵ３０４に適合する任意のものでよい。
【００３８】
ハードティスク・ドライブ３０８にはさらに、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=t0での初期値y(t0)を与える。これは実際は、実行させるシミュレーションに応じて、オペレータが事前にセットしておくものである。
【００６０】
ステップ５０４では、管理用の実行コードは、個々のストランドの常微分方程式を解くストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎを一斉にスタートさせ、並列に動作させる。
【００６１】
起動されると、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは各々、本発明の特徴を含むＯＤＥ４５のスキームにより、常微分方程式を解く計算を行う。各々のストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの処理の詳細は、図６のフローチャートで詳しく説明する。
【００６２】
管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎが全て終了するのを待つ。というのは、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎは、計算量が同等とは限らず、また、本発明の特徴によれば、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎのどのストランドも、エラーが検出されたことに応答して再計算を行うことがあるので、その分、遅延する可能性があるからである。
【００６３】
管理用の実行コードは、ステップ５０８で、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎの１ステップ分の計算ステップが終了したのを確認すると、ステップ５１０で、シミュレーションの終了かどうかの判断を行う。
【００６４】
シミュレーションの終了は、オペレータの操作、予定したシナリオの終了、または予定時間経過などによって決定される。管理用の実行コードは、シミュレーションの終了でないと判断すると、ステップ５０４に戻り、次のステップの計算に進む。もしシミュレーションの終了であるなら、管理用の実行コードは、シミュレーションを終了させる。
【００６５】
次に、図６のフローチャートを参照して、ストランド５０６＿１、５０６＿２、５０６＿３、・・・、５０６＿Ｎで個々に実行される計算処理について説明する。
【００６６】
背景技術でも説明したが、この実施例では、ルンゲ・クッタ・フェールベルク法において、ＯＤＥ４５のスキームで、連立常微分方程式をコンピュータによる数値計算により解くものとする。なお、本発明は、ルンゲ・クッタ・フェールベルク法の特定の次数に限定されず、ＯＤＥ３４、ＯＤＥ５６などのその他のスキームも使用可能である。なお、ルンゲ・クッタ・フェールベルク法のより詳しい説明は、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))に着目するものとする。
【００７０】
ここでは、x0 = xnとして、処理を開始する。
xnを時間、
ki(但し、i = 1,...,6)を中間変数、
hをステップ・サイズ（但し、これは定数でなく変数）、
ai, bij, ci,c*i (ここで、i,jは添字)を所定の定数とすると、
k1 = hf[j](xn,yn)
k2 = hf[j](xn+a2h,yn+b21k1)
k3 = hf[j](xn+a3h,yn+b31k1+b32k2)
k4 = hf[j](xn+a4h,yn+b41k1+b42k2+b43k3)
k5 = hf[j](xn+a5h,yn+b51k1+b52k2+b53k3+b54k4)
k6 = hf[j](xn+a6h,yn+b61k1+b62k2+b63k3+b64k4+b65k5)
【数２】
【００７１】
ステップ６０８では、次の式により誤差Δが計算される。
【数３】
【００７２】
ステップ６１０では、計算した誤差Δが、予め定めた閾値Δ0より大きいかどうか、判断される。
【００７３】
ステップ６１０で、Δ > Δ0なら、ステップ６１２で、次の式により、hpを計算する。ただし、計算されたhpとx0を加えた値がxn+1を超える場合には、x0 + hp = xn+1とする。
【数４】
【００７４】
次に、ステップ６１４では、次の式に従い、i = 1, .., Nの５次エルミート補間値U[i]が計算される。
【数５】
但し、
【数６】
【００７５】
ここで、f[i]lは、次の式であらわされる。
【数７】
【００７６】
さらに、
【数８】
【００７７】
ところで、k1〜k6の右辺にあらわれるynは実際は、
yn ≡ yn(y[1],y[2],...,y[N])のように、
y[1],y[2],...,y[N]の関数である。
【００７８】
そこで、数５で決定される補間式Ｕ[i](xn+hp)を、y[i]に置き換えて(但し、i = 1, ,,,, Nの各々に亘る)、ステップ６１６では、x0+hpから出発して、ステップ６０６と同様のＯＤＥ４５の計算が行われる。なお、補間式Ｕ[i](xn+hp)のより詳しい説明については、前述の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は除く)にxn+hpを渡し、ストランドiが補間式Ｕ[i](xn+hp)を計算して、その値を、ストランドｊに送り返してくるという処理により、ストランドｊは、Ｕ[i](xn+hp)(但し、i = 1, ,,,, N)の値を集める。この処理は、図７に関連して、後で説明する。
【００８０】
ステップ６１８では、ステップ６０８と同様の処理で誤差が計算され、ステップ６２０では、ステップ６１０と同様の処理で、誤差と閾値が比較され、誤差が閾値より大きいと判断されると、ステップ６１２に戻る。
【００８１】
ステップ６２０で、誤差が閾値より小さいか閾値と等しいと判断されると、ステップ６２２に進む。
【００８２】
ステップ６２２では、x0+hp < xn+1かどうかが判断され、そうであれれば、ステップ６１２に戻って、hpが再計算され、その再計算されたhpに基づき、ステップ６１４で補間値Ｕが計算され、ステップ６１６でＯＤＥ４５の計算が行われる。
【００８３】
ルンゲ・クッタ・フェールベルク法の誤差式の性質により、ステップ６１０からステップ６１２に分岐した直後は必ず、ステップ６２２での判断はx0+hp < xn+1であり、その後、ステップ６１２、ステップ６１４、ステップ６１６、ステップ６１８、ステップ６２０、ステップ６２２のループが廻る毎にhpの値は次第に増加する。
【００８４】
そうして、x0+hpがxn+1に等しくなると、ステップ６２４に進み、下記の式で次ステップhの計算をする。ここで、式右辺のhは、hpで読み替えるものとする。
【数９】
【００８５】
ステップ６２６では、x0 = xn+1によって時間が進められ、ステップ６２８では、
状態y0が、yn+1で更新される。そして処理は、判断ステップ６０４に戻る。
【００８６】
判断ステップ６１０に戻って、誤差Δが閾値Δ0に等しいか、閾値Δ0より小さい場合は、既に説明したステップ６２４、６２６、及び６２８を経て、判断ステップ６０４に戻る。
【００８７】
図７は、ストランド１〜Ｎの動作を示す概要タイミング図である。ストランド１〜Ｎに対応する実行コードによる計算の結果、ストランドｊだけが、エラー、すなわち、ステップ６１０で誤差Δ > 閾値Δ0となったとする。すると、ストランド１〜ｊ−１と、ストランドｊ＋１〜Ｎが、ステップ６２４、ステップ処理６２６、及びステップ６２８を経てxn+1に到達している一方で、ストランドｊは、ステップ６１２から再計算の処理に入る。ストランドｊは、ステップ６１４で補間値を計算するが、その際、自己充足的に計算できるのは自ストランドの補間値だけなので、その常微分方程式の右辺に現れるy[i]に対応するそれぞれのストランドに、x0+hpの値を渡して、それらのストランドから、計算された補間値を受け取る。それぞれのストランドでは、ステップ６１４で説明した補間の計算式が用いられる。
【００８８】
図７では、そのようにして、対応するそれぞれのストランドで計算された補間値は、それぞれ、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法にも適用可能である。
【００９３】
さらに、補間の計算式は、実施例に示されている５次のエルミート補間式に限定されず、必要な精度を満たすものであるなら、線形補間など、ニュートンの補間式など、任意の補間式を用いることができる。
【００９４】
また、上記実施例は、シングル・プロセッサを例として説明されたが、これには限定されず、本発明の技法は、マルチプロセッサ、マルチコアなどの複数のＣＰＵをもつシステムにでも同様に適用することができる。
【符号の説明】
【００９５】
３１６ ハードディスク・ドライブ
４０２ シミュレーション・モデリング・ツール
４０４ ソース・コード
４０６ 解析モジュール
４０４ ソースコード
４０８ ストランド作成モジュール
４１０ コード生成モジュール
４１２ コンパイラ
４１４ 実行環境
【特許請求の範囲】
【請求項１】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くための方法であって、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップとを有する、
常微分方程式を解くための方法。
【請求項２】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取るステップをさらに有する、
請求項１の方法。
【請求項３】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項２の方法。
【請求項４】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項１の方法。
【請求項５】
前記補間式が、５次エルミート補間式である。請求項３の方法。
【請求項６】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くためのプログラムあって、
前記コンピュータをして、
前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップ、
とを実行させる、
常微分方程式を解くためのプログラム。
【請求項７】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取るステップをさらに有する、
請求項６のプログラム。
【請求項８】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項７のプログラム。
【請求項９】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項６のプログラム。
【請求項１０】
前記補間式が、５次エルミート補間式である。請求項８のプログラム。
【請求項１１】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くためのシステムあって、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップを実行させる手段とを有する、
常微分方程式を解くためのシステム。
【請求項１２】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取る手段をさらに有する、
請求項１１のシステム。
【請求項１３】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項１２のシステム。
【請求項１４】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項１１のシステム。
【請求項１５】
前記補間式が、５次エルミート補間式である。請求項１３のシステム。
【請求項１】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くための方法であって、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算するステップと、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップとを有する、
常微分方程式を解くための方法。
【請求項２】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取るステップをさらに有する、
請求項１の方法。
【請求項３】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項２の方法。
【請求項４】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項１の方法。
【請求項５】
前記補間式が、５次エルミート補間式である。請求項３の方法。
【請求項６】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くためのプログラムあって、
前記コンピュータをして、
前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算するステップと、
前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップ、
とを実行させる、
常微分方程式を解くためのプログラム。
【請求項７】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取るステップをさらに有する、
請求項６のプログラム。
【請求項８】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項７のプログラム。
【請求項９】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項６のプログラム。
【請求項１０】
前記補間式が、５次エルミート補間式である。請求項８のプログラム。
【請求項１１】
コンピュータの処理により、埋め込み型ルンゲ・クッタ法を用いて連立常微分方程式を解くためのシステムあって、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、次数Ｎ（Ｎは所定の整数）の状態値と次数Ｎ＋１の状態値の間の誤差を値を計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記誤差が所定の閾値を超えたことに応答して、前記誤差と前記閾値とから、ステップ・サイズを計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズに対応する補間値を計算する手段と、
前記コンピュータの処理により、前記連立常微分方程式の各常微分方程式毎に、前記計算されたステップ・サイズと、前記補間値を用いて、再計算するステップを実行させるステップを実行させる手段とを有する、
常微分方程式を解くためのシステム。
【請求項１２】
前記前記誤差が所定の閾値を超えたことに応答して、前記計算されたステップ・サイズを以って、前記連立常微分方程式の他の常微分方程式の計算式から前記補間値の計算値を受け取る手段をさらに有する、
請求項１１のシステム。
【請求項１３】
前記埋め込み型ルンゲ・クッタ法が、ＯＤＥ４５である、請求項１２のシステム。
【請求項１４】
前記埋め込み型ルンゲ・クッタ法が、ルンゲ・クッタ・フェールベルク法である、請求項１１のシステム。
【請求項１５】
前記補間式が、５次エルミート補間式である。請求項１３のシステム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【公開番号】特開２０１１−１８６９９１（Ｐ２０１１−１８６９９１Ａ）
【公開日】平成２３年９月２２日（２０１１．９．２２）
【国際特許分類】
【出願番号】特願２０１０−５４２５１（Ｐ２０１０−５４２５１）
【出願日】平成２２年３月１１日（２０１０．３．１１）
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
【公開日】平成２３年９月２２日（２０１１．９．２２）
【国際特許分類】
【出願日】平成２２年３月１１日（２０１０．３．１１）
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
[ Back to top ]