"格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、前記複数の粒子の前記擾乱項の値を前記格納部に格納し、前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する処理をコンピュータに実行させるプログラム。前記コンピュータは、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項1記載のプログラム。前記コンピュータは、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項1記載のプログラム。前記コンピュータは、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項1、2、又は3記載のプログラム。現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項を含み、前の時刻における該複数の粒子の位置、速度、及び密度から該現在時刻における該複数の粒子の位置、速度、及び密度を求める計算式の情報を格納する格納部と、前記複数の粒子の前記擾乱項の値を求め、該複数の粒子の該擾乱項の値を前記格納部に格納し、初期時刻における該複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用い、該複数の粒子の該擾乱項の値を前記計算式に適用することで、前記現在時刻における該複数の粒子の位置、速度、及び密度を求め、得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の該擾乱項の値と該現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、該初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求める処理部と、前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力する出力部とを備えることを特徴とするシミュレーション装置。前記処理部は、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項5記載のシミュレーション装置。前記処理部は、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項5記載のシミュレーション装置。前記処理部は、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項5、6、又は7記載のシミュレーション装置。コンピュータにより実行されるシミュレーション方法であって、格納部に格納された、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱項の情報を用いて、該擾乱項の値を求め、前記複数の粒子の前記擾乱項の値を前記格納部に格納し、前記格納部に格納された、前記複数の粒子の前記擾乱項を含む計算式であって、前の時刻における該複数の粒子の位置、速度、及び密度から前記現在時刻における該複数の粒子の位置、速度、及び密度を求める該計算式の情報を用い、初期時刻における該複数の粒子の位置、速度、及び密度を該前の時刻における該複数の粒子の位置、速度、及び密度として用い、前記複数の粒子の前記擾乱項の値を該計算式に適用することで、該現在時刻における該複数の粒子の位置、速度、及び密度を求め、得られた前記複数の粒子の位置、速度、及び密度を前記前の時刻における前記複数の粒子の位置、速度、及び密度として用いて該複数の粒子の前記擾乱項の値と前記現在時刻における該複数の粒子の位置、速度、及び密度を求める処理を繰り返すことで、前記初期時刻より後の複数の時刻における該複数の粒子の位置、速度、及び密度を順次求め、前記複数の時刻における前記複数の粒子の位置をシミュレーション結果として出力することを特徴とするシミュレーション方法。前記コンピュータは、前記現在時刻における各粒子の位置又は速度に付加された擾乱の向きを示す乱数を用いて、前記複数の粒子の位置又は速度に付加された前記擾乱項の値を求めることを特徴とする請求項9記載のシミュレーション方法。前記コンピュータは、前記現在時刻における各粒子の加速度、圧力、又は密度に付加された擾乱の大きさを示す乱数を用いて、前記複数の粒子の加速度、圧力、又は密度に付加された前記擾乱項の値を求めることを特徴とする請求項9記載のシミュレーション方法。前記コンピュータは、前記前の時刻における各粒子の位置、速度、及び密度を用いて該前の時刻における各粒子のエネルギーを求め、得られたエネルギーが前記初期時刻における該粒子のエネルギーより大きい場合、該粒子の擾乱項の値を零に設定することを特徴とする請求項9、10、又は11記載のシミュレーション方法。上述したGodunovSPH法等を用いても、粒子群のような粒子直径の10倍程度以上の規模の大きなスケールの運動から、各粒子の乱雑な運動が生成されてしまい、大きなスケールの整った流れが散逸してしまう結果、シミュレーションの精度が悪化してしまうという問題は、依然として残されている。例えば、水面波の伝播等の波動解析において、粒子法を用いて数値計算を進めて行くと、整った流れから各粒子の乱雑な運動が生成されてしまい、波高が現実よりも過剰に減衰してしまうことがある。図1は、このような整った流れを示している。図1の破線101~103は、流体中の整った流れの流線を表し、ベクトル121~125は、それぞれ粒子111~115の速度ベクトルを表す。このとき、粒子111及び113は、粒子111に働くベクトル131及び粒子113に働くベクトル132が示すように、流線101~103から外れる方向の力を粒子114から受けてしまう。このように、粒子直径の10倍程度以上の規模の大きなスケールの整った流れからは、剛体球の玉突き衝突のような現象が起こり、各粒子の乱雑な運動が生成されるとともに、流線101~103が示す整った運動が壊されてしまうものと考えられる。この問題は、系の運動エネルギーが保存しても起こる、粒子法に顕著な問題である。粒子法による水面波の散逸は、粒子直径の10倍程度以上の規模の大きなスケールの整った流れから各粒子の乱雑な運動が生成されるプロセスと、生成された乱雑な運動が人工粘性項により散逸されるプロセスの2つのプロセスによるものと考えられる。GodunovSPH法は、人工粘性項の効果を小さくするように開発された方法であるが、それでも小さなスケールでの乱雑な運動には散逸効果が現れ、運動を散逸させている。本発明の課題は、粒子法を用いた運動シミュレーションにおいて、粒子群による大きなスケールの運動から各粒子の乱雑な運動が生成されることを防ぐことである。シミュレーション装置は、格納部、処理部、及び出力部を備える。格納部は、現在時刻における複数の粒子の位置、速度、加速度、圧力、又は密度に付加された擾乱じょうらん項を含み、前の時刻におけるそれらの粒子の位置、速度、及び密度から現在時刻におけるそれらの粒子の位置、速度、及び密度を求める計算式の情報を格納する。処理部は、上記複数の粒子の擾乱項の値を求め、それらの粒子の擾乱項の値を格納部に格納する。次に、初期時刻におけるそれらの粒子の位置、速度、及び密度を前の時刻におけるそれらの粒子の位置、速度、及び密度として用い、それらの粒子の擾乱項の値を上記計算式に適用することで、現在時刻におけるそれらの粒子の位置、速度、及び密度を求める。そして、得られた上記複数の粒子の位置、速度、及び密度を前の時刻におけるそれらの粒子の位置、速度、及び密度として用いてそれらの粒子の擾乱項の値と現在時刻におけるそれらの粒子の位置、速度、及び密度を求める処理を繰り返す。これにより、初期時刻より後の複数の時刻におけるそれらの粒子の位置、速度、及び密度を順次求める。出力部は、上記複数の時刻における上記複数の粒子の位置をシミュレーション結果として出力する。開示のシミュレーション装置によれば、粒子法を用いた運動シミュレーションにおいて、各粒子の位置、速度、加速度、圧力、又は密度に擾乱項を付加することで、粒子群による大きなスケールの運動から各粒子の乱雑な運動が生成されることが防止される。例えば、水面波等の問題に対するシミュレーションにおいて、整った流れから生成される乱雑な運動による波高の過剰な減衰を防止することができ、その結果、シミュレーション精度を向上させることが可能となる。以下、図面を参照しながら、実施形態を詳細に説明する。粒子法の数値計算アルゴリズムには、各粒子の位置及び速度を更新する手順が含まれる。この手順としては、以下のような形式が採用されることが多い。式11及び12において、下付添え字a及びbはそれぞれ粒子a及びbを表し、上付添え字n及びn1はシミュレーションにおける時刻時間ステップを表す。dtは時間について離散化した際の時間刻み幅を表し、は、粒子aが粒子bから受ける力を表す。ここで、各粒子の速度に擾乱を与える場合、δを小さな実数又は乱数とし、Γを実数として、を単位ベクトルとして、式11を以下のような形式に拡張する。の向きは、粒子群の速度分布や位置関係等、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。この場合、図2に示すように、各粒子の速度ベクトル201に擾乱ベクトル203が外力として加えられ、速度ベクトル202が生成される。また、各粒子の位置に擾乱を与える場合、式12を以下のような形式に拡張する。この場合、図3に示すように、各粒子の位置ベクトルに擾乱ベクトル301が変位として加えられる。擾乱ベクトル301の向きは、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。また、各粒子の加速度に擾乱を与える場合、式11を以下のような形式に拡張する。この場合、粒子aが粒子bから受ける力に擾乱が加えられる。さらに、流体計算を行なう際には各粒子の速度、位置、又は加速度に擾乱を与える代わりに、各粒子の圧力もしくは密度に擾乱を与えても良い。この場合、各粒子の擾乱を加えた圧力もしくは密度を以下のように決め、粒子aの圧力もしくは密度として用いる。式16について圧力の擾乱として付加されている実数δpや式17について密度の擾乱として付加されている実数δρは、粒子群の速度分布や位置関係等、計算中の粒子の情報から確定的に決めても良く、乱数を用いて決めても良い。このように、各粒子の位置、速度、加速度、圧力、又は密度に擾乱を与えることで、人工粘性項により散逸してしまう各粒子の乱雑な運動を付加することができ、散逸された運動エネルギーを注入することができる。これにより、大きなスケールの運動から小さなスケールのランダムな運動へと運動エネルギーが流入してしまう現象を防ぐことが可能になる結果、シミュレーションの精度を向上させることができる。粒子のランダムな運動では、人工粘性項による効果により運動エネルギーが散逸してしまうが、粒子のランダムな運動以外の大きなスケールの運動では、人工粘性項による効果は小さいと考えられる。さらに、微小なスケールのランダムな運動がある程度運動エネルギーを保有している場合は、粒子直径の10倍程度以上の規模の大きなスケールの運動から運動エネルギーが流入しなくなると考えられる。したがって、各粒子の物理量に擾乱を与えて粒子のランダムな運動エネルギーを注入すれば、粒子直径の10倍程度以上の規模の大きなスケールの運動を散逸させることなく、各粒子の位置を計算することができる。図4は、実施形態のシミュレーション装置の構成例を示している。図4のシミュレーション装置は、処理部401、格納部402、及び出力部403を備え、初期条件411に基づいて粒子法による数値計算を行って、シミュレーション結果412を出力する。格納部402は、現在時刻における複数の粒子についての物理量に付加された擾乱項を含み、前の時刻におけるそれらの粒子の位置、速度、及び密度から現在時刻におけるそれらの粒子の位置、速度、及び密度を求める計算式の情報を格納する。なお、本実施形態では、各粒子の速度の物理量に対して擾乱を加えるシミュレーションについて説明をおこなうが、速度以外の物理量、例えば、位置、加速度、密度や圧力等に対して擾乱を加えてシミュレーションを行っても良い。処理部401は、初期条件411に含まれる、初期時刻におけるそれぞれの粒子の位置、速度、及び密度を用いて、それぞれの粒子の速度についての擾乱項の値を求め、得られた擾乱項の値を格納部402に格納する。次に、初期時刻におけるそれぞれの粒子の位置、速度、及び密度を前の時刻におけるそれぞれの粒子の位置、速度、及び密度として用い、それぞれの粒子の速度についての擾乱項の値を計算式に適用することで、現在時刻におけるそれぞれの粒子の位置、速度、及び密度を求める。そして、処理部401は、得られたそれぞれの粒子の位置、速度、及び密度を前の時刻におけるそれぞれの粒子の位置、速度、及び密度として用いて、それぞれの粒子の速度についての擾乱項の値と現在時刻におけるそれぞれの粒子の位置、速度、及び密度を求める処理を繰り返す。これにより、初期時刻より後の複数の時刻におけるそれぞれの粒子の位置、速度、及び密度が、順次求められる。出力部403は、得られた複数の時刻におけるそれぞれの粒子の位置をシミュレーション結果412として出力する。まず、以下の3次元流体の運動方程式を用いて、粒子法によりシミュレーションを行う場合を考える。式21~23において、ρx,y,z,t、px,y,z,t、及びは、それぞれ、流体の密度場、圧力場、及び速度場を現し、cは音速を表す。はベクトル微分演算子を表し、x、y、及びzは空間座標を表し、tは時刻を表す。は、z方向鉛直方向の単位ベクトルを表し、gは重力加速度を表す。ρ0は、圧力が0のときの密度を表す。式21は質量保存則を表し、左辺は流体粒子の密度時間変化率、右辺は流体粒子の体積変化率を表す。すなわち、式21は、流体粒子の密度変化率と体積変化率との関係を表している。式22は、流体粒子についての運動方程式である。式22の左辺は流体粒子の慣性項、右辺第1項は流体の圧力による力、右辺第2項は重力を表している。式21~23をGodunovSPH法に基づいて離散化し、各粒子の速度に擾乱項を付加すると、以下のような計算式が得られる。式24~27において、は、それぞれ時刻nにおける粒子aの位置ベクトル、速度ベクトル、及び密度を表す。粒子bについても同様の表記が用いられる。dtは時間刻み幅を表し、前の時刻nにおける各粒子の位置ベクトル、速度ベクトル、及び密度が既知であるとき、式24~27により、現在時刻n1における粒子aの位置ベクトル、速度ベクトル、及び密度が求められる。ここで、Wはカーネル関数であり、は再標準化行列である。3次元の場合、再標準化行列は以下のように定義される。式24~28において、は、以下のように定義される。式25及び26において、は、粒子aと粒子bの間での1次元リーマン問題を解いて得られた時空間の中間値を表し、具体的には以下のように決定される。まず、粒子a及びbに関して、以下の特性関数を定める。さらに、勾配を以下のように計算する。式46は、速度勾配テンソルを表している。式41~50を用いて、上記時空間の中間値を以下のように決定する。ところで、式25の右辺第4項が速度に付加された擾乱項であり、以下のように決定される。時刻nにおける粒子aのエネルギーは、以下のように書くことができる。式61のエネルギーを用いて、擾乱項に含まれるを以下のように定義する。次に、乱数γ1及びγ2を用いて、乱数ベクトルを以下のように決定する。乱数γ1及びγ2としては、例えば、-1から1までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、粒子aの速度ベクトルを基に、擾乱の向きを表す単位ベクトルを以下のように決定する。図5は、図4のシミュレーション装置が式24~64を用いて3次元粒子群の位置を解析するシミュレーション処理の例を示すフローチャートである。この場合、式24~64に相当する計算式の情報が格納部402に格納されている。計算式の情報には、それぞれの式に含まれる変数、定数等の情報と、それらの変数及び定数を用いた演算の情報が含まれる。処理部401は、この計算式の情報を参照しながら、初期時刻0から時刻NNは1以上の整数までの各時刻における各粒子の位置ベクトル、速度ベクトル、及び密度を計算する。まず、処理部401は、時刻を表す変数nに0を設定しステップ501、初期条件411を取得するステップ502。初期条件411はユーザから与えられる場合もあり、処理部401が所定のアルゴリズムで生成する場合もある。この初期条件411には、初期時刻0における各粒子の位置ベクトル、速度ベクトル、密度、及びエネルギーの値と、位置ベクトルが取り得る値の範囲等が含まれる。次に、処理部401は、式61により、粒子群中の最初の粒子を処理対象として選択し、時刻nにおける処理対象粒子のエネルギーを計算するステップ503。そして、時刻nにおけるエネルギーを初期時刻0におけるエネルギー初期エネルギーと比較するステップ504。時刻nにおけるエネルギーが初期エネルギー以下であればステップ504,YES、式62によりを計算するステップ506。時刻nにおけるエネルギーが初期エネルギーより大きければステップ504,NO、式62に従ってを0に設定するステップ505。ステップ505の処理を行うことで、処理対象粒子のエネルギーが初期エネルギーを超えた場合に速度の擾乱項の効果をなくすことができ、ランダムな運動エネルギーの過剰な注入が防止される。次に、処理部401は、乱数γ1及びγ2を発生し、式63及び64を用いて、それぞれの粒子の速度についての擾乱項式25の右辺第4項の値を計算し、その値を格納部402に格納するステップ507。そして、時刻nにおける各粒子の位置ベクトル、速度ベクトル、及び密度と速度の擾乱項の値を用いて、式24~55により、時刻n1における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算するステップ508。次に、処理部401は、粒子群中のすべての粒子について処理が終了したか否かを判定するステップ509。未処理の粒子が残っていればステップ509,NO、次の粒子を処理対象として選択して、ステップ503~508の処理を繰り返す。そして、すべての粒子について処理が終了すればステップ509,YES、時刻n1における各粒子の位置ベクトルを出力部403に転送する。出力部403は、転送された各粒子の位置ベクトルを、時刻n1におけるシミュレーション結果412として所定の形式で出力するステップ510。シミュレーション結果412の出力方法としては、表示装置の画面に粒子群の画像情報として表示したり、有線又は無線の通信回線を介してプリンタや情報処理装置に転送したりする方法が考えられる。シミュレーション結果412として、各粒子の速度ベクトル、密度等を含めて出力しても良い。次に、処理部401は、nがN-1に達したか否かを判定するステップ511。nがN-1より小さければステップ511,NO、n=n1と設定してステップ512、ステップ503~510の処理を繰り返す。そして、nがN-1に達するとステップ511,YES、シミュレーション処理を終了する。なお、図5のシミュレーション処理において、必ずしもすべてのステップを実行する必要はなく、シミュレーション対象や初期条件に応じて一部のステップを省略することも可能である。例えば、ランダムな運動エネルギーが過剰に注入される可能性が低い場合は、ステップ504及び505の処理を省略しても良い。3次元粒子群の代わりに2次元粒子群の位置を解析する場合は、3次元ベクトルを2次元ベクトルに変更し、計算式の一部を変更することで、図5と同様のシミュレーション処理を行うことができる。この場合、式28の再標準化行列の代わりに、以下のような再標準化行列が用いられる。また、式46の速度勾配テンソルの代わりに、以下のような速度勾配テンソルが用いられる。さらに、乱数γ1及びγ2の代わりに乱数γを用いて、以下のように回転行列Aが決定される。乱数γとしては、例えば、-1から1までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、式64の代わりに次式を用いて、速度の擾乱の向きを表す単位ベクトルを決定する。ここで、自由表面流での一般的なベンチマークテストとして、2次元流体のダム崩壊問題に対するシミュレーションの結果を説明する。ダム崩壊問題は、初期時刻0で壁に接して静止している流体の塊が、一定の重力の元で時間の経過とともに崩れていく現象である。図6は、初期時刻0における流体の2次元形状を示している。図6において、x軸、y軸、及びx=0.584mの直線602がダムの左壁、底、及び右壁にそれぞれ対応し、x=0~0.146m及びy=0~0.292mの範囲に長方形601の流体が静止している。時刻0t=0における初期条件としては、長方形601内に格子間隔0.00365mの一様格子を設けて粒子群を配置し、各粒子の初期速度をm/sとし、初期密度をとした。密度ρ0は1000kg/m3、重力加速度gは9.8m/s2、音速cは23.92m/sとした。時間刻み幅dtとしては1.526×105sを用いた。図7及び図8は、式25に速度についての擾乱項を付加しない標準的なGodunovSPH法によるシミュレーション結果の画像を示しており、t=0~9.99458sの12個の時刻における流体の形状を表している。一方、図9及び図10は、式25の速度についての擾乱項を付加した場合のシミュレーション結果の画像を示しており、t=0~9.99898sの12個の時刻における流体の形状を表している。速度の擾乱項を付加しない場合には、図8のt=6.87602s以降において波がほとんど消えてしまっている。これに対して、速度の擾乱項を付加した場合には、図10のt=6.87432s以降においても波が残っており、速度の擾乱項を付加しない場合の図8のt=6.87602s以降において見られた波の過剰な減衰が抑止されることにより、波の減衰性が改善されていることが分かる。上述したように、式25では各粒子の速度に擾乱を与えているが、速度の代わりに位置、加速度、圧力、又は密度に擾乱を与えても良い。例えば、各粒子の位置に擾乱を与える場合は、式25及び27を以下のように書き換えて、図5のシミュレーション処理を実行することができる。式76は式25から右辺第4項の速度についての擾乱項を取り除いた計算式に相当し、式77の右辺第3項が位置に付加された擾乱項である。この場合、処理部401は、図5のステップ507において、式25の右辺第4項を計算する代わりに、式77の右辺第3項を計算する。乱数を用いて擾乱の向きを表す単位ベクトルを決定する方法は、速度に擾乱を与える場合と同様である。そして、図5のステップ508において、式25及び27の代わりに式76及び77を用いて、時刻n1における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。次に、各粒子の加速度に擾乱を与える場合は、式25を以下のように書き換えて、図5のシミュレーション処理を実行することができる。式78は式25から右辺第4項の速度についての擾乱項を取り除き、右辺第2項のをに置き換えた計算式に相当する。式78の右辺第2項のうち、を含む部分が加速度に付加された擾乱項である。この場合、処理部401は、図5のステップ507において、式25の右辺第4項を計算する代わりに、式78の擾乱項を計算する。このとき、乱数γ1及びγ2の代わりに乱数γを用いて擾乱項を計算する。そして、図5のステップ508において、式25の代わりに式78を用いて、時刻n1における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。次に、各粒子の密度に擾乱を与える場合は、式25を式76に置き換え、式26を以下のように書き換えることができる。式79の右辺第3項の乱数δが密度に付加された擾乱項である。この場合、処理部401は、図5のシミュレーション処理の代わりに、図11に示すようなシミュレーション処理を実行する。図11のシミュレーション処理は、図5のシミュレーション処理からステップ503~506を除いた処理に相当し、ステップ1101~1108の処理は、図5のステップ501、502、及び507~512の処理に対応する。ただし、処理部401は、ステップ1103において、式25の右辺第4項を計算する代わりに、乱数δを発生し、その値を密度についての擾乱項の値として格納部402に格納する。乱数δとしては、例えば、-0.00002から0.00002までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、ステップ1104において、時刻nにおける各粒子の位置ベクトル、速度ベクトル、及び密度と密度の擾乱項の値を用いて、時刻n1における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。この場合、式24、76、79、及び27~55が計算式として用いられる。次に、各粒子の圧力に擾乱を与える場合は、式25を式76に置き換え、圧力の中間値を求める式55を以下のように書き換えて、図11のシミュレーション処理を実行することができる。式80の右辺第2項の乱数δpが圧力に付加された擾乱項である。この場合、処理部401は、ステップ1103において、乱数δpを発生し、その値を圧力についての擾乱項の値として格納部402に格納する。乱数δpとしては、例えば、-0.00002から0.00002までの一様乱数が用いられるが、これ以外の乱数を用いても構わない。そして、ステップ1104において、時刻nにおける各粒子の位置ベクトル、速度ベクトル、及び密度と圧力の擾乱項の値を用いて、時刻n1における処理対象粒子の位置ベクトル、速度ベクトル、及び密度を計算する。この場合、式24、76、26~54、及び80が計算式として用いられる。図4のシミュレーション装置は、水面波の伝播等の波動解析に限らず、流体、弾性体等の連続体や粉粒体の運動を解析する他のシミュレーション処理にも利用できる。また、シミュレーション処理で用いられる計算式は式24~80に限られるものではなく、各粒子の物理量に擾乱を与える他の計算式を用いても構わない。図4のシミュレーション装置は、例えば、図12に示すような情報処理装置コンピュータを用いて実現することが可能である。図12の情報処理装置は、CentralProcessingUnitCPU1201、メモリ1202、入力装置1203、出力装置1204、外部記憶装置1205、媒体駆動装置1206、及びネットワーク接続装置1207を備える。これらはバス1208により互いに接続されている。メモリ1202は、例えば、ReadOnlyMemoryROM、RandomAccessMemoryRAM、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、CPU1201は、メモリ1202を利用してプログラムを実行することにより、図5又は図11のシミュレーション処理を行う。メモリ1202は、図4の格納部402としても使用できる。入力装置1203は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置1204は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置1204は、図4の出力部403としても使用できる。外部記憶装置1205は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置1205には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置1205にプログラム及びデータを格納しておき、それらをメモリ1202にロードして使用することができる。媒体駆動装置1206は、可搬型記録媒体1209を駆動し、その記録内容にアクセスする。可搬型記録媒体1209は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体1209には、CompactDiskReadOnlyMemoryCD-ROM、DigitalVersatileDiskDVD、UniversalSerialBusUSBメモリ等も含まれる。オペレータは、この可搬型記録媒体1209にプログラム及びデータを格納しておき、それらをメモリ1202にロードして使用することができる。このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ1202、外部記憶装置1205、及び可搬型記録媒体1209のような、物理的な非一時的な記録媒体が含まれる。ネットワーク接続装置1207は、通信ネットワーク103に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置1207を介して受け取り、それらをメモリ1202にロードして使用することができる。ネットワーク接続装置1207は、図4の出力部403としても使用できる。開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。"