JP3735128B2

JP3735128B2 - シミュレーション装置

Info

Publication number: JP3735128B2
Application number: JP28152493A
Authority: JP
Inventors: 純司金児
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 1993-11-10
Filing date: 1993-11-10
Publication date: 2006-01-18
Anticipated expiration: 2021-01-18
Also published as: JPH07134704A; US5646869A

Description

【０００１】
【産業上の利用分野】
本発明は不規則外乱を受ける系の数値シミュレーションを行うシミュレーション装置に関する。
【０００２】
【従来の技術】
自然科学の分野において、一般に物理系を記述する数学モデルの解を解析的に求めるのは難しく、その数学モデルの初期値問題は数値処理により、近似的に解かれる場合が多い。近年の情報処理装置の発達を背景にして、こうした数値処理を行う種々のシミュレーション装置が考案されているが、シミュレーション精度を上げるためには膨大な計算量が必要となり、汎用性が損なわれることがある。そこでシミュレーションの信頼性と汎用性を考慮しつつ、いかにシミュレーション装置を構築するかが重要な問題である。
【０００３】
物理系として、工業プラントやロボットなどの挙動を扱う工学の分野においては、センサからの観測データをもとに物理系の挙動を推定し、この物理系を制御することが重要な課題である。この推定処理は一般にシミュレーション装置を用いた数値シミュレーションにより行われる。
【０００４】
シミュレーション装置の性能は数値シミュレーションの結果の信頼性を左右し、この結果を用いて物理系を制御する制御装置の性能を決定付けることになる。このため、数値シミュレーションの処理対象である物理系や、制御装置に要求される信頼性に応じて、多種多様なシミュレーション装置が開発されている。
【０００５】
処理対象となる物理系のモデルは自然法則に基づいて導出されるが、その際に、現実の多くの物理系に混在している雑音などの不規則外乱の影響が無視され、理想化あるいは単純化が施されている場合が多い。しかし対象とする物理系によっては、不規則外乱の影響が無視できないものもある。
【０００６】
例えば人工衛星の制御装置を設計する際には、太陽の放射圧、重力、宇宙空間に浮遊する他の衛星の残骸との衝突などの不規則外乱を無視することができない。また工業プラントやロボットなどの制御装置においては、制御対象を観測するセンサからの信号に混入する雑音などの不規則外乱が無視できないことがある。
【０００７】
このように、不規則外乱の影響を受ける物理系の解析や制御を目的とするシミュレーションにおいては、不規則外乱を反映した物理モデルを与え、これを数値処理する必要がある。
【０００８】
一般に任意の物理系のモデルは次のような常微分方程式により表現できる。
【０００９】
【数１】

【００１０】
ここで、ｘ＝ｘ（ｔ）は物理系における変位や温度等の、シミュレーションが対象とする物理量であり、時刻ｔに依存して変動する。外１はｘのｔに関する
【００１１】
【外１】

【００１２】
微分を表す。ｆ（ｘ，ｔ）はｘとｔの関数であり、ｘ₀は時刻ｔ₀において与えられるｘの初期値である。
（１）式の物理モデルに対する数値シミュレーションにおいては、時刻ｔ_n＝ｔ₀＋ｎｈにおける物理量ｘ（ｔ_n）の近似値をｘ_nとし、次式により次の時刻ｔ_n+1における物理量ｘ（ｔ_n+1）の近似値ｘ_n+1を計算する。
【００１３】
【数２】

【００１４】
ただし、ｘ₀は近似値ではなく、時刻ｔ₀におけるｘの初期値である。ここでｈはシミュレーションの時間ステップの刻み幅であり、Ψ（ｘ_n，ｔ_n，ｈ）は時刻ｔ_nにおける近似値ｘ_nの時間増分を与えている。このΨ（ｘ_n，ｔ_n，ｈ）をｘ_n、ｔ_n、ｈを用いてどのように構築するかが重要であり、数値シミュレーションの結果の信頼性を左右することになる。
【００１５】
（２）式のΨ（ｘ_n，ｔ_n，ｈ）の構築方法はテーラ級数項の方法とルンゲ−クッタ法の２種類に大別される。このうちテーラ級数項の方法は原理的には単純で、Ψ（ｘ_n，ｔ_n，ｈ）をシミュレーションに必要な次数までのテーラ級数項により与える方法である。しかしこの方法では一般に、ｐ次の項の係数に関数ｆ（ｘ，ｔ）のｐ次の導関数∂^pｆ（ｘ，ｔ）／∂ｘ^pが含まれる。関数ｆ（ｘ，ｔ）がｘの線型関数等の簡単な関数形で与えられる場合は、導関数は簡単に求められるが、ｆ（ｘ，ｔ）が複雑な非線型関数の場合、その高階の導関数をｘ、ｔの関数として与えることは困難である。
【００１６】
またこのテーラ級数項の方法ではオペレータがｆ（ｘ，ｔ）の導関数を求めて、シミュレーション装置に入力しなければならず、さらに導関数の値の計算はシミュレーション装置にとって多大な負担となるので汎用性に乏しい。
【００１７】
そこで数値シミュレーションでは、もう一つの方法として、次式によりΨ（ｘ_n，ｔ_n，ｈ）を構築するルンゲ−クッタ法が多く用いられている。
【００１８】
【数３】

【００１９】
ここでｓは自然数であり、ｋパラメータｋ_iは適切に選ばれたｘ、ｔに対応する関数ｆ（ｘ，ｔ）の値で、次式のような計算処理により生成される。
【００２０】
【数４】

【００２１】
（３）、（４）式でｂパラメータｂ_iやａパラメータａ_ijはｋパラメータに対する重みを表すパラメータであり、シミュレーション装置内のメモリに予め格納されている。ｂパラメータとａパラメータの選び方によって近似値ｘ_nに対する増分Ψ（ｘ_n，ｔ_n，ｈ）の値が異なり、近似値ｘ_n+1の精度が左右されるため、結果としてシミュレーション装置の性能が決められる。
【００２２】
このｋパラメータを用いた方法は、（２）式の数値処理にあたって係数関数ｆ（ｘ，ｔ）の微分形を与える必要がないため、テーラ級数の方法よりも遙かに汎用的である。
【００２３】
時刻ｔ_nにおける（１）式の解をｘ（ｔ_n）とすると、数値シミュレーションによる近似値ｘ_nのｘ（ｔ_n）に対する誤差はｘ（ｔ_n）−ｘ_nで与えられる。シミュレーション結果の信頼性は誤差ｘ（ｔ_n）−ｘ_nが小さいほど高くなる。この誤差の大きさを評価する尺度として、刻み幅ｈが用いられる。誤差がｈ^r+1程度であることが保証されるとき、そのシミュレーション装置あるいはシミュレーションの結果はｒ次の近似精度を持つといわれ、ｒの値が大きいほど、信頼性の高いシミュレーション結果が得られる。
【００２４】
対象とする物理系が雑音等の不規則外乱の影響を受けている場合は、この系の物理モデルは、（１）式の代わりに次のような伊藤型確率微分方程式を用いて記述される。
【００２５】
【数５】

