第2章　数値シミュレーション

　数値シミュレーションと言えば，構造解析，振動解析，機構解析，気流解析など市販の専用ソフトを想像される方が多いようである． しかしそれだけではない．極端に思われるかもしれないが，電卓で検算することもシミュレーションである． 本章で理解して頂きたいことは「計算機で簡単に計算できることは面倒がらずにやろう」と言うことである． 「難しい．分からない．教えて欲しい」などと言っていないで．数値シミュレーションを試してみよう．やってみれば意外と簡単である．
　ここでは，簡単な数値シミュレーションの例を幾つか示す．数値シミュレーションが強力な助っ人であることを理解し，数値シミュレーションと仲良くなって欲しい．

2.1　確率的問題のシミュレーション利用

　統計や確率の問題は，乱数を使った数値シミュレーションに適している． 統計学や確率論を知らなくても，殆どのことはシミュレーションで予測可能である．

例題：サイコロ実験結果の確率分布

　普通のサイコロをm回振って出た目の数の総和を得点yとする．この試行をn回実施したとき，得点yの確率分布を求めてみよう． 以下に，そのoctaveプログラムとその実行結果を示す．

% Prog_Simu1.m
clear;
close('all');
m = 10;　　　　　　　　　　　% 試験1回で得られるデータ数は10個とします
n = 100000;　　　　　　　　 % 試験回数は100000回とします
x = ceil(rand(n,m)*6);　　% 一様乱数で作ったサイの目の行列
sx = sum(x,2);　　　　　　 % 各列の和（評価量）
y = sort(sx);　　　　　　　 % 評価量を昇順に並べ替え
yi = unique(y);　　　　　　% 評価量の種類
for ix=1:length(yi)
　ygroup=find(y==yi(ix));
　ni(ix,1)=length(ygroup);% 各評価量毎の出現数
end
bar(yi,ni/n) % 確率密度分布

2.2　動的問題のシミュレーション利用

　時間の経過とともに変化する動的な問題は, ラプラス変換や線形代数を使って状態空間モデルに表してから， シミュレーションするのがスマートである．
　しかし，状態空間モデルやラプラス変換，線形代数の知識が必要なので，簡単には計算できない． でも，微分方程式に変位や速度，加速度を逐次代入して状態の変化を追ってゆくシミュレーションなら 簡単である．先に状態空間モデルを作る方法と比べたら，かなり原始的であるが， それなりの結果を示してくれる．振動解析を知らなくてもある程度のことはできるのだ．

例題：ばねの自由振動

　剛性が$\ k=10\ \rm{(N/m)}\ $で粘性が$\ c=1\ \rm{(Ns/m)}\ $のばねの一方に質量$\ m=100\ \rm{kg}\ $の錘が取り付けられ， 他方が天井に固定されているとしよう． 今，つり合って静止した錘を鉛直下方に$\ 0.1\ \rm{m}\ $引っ張って離した後の錘の位置，速度，加速度 をシミュレーションで示す．
この場合，初期の加速度は$\ \ddot{x}(0)\ $は未知で，速度は$\ \dot{x}(0)=0\ \rm{m/s}\ $，位置は$\ x(0)=-0.1\ \rm{m}\ $である． また，錘の重心に作用する慣性力$\ m\ddot{x}\ $とばねの粘性力$\ c\dot{x}\ $，剛性力$\ kx\ $の和は0である． よって，初期の加速度は$\ \ddot{x}(0)=-kx/m\ $となる． 次に，微少時間$\ dt\ $だけ経過したときの，速度は$\ \dot{x}(dt)=\dot{x}(0)+\ddot{x}(0)dt\ $， 位置は$\ x(dt)=x(0)+\dot{x}(0)dt\ $，加速度は$\ \ddot{x}(dt)=-(c\dot{x}(dt)+kx(dt))/m\ $として求められる． これを繰り返してゆけば，錘の位置をシミュレーションできる．

% Prog_Simu_2.m
clear;
close;
n=1000;　　　　　　　　　　% 試行回数
k=10;　　　　　　　　　　　% ばねの剛性 10 N/m
c=1;　　　　　　　　　　　　% ばねの粘性 1 Ns/m
m=100;　　　　　　　　　　　% 質量 100 kg
v0=0;　　　　　　　　　　　% 初期速度 0 m/s
p0=-0.1;　　　　　　　　　% 初期位置 -0.1 m
dt=0.5;　　　　　　　　　　% サンプリング周期 0.05 sec
T=[];P=[];V=[];A=[];
for ix=1:n
a0=-(c*v0+k*p0)/m;　　　 % 現時点の加速度
v1=v0+a0*dt;　　　　　　　% 次の時点の速度
p1=p0+v1*dt;　　　　　　　% 次の時点の位置
t=dt*ix;　　　　　　　　　　% 次の時点の時刻
v0=v1;　　　　　　　　　　　% 次の速度
p0=p1;　　　　　　　　　　　% 次の位置
A=[A;a0];　　　　　　　　　% 加速度ベクトル
V=[V;v0];　　　　　　　　　% 速度ベクトル
P=[P;p0];　　　　　　　　　% 位置ベクトル
T=[T;t];　　　　　　　　　% 時刻ベクトル
end
subplot(311),plot(T,A),grid,title('Accerreration')
subplot(312),plot(T,V),grid,title('Verocity')
subplot(313),plot(T,P),grid,title('Position')

2.3　その他の問題のシミュレーション利用

　前節までに，確率的な問題と動的問題の数値シミュレーション例を示したが， それ以外にも沢山のシミュレーション事例を考えることができる． ここでは，円の面積が$\pi r^2$であることを利用し，数値シミュレーションで円周率$\pi$を求めている．

例題：円周率を求める

モンテカルロシミュレーション（乱数を使ったシミュレーション）の例題として紹介されていた記憶がある． 0～1の範囲の一様乱数を2つずつ発生させて，$xy$座標上の点$(x_{i},y_{i})$を発生させ， $x_{i}^{2}+y_{i}^{2} \leq 1/$の点の割合から円周率を求めるというものだった．

% Prog_Simu_3.m
clear;
n=10000;　　　　　　　　　　% 試行回数
x=rand(n,1);　　　　　　　　% 一様乱数[0,1]からなるx座標のベクトル
y=rand(n,1);　　　　　　　　% 一様乱数[0,1]からなるy座標のベクトル
r=sqrt(x.^2+y.^2);　　　　　% 原点から各点までの距離
[N,m]=find(r<1);　　　　　　% 距離が1未満の点を数える
pi_hat=4*length(N)/n;　　　% 全体に占める距離が1未満の点の割合

2.4　実験準備としてのシミュレーション利用

　実験条件の設定をミスしたために，無駄になった実験は相変わらず多い． この手の失敗は，実験を模擬したシミュレーションでかなり防止できる． 知りたいことを未知数にした実験モデル式を作り，結果に影響する実験条件の設定を色々変えて シミュレーションを行う．
　次に，得られた模擬実験データを解析して未知数の値を導き出せることを確認する． この作業だけで，実験もデータ解析プログラムも正しいことを確認できる． つまり，実験を行う前に，シミュレーションとデータ処理のプログラムを完成させておくのである．
　この作業を忘れずに行えば， いざ実験を完了した後に「実験計画が不備で因子Aの影響がわかりませんでした」なんていう お粗末な報告は無くなるはずである．

問題2.1

　半径が1mの球の体積をシミュレーションで求めよ． ただし，円周率$\pi$を用いないこと． 円周率をモンテカルロ法で求めるシミュレーションを拡張すればよい．
【狙い】シミュレーションの練習問題

問題2.2

　寿命分布が標準正規分布に従う母集団から$k$個のデータを取り出し， 標本平均$\hat{\mu}$と標本分散$\hat{\sigma^{2}}$から母集団の寿命分布を推定する． $k=3, 10, 30$の場合について予想される寿命分布を図示しなさい． また，真の寿命分布も図示して，予想した分布との違いや共通点を述べなさい．
【狙い】 たった3個のデータから予想した確率分布で，母集団の良否を判定するにはどうしたら良いか？ このような現実的な問題を，簡単なシミュレーションで検討できることに感動して欲しい．

問題1.3

　地球の中心を通る真直ぐなトンネルを掘り，地表から鉄球を落下させたとする． このとき，トンネルを通過して地球の裏側に到着する間の鉄球の位置と速度と加速度の変化を図示し， 地球の裏側に到着するまでの時間を求めよ． ただし，空気抵抗は無視できるとし，地球に関する情報は以下で与えられているとする．

• 地球の半径：6378km
• 地球の密度：5.52g/$\rm{cm}^2$
• 地表での重力加速度：9.8m/sec$^2$

【狙い】webで「地球　トンネル」のキーワード検索すれば，この問題の厳密解法が見つかるが， とても面倒な計算である． しかし，シミュレーションを使えば，正解に近い答えを容易に得られることを体感する．