"連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション方法において、入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出し、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対する時間積分を行い、単位時間経過後の各粒子の速度を算出し、前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出し、前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出し、前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行い、前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出し、前記単位時間経過後の各粒子の速度を基に時間による積分を行い、前記単位時間経過後の各粒子の位置を算出し、前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記粒子の全てについて前記初期状態から所定時間が経過するまで繰返し、各単位時間経過後における速度、密度、圧力及び位置を取得することを特徴とするシミュレーション方法。前記シミュレーション方法において、前記カーネル関数は、全空間で積分すると1になり、前記粒子のうち、一つの粒子と影響を受ける他の粒子との間の距離が0に近づくと極限値がディラックのデルタ関数になり、さらに1階以上微分可能であることを特徴とする請求項1記載のシミュレーション方法。前記シミュレーション方法において、前記カーネル関数は、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを正規分布として表現するガウシアン関数であることを特徴とする請求項1又は請求項2記載のシミュレーション方法。前記シミュレーション方法において、前記カーネル関数は、有限の領域外で0になる5次のスプライン関数であることを特徴とする請求項1又は請求項2記載のシミュレーション方法。前記シミュレーション方法において、前記離散化された加速度運動方程式は、圧力勾配に応じた加速度を表すとともに、密度が大きいほど加速度が小さくなるように加速度を表すことを特徴とする請求項1~請求項4のいずれか一項に記載のシミュレーション方法。前記シミュレーション方法において、前記離散化された加速度運動方程式は、ハミルトン方程式を満たすことを特徴とする請求項5記載のシミュレーション方法。前記シミュレーション方法において、前記離散化された加速度運動方程式は、各粒子について対称形となるように調整用パラメータを用いて前記加速度運動方程式を変形したことを特徴とする請求項5又は請求項6記載のシミュレーション方法。前記シミュレーション方法において、前記調整用パラメータが1であることを特徴とする請求項7記載のシミュレーション方法。前記シミュレーション方法において、前記調整用パラメータが0であることを特徴とする請求項7記載のシミュレーション方法。連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション装置において、前記各粒子のうち、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された加速度運動方程式を用いて、各粒子の加速度と境界面からの反発力とを算出し、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対する時間積分を行うことで、単位時間経過後の各粒子の速度を算出する速度算出部と、前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて求めた各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出する密度算出部と、前記密度算出部が各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行う前記平滑化部と、前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出する圧力算出部と、前記単位時間経過後の各粒子の速度を基に時間による積分を行うことで、前記単位時間経過後の各粒子の位置を算出する位置算出部と、前記速度算出部、前記密度算出部、前記圧力算出部、前記平滑化部及び前記位置算出部に対して、初期状態の速度、密度、圧力及び位置の初期値を与え、前記初期状態から所定時間経過後までの各単位時間経過後における速度、密度、圧力及び位置の算出を行わせ、算出された各単位時間経過後における速度、密度、圧力及び位置を取得する制御部と、を備えたことを特徴とするシミュレーション装置。連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーションプログラムにおいて、コンピュータに、入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、前記各粒子のうち、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出させ、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対する時間積分を行うことで、単位時間経過後の各粒子の速度を算出させ、前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出させ、前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出させ、前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行わせ、前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出させ、前記単位時間経過後の各粒子の速度を基に時間による積分を行うことで、前記単位時間経過後の各粒子の位置を算出させ、前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記各粒子の全てについて前記初期状態から所定時間が経過するまで繰り返させ、各単位時間経過後における速度、密度、圧力及び位置を取得させることを特徴とするシミュレーションプログラム。しかしながら、従来の粒子法を用いたシミュレーションの場合、表面波の伝播を長時間扱うと、表面波が減衰してしまうという問題点がある。これについては、実際に標準的なSPH法で波の伝播解析を行ったところ、一定の時間を過ぎると波高の大きな減衰や発散が生じてしまい、造波試験を実施できないことが分かった。このように、従来の粒子法を用いたシミュレーションでは、現実との乖離が大きくなってしまうため、現実の水槽内に設置した造波板により造波を行う水理実験を減らすことができなかった。開示の技術は、上記に鑑みてなされたものであって、波高が減衰しない現実に近い解析を行うことができ、水理実験の実施回数を減らすことができるシミュレーション方法、シミュレーション装置及びシミュレーションプログラムを提供することを目的とする。本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムは、一つの態様において、連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション方法において、入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出し、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対する時間積分を行い、単位時間経過後の各粒子の速度を算出し、前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出し、前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出し、前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行い、前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出し、前記単位時間経過後の各粒子の速度を基に時間による積分を行い、前記単位時間経過後の各粒子の位置を算出し、前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記粒子の全てについて前記初期状態から所定時間が経過するまで繰返し、各単位時間経過後における速度、密度、圧力及び位置を取得する。本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムの一つの態様によれば、波高が減衰しない現実に近い解析を行うことができ、水理実験の実施回数を減らすことができるという効果を奏する。以下に、本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムの実施例を図面に基づいて詳細に説明する。なお、以下の実施例により本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムが限定されるものではない。図1は、実施例に係るシミュレーション装置のブロック図である。図1に示すように、本実施例に係るシミュレーション装置は、ユーザインタフェース1、制御部2、結果記憶部3、速度算出部4、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8を有している。ユーザインタフェース1は、シミュレーション装置を使用するユーザからの数値の入力、ユーザへのシミュレーション結果の出力を行なうインタフェース装置であり、具体的にはキーボード等の入力装置及び表示装置等の出力装置である。制御部2は、初期状態における各粒子の質量、速度、密度、圧力及び位置、すなわち各粒子の質力、並びに、速度、密度、圧力及び位置の初期値の入力をユーザインタフェース1から受ける。そして、制御部2は、受信した各粒子の質量、速度、密度、圧力及び位置の初期値を速度算出部4へ出力する。また、制御部2は、受信した各粒子の質量及び密度の初期値を密度算出部5へ出力する。また、制御部2は、受信した各粒子の圧力の初期値を圧力算出部6へ出力する。また、制御部2は、受信した各粒子の質量を平滑化部8へ出力する。さらに、制御部2は、受信した各粒子の位置の初期値を位置算出部7へ出力する。また、制御部2は、シミュレーションの終了のタイミングの指定をユーザインタフェース1から受ける。本実施例では、シミュレーションの終了のタイミングの指定は、ステップ単位で行うものとし、ここではシミュレーションがLステップ後に終了するよう指定されたものとする。ここで、ステップについて説明する。1ステップとは、本実施例ではΔt時間経過した後の状態をいう。すなわち、n1≦n≦Lステップ後とは、nΔt時間経過した後のことを指す。以下では、初期状態を1ステップ目とし、nステップ後のステップのことをステップnと言う。すなわち、初期状態は1回目のステップといえる。以下では、既に各物理量の算出が終わった最新のステップをステップnとし、次に物理量を計算するステップをステップn1として説明する。ここで、本実施例では、ステップ単位で終了のタイミングを指定したが、これは他の方法でも良く、例えば、初期状態から終了までの時間を指定してもよい。また、制御部2は、ユーザインタフェース1より入力を受けた各粒子の位置の初期値から、各粒子の影響範囲内にある他の粒子のリストを作成する。以下では、各粒子の影響範囲内にある他の粒子を近傍粒子と呼ぶ。また、近傍粒子のリストを近傍粒子リストと呼ぶ。図2は、影響範囲及び近傍粒子について説明するための図である。ここで、図2を参照して、各粒子の影響範囲及び近傍粒子について説明する。ここでは、着目粒子と他の粒子との距離が2h以内の場合に、相互に影響を与え合うものとして説明する。ここで、影響とは、例えば、粒子が移動する場合に、他の粒子へ力を加えるなどのことを言う。ここでは、図2の粒子101を例に説明する。粒子101の影響範囲102は、粒子101を中心として半径2hの円の中に含まれる範囲である。この半径2hは、粒子101の影響半径である。そして、影響範囲102の中に含まれる他の粒子が、粒子101の近傍粒子である。図2では、例えば、粒子103のように、斜線で表されている粒子が近傍粒子である。すなわち、粒子101が影響を与えることになる他の粒子は、斜線で表される近傍粒子ということになる。逆に、斜線で表される近傍粒子からのみ、粒子101は影響を受けることになる。ここで、近傍粒子が数十個程度含まれるように半径2hを設定すれば、着目した粒子と他の粒子間の相互作用を解析するには十分であることが経験的に知られている。また、粒子101の位置を表す位置ベクトルをraとし、粒子103の位置を表す位置ベクトルをrbとすると、粒子101と粒子103との間の距離は|ra-rb|となる。すなわち、|ra-rb|≦2hを満たす場合に、粒子103は粒子101の近傍粒子とされる。制御部2は、例えば図2の粒子101に対する近傍リストを作る場合、まず、影響範囲102に含まれる粒子103のような斜線で表される各粒子を粒子101の近傍粒子として抽出する。例えば、制御部2は、着目粒子である粒子101の位置と他の粒子の位置から粒子101と各他の粒子との距離を求め、その距離が2h以下となる粒子を近傍粒子として抽出する。図1に戻って、制御部2は、近傍粒子の抽出を全ての粒子に対して行う。そして、制御部2は、各着目粒子と抽出した各着目粒子の近傍粒子とが対応する近傍粒子リストを作成する。そして、制御部2は、作成した近傍粒子リストを速度算出部4、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8へ出力する。また、制御部2は、ステップn1における各着目粒子の速度の入力を速度算出部4から受ける。また、制御部2は、ステップn1における各着目粒子の密度及び圧力の入力を平滑化部8から受ける。さらに、制御部2は、ステップn1における各着目粒子の位置の入力を位置算出部7から受ける。ここで、1ステップとは、上述したように、Δt時間経過した後の状態である。そこで、例えば、ステップn1における速度とは、ステップnにおける速度からΔt時間後の速度である。すなわち、制御部2は、初期状態からΔt時間経過した後の、各着目粒子の速度、密度、圧力及び位置の入力を受ける。次に、制御部2は、初期状態から2Δt時間経過した後の、各着目粒子の速度、密度、圧力及び位置の入力を受ける。このように、制御部2は、既にステップnにおける速度、密度、圧力及び位置を受信している場合には、次にステップn1における各着目粒子の速度、密度、圧力及び位置を受信する。そして、制御部2は、指定されたシミュレーションの終了のタイミングまで、各ステップにおける各着目粒子の速度、密度、圧力及び位置を順次受信する。さらに、制御部2は、ステップn1における各着目粒子の最新の速度、密度、圧力及び位置から1ステップ後の情報の入力を受けると、受信したステップn1における各情報を結果記憶部3へ記憶させる。また、制御部2は、ステップn1における各着目粒子の密度、圧力及び位置を、次の計算に用いる各着目粒子の密度、圧力及び位置として、速度算出部4へ出力する。また、制御部2は、ステップn1における各着目粒子の密度を次の計算に用いる密度として密度算出部5へ出力する。また、制御部2は、ステップn1における各着目粒子の圧力を次の計算に用いる現在の圧力として圧力算出部6へ出力する。また、制御部2は、ステップn1における各着目粒子の速度を位置算出部7へ出力する。さらに、制御部2は、受信したステップn1における各着目粒子の位置を用いてステップn1における各着目粒子の近傍粒子を求め、近傍粒子リストを作成する。そして、制御部2は、作成した近傍粒子リストを速度算出部4、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8へ出力する。このように、制御部2は、ステップ毎の速度、密度、圧力及び位置、並びに近傍粒子リストを各部へ出力する。そして、制御部2は、現在のデータを用いて1ステップ後のデータを算出するよう、速度算出部4、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8に対して指示し、ステップLまで繰り返させる。結果記憶部3は、メモリやハードディスクなどの記憶装置である。速度算出部4は、例えば、カーネル関数として、次の式1で表される5次のスプライン関数を記憶している。ここで、αDは、2次元の計算を行なう場合にはαD=7/4πh2という値であり、3次元の計算を行なう場合にはαD=21/16πh3という値である。そして、rは、カーネル関数を用いた計算の対象としている2つの粒子の距離を表し、hは影響範囲を表しており、qはq=r/hと定義される。式1は、ある粒子における、他の粒子へ影響を及ぼす寄与の大きさを表している。そして、式1は、粒子同士が離れると一方の粒子が他方の粒子に影響を及ぼす寄与の大きさが小さくなることを表わしている。そして、式1は、イ全空間で積分すると1になる、ロ影響範囲が0に近づくと、すなわち、一方の粒子と影響を受ける他の粒子との間の距離が0に近づくと極限値がディラックのデルタ関数になる、ハ1階以上微分可能である、という3つの条件を満たしている。さらに、式1は、qが0より小さい又は2より大きい場合の値が0になることから、有限の領域の外で値が0となる性質を有している。この性質は、コンパクトサポートと呼ばれる場合がある。ここで、圧縮場では、d2W/dq2が負の領域でtensileinstabilityと呼ばれる数値不安定が生じやすい。そこで、カーネル関数としては、d2W/dq2が負となる領域が比較的狭く、不安定が発生しにくいと考えられる式1で表される5次のスプライン関数を用いることが特に好ましい。ただし、数値的不安定さを抑える要求に応じて、カーネル関数は他のものを用いることもできる。例えば、他のカーネル関数としては、次の式2~4のようなものを用いてもよい。ここで、q=r/hである。また、αDは、2次元の計算を行なう場合にはαD=1/πh2という値であり、3次元の計算を行なう場合にはαD=1/π3/2h3という値である。式2は、ガウシアン関数である。ここで、αDは、2次元の計算を行なう場合にはαD=2/πh2という値であり、3次元の計算を行う場合にはαD=5/4πh3という値である。ここで、αDは、2次元の計算を行なう場合にはαD=10/7πh2という値であり、3次元の計算を行なう場合にはαD=1/πh3という値である。式2~4はいずれも、ある粒子における、他の粒子へ影響を及ぼす寄与の大きさを表している。そして、式2~4はいずれも、粒子同士が離れると一方の粒子が他方の粒子に影響を及ぼす寄与の大きさが小さくなることを表わしている。そして、式2~4はいずれも、イ、ロ及びハの条件を満たしている。特に、式2のガウシアン関数は、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを正規分布として表現している。そして、式2を用いた場合、後で詳述する、速度、密度、圧力及び位置の算出において、数学的には厳密な値が求められる。しかし、式2を用いた場合、Wr,h=0になる点を求めることができない。この場合、影響の範囲がはっきりしなくなることから、影響半径をどの値に取ったらいいのかが不明確になってしまう。そこで、本実施例では、精度とのバランスを考慮してカーネル関数として式1を用いている。さらに、速度算出部4は、例えば、粒子の加速度を表す離散化された加速度運動方程式として次の式5を記憶している。また、添字aは、粒子aについての情報であることを表しており、添字bは粒子bについての情報であることを表している以下、添字のa及びbは同様であるものとする。。すなわち、式5の運動方程式は、粒子aの加速度を表している。ここで、粒子aとは、全ての粒子の中から任意に抽出した粒子である。そして、粒子bとは、粒子aの近傍粒子である。すなわち、和の範囲は粒子aの近傍粒子について亘ることになる。そして、ρaは粒子aの密度であり、ρbは粒子bの密度である。そして、ρabここでは、ρabの上部の横棒の記号を省いて記載している。はρaとρbとの算術平均を表している。また、paは粒子aにかかる圧力であり、pbは粒子bにかかる圧力である。また、Wabは粒子aと粒子bにおけるカーネル関数であり、本実施例では式1で表される。さらに、∇aは、粒子aの位置におけるベクトル微分演算子を表している。また、mbは粒子bの質量である。また、gの項ここで、gはベクトルを表しているものとするは外力項を表している。式5の右辺の第1項は圧力勾配を表している。また、式5の右辺の第2項は人工粘性を表している。さらに、式5の第3項は重力を表している。さらに、式5の右辺の第2項の人工粘性項は、粒子aと粒子bとが互いに近づく場合に限り粒子同士の衝突を妨げる向きに力を及ぼすための項である。人工粘性項におけるβavは、人工粘性の強さを調整するパラメータである。また、人工粘性項におけるφは次の式6で表される。ここで、rba=rb-raここで、rba、rb及びraはそれぞれベクトルを表しているものとする。であり、ra、rbはそれぞれ、粒子aの位置ベクトルと粒子bの位置ベクトルとを表している。すなわち、rbaは粒子aと粒子bとの距離である。vba=vb-vaここで、vba、vb及びvaはそれぞれベクトルを表しているものとする。である。また、ηは、分母が発散しないためのパラメータである。ここで、式5の導出について説明する。流体の運動方程式は次の式7で表される。ここで、vここで、vはベクトルを表しているものとする。は速度ベクトルを表し、ρは流体の密度を表しており、pは流体の圧力を表しており、gの項ここで、gはベクトルを表しているものとする。は外力項を表している。また、Δは、ラプラス演算子を表している。そして、式7の右辺の第1項は圧力勾配を表し、第2項は粘性応力による力を表し、第3項は外力を表している。ここで、圧力勾配を表す-∇P/ρは、圧力Pの大きいところから小さなところに向かって加速度が働き、さらに密度ρが大きいほど加速度が小さくなることを表している。さらに、∇P/ρは、次の式8のように変形できる。ここで、式8におけるσは、次に行うように式7に対して粒子法を用いて離散化を行った場合に、各粒子について対称な形にする場合を一般化するための調整用のパラメータである。次に、式7の圧力勾配を、式8を用いて変形したうえで、SPH法を用いて離散化を行う。ここで、SPH法では、粒子aの物理量Aは、次の式9のように離散化表現される。物理量Aは例えば、圧力や密度などである。そこで、式7の圧力勾配を式8を用いて変形した式を、式9を用いて離散化表現することで、次の式10のように表すことができる。ここで、式10では、式7における粘性応力の項が省かれている。これは、本実施例で解析の対象としている波などの粘性が低い流体においては、粘性応力が無視してよいほどの力であるためであり、離散化後の式から除いている。この式10で分かるように、運動方程式が、粒子a及び粒子bのそれぞれに対して対称な式として表されている。すなわち、式10を用いることで、粒子aと粒子bとの間で相互に働く力が等しくなる。このため、式10を用いた解析においては、運動量保存側やエネルギー保存則を守ることが容易となる。そして、調整用のパラメータをσ=1とすると、式5が求まる。ここで、本実施例では、調整用パラメータをσ=1とした場合としたが、σの値は他の値をとってもよく、例えば、従来粒子法において一般的に用いられてきたσ=0を用いてもよい。さらに、後述するように、本実施例では時間積分をするにあたり、エネルギー散逸が比較的生じにくいSymplecticEuler法を用いる。そして、SymplecticEuler法では、運動方程式は次の式11及び式12のHamilton方程式で記述される。ここで、ea=mavaここで、eaはベクトルを表しているものとする。は、粒子aの運動量ベクトルである。また、∇rは位置に関する微分演算子であり、∇eは運動量に関する微分演算子である。また、Hは全エネルギHamiltonianである。SymplecticEuler法は、Hamilton方程式で記述される系を時間積分するスキームであり、次の式13及び式14で表されるステップで実行される。ここで、添え字nは時間ステップ数を表し、Δtは時間増分を表している。式13及び式14は両辺にn1ステップの値を含むので、一般には陰的公式である。しかし、HamiltonianがHe,r=TeUrとeのみに依存する項と、rのみに依存する項の和として表されるときには、式13は次の式15として表せ、式14は次の式16として表される。式15及び式16は、右辺に未知の値を含まないため、陽的に解くことができる。これを、粒子法に適用すると、Hamiltonianは、次の式17のように定義できる。ここで、Nは総粒子数である。また、式17の第1項は運動エネルギーを表し、第2項は粒子間ポテンシャルを表し、第3項は外力ポテンシャルを表している。ただし、∇Ga=gここで、gはベクトルを表しているものとする。である。そして、式17と定義した場合、式12から式5の運動方程式を得ることができる。さらに、式13からdra/dt=vaとなる。すなわち、式5は、Hamiltonsystemを形成していることがわかり、式5は、Symplectic法を用いる運動方程式として適切であるといえる。ここで、従来は運動方程式として、例えば、次の式18などが用いられてきた。ここで、Πabは、粘性応力を表す粘性項であり、次の式19のように表される。ここで、μは粘性係数であり、水の場合はμ=10-3Pasである。また、ξは物理粘性を調整するパラメータであり、例えば、ξ=4.96333などである。また、式19の右辺の第2項は物理量の不連続に起因する数値誤差を是正する人工粘性項である。そして、式18で表される運動方程式を用いてシミュレーションを行った場合、波高に現実よりも大きな減衰などが生じてしまう。これに対して、本実施例における運動方程式である式5から人工粘性を省いた運動方程式は、Hamilton方程式で表されるエネルギー保存則から導出される。すなわち、式5の運動方程式はエネルギー保存則や運動量保存則を的確に反映することができ、波高の減衰などを抑えることができる。したがって、本実施例における運動方程式である式5を用いて解析を行うことで、式18のような従来の式を用いた場合と比較して、波動伝播のような安定性が要求される解析をより適切に行うことができる。図1に戻って、速度算出部4の説明を続ける。速度算出部4は、運動方程式5を用いて加速度dv/dtを算出する。次に、速度算出部4は、境界からの反発力を次の式20を用いて求める。ここで、図3を参照して、式20について説明する。図3は、実施例で用いた境界からの反発力を説明するための図である。ここでは、図3における粒子200に加えられる境界からの反発力について説明する。図3に示すように、dは境界面201からの所定の距離である。そして、境界層202、境界面からdの幅の層を示している。また、sは粒子200の境界層202への侵入距離である。本実施例では、境界面201から幅dの境界層202を考え、境界層202内に粒子200が突入した場合に、粒子200に対して境界面201から粒子200の向きに反発加速度が加わるものとする。反発力fnは、弾性ばねによるものと同じ形式でfn=ksとして与える。ばね剛性係数kは、音速c0の粒子が有する運動エネルギーmc02/2となり、粒子が境界層202へdの長さだけ進入したときの弾性歪エネルギーkd2/2が同程度となるように定める。すなわち、連続体の運動の解析を行なうシミュレーションにおいては、粒子200は音速を超えることはないと考えられるので、音速で粒子200が進入した場合に境界面201においてエネルギーがなくなるように設定しておけば、粒子200は境界面を突き抜けることはなくなる。このため、βを1程度のパラメータとする。このようにβを設定することで、式20が導かれる。式20で表される反発力を用いることで、境界壁からの粒子のすり抜けを防止し、導入する摩擦やfreeslip条件の正確性を向上でき、さらに、境界で物理的に無意味なエネルギー消費の発生を抑えることができる。また、式20は簡単な式で表されているので、式20で表される反発力を用いた場合には、計算負荷を抑えることができる。ここで、波高の現実よりも大きな減衰などの各種問題が発生しにくい式として本実施例では式20を用いたが、この反発力は他の式を用いてもよい。例えば、反発力が過大になる可能性を考慮して、粒子が境界面に近づいている場合には、弾性ばねを考慮し、逆に粒子が境界面から遠ざかる場合には、弾性ばねを考慮しないように反発力を設定してもよい。このような反発力は一方向弾性ばね法と呼ばれることがある。エネルギー保存の観点からは粒子と境界面の相対速度に関係なく弾性ばねを考慮する式20のような反発力が好ましいが、境界での不自然な反発を避けるために一方向弾性ばね法をもちいてもよい。そして、速度算出部4は、運動方程式5を用いて算出した加速度dv/dtと式20を用いて求めた反発力とを加えてステップnにおける加速度Fnaここで、Fnaは加速度ベクトルを表しているものとする。を算出する。次に、速度算出部4は、算出した加速度Fnaを用いて次の式21で表される時間積分を行い、ステップn1における速度を求める。ここで、vnはステップnにおける速度を表している。また、時間増分Δtは、本実施例では次の式22とした。以下の説明での時間増分Δtは全て式22とした。ここで、αΔtは発散防止用の係数で、本実施例では一例として0.1とした。速度算出部4は、以上に説明した、運動方程式を用いて算出した加速度に境界面からの反発力を加えてステップnにおける加速度を算出し、この加速度を用いてステップn1の速度を算出するという処理を、全ての粒子に対して行う。そして、速度算出部4は、算出したステップn1における各粒子の速度を制御部2、密度算出部5及び位置算出部7へ出力する。密度算出部5は、式1で表されるカーネル関数を記憶している。さらに、密度算出部5は、次の式23で表される粒子の密度の時間変化を表す密度変動の連続方程式を離散化した式を記憶している。ここで、式23の導出について説明する。離散化しない場合の流体の密度の時間変化の連続方程式は次の式24で表される。この式24に対して、上述した式9で表されるSPH法における粒子aの物理量Aを用いて離散化することで式23が得られる。この連続方程式は、各粒子の質量を重みとして、重みを付けたカーネル関数で表される各粒子が有する密度を重ね合わせたものである。図1に戻って、密度算出部5は、各粒子の質量、ステップnにおける各粒子の位置及びステップnにおける近傍粒子リストの入力を制御部2から受ける。また、密度算出部5は、ステップn1における速度を速度算出部4から受ける。そして、密度算出部5は、受信した各粒子の質量、ステップn1における各粒子の速度、ステップnにおける各粒子の位置及びステップnにおける近傍粒子リストを式23に代入する。そして、密度算出部5は、ステップn1における密度変動である密度時間微分dρa/dt=Dan1を算出する。次に、密度算出部5は、前回のステップで算出したステップnにおける密度及び算出したステップn1における密度変動を用いて次の式25で表される時間積分を行い、ステップn1の密度を求める。密度算出部5は、密度変動の算出及び算出した密度変動を用いた密度の算出を全ての粒子に対して行う。そして、密度算出部5は、算出したステップn1における各粒子の密度を圧力算出部6及び平滑化部8へ出力する。平滑化部8は、ステップn1における各粒子の密度の入力を密度算出部5から受ける。平滑化部8は、カウンタを有している。そして、平滑化部8は、ステップn1における各粒子の密度及び圧力の入力の回数をカウンタによりカウントし、20回に一回例えば、20の倍数のステップ毎、受信したステップn1における各粒子の密度及び圧力に対して以下の平滑化を行う。この20回に一回のステップが所定ステップの一例にあたる。ここで、本実施例では、20回に一回の割合としたが、これは他の割合としてもよく、さらに、定期的でなくてもよい。この平滑化を行うステップは、密度及び圧力の算出による表面における粒子のみだれの状態を考慮して決定されることが好ましい。平滑化部8は、受信したステップn1における密度に対して、次の式26を用いて圧力の平滑化を行う。ここで、ρanewは、平滑化処理後の密度である。そして、平滑化部8は、平滑化処理後の密度をステップn1における密度として圧力算出部6へ出力する。また、平滑化部8は、平滑化を行わなかった場合は、密度算出部5から受信したステップn1における各粒子の密度を制御部2に出力し、平滑化を行った場合には、平滑化後のステップn1の各粒子の密度を制御部2に出力する。後述するように、密度の平滑化に伴い圧力に対しても平滑化が行われる。このように、密度及び圧力の平滑化を行うことで、圧力振動を抑制することができる。ここで、本実施例では、密度算出部5が算出したステップn1の各粒子の密度は常に平滑化部8に渡されているが、これは次のような方法でもよい。例えば、密度算出部5は、今回処理を行なうステップが20回に一回のステップに該当するか否かを判定し、今回処理を行なうステップが20回に一回のステップの場合には、密度算出部5は、平滑化部8へステップn1の各粒子の密度を出力する。これに対して、今回処理を行うステップが20回に一回のステップ以外の場合には、密度算出部5は、制御部2へステップn1の各粒子の密度を出力する。このような方法を取ることで、平滑化部8における処理を軽減できる。圧力算出部6は、ステップn1における各粒子の密度の入力を密度算出部5から受ける。また、20回に一回のステップにおいて、平滑化されたステップn1における各粒子の密度の入力を平滑化部8から受ける。そして、圧力算出部6は、平滑化部8からステップn1における各粒子の密度の入力を受けた場合、圧力算出部6から受信したステップn1における各粒子の密度を平滑化部8からステップn1における各粒子の密度に更新する。そして、圧力算出部6は、次の式27で表される状態方程式にステップn1における密度を用いて、ステップn1における圧力を算出する。ここで、ρ0は平均密度を表す。本実施例では、一例として、ρ0=1g/cm3、c0=50m/secとして計算を行った。ここで、状態方程式により、密度の大きな領域は圧力が大きくなり、その周囲へ向かって圧力勾配が働くため、式27を用いた結果として密度の凹凸がならされる。このように、密度の凹凸がならされることで、密度が均一に近くなり、より数値計算の正確性を向上することができる。また、密度の揺らぎの程度δp=ρ/ρ0-1は、流体の典型的な流速v0と音速を用いて表されるv02/c02以下と見積もることができる。そして、本実施例における流速は、音速50m/secに比べて、十分小さい値であるので、密度の揺らぎも十分小さい値となる。そして、本実施例では流体の体積変化が0.1程度であったことから、密度変化も同様に0.1程度であり、数値計算上問題となる値ではない。圧力算出部6は、ステップn1における圧力の算出を全ての粒子について行う。そして、圧力算出部6は、算出したステップn1における各粒子の圧力を制御部2へ出力する。ここで、本実施例では、圧力算出部6は、常に密度算出部5からステップn1における各粒子の密度を受信し、平滑化部8から平滑化が施された密度を受信した場合に、ステップn1における各粒子の密度を平滑化後の値に更新している。しかし、次のような方法でもよい。すなわち、密度算出部5が、今回処理を行なうステップが20回に一回のステップに該当するか否かを判定し、今回処理を行なうステップが20回に一回のステップであれば、密度算出部5は、圧力算出部6にはステップn1における各粒子の密度を出力せず、平滑化部8にのみ出力する。そして、圧力算出部6は、平滑化部8からのみステップn1における各粒子の密度を受信する。一方、密度算出部5が、今回処理を行なうステップが20回に一回のステップ以外のステップと判定した場合、密度算出部5は、圧力算出部6と平滑化部8の双方にステップn1における各粒子の密度を出力する。そして、圧力算出部6は、密度算出部5及び平滑化部8の双方からステップn1における各粒子の密度を受信する。このような方法を用いることで、圧力算出部6の処理を軽減することができる。位置算出部7は、各粒子の位置の初期値の入力を制御部2から受ける。さらに、位置算出部7は、n1ステップにおける各粒子の速度の入力を速度算出部4から受ける。そして、位置算出部7は、受信したn1ステップにおける速度を用いて、次の式28で表される時間積分を行い、ステップn1における粒子の位置を求める。位置算出部7は、ステップn1における粒子の位置の算出を全ての粒子に対して行う。そして、位置算出部7は、算出したステップn1における各粒子の位置を制御部2へ出力する。このように、時間積分により次のステップの速度をまず求め、次に、求めた次のステップの速度を用いて時間積分を行い、次のステップの密度、圧力及び位置を算出する時間積分法はSymplecticEuler法と呼ばれることがある。シミュレーションにおいて時間積分を行なう手順を時間積分スキームと称することがあり、このSymplecticEuler法は、Hamilton方程式で記述される系に対する時間積分スキームの一種である。以上に説明したように、本実施例にかかるシミュレーションでは、まず速度を時間微分し、次のステップの速度を求め、その後に次のステップの密度及び圧力を求め、最後に次のステップの位置を求めることになる。すなわち、本実施例に係るシミュレーションにおいては、ステップn1における密度及びステップn1の位置のいずれにおいても、それらを求める場合にステップnの位置を用いている。すなわち、密度及び位置を求めるタイミングを一致することができる。これに対して、従来一般的に用いられている、leapfrog法は、計算の最初に位置座標を半ステップ進め、以降、位置座標と速度を半ステップずつずらしたまま時間積分を行う方法である。この方法では、密度と位置を求めるタイミングがずれてしまう。この点、密度の算出においては、位置が大きく影響する。そのため、密度と位置とは同じタイミングで算出することが数値計算を適正に行ううえで好ましい。すなわち、本実施例で用いたSymplecticEuler法による時間微分を行うことで、エネルギーの保存を数値計算上正確に保存することができる。次に、図4を参照して、本実施例に係るシミュレーションの処理の流れについて説明する。図4は、実施例に係るシミュレーションの処理のフローチャートである。ここでは、平滑化部8が、平滑化処理を施す20回に1回のステップとして、20の倍数のステップ毎に平滑化を施す場合を例に説明する。制御部2は、各粒子の質量、並びに各粒子の速度、密度、圧力及び位置の初期値を含む入力データをユーザインタフェース1から取得するステップS1。そして、制御部2は、受信した各粒子の質量、速度、密度、圧力及び位置の初期値を速度算出部4へ出力する。また、制御部2は、受信した各粒子の質量及び密度の初期値を密度算出部5へ出力する。また、制御部2は、受信した各粒子の圧力の初期値を圧力算出部6へ出力する。また、制御部2は、受信した各粒子の質量を平滑化部8へ出力する。さらに、制御部2は、受信した各粒子の位置の初期値を位置算出部7へ出力する。制御部2は、受信した各粒子の位置を用いて近傍粒子リストを作成するステップS2。速度算出部4は、式5の運動方程式を用いて粒子の現在のステップにおける加速度を算出するステップS3。次に、速度算出部4は、粒子が境界近傍層に含まれる位置にある境界近傍の粒子か否かを判定するステップS4。速度算出部4が境界近傍の粒子で無いと判定した場合ステップS4否定、ステップS6へ進む。これに対して、境界近傍の粒子である場合ステップS4肯定、速度算出部4は、式20を用いて境界からの反発力を算出するステップS5。そして、速度算出部4は、現在のステップにおける加速度に反発力を加えた粒子の加速度を求めるステップS6。次に、速度算出部4は、現在のステップにおける速度及び算出した現在のステップにおける加速度を用いて式21の時間積分を行い、次のステップにおける粒子の速度を算出する。そして、制御部2は、速度算出部4によって次のステップにおける速度の算出が全ての粒子について行われたか否かを判定するステップS7。制御部2は、次のステップにおける速度の算出を行っていない粒子があると判定した場合ステップS7否定、速度算出部4に対してステップS3からステップS6までを繰り返させる。これに対して、全ての粒子について次のステップにおける速度の算出が完了している場合ステップS7肯定、密度算出部5は、式23の連続方程式を用いて次のステップにおける密度時間微分を算出するステップS8。そして、速度算出部4は、算出した次のステップにおける各粒子の速度を、制御部2、密度算出部5及び位置算出部7へ出力する。次に、密度算出部5は、現在のステップにおける密度及び次のステップにおける密度時間微分を用いて式25の時間積分を行い、次のステップにおける密度を算出するステップS9。そして、密度算出部5は、算出した密度を圧力算出部6及び平滑化部8へ出力する。圧力算出部6は、密度算出部5から受信した次のステップにおける密度を式27の状態方程式に用いて次のステップにおける圧力を算出するステップS10。位置算出部7は、速度算出部4から受信した次のステップにおける速度を用いて式28の時間積分を行い、次のステップの粒子の位置を算出するステップS11。平滑化部8は、次のステップにおける密度の入力を密度算出部5から受ける。そして、平滑化部8は、現在のステップが20の倍数のステップか否かを判定するステップS12。平滑化部8は、現在のステップが20の倍数のステップで無いと判定した場合ステップS12否定、密度算出部5から受信した次のステップを制御部2へ出力し、ステップS15へ進む。これに対して、現在のステップが20の倍数のステップの場合ステップS12肯定、平滑化部8は、受信した次のステップにおける密度に式26を用いて平滑化を行うステップS13。そして、平滑化部8は、平滑化を施した次のステップの密度を制御部2及び圧力算出部6へ出力する。そして、圧力算出部6は、平滑化部8から受信した平滑化が施された次のステップの密度を式27の状態方程式に用いて次のステップにおける平滑化された圧力を算出するステップS14。そして、圧力算出部6は、平滑化された圧力に次のステップの圧力を更新する。その後、圧力算出部6は、算出した次のステップにおける圧力を制御部2へ出力する。そして、制御部2は、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8によって次のステップにおける密度、圧力及び位置の算出が全ての粒子について行われたか否かを判定するステップS15。制御部2は、次のステップにおける密度、圧力及び位置の算出を行っていない粒子があると判定した場合ステップS15否定、密度算出部5、圧力算出部6、位置算出部7及び平滑化部8に対してステップS8からステップS14までを繰り返させる。これに対して、全ての粒子について次のステップにおける密度、圧力及び位置の算出が完了している場合ステップS15肯定、制御部2は、今回のステップの計算を完了し、計算対象の時刻を1ステップ進めるステップS16。そして、制御部2は、今回のステップの各粒子の速度、密度、圧力及び位置を含む計算結果データをユーザインタフェース1に出力するステップS17。さらに、制御部2は、今回のステップの結果データを結果記憶部3に格納する。そして、制御部2は、指定された最終ステップが完了したか否かを判定するステップS18。制御部2は、全ステップが完了していないと判定した場合ステップS18否定、ステップS2へ戻る。これに対して、最終ステップが完了している場合ステップS18肯定、制御部2は、シミュレーションの処理を終了する。水理実験との比較次に、本実施例に係るシミュレーションの正確さを調べるため、水理実験との比較を行った場合の結果について説明する。図5は、水理実験に用いた造波水槽の概念図である。この水理実験では、垂直護岸概要条件を用いて実験を行った。そして、この水理実験では、図5に示す位置に、波高計301~306を配置し、波高記録を得た。図5の護岸401は模擬的な防波堤であり、斜面402は模擬的な水底から護岸401に向かう模擬的な斜面である。そして、造波板403が振幅運動を行うことによって、造波がなされる。そして、水槽の波高計301から護岸401までの距離を18.5mとし、水深を0.679mとした。また、斜面402は、10分の1の勾配を有している。そして、この検証では、図5と同様の水槽の波高計301から護岸401に亘る解析モデルを用いて、上述したシミュレーション方法により解析を行い、波高履歴及び越波量を図5で示される水槽を用いた水理実験の結果と比較した。ここで、シミュレーションに用いた造波板は、図5の水槽の造波板と同じ位置にあると想定して解析を行なっている。そして、波高計306において、水槽実験によって得られた波高履歴を再現するよう、造波板403の移動速度を無反射造波条件式の定常波項を省略した次の式29を用いて計算を行った。ここで、ηDは目標波形、η0は造波板前面の水位変動、ωは角周波数であり、xwallは造波板403の位置座標である。また、Bは次の式30を満たす。さらに、Hは静水面高さ、Kは波数である。Kは次の式31から求めることができる。ここで、gは重力加速度である。表面波伝播の過程のみを調べるためにまず、護岸を設置しない状態での水理実験、本実施例に係るシミュレーション方法の数値解析結果、及び従来手法による数値解析結果を比較する。実験条件は、波高20cm、周期2secである。数値解析における粒子数は12万1980体、初期粒子間隔は1cm、影響半径=2hは3cmである。従来手法による数値解析としては、次の点が本実施例に係るシミュレーション方法と異なるものを採用した。すなわち、従来手法による数値解析では、Aカーネル関数として、数式4で表されるcubicspline関数を用いる。B運動方程式として式18を用いる。ここで、式19において、βav=0.1、η=0.01hとした。C時間積分スキームとしてはleapfrog法を用いる。ここで、leapfrog法とは、計算の最初に位置座標を半ステップ進め、以降、位置座標と速度を半ステップずつずらしたまま時間積分を行う方法である。図6は、波高計302の位置における波高の解析結果と水理実験との比較を表す図である。図7は、波高計303の位置における波高の解析結果と水理実験との比較を表す図である。図8は、波高計304の位置における波高の解析結果と水理実験との比較を表す図である。図9は、波高計305の位置における波高の解析結果と水理実験との比較を表す図である。図10は、波高計306の位置における波高の解析結果と水理実験との比較を表す図である。図6~図10における実線511~551で表されるグラフが本実施例に係るシミュレーション方法を用いた数値解析結果である。図6~図10は、いずれも、横軸が時刻であり、縦軸が波高である。また、図6~図10における点線512~552で表されるグラフが水理実験による結果である。また、図6~図10における一点鎖線513~553で表されるグラフが本実施例に係るシミュレーション方法を用いた数値解析結果である。図6~図10から分かるように、従来手法においては、造波板403に近い波高計302では水理実験と波高がほぼ一致しているが、波の伝播距離が増すにつれて波高が大きく減衰している。そして、波高計306の位置では、従来手法で求められた波高が水理実験と大きく乖離している。例えば、図6の点線512と一点鎖線513とは近似している。しかし、図10の点線552と一点鎖線553とは大きなズレが生じている。これに対して、本実施例に係るシミュレーション方法を用いた数値解析の場合には、図6~図10に示すように、造波板403から遠くなり波の伝播距離が長くなっても、解析結果が水理実験とほぼ一致している。例えば、図6の点線512と実線511とは近似しており、さらに図10の点線552と実線551とも近似している。すなわち、従来手法と比較して本実施例に係るシミュレーション方法は、より水理実験に近い結果を得ることができる。次に、カーネル関数、運動方程式又は時間積分スキームによる、数値解析の正確性の向上に対する寄与の度合いを測るため、カーネル関数、運動方程式又は時間積分スキームのそれぞれを変更した場合の計算結果について説明する。図11は、実施例に係るシミュレーション方法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。図11では、カーネル関数、運動方程式又は時間積分スキームをそれぞれ前述した従来手法におけるA~Cに切り替えている。図11は、横軸が時刻であり、縦軸が波高である。ここで、太実線601が水理実験を示し、細実線602が本発明を示し、点線603がカーネル関数を切り替えた場合を示し、一点鎖線604が運動方程式を切り替えた場合を示し、二点鎖線605が時間積分スキームを切り替えた場合を示している。図11から分かるように、太実線601と細実線602とはもっとも近似している。そして、カーネル関数、運動方程式又は時間積分スキームのそれぞれを切り替えた場合、細実線602と比べて、太実線601から離れていく。その中で、カーネル関数を切り替えた点線603が最も大きく太実線601から乖離している。また、図12は、従来手法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。図12では、従来手法としてA~Cを満たす式を用いた場合とし、それらA~Cを、本実施例で用いた式1のカーネル関数、式5の運動方程式及びSymplecticEuler法を用いた時間積分に切り替えている。図12は、横軸が時刻であり、縦軸が波高である。ここで、太実線611が水理実験を示し、細実線612が本発明を示し、点線613がカーネル関数を切り替えた場合を示し、一点鎖線614が運動方程式を切り替えた場合を示し、二点鎖線615が時間積分スキームを切り替えた場合を示している。図12から分かるように、太実線611と細実線612とはもっとも乖離している。そして、カーネル関数、運動方程式又は時間積分スキームのそれぞれを切り替えた場合、細実線612と比べて、太実線611に近づいていく。その中で、カーネル関数を切り替えた点線613が最も太実線601に近似している。すなわち、カーネル関数として、式1で表される5次のスプライン関数を用いることが、シミュレーションの正確性の向上に最も効果が大きいといえる。さらに、護岸からの反射波の影響について調べるために直立護岸を設置した状態での水理実験と本実施例に係るシミュレーション方法による数値解析結果との比較を行う。実験条件として、瀬戸内条件及び外洋条件の2種類の条件を設定した。図13は、垂直護岸を設置した実験条件を示す図である。ここでは、瀬戸内条件及び外洋条件の条件として、図13に示すそれぞれの値を用いた。図14は、外洋条件を用いた場合の波高履歴を示す図である。また、図15は、瀬戸内条件を用いた場合の波高履歴を示す図である。図14及び図15のいずれも横軸が時刻であり、縦軸が波高である。そして、実線701及び実線711が本実施例に係るシミュレーション方法を用いた数値解析結果である。また、点線702及び点線712が水理実験による結果である。図14において、実線701と点線702とは波高、周期及び波形のいずれも十分に近似している。また、図15において、実線711と点線712とは波高、周期及び波形のいずれも十分に近似している。すなわち、本実施例に係るシミュレーション方法を用いた数値解析では、波高、周期及び波形とも、水理実験が十分に再現できているといえる。以上に説明したように、本実施例に係るシミュレーション方法を用いた場合と従来の粒子法を用いた数値解析とを比較した場合、後者では水理実験の結果と比較して波高が大きく減衰するが、前者では波高の減衰がほとんど無く、水理実験の波高とよく一致する。また、護岸前面の砕波を生じる非線形性の高い領域においても、本実施例に係るシミュレーション方法を用いた場合、波高が実測と良く一致する。これにより、護岸砕波過程における越波量や衝撃圧については造波水槽を用いた水理実験が行われるが、本実施例に係るシミュレーション方法による数値解析を行えば、水理実験をシミュレーションに置き換えることにより水理実験の回数の低減ができ、さらに、実海域の波浪に対する応答も直接評価できる。したがって、護岸設計コストの低減及び合理化に大きく寄与することができる。また、上記実施例で説明した各種の処理は、あらかじめ用意されたプログラムをコンピュータで実行することによって実現することができる。そこで、以下では、図16を用いて、図1に示したシミュレーション装置と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。図16は、シミュレーションプログラムを実行するコンピュータを示す図である。図16に示すように、コンピュータ1000は、キャッシュ1010と、CPUCentralProcessingUnit1020と、HDDHardDiskDrive1030と、RAMRandomAccessMemory1040及びROMReadOnlyMemory1050を有する。キャッシュ1010、CPU1020、HDD1030、RAM1040及びROM1050は、それぞれバスによって接続されている。HDD1040には、図1に示したシミュレーション装置と同様の機能を発揮する各種のシミュレーションプログラム1031が予め記憶されている。そして、CPU1020は、これらの携帯端末制御プログラム1031を読み出して実行する。これにより、図16に示すように、シミュレーションプログラム1031は、シミュレーションプロセス1021になる。なお、上記したシミュレーションプログラム1031については、必ずしもHDD1030に記憶させなくてもよい。例えば、コンピュータ1000に挿入されるフレキシブルディスク、CDCompactDisk-ROM、DVDDigitalVersatileDisc、光磁気ディスク、ICIntegratedCircuitカードなどの可搬用の物理媒体にシミュレーションプログラム1031を記憶させてもよい。または、コンピュータ1000の内外に備えられるハードディスクドライブHDDなどの固定用の物理媒体にシミュレーションプログラム1031を記憶させてもよい。または、公衆回線、インターネット、LANLocalAreaNetwork、WANWideAreaNetworkなどを介してコンピュータ1000に接続される他のコンピュータまたはサーバにシミュレーションプログラム1031を記憶させてもよい。そして、コンピュータ1000は、上述したフレキシブルディスク等から各プログラムを読み出して実行するようにしてもよい。"