【００２６】
ここで、ｘ＝ｘ（ｔ）は対象とする物理系における状態量であり、ｆ（ｘ）、ｇ（ｘ）はｘの関数である。またｗ（ｔ）は雑音モデルを表し、その値はブラウン運動モデルの確率過程に依存する。このモデルは（１）式の右辺に、不規則外乱による第２項を付加した形になっている。ｘ₀は時刻ｔ₀において与えられるｘの初期値である。
【００２７】
不規則外乱を受ける物理モデルを表す（５）式は、（１）式と異なって不規則外乱による付加項があるため、その数値シミュレーションを行うためには、（３）、（４）式のΨ（ｘ_n，ｔ_n，ｈ）をそのまま用いることはできない。この場合には、確率過程に依存する不規則外乱の項を扱うために、（３）、（４）式のΨ（ｘ_n，ｔ_n，ｈ）を確率微分方程式の場合に拡張する必要がある。
【００２８】
この拡張された時間増分Ψ（ｘ_n，ｔ_n，ｈ）をｋパラメータを用いて構築する従来の方法として、確率的Heun法、１．５次精度を保証する３段陽的解法、漸近的に有効な解法がある。
【００２９】
確率的Heun法（RUMELIN, SIAM J. Vol.19, No.3, June 1982）は（５）式の物理モデルに対するΨ（ｘ_n，ｔ_n，ｈ）として、（３）式の右辺に確率過程に依存する項を付加したものを用いている。この方法は確率過程を扱う不規則外乱を受ける物理モデルに対する（３）式の素直な拡張であり、１次の近似精度を保証している。また（３）式と異なってｋパラメータが関数ｇ（ｘ）とｇ（ｘ）の１階の導関数の値を含むことを特徴とする。
【００３０】
１．５次精度を保証する３段陽的解法（斉藤、三井日本応用数理学会論文誌 Vol.2, No.1, 1992）は２種類の確率過程を用いて、確率的Heun法のΨ（ｘ_n，ｔ_n，ｈ）をさらに拡張している。この方法は一般に１．５次の高い近似精度を保証しており、ｋパラメータは関数ｆ（ｘ）とｇ（ｘ）のそれぞれの１階の導関数ならびに、ｇ（ｘ）の２階の導関数の値を含んでいる。
【００３１】
漸近的に有効な解法（NEWTON, SIAM J. Vol.51, No.2, April 1991）は確率的Heun法のΨ（ｘ_n，ｔ_n，ｈ）にさらに（ｈ）^1/2に関する項を付加して、ｋパラメータの計算処理を簡単化している。この方法では、ｋパラメータが関数ｇ（ｘ）の１階の導関数の値を含まず、関数ｆ（ｘ）とｇ（ｘ）の値を用いて計算できるので、汎用性に優れたシミュレーション方法を与えている。この方法により保証される近似精度は確率的Heun法と同様に１次である。
【００３２】
【発明が解決しようとする課題】
しかしながら上述のような従来の時間増分の構築方法による数値シミュレーションには以下のような問題がある。
【００３３】
確率的Heun法では、近似値ｘ_nに対する時間増分Ψ（ｘ_n，ｔ_n，ｈ）を計算する時、既に計算されているｘ_nの値とｔ_n、ｈを用いて２組のｋパラメータを計算し、ｋパラメータに適当な重みを付けて加算する。このとき一方の組のｋパラメータの計算式が、不規則外乱を記述している関数ｇ（ｘ）の１階の導関数を含むため、数値シミュレーションを行う際に、関数ｇ（ｘ）の１階の導関数の値を計算しなければならない。このためテーラ級数項の方法と同様に計算量が増大して、シミュレーション装置の負荷が大きくなるという欠点をもち、汎用的でない。また保証されるシミュレーションの近似精度は１次に留まっており、必ずしもその信頼性は高くない。
【００３４】
１．５次精度を保証する３段陽的解法は、確率的Heun法よりもシミュレーション結果の信頼性は高いが、ｋパラメータの計算式が関数ｆ（ｘ）とｇ（ｘ）の１階の導関数ならびに、ｇ（ｘ）の２階の導関数を含むため、これらの導関数の関数形を与える必要がある。関数ｇ（ｘ）だけでなく、関数ｆ（ｘ）の導関数の値が必要な上、さらに関数ｇ（ｘ）の２階の導関数の値も必要とするため、確率的Heun法よりも著しく計算量が増大し、実際にはシミュレーション装置として実用化するのは困難である。
【００３５】
漸近的に有効な解法は、確率的Heun法と１．５次精度を保証する３段陽的解法と比べて、数値シミュレーションを行う際に、関数ｇ（ｘ）等の導関数を与える必要がないという利点を有する。この方法の場合、ｋパラメータの計算式は関数ｆ（ｘ）とｇ（ｘ）により記述されており、数値処理が容易で、汎用性に優れている。しかし保証されるシミュレーションの近似精度は１次に留まっており、その信頼性は低いという問題がある。
【００３６】
不規則外乱の影響を受ける系のシミュレーションを行うには、（５）式に示される確率微分方程式の解を数値的に求めるシミュレーション装置の開発が不可欠である。それにも係わらず、従来のシミュレーション方法が信頼性と汎用性の面で上述のような問題を抱えているため、実用的なシミュレーション装置はほとんど開発されていない。
【００３７】
本発明は、不規則外乱を受ける系の数値シミュレーションにおいて、計算量が少なくて汎用性に優れ、かつ、高い近似精度を持つシミュレーション結果が得られるシミュレーション装置を提供することを目的とする。
【００３８】
【課題を解決するための手段】
図１は本発明の原理図である。
本発明は、不規則外乱を受ける物理系を、状態量の時間変化を表す関数と乱数によって記述し、時間刻み毎に、状態量の初期値から、状態量の近似値を陽的に求め、シミュレーション結果を出力する情報処理装置におけるシミュレーション装置１０である。
【００３９】
本シミュレーション装置１０は、初期設定入力手段１１と、増分生成手段１２、近似値更新手段１３、出力手段１４を有する。
初期設定入力手段１１は、初期データである上記関数と初期値を増分生成手段１２に入力する。
【００４０】
増分生成手段１２は、上記乱数に対応する第１の乱数と第２の乱数を発生し、これらの乱数と時間刻みを用いて、初期設定入力手段１１より与えられる関数の関数値として、３種類以上の異なるｋパラメータを状態量の初期値から生成する。そして、時間刻みの平方根に比例し、かつ、第１の乱数と第２の乱数に依存しない項を用いて、生成したｋパラメータからルンゲ−クッタ型の状態量の時間増分を生成する。
【００４１】
近似値更新手段１３は、増分生成手段１２が生成した時間増分を用いて状態量の近似値を更新し、更新した近似値を出力手段１４と増分生成手段１２に与える。
【００４２】
増分生成手段１２は、上記の処理に従い、さらに近似値更新手段１３から与えられる更新された近似値に対する時間増分を生成し、近似値更新手段１３に与える。
【００４３】
出力手段１４は、近似値更新手段１３から与えられる時間刻み毎の状態量の各近似値を、シミュレーション結果として外部に出力する。
上記増分生成手段１２は、さらに数値処理を簡単化するために、不図示のパラメータ生成手段と時間増分生成手段、第１のパラメータ格納手段、第２のパラメータ格納手段を備えるものとする。
【００４４】
第１のパラメータ格納手段は、パラメータ生成手段が用いる第１のパラメータを格納し、第２のパラメータ格納手段は、時間増分生成手段が用いる第２のパラメータを格納する。
【００４５】
パラメータ生成手段は、時間刻み、第１の乱数、第２の乱数のうちの一つと、第１のパラメータおよび状態量の近似値とを用いて、状態量に対応する標本状態量を生成し、与えられた関数を用いて、生成した標本状態量から上記ｋパラメータを生成し、時間増分生成手段に与える。
【００４６】
時間増分生成手段は、パラメータ生成手段から与えられるｋパラメータと、第１の乱数、第２の乱数、時間刻みおよび第２のパラメータとを用いて、状態量の時間増分を生成する。
【００４７】
【作用】
増分生成手段１２が、２種類の乱数を用いて時間増分を生成するので、時間増分を時間刻みの巾について展開した時に、１種類の乱数のみを用いて時間増分を生成するシミュレーション装置よりも、高い次数の項を含むことが可能になる。
【００４８】
また時間増分が、時間刻みの平方根に比例し、かつ、上記２種類の乱数の値に依存しない項を含むため、（５）式における関数ｆ（ｘ）やｇ（ｘ）の導関数の値を計算する必要がなくなり、このような項を用いないシミュレーション装置よりも数値処理が容易になる。
【００４９】
また状態量の近似値から標本状態量を生成し、生成した標本状態量からｋパラメータを生成するパラメータ生成手段と、ｋパラメータから時間増分を生成する時間増分生成手段を備えることにより、数値処理を段階的に行うことができる。
【００５０】
【実施例】
以下、図面を参照しながら、本発明の実施例を説明する。
本実施例においては、不規則外乱の影響を受ける系として、（５）式の形で記述される物理モデルのシミュレーションを行う場合を考える。ここで状態量ｘ（ｔ）は一般にｄ次元の状態量であり、（５）式はｄ次元の確率微分方程式を記述する数学モデルに対応する。
【００５１】
本実施例のシミュレーション装置では、時刻ｔ_n＝ｔ₀＋ｎｈにおける状態量ｘ（ｔ）の近似値外２（以後ｘ_nのバーと記す。）と、時刻ｔ_n+1＝ｔ₀＋
【００５２】
【外２】

【００５３】
（ｎ＋１）ｈにおける状態量ｘ（ｔ）の近似値外３（以後ｘ_n+1のバーと記
【００５４】
【外３】

【００５５】
す。）の関係を次式により記述する。
【００５６】
【数６】

【００５７】
Ψ（ｘ_nのバー，ｔ_n，ｈ）は時刻ｔ_nにおける状態量ｘ（ｔ）の近似値ｘ_nのバーに対する時間増分を表している。近似値ｘ_nのバー、ｘ_n+1のバー、時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）もまたｄ次元の値を持つ。
【００５８】
図２は本実施例のシミュレーション装置の構成図である。図２において初期設定入力部２１は図１の初期設定入力手段１１に相当し、外部から与えられる初期データである（５）式における関数ｆ（ｘ）とｇ（ｘ）の具体的な関数形と、時刻ｔ₀における状態量ｘの初期値ｘ₀＝ｘ（ｔ₀）とを初期設定としてｋパラメータ生成装置２２に入力する。
【００５９】
ｋパラメータ生成装置２２、時間増分生成装置２４、メモリ２３─１〜２３─４、メモリ２５─１〜２５─４は、それぞれパラメータ生成手段、時間増分生成手段、第１のパラメータ格納手段、第２のパラメータ格納手段に相当する。そして、これらは乱数発生装置２６─１、２６─２と共に、図１の増分生成手段１２を形成している。
【００６０】
乱数発生装置２６─１、２６─２はそれぞれ、第１の乱数、第２の乱数を発生する。
ｋパラメータ生成装置２２と時間増分生成装置２４は、初期設定入力手段１１または近似値更新装置２７からの入力に基づいて、時刻ｔ_nにおける状態量ｘ（ｔ）の近似値ｘ_nのバーに対する時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を生成し、近似値更新装置２７に入力する。
【００６１】
近似値更新装置２７は図１の近似値更新手段１３に相当し、時間増分生成装置２４から入力される時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を近似値ｘ_nのバーに加算してこれをｘ_n+1のバーとし、ｋパラメータ生成装置２２と出力部２８に与える。
【００６２】
出力部２８は図１の出力手段１４に相当し、近似値更新装置２７から逐次与えられる状態量の近似値をまとめてシミュレーション結果とし、外部に出力する。図２のシミュレーション装置では、（５）式の右辺第２項の不規則外乱を記述する項ｇ（ｘ）ｗ（ｔ）を数値処理して、近似値ｘ_nのバーに反映させるため、４種類のｋパラメータを導入し、時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を次式に基づいて生成する。
【００６３】
【数７】

【００６４】
（７）式の右辺第１項は常微分方程式に関するルンゲ−クッタ法の場合に類似した項であり、第３、第４項は（５）式の解を１．５次以上の近似精度で求めるために必要な項である。ここでパラメータΔＷ_n+1、外４（以後ΔＷ_n+1の
【００６５】
【外４】

【００６６】
波と記す。）はそれぞれ、時刻ｔ_nにおける平均０、分散１の正規乱数のサンプルξ_n+1、外５（以後ξ_n+1の波と記す。）を用いてΔＷ_n+1＝ξ_n+1（ｈ
【００６７】
【外５】

【００６８】
）^1/2、ΔＷ_n+1の波＝ξ_n+1の波（ｈ）^1/2のように与えられる。サンプルξ_n+1、ξ_n+1の波は、それぞれ第１の乱数、第２の乱数に相当し、近似値ｘ_nのバーを更新する度に、乱数発生装置２６─１、２６─２より独立に生成される。
【００６９】
（７）式の右辺第４項は、サンプルξ_n+1、ξ_n+1の波の値に依存せず、（ｈ）^1/2に比例する。第３、第４項があるために、関数ｆやｇの導関数の値を用いなくても、高い近似精度を実現できる。従ってこれらの項は、数値処理を容易にし、シミュレーション装置の汎用性を高めるために必要である。ここでνは適切な正の実数である。
【００７０】
（７）式においてｂパラメータｂ_i、外６（以後ｂ_iのバーと記す。）、
【００７１】
【外６】

【００７２】
外７（以後ｂ_iの波と記す。）、外８（以後ｂ_iの山と記す。）、（ｉ
【００７３】
【外７】

【００７４】
【外８】

【００７５】
＝１，・・・，ｓ）はそれぞれ、時間増分を生成する際のｋパラメータｋ_i、外９（以後ｋ_iのバーと記す。）、外１０（以後ｋ_iの波と記す。）、
【００７６】
【外９】

【００７７】
【外１０】

【００７８】
外１１（以後ｋ_iの山と記す。）、（ｉ＝１，・・・，ｓ）に対する重み付け
【００７９】
【外１１】

【００８０】
のパラメータで、予めメモリ２５─１、２５─２、２５─３、２５─４、に格納されている。これらのｂパラメータは第２のパラメータに相当する。
ｋパラメータｋ_i（ｉ＝１，・・・，ｓ）は近似値ｘ_nのバーに対して、（８）式により与えられる関数ｆ（ｘ）の値で、ｋパラメータｋ_iのバー、ｋ_iの波、ｋ_iの山（ｉ＝１，・・・，ｓ）はそれぞれ、（９）、（１０）、（１１）式により与えられる関数ｇ（ｘ）の値である。これらの４種類のｋパラメータはｋパラメータ生成装置２２により生成される。
【００８１】
（８）、（９）、（１０）、（１１）式におけるａパラメータａ_ij、外１２
【００８２】
【外１２】

【００８３】
（以後ａ_ijのバーと記す。）、外１３（以後ａ_ijの波と記す。）、外１
【００８４】
【外１３】

【００８５】
４（以後ａ_ijの山と記す。）、（ｉ＝１，・・・，ｓ）はそれぞれ、ｋパラメ
【００８６】
【外１４】

【００８７】
ータを生成する際のｋパラメータｋ_i、ｋ_iのバー、ｋ_iの波、ｋ_iの山（ｉ＝１，・・・，ｓ）に対する重み付けのパラメータで、予めメモリ２３─１、２３─２、２３─３、２３─４、に格納されている。これらのａパラメータは第１のパラメータに相当する。
【００８８】
次にこれらのａパラメータとｂパラメータの決定方法について説明する。
一般にルンゲ−クッタ法による微分方程式の近似解ｘ_nのバーに対して、ｒ次の近似精度を保証するためには、時刻ｔ_n+1における解ｘ（ｔ_n+1）を、時刻ｔ_nにおける解ｘ（ｔ_n）のまわりでｈに関してテーラ展開し、ｈのｒ次の次数の項までが一致するように時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を決める必要がある。
【００８９】
本シミュレーション装置においては、（５）式の近似解ｘ_nのバーに対して、１．５次の近似精度を保証するために、ｈに関して１．５次の項までのテーラ展開を用いることにする。時刻ｔ_n+1における（５）式の解ｘ（ｔ_n+1）のｘ（ｔ_n）のまわりでのテーラ展開は、次式で記述されることが知られている（Ito, Applied Math. and Optimization, Vol.1, 1975）。
【００９０】
【数８】

【００９１】
（１２）式におけるＦとＧはｄ次元の微分作用素で、次式により定義される。
【００９２】
【数９】

【００９３】
また（１２）式のΔＷ_n+1、ΔＷ_n+1の波はそれぞれ独立なブラウン運動モデルの増分を表し、（７）式のΔＷ_n+1、ΔＷ_n+1の波と同様に正規乱数のサンプルξ_n+1、ξ_n+1の波を用いてΔＷ_n+1＝ξ_n+1（ｈ）^1/2、ΔＷ_n+1の波＝ξ_n+1の波（ｈ）^1/2のように定められる。Ｏ（ｈ²）はｈに関して２次以上の次数の項を表す。
【００９４】
ここで時刻ｔ_n+1における近似値ｘ_n+1のバーを与える（６）式の右辺を、近似値ｘ_nのバーのまわりでｈに関してテーラ展開し、ｈに関して１．５次の項までの各項の係数を（１２）の各項の係数と比較すれば、ａパラメータとｂパラメータに対する拘束条件を得ることができる。このとき、（７）式右辺の第３、第４項があるために、関数ｆ、ｇの導関数を用いることなく、（１２）式右辺のＯ（ｈ²）以外の項を得ることができる。こうして求められた拘束条件を次に示す。
【００９５】
【数１０】

【００９６】
【数１１】

【００９７】
ただし、ｃ_i、外１５（以後ｃ_iのバーと記す。）、外１６（以後ｃ_i
【００９８】
【外１５】

【００９９】
【外１６】

【０１００】
の波と記す。）、外１７（以後ｃ_iの山と記す。）は各ａパラメータを用い
【０１０１】
【外１７】

【０１０２】
て次式により定義する。
【０１０３】
【数１２】

【０１０４】
ａパラメータとｂパラメータに対する拘束条件が（１４．１）〜（１４．２２）式の２２個であることを考慮して、ｓ＝４と選べば、シミュレーションに必要なａパラメータとｂパラメータの総数は４０個になる。ここで、簡単のためｂ_i＝ｂ_iのバー、ａ_ij＝ａ_ijのバー、ｃ_i＝ｃ_iの波、ｂ₃の波＝ｂ₃の山、ｂ₄の波＝ｂ₄の山と置くと、パラメータの総数は２５個になる。またこのとき、（１４．１）、（１４．４）の２式が同値となり、（１４．２）、（１４．９）、（１４．１０）、（１４．１１）の４式が同値となり、（１４．１９）、（１４．２１）の２式が同値となる。この結果２２個の拘束条件は１７個に縮退する。従ってａパラメータとｂパラメータを決定するにあたっては、８個のパラメータを自由パラメータとして選び、さらにνの値を適当に与える必要がある。
【０１０５】
一例としてａ₃₁の波＝ａ₄₁の波＝ａ₄₂の波＝０、ａ₃₂＝１／４、ａ₄₃の山＝１／６、ｃ₂＝１／２、ｃ₄＝１、ｃ₄の山＝０、ν＝３と選び、これらを１７個の拘束条件を表す代数方程式に代入して、他の独立パラメータを求めると、各メモリに格納される４０個のａパラメータ、ｂパラメータが次のように決定される。
【０１０６】
【数１３】

【０１０７】
【数１４】

【０１０８】
（１６．１）式のｂパラメータの値およびν＝３を（７）式に代入すると、時間増分を与える式は次式のようになる。
【０１０９】
【数１５】

【０１１０】
また（１６．２）式のａパラメータの値およびν＝３を（８）〜（１１）式に代入すると、ｋパラメータを与える式は次式のようになる。
【０１１１】
【数１６】

【０１１２】
【数１７】

【０１１３】
（６）、（１７）、（１８．１）、（１８．２）式により生成される近似値ｘ_n+1のバーは、ｈに関して１．５次の次数の項まで、（１２）式の右辺のテーラ展開と一致するため、１．５次の近似精度が保証されている。また関数ｆやｇの導関数を必要としないため、オペレータがこれらの導関数の関数形を求める必要がなく、近似値の数値計算における計算量も少なくて済む。
【０１１４】
図３は図２のシミュレーション装置の動作の概要を示したフローチャートである。
シミュレーションの考察対象の物理系が与えられると、オペレータは前処理Ｐ１において対象とする物理系を記述する物理モデルとその初期状態を決める。この物理モデルと初期状態は（５）式に対応している。そしてここで決めた関数ｆ（ｘ）、ｇ（ｘ）の関数形と、状態量ｘの初期値ｘ₀＝ｘ（ｔ₀）を初期データとしてシミュレーション装置に入力する。
【０１１５】
次にシミュレーション装置の内部では、初期設定入力部２１が与えられた初期データをｋパラメータ生成装置２２に入力する。ただし、このときｘ₀のバー＝ｘ₀と設定する。
【０１１６】
ステップＱ１でｋパラメータ生成装置２２は、関数ｆ（ｘ）、ｇ（ｘ）の関数形を格納した後、メモリ２３─１〜２３─４に格納されたａパラメータと乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルξ_n+1、ξ_n+1の波を用いて、初期値ｘ₀より各ｋパラメータを逐次生成する。このｋパラメータ生成処理はｋパラメータ生成装置２２内に格納された（８）〜（１１）式を用いて行う。ただし、このときｎ＝０、ｓ＝４とする。
【０１１７】
ステップＱ２で時間増分生成装置２４は、メモリ２５─１〜２５─４に格納されたｂパラメータと乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルξ_n+1、ξ_n+1の波を用いて、生成されたｋパラメータから時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を生成する。この時間増分生成処理は時間増分生成装置２４内に格納された（７）式を用いて行う。ただし、このときｎ＝０、ｓ＝４、ν＝３とする。
【０１１８】
ステップＱ３で近似値更新装置２７は、格納している（６）式を用いて、生成された時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）とｘ_nのバーから、次のシミュレーションの参照時刻ｔ_n+1における近似値ｘ_n+1のバーを計算し、これをｋパラメータ生成装置２２と出力部２８に出力する。ただし、このときｎ＝０とする。
【０１１９】
ステップＱ４で近似値更新装置２７は、ｎがシミュレーションの参照終了時刻に対応するかどうかを判定する。
まだ参照終了点の近似値が得られていない場合はｎ＝ｎ＋１とし、ステップＱ１に戻って新たにｋパラメータを生成し、ステップＱ２で時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）を生成し、ステップＳ３で次の参照点の近似値ｘ_n+1のバーを求める。以後ステップＱ１〜Ｑ４を繰り返し、ステップＱ４で参照終了時刻に対応した最終参照点に到達したことを確認すると処理を終える。
【０１２０】
図４〜７は図３のステップＱ１のｋパラメータ生成装置２２の詳細な動作を示すフローチャートである。図５〜７におけるｘ_n ^(m)、外１８（以後ｘ_n ^(m
【０１２１】
【外１８】

【０１２２】
⁾のバーと記す。）、外１９（以後ｘ_n ^(m)の波と記す。）、外２０（
【０１２３】
【外１９】

【０１２４】
【外２０】

【０１２５】
以後ｘ_n ^(m)の山と記す。）、（ｍ≧１）はそれぞれ、ｉ≧２の場合の（８）、（９）、（１０）、（１１）式において、関数ｆまたはｇによりｋパラメータに変換される（）内の数値を表し、標本状態量に相当する。ただし、ｍ＝ｉ−１とする。ｉ＝１の場合はｘ_nのバーが標本状態量に相当する。
【０１２６】
図４のステップＳ１においてまず、初期設定入力部２１または近似値更新装置２７から与えられたｘ_nのバーを関数ｆと関数ｇを用いて変換し、パラメータｋ₁、ｋ₁のバー、ｋ₁の波、ｋ₁の山を生成する。ただしｋ₁のバー、ｋ₁の波、ｋ₁の山は同じ値なので、これらの値の生成処理は一回でよい。
【０１２７】
次に図５のステップＳ２において、ステップＳ１で生成されたｋ₁、ｋ₁のバー、ｋ₁の波、ｋ₁の山の値から、パラメータｋ₂、ｋ₂のバー、ｋ₂の波、ｋ₂の山を生成する。
【０１２８】
パラメータｋ₂を生成するには、まずｋ₁、ｋ₁のバー、ｋ₁の波の値にメモリ２３─１、２３─２、２３─３に格納されているａ₁₁、ａ₁₁のバー、ａ₁₁の波の値を乗算する。さらに乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルを用いて（８）式の（）内の第２、第３、第４項を生成し、それらをｘ_nのバーに順次加算して、ｘ_n ⁽¹⁾を生成する。次に生成されたｘ_n ⁽¹⁾の値を関数ｆを用いて変換し、ｋ₂の値を得る。
【０１２９】
パラメータｋ₂のバーを生成するには、まずｋ₁、ｋ₁のバー、ｋ₁の山の値にメモリ２３─１、２３─２、２３─４に格納されているａ₁₁、ａ₁₁のバー、ａ₁₁の山の値を乗算する。さらに乱数発生装置２６─１が発生する正規乱数のサンプルを用いて（９）式の（）内の第２、第３、第４項を生成し、それらをｘ_nのバーに順次加算して、ｘ_n ⁽¹⁾のバーを生成する。次に生成されたｘ_n ⁽¹⁾のバーの値を関数ｇを用いて変換し、ｋ₂のバーの値を得る。
【０１３０】
パラメータｋ₂の波を生成するには、まずｋ₁、ｋ₁の山の値にメモリ２３─１、２３─４に格納されているａ₁₁、ａ₁₁の山の値を乗算して（１０）式の（）内の第２、第３項を生成し、それらをｘ_nのバーに順次加算して、ｘ_n ⁽¹⁾の波を生成する。次に生成されたｘ_n ⁽¹⁾の波の値を関数ｇを用いて変換し、ｋ₂の波の値を得る。
【０１３１】
パラメータｋ₂の山を生成するには、まずｋ₁の山の値にメモリ２３─４に格納されているａ₁₁の山の値を乗算して（１１）式の（）内の第２項を生成し、ｘ_nのバーに加算して、ｘ_n ⁽¹⁾の山を生成する。次に生成されたｘ_n ⁽¹⁾の山の値を関数ｇを用いて変換し、ｋ₂の山の値を得る。
【０１３２】
図６、図７のパラメータｋ₃、ｋ₃のバー、ｋ₃の波、ｋ₃の山およびｋ₄、ｋ₄のバー、ｋ₄の波、ｋ₄の山の生成処理についても、基本的に図５の処理と同様である。
【０１３３】
図６のステップＳ３では、ステップＳ１で生成されたｋ₁、ｋ₁のバー、ｋ₁の波、ｋ₁の山の値と、ステップＳ２で生成されたｋ₂、ｋ₂のバー、ｋ₂の波、ｋ₂の山の値から、ｘ_n ⁽²⁾、ｘ_n ⁽²⁾のバー、ｘ_n ⁽²⁾の波、ｘ_n ⁽²⁾の山の値を求め、パラメータｋ₃、ｋ₃のバー、ｋ₃の波、ｋ₃の山を生成する。
【０１３４】
図７のステップＳ４では、ステップＳ１で生成されたｋ₁、ｋ₁のバー、ｋ₁の波、ｋ₁の山の値と、ステップＳ２で生成されたｋ₂、ｋ₂のバー、ｋ₂の波、ｋ₂の山の値と、ステップＳ３で生成されたｋ₃、ｋ₃のバー、ｋ₃の波、ｋ₃の山の値とから、ｘ_n ⁽³⁾、ｘ_n ⁽³⁾のバー、ｘ_n ⁽³⁾の波、ｘ_n ⁽³⁾の山の値を求め、パラメータｋ₄、ｋ₄のバー、ｋ₄の波、ｋ₄の山を生成する。
【０１３５】
図８は時間増分生成装置２４の詳細な動作を示すフローチャートである。
時間増分生成装置２４は、まず時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）の値をゼロに設定する。
【０１３６】
次にｋパラメータ生成装置２２から与えられるｋ₁、ｋ₂、ｋ₃、ｋ₄の値にメモリ２５─１に格納されたｂパラメータｂ₁、ｂ₂、ｂ₃、ｂ₄の値をそれぞれ乗算して、（７）式右辺の第１項の値を生成する。そして生成された第１項の値をΨ（ｘ_nのバー，ｔ_n，ｈ）に加算する。
【０１３７】
次にｋパラメータ生成装置２２から与えられるｋ₁のバー、ｋ₂のバー、ｋ₃のバー、ｋ₄のバーの値にメモリ２５─２に格納されたｂパラメータｂ₁のバー、ｂ₂のバー、ｂ₃のバー、ｂ₄のバーの値をそれぞれ乗算する。さらに乱数発生装置２６─１が発生する正規乱数のサンプルを用いて（７）式右辺の第２項の値を生成する。そして生成された第２項の値をΨ（ｘ_nのバー，ｔ_n，ｈ）に加算する。
【０１３８】
次にｋパラメータ生成装置２２から与えられるｋ₁の波、ｋ₂の波、ｋ₃の波、ｋ₄の波の値にメモリ２５─３に格納されたｂパラメータｂ₁の波、ｂ₂の波、ｂ₃の波、ｂ₄の波の値をそれぞれ乗算する。さらに乱数発生装置２６─２が発生する正規乱数のサンプルを用いて（７）式右辺の第３項の値を生成する。そして生成された第３項の値をΨ（ｘ_nのバー，ｔ_n，ｈ）に加算する。
【０１３９】
次にｋパラメータ生成装置２２から与えられるｋ₁の山、ｋ₂の山、ｋ₃の山、ｋ₄の山の値にメモリ２５─４に格納されたｂパラメータｂ₁の山、ｂ₂の山、ｂ₃の山、ｂ₄の山の値をそれぞれ乗算して、（７）式右辺の第４項の値を生成する。そして生成された第３項の値をΨ（ｘ_nのバー，ｔ_n，ｈ）に加算する。
【０１４０】
こうして得られたΨ（ｘ_nのバー，ｔ_n，ｈ）の値が時刻ｔ_nにおける近似値ｘ_nのバーに対する時間増分に相当する。
近似値更新装置２７は時間増分生成装置２４から与えられる時間増分Ψ（ｘ_nのバー，ｔ_n，ｈ）の値をｘ_nのバーの値に加算し、次の参照時刻ｔ_n+1における状態量ｘ（ｔ_n+1）の近似値ｘ_n+1のバーを生成する。そして生成されたｘ_n+1のバーの値をｋパラメータ生成装置２２と出力部２８に与える。
【０１４１】
出力部２８は近似値更新装置２７から与えられる各参照時刻の近似値を格納しておき、近似値の更新が終了すると、これらの近似値をシミュレーション結果として、外部のディスプレイ装置や、シミュレーション対象の物理系を制御する制御装置等に出力する。
【０１４２】
次に本実施例のシミュレーション装置を具体的に、不規則外乱を受ける物理系に適用した結果を説明する。
図９は地球を周回する人工衛星の衛星軌道を示している。衛星軌道は、地球の大気密度の急速な変動と上層大気中の外乱等の影響を受けて変化し、これらが不規則外乱として作用する。
【０１４３】
図９において、状態量ｘは衛星軌道を含む平面内で、地球の中心から人工衛星３０に向かう放射方向における、与えられた衛星軌道からの摂動を表す。ここでｘ₁＝ｘ、ｘ₂＝ｄｘ／ｄｔと置くと、図９の衛星軌道は次のような２次元のモデルにより記述される。
【０１４４】
【数１８】

【０１４５】
ここで、Ｘは（５）式の状態量に相当する２次元の状態量を表し、関数ｆとｇは（５）式の右辺の関数に相当し、２次元の関数値をもつ。またａ、ｂ、ｃは定数である。
【０１４６】
図１０は（１９）式で記述される衛星軌道のモデルを対象にして、図２のシミュレーション装置により数値シミュレーションを行った結果を示す。図１０のシミュレーションでは、ａ＝０．７５、ｂ＝１．５、ｃ＝０．５とした。またシミュレーションの時間ステップの刻み幅ｈはｈ＝２^-8（秒）とし、時刻ｔ₀＝０における初期条件として、ｘ₁（０）＝１．５（ｍ）、ｘ₂（０）＝−０．９（ｍ／ｓ）を与えた。またメモリ２３─１〜２３─４、２５─１〜２５─４に格納するａパラメータ、ｂパラメータとしては、（１６）式の値を用いた。
【０１４７】
図１１は（１９）式で不規則外乱の影響を無視した理想的な衛星軌道のモデルに対するシミュレーション結果を示す。図１０を図１１と比べると、大気密度の変動等の不規則外乱が人工衛星の軌道に与える影響が無視できないことがわかる。
【０１４８】
次に本シミュレーション装置の精度を他のシミュレーション装置と比較するために、次式で与えられる１次元のモデルを考える。
【０１４９】
【数１９】

【０１５０】
これは、（５）式においてｘを１次元の状態量とし、関数ｆ、ｇの関数形をｆ（ｘ）＝ａｘ、ｇ（ｘ）＝ｂｘと置いたモデルに相当する。ただし、ａ、ｂは定数であり、ｘ₀は時刻ｔ₀におけるｘの初期値である。
【０１５１】
（２０）式の確率微分方程式は関数ｆ、ｇの関数形が簡単なため、その厳密解を解析的に求めることができる。この厳密解は次式で与えられる。
【０１５２】
【数２０】

【０１５３】
ここでＷ（ｔ）はブラウン運動モデルの確率過程で、雑音モデルｗ（ｔ）に対応している。またｅｘｐの中の−（１／２）ｂ²ｔの項は、確率積分に現れる特徴的な項である。
【０１５４】
図１２は（２１）式でａ＝１、ｂ＝１、ｘ₀＝１と置き、２^-9秒の時間刻みでサンプリングした結果を示す。また図１３、図１４はそれぞれ０．５次の近似精度のオイラ法を用いたシミュレーション装置、本実施例のシミュレーション装置によるシミュレーション結果を示す。図１３、図１４のシミュレーションにおいてもａ＝１、ｂ＝１、ｘ₀＝１と置き、シミュレーションの時間刻みはｈ＝２^-6（秒）とした。
【０１５５】
図１２、図１３、図１４を比較すると、０．５次の近似精度のオイラ法による結果よりも、１．５の近似精度を保証する本実施例のシミュレーション結果の方が厳密解のサンプリング結果に近いことがわかる。
【０１５６】
本実施例では、ａパラメータとｂパラメータの値の一例として（１６．１）、（１６．２）式の値を採用したが、本発明はこれに限定されるものではなく、（１４．１）〜（１４．２２）および（１５）式の拘束条件を満足する他の値も用いることができる。
【０１５７】
また（７）〜（１１）式の形の時間増分を用いれば、１．５次より高い近似精度を保証するシミュレーション装置の実現も可能である。
【０１５８】
【発明の効果】
本発明のシミュレーション装置によれば、不規則外乱の影響を受ける系に対して、１．５次の高精度を保証するシミュレーション結果を得ることができる。
【０１５９】
また数学モデルで与えられる関数の導関数を求める必要がなく、与えられた関数形のみを入力すればよいので、オペレータの負担が軽減される。さらに導関数の値を計算しなくて済むので、数値処理における計算量が少なく、汎用性に富むシミュレーション装置を提供する。
【０１６０】
このように、不規則外乱の影響を受ける系に対して、高い信頼性を持つシミュレーション結果を容易に得ることができるので、対象とする系の挙動の精密な解析と、その系を制御する高性能な制御装置の実現が可能になる。
【図面の簡単な説明】
【図１】本発明の原理図である。
【図２】本発明の一実施例のシミュレーション装置の構成図である。
【図３】一実施例のシミュレーション装置の動作を示すフローチャートである。
【図４】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図５】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図６】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図７】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図８】一実施例の時間増分生成装置の動作を示すフローチャートである。
【図９】人工衛星の衛星軌道を示す図である。
【図１０】図９の衛星軌道のモデルに対する本実施例のシミュレーション装置によるシミュレーション結果を示す図である。
【図１１】不規則外乱を受けない衛星軌道のモデルに対するシミュレーション結果を示す図である。
【図１２】一具体例のモデルの厳密解を示す図である。
【図１３】図１２のモデルに対する０．５次のオイラ法によるシミュレーション結果を示す図である。
【図１４】図１２のモデルに対する本実施例のシミュレーション装置によるシミュレーション結果を示す図である。
【符号の説明】
１０シミュレーション装置
１１初期設定入力手段
１２増分生成手段
１３近似値更新手段
１４出力手段
２１初期設定入力部
２２ｋパラメータ生成装置
２３─１、２、３、４メモリ
２５─１、２、３、４メモリ
２４時間増分生成装置
２６─１、２乱数発生装置
２７近似値更新手段
２８出力部
３０人工衛星

Claims

不規則外乱を受ける物理系の状態量ｘの時間変化を、関数と乱数によって記述し、時間刻みｈ毎に、状態量ｘの初期値から、状態量ｘの近似値を陽的に求め、シミュレーション結果を出力するシミュレーション装置であって、
関数ｆ（ｘ）およびｇ（ｘ）と時刻ｔ_０における状態量ｘの

を入力する初期設定入力手段と、

を格納する第１のパラメータ格納手段と、

を格納する第２のパラメータ格納手段と、
時刻ｔ _ｎにおける状態量ｘの

を格納する近似値格納手段と、
ｎ番目の時刻ｔ_ｎ＝ｔ_０＋ｎｈにおける第１の乱数ξ、時刻ｔ_ｎにおける

、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、および前記ａパラメータを用いて、状態量ｘの

から

を計算するための

というｋパラメータ計算式を格納する第１の計算式格納手段と、
前記ｋパラメータ、時間刻みｈ、正の実数ν、ΔＷ、

、および前記ｂパラメータを用いて状態量ｘの時間増分Ψを計算するための

という時間増分計算式を格納する第２の計算式格納手段と、
時間増分Ψを用いて状態量ｘの

から時刻ｔ _ｎ＋１における状態量ｘの

を計算するための

という近似値計算式を格納する第３の計算式格納手段と、
第１の乱数ξ、

、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、前記ａパラメータ、および前記ｋパラメータ計算式を用いて、前記ｋパラメータを生成するｋパラメータ生成手段と、
時間増分Ψの値を格納する時間増分格納手段と、
前記時間増分格納手段に０を設定した後、時間刻みｈ、ｂパラメータｂ _ｉ、ｋパラメータｋ _ｉ、および前記時間増分計算式を用いて該時間増分計算式の第１項の値を生成して該時間増分格納手段に加算し、ΔＷ、

、および該時間増分計算式を用いて該時間増分計算式の第２項の値を生成して該時間増分格納手段に加算し、

、および該時間増分計算式を用いて該時間増分計算式の第３項の値を生成して該時間増分格納手段に加算し、時間刻みｈ、正の実数ν、

、および該時間増分計算式を用いて該時間増分計算式の第４項の値を生成して該時間増分格納手段に加算することで、時間増分Ψを生成する増分生成手段と、
生成された時間増分Ψと前記近似値計算式を用いて状態量ｘの

を生成し、前記近似値格納手段の

を

に更新する近似値更新手段と、
前記近似値更新手段から与えられる時間刻みｈ毎の状態量ｘの各近似値を、前記物理系に対する前記シミュレーション結果として出力する出力手段と
を有することを特徴とするシミュレーション装置。
前記近似値更新手段は更新された近似値を前記増分生成手段に入力し、
前記増分生成手段は前記近似値更新手段より与えられる前記更新された近似値から、前記時間増分を生成することを特徴とする請求項１記載のシミュレーション装置。
前記増分生成手段は前記第１の乱数を発生する第１の乱数発生手段と、前記第１の乱数とは独立な前記第２の乱数を発生する第２の乱数発生手段を有することを特徴とする請求項１または２記載のシミュレーション装置。
前記第１の乱数と前記第２の乱数は、互いに独立な、平均０、分散１の正規乱数の系列に属することを特徴とする請求項３記載のシミュレーション装置。
前記第１のパラメータ格納手段および第２のパラメータ格納手段は、前記ａパラメータおよびｂパラメータとして、前記時間増分を前記時間刻みの巾について展開した時の、前記時間刻みについて１．５次までの各項の係数が、前記物理系を記述する前記状態量に関する微分方程式の解を前記時間刻みの巾について展開した時の、前記時間刻みについて１．５次までの各項の係数に一致するような値を格納することを特徴とする請求項１または２記載のシミュレーション装置。
前記第１のパラメータ格納手段および第２のパラメータ格納手段は、前記ａパラメータおよびｂパラメータとして、

という条件を満たす値を格納することを特徴とする請求項５記載のシミュレーション装置。
不規則外乱を受ける物理系の状態量ｘの時間変化を、関数と乱数によって記述し、時間刻みｈ毎に、状態量ｘの初期値から、状態量ｘの近似値を陽的に求め、シミュレーション結果を出力するシミュレーション方法であって、
初期設定入力手段が、関数ｆ（ｘ）およびｇ（ｘ）と時刻ｔ_０における状態量ｘの

を入力し、
ｋパラメータ生成手段が、ｎ番目の時刻ｔ_ｎ＝ｔ_０＋ｎｈにおける第１の乱数ξ、時刻ｔ_ｎにおける

、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、第１のパラメータ格納手段に格納された

、および第１の計算式格納手段に格納された、

というｋパラメータ計算式を用いて、近似値格納手段に格納された時刻ｔ _ｎにおける状態量ｘの

から

を生成し、
増分生成手段が、状態量ｘの時間増分Ψの値を格納する時間増分格納手段に０を設定した後、時間刻みｈ、第２のパラメータ格納手段に格納されたｂパラメータｂ _ｉ（ｉ＝１，．．．，ｓ）、生成されたｋパラメータｋ _ｉ、および第２の計算式格納手段に格納された、

という時間増分計算式を用いて、該時間増分計算式の第１項の値を生成して該時間増分格納手段に加算し、ΔＷ、前記第２のパラメータ格納手段に格納された

、生成された

、および該時間増分計算式を用いて該時間増分計算式の第２項の値を生成して該時間増分格納手段に加算し、

、前記第２のパラメータ格納手段に格納された

、生成された

、および該時間増分計算式を用いて該時間増分計算式の第３項の値を生成して該時間増分格納手段に加算し、時間刻みｈ、正の実数ν、前記第２のパラメータ格納手段に格納された

、生成された

、および該時間増分計算式を用いて該時間増分計算式の第４項の値を生成して該時間増分格納手段に加算することで、時間増分Ψを生成し、
近似値更新手段が、生成された時間増分Ψと第３の計算式格納手段に格納された、

という近似値計算式を用いて、状態量ｘの

から時刻ｔ_ｎ＋１における状態量ｘの

を生成して、前記近似値格納手段の

を

に更新し、
出力手段が、前記近似値更新手段から与えられる時間刻みｈ毎の状態量ｘの各近似値を、前記物理系に対する前記シミュレーション結果として出力する
ことを特徴とするシミュレーション方法。