JPH07134704A

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

Info

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

Abstract

(57)【要約】【目的】本発明は不規則外乱を受ける系の数値シミュ
レーションに関し、高い信頼性と汎用性を持つシミュレ
ーション装置を提供することを目的とする。【構成】本発明のシミュレーション装置１０は、初期
設定入力手段１１、増分生成手段１２、近似値更新手段
１３、出力手段１４を有する。初期設定入力手段１１は
状態量の初期値と、不規則外乱を受ける系を記述する関
数を増分生成手段１２に入力し、増分生成手段１２は２
種類の乱数系列を用いて、上記関数の導関数を用いるこ
となく、近似値の時間増分を生成する。近似値更新手段
１３は生成された時間増分により、近似値を更新し、出
力手段１４はシミュレーション結果を出力する。得られ
るシミュレーション結果は１．５次の近似精度を保証し
ている。

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, 199
2）は２種類の確率過程を用いて、確率的Heun法のΨ
（ｘ_n，ｔ_n，ｈ）をさらに拡張している。この方法は
一般に１．５次の高い近似精度を保証しており、ｋパラ
メータは関数ｆ（ｘ）とｇ（ｘ）のそれぞれの１階の導
関数ならびに、ｇ（ｘ）の２階の導関数の値を含んでい
る。

【００３１】漸近的に有効な解法（NEWTON, SIAM J. Vo
l.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）のまわりでのテ
ーラ展開は、次式で記述されることが知られている（It
o, 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/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+ ₁のバーの値をｋ
パラメータ生成装置２２と出力部２８に与える。

【０１４１】出力部２８は近似値更新装置２７から与え
られる各参照時刻の近似値を格納しておき、近似値の更
新が終了すると、これらの近似値をシミュレーション結
果として、外部のディスプレイ装置や、シミュレーショ
ン対象の物理系を制御する制御装置等に出力する。

【０１４２】次に本実施例のシミュレーション装置を具
体的に、不規則外乱を受ける物理系に適用した結果を説
明する。図９は地球を周回する人工衛星の衛星軌道を示
している。衛星軌道は、地球の大気密度の急速な変動と
上層大気中の外乱等の影響を受けて変化し、これらが不
規則外乱として作用する。

【０１４３】図９において、状態量ｘは衛星軌道を含む
平面内で、地球の中心から人工衛星３０に向かう放射方
向における、与えられた衛星軌道からの摂動を表す。こ
こでｘ₁＝ｘ、ｘ₂＝ｄｘ／ｄｔと置くと、図９の衛星
軌道は次のような２次元のモデルにより記述される。

【０１４４】

【数１８】

【０１４５】ここで、Ｘは（５）式の状態量に相当する
２次元の状態量を表し、関数ｆとｇは（５）式の右辺の
関数に相当し、２次元の関数値をもつ。またａ、ｂ、ｃ
は定数である。

【０１４６】図１０は（１９）式で記述される衛星軌道
のモデルを対象にして、図２のシミュレーション装置に
より数値シミュレーションを行った結果を示す。図１０
のシミュレーションでは、ａ＝０．７５、ｂ＝１．５、
ｃ＝０．５とした。またシミュレーションの時間ステッ
プの刻み幅ｈはｈ＝２^-8（秒）とし、時刻ｔ₀＝０にお
ける初期条件として、ｘ₁（０）＝１．５（ｍ）、ｘ₂
（０）＝−０．９（ｍ／ｓ）を与えた。またメモリ２３
─１〜２３─４、２５─１〜２５─４に格納するａパラ
メータ、ｂパラメータとしては、（１６）式の値を用い
た。

【０１４７】図１１は（１９）式で不規則外乱の影響を
無視した理想的な衛星軌道のモデルに対するシミュレー
ション結果を示す。図１０を図１１と比べると、大気密
度の変動等の不規則外乱が人工衛星の軌道に与える影響
が無視できないことがわかる。

【０１４８】次に本シミュレーション装置の精度を他の
シミュレーション装置と比較するために、次式で与えら
れる１次元のモデルを考える。

【０１４９】

【数１９】

【０１５０】これは、（５）式においてｘを１次元の状
態量とし、関数ｆ、ｇの関数形をｆ（ｘ）＝ａｘ、ｇ
（ｘ）＝ｂｘと置いたモデルに相当する。ただし、ａ、
ｂは定数であり、ｘ₀は時刻ｔ₀におけるｘの初期値で
ある。

【０１５１】（２０）式の確率微分方程式は関数ｆ、ｇ
の関数形が簡単なため、その厳密解を解析的に求めるこ
とができる。この厳密解は次式で与えられる。

【０１５２】

【数２０】

【０１５３】ここでＷ（ｔ）はブラウン運動モデルの確
率過程で、雑音モデルｗ（ｔ）に対応している。またｅ
ｘｐの中の−（１／２）ｂ²ｔの項は、確率積分に現れ
る特徴的な項である。

【０１５４】図１２は（２１）式でａ＝１、ｂ＝１、ｘ
₀＝１と置き、２^-9秒の時間刻みでサンプリングした結
果を示す。また図１３、図１４はそれぞれ０．５次の近
似精度のオイラ法を用いたシミュレーション装置、本実
施例のシミュレーション装置によるシミュレーション結
果を示す。図１３、図１４のシミュレーションにおいて
もａ＝１、ｂ＝１、ｘ₀＝１と置き、シミュレーション
の時間刻みはｈ＝２^-6（秒）とした。

【０１５５】図１２、図１３、図１４を比較すると、
０．５次の近似精度のオイラ法による結果よりも、１．
５の近似精度を保証する本実施例のシミュレーション結
果の方が厳密解のサンプリング結果に近いことがわか
る。

【０１５６】本実施例では、ａパラメータとｂパラメー
タの値の一例として（１６．１）、（１６．２）式の値
を採用したが、本発明はこれに限定されるものではな
く、（１４．１）〜（１４．２２）および（１５）式の
拘束条件を満足する他の値も用いることができる。

【０１５７】また（７）〜（１１）式の形の時間増分を
用いれば、１．５次より高い近似精度を保証するシミュ
レーション装置の実現も可能である。

【０１５８】

【発明の効果】本発明のシミュレーション装置によれ
ば、不規則外乱の影響を受ける系に対して、１．５次の
高精度を保証するシミュレーション結果を得ることがで
きる。

【０１５９】また数学モデルで与えられる関数の導関数
を求める必要がなく、与えられた関数形のみを入力すれ
ばよいので、オペレータの負担が軽減される。さらに導
関数の値を計算しなくて済むので、数値処理における計
算量が少なく、汎用性に富むシミュレーション装置を提
供する。

【０１６０】このように、不規則外乱の影響を受ける系
に対して、高い信頼性を持つシミュレーション結果を容
易に得ることができるので、対象とする系の挙動の精密
な解析と、その系を制御する高性能な制御装置の実現が
可能になる。

【図面の簡単な説明】

【図１】本発明の原理図である。

【図２】本発明の一実施例のシミュレーション装置の構
成図である。

【図３】一実施例のシミュレーション装置の動作を示す
フローチャートである。

【図４】一実施例のｋパラメータ生成装置の動作を示す
フローチャートである。

【図５】一実施例のｋパラメータ生成装置の動作を示す
フローチャートである。

【図６】一実施例のｋパラメータ生成装置の動作を示す
フローチャートである。

【図７】一実施例のｋパラメータ生成装置の動作を示す
フローチャートである。

【図８】一実施例の時間増分生成装置の動作を示すフロ
ーチャートである。

【図９】人工衛星の衛星軌道を示す図である。

【図１０】図９の衛星軌道のモデルに対する本実施例の
シミュレーション装置によるシミュレーション結果を示
す図である。

【図１１】不規則外乱を受けない衛星軌道のモデルに対
するシミュレーション結果を示す図である。

【図１２】一具体例のモデルの厳密解を示す図である。

【図１３】図１２のモデルに対する０．５次のオイラ法
によるシミュレーション結果を示す図である。

【図１４】図１２のモデルに対する本実施例のシミュレ
ーション装置によるシミュレーション結果を示す図であ
る。

【符号の説明】

１０シミュレーション装置１１初期設定入力手段１２増分生成手段１３近似値更新手段１４出力手段２１初期設定入力部２２ｋパラメータ生成装置２３─１、２、３、４メモリ２５─１、２、３、４メモリ２４時間増分生成装置２６─１、２乱数発生装置２７近似値更新手段２８出力部３０人工衛星

Claims

【特許請求の範囲】

【請求項１】不規則外乱を受ける物理系の状態量の時
間変化を、関数と乱数によって記述し、時間刻み毎に、
該状態量の初期値から、前記状態量の近似値を陽的に求
め、シミュレーション結果を出力する情報処理装置にお
いて、前記関数と前記初期値を入力する初期設定入力手段（１
１）と、前記乱数に対応する第１の乱数と第２の乱数および前記
時間刻みとを用いて、前記初期値から、前記初期設定入
力手段（１１）より与えられる前記関数の関数値とし
て、３種類以上の異なるｋパラメータを生成し、生成し
た前記ｋパラメータを用いて、前記状態量の時間増分を
生成する増分生成手段（１２）と、前記時間増分を用いて前記近似値を更新する近似値更新
手段（１３）と、前記近似値更新手段（１３）から与えられる前記時間刻
み毎の前記状態量の各近似値を、前記物理系に対する前
記シミュレーション結果として出力する出力手段（１
４）とを、有することを特徴とするシミュレーション装
置（１０）。
【請求項２】前記増分生成手段（１２）は、前記第
１の乱数と、前記第２の乱数と、前記時間刻みの平方根
に比例しかつ前記第１の乱数と第２の乱数に依存しない
項とを用いて、前記ｋパラメータから前記状態量の時間
増分を生成することを特徴とする請求項１記載のシミュ
レーション装置（１０）。
【請求項３】前記近似値更新手段（１３）は更新した
近似値を前記増分生成手段（１２）に入力し、前記増分生成手段（１２）は前記近似値更新手段（１
３）より与えられる前記更新された近似値から、前記時
間増分を生成することを特徴とする請求項２記載のシミ
ュレーション装置（１０）。
【請求項４】前記増分生成手段（１２）は前記第１の
乱数を発生する第１の乱数発生手段（２６─１）と、前
記第１の乱数とは独立な前記第２の乱数を発生する第２
の乱数発生手段（２６─２）を有し、前記時間増分は、前記時間刻みに比例する項と、前記時
間刻みの平方根と前記第１の乱数の積に比例しかつ前記
第２の乱数に依存しない項と、前記時間刻みの平方根と
前記第２の乱数の積に比例しかつ前記第１の乱数に依存
しない項と、前記時間刻みの平方根に比例しかつ前記第
１の乱数と第２の乱数に依存しない項とを含むことを特
徴とする請求項１記載のシミュレーション装置（１
０）。
【請求項５】前記第１の乱数と前記第２の乱数は、互
いに独立な、平均０、分散１の正規乱数の系列に属する
ことを特徴とする請求項４記載のシミュレーション装置
（１０）。
【請求項６】前記増分生成手段（１２）はさらに、前
記時間刻み、前記第１の乱数、前記第２の乱数のうちの
少なくとも一つと、前記状態量の近似値とを用いて、前
記状態量に対応する標本状態量を生成し、前記関数を用
いて、生成した前記標本状態量から前記ｋパラメータを
生成するパラメータ生成手段（２２）と、前記パラメータ生成手段（２２）から与えられる前記ｋ
パラメータと、前記第１の乱数、前記第２の乱数および
前記時間刻みとを用いて、前記時間増分を生成する時間
増分生成手段（２４）とを有することを特徴とする請求
項３記載のシミュレーション装置（１０）。
【請求項７】前記増分生成手段（１２）はさらに、第
１のパラメータ格納手段（２３─１〜２３─４）と第２
のパラメータ格納手段（２５─１〜２５─４）を有し、前記パラメータ生成手段（２２）は、前記第１のパラメ
ータ格納手段（２３─１〜２３─４）に格納された第１
のパラメータを用いて、前記ｋパラメータを生成し、前記時間増分生成手段（２４）は、前記第２のパラメー
タ格納手段（２５─１〜２５─４）に格納された第２の
パラメータを用いて、前記ｋパラメータから、前記時間
増分を生成することを特徴とする請求項６記載のシミュ
レーション装置（１０）。
【請求項８】前記第１のパラメータと前記第２のパラ
メータは、前記時間増分を前記時間刻みの巾について展
開した時の、前記時間刻みについて１．５次までの各項
の係数が、前記物理系を記述する前記状態量に関する微
分方程式の解を前記時間刻みの巾について展開した時
の、前記時間刻みについて１．５次までの各項の係数に
一致するような値であることを特徴とする請求項７記載
のシミュレーション装置（１０）。
【請求項９】前記パラメータ生成手段（２２）は一つ
の前記近似値に対して、種類毎に、複数個のｋパラメー
タを生成し、前記近似値と既に生成された前記ｋパラメ
ータの線型結合により、前記標本状態量を生成すること
を特徴とする請求項６記載のシミュレーション装置（１
０）。
【請求項１０】前記パラメータ生成手段（２２）は、
前記関数を用いて、前記標本状態量を前記ｋパラメータ
に変換することを特徴とする請求項９記載のシミュレー
ション装置（１０）。
【請求項１１】前記時間増分生成手段（２４）は、前
記ｋパラメータの線型結合により、前記時間増分を生成
することを特徴とする請求項１０記載のシミュレーショ
ン装置（１０）。
【請求項１２】不規則外乱を受ける物理系の状態量の
時間変化を、関数と乱数によって記述し、時間刻み毎
に、該状態量の初期値から、前記状態量の近似値を陽的
に求め、シミュレーション結果を出力する情報処理装置
において、前記関数と前記初期値を入力し、前記乱数に対応する第１の乱数と第２の乱数および前記
時間刻みとを用いて、前記初期値から、前記関数の関数
値として、３種類以上の異なるｋパラメータを生成し、生成した前記ｋパラメータを用いて、前記状態量の時間
増分を生成し、前記時間増分を用いて前記近似値を更新し、前記時間刻み毎の前記状態量の各近似値を前記シミュレ
ーション結果として出力する各ステップから成ることを
特徴とするシミュレーション方法。
【請求項１３】前記更新した近似値から、さらに前記
関数に基づき、前記第１の乱数と第２の乱数および、前
記時間刻みとを用いて、前記状態量の時間増分を生成す
るステップを有することを特徴とする請求項１２記載の
シミュレーション方法。
【請求項１４】前記第１の乱数と、前記第２の乱数
と、前記時間刻みの平方根に比例しかつ前記第１の乱数
と第２の乱数に依存しない項とを用いて、前記ｋパラメ
ータから前記状態量の時間増分を生成することを特徴と
する請求項１２記載のシミュレーション方法。
【請求項１５】前記時間刻みに比例する項と、前記
時間刻みの平方根と前記第１の乱数の積に比例しかつ前
記第２の乱数に依存しない項と、前記時間刻みの平方根
と前記第２の乱数の積に比例しかつ前記第１の乱数に依
存しない項と、前記時間刻みの平方根に比例しかつ前記
第１の乱数と第２の乱数に依存しない項とを用いて、前
記ｋパラメータから前記状態量の時間増分を生成するこ
とを特徴とする請求項１２記載のシミュレーション方
法。