前々回、前回と質点系の振動の話を書きました。
詳細は以前の記事を参照してください。
前回の「多質点振動の弾性振動」の場合は解こうとする対象がマトリクスなので、fortranのプログラミングに慣れた人でも作成までに３日程度かかることを紹介しました。
今回はこの問題をMathematicaで解こうとするとどういうアプローチになるのかを紹介します。
上の図の振動モデルは最終的に下記の式として整理することができます。
この式を導く流れは前回の記事を参照してください。
さて、上記の式は２階の常微分方程式を解くという問題に帰着します。
しかもその対象はマトリクス。
そう、そんなときこそMathematica！
Mathematicaならマトリクスを対象とした２階の常微分方程式を解くことなんて一発です。
実際、ヒューリンクスブログの執筆者で、主にMathematicaについて書いている守谷くんに「振動方程式って２階の常微分方程式でマトリクスはこうなるんだよ」って伝えたら３時間ほどでプログラミングしてくれました。
そのサンプルが下記です。
先ず、入力する波形を定義します。
g[t_] := d Exp[-h (2 Pi)/(T 5) t^2] t^2 (Sin[(2 Pi)/T t] + Sin[(2 Pi)/(0.7 T) t]) /. {d -> 10., h -> 0.1, T -> 0.2};
Plot[g[t], {t, 0, 10}, PlotRange -> All]
次に各マトリクスを定義。
n = 5 - 1;
mM = SparseArray[{{i_, i_} -> 1.5}, {n, n}];
mM // Normal // MatrixForm
mC = SparseArray[{{i_, i_} :> (2 i - 1) 500., {i_, j_} /;
j == i + 1 :> -500. i, {i_, j_} /; j == i - 1 :> -500. j}, {n, n}];
mC // Normal // MatrixForm
mK = SparseArray[{{i_, i_} :> (2 i - 1) 500., {i_, j_} /;
j == i + 1 :> -500. i, {i_, j_} /; j == i - 1 :> -500. j}, {n, n}];
mK // Normal // MatrixForm;
そして、初期条件を与えてNDSolve一発！
v[x_] := Table[x, {n}];
sol = NDSolve[{mM.y''[t] + mC.y'[t] + mK.y[t] == -mM.v[1.] g''[t],
y[0] == v[0], y'[0] == v[0]}, y, {t, 0, 10}, SolveDelayed -> True]
これだけではどんな解析なのか不明ですね。
Manipulateを使って実際にモデルが揺れる様子も分かるようにしてくれました。
fy[i_, t_?NumberQ] := Piecewise[{{g[t], i == n}}, (y /. sol[[1]])[t][[i]] + g[t]];
cvs = fy[#, t] & /@ Range[n];
Plot[cvs, {t, 0.5, 2.5}, PlotRange -> All]
Manipulate[
Grid[
{{ListPlot[{tmp = {Append[((y /. sol[[1]])[t][[#]] & /@Range[n]) + g[t], g[t]], Range[n, 0, -1]} // Transpose, tmp},
Joined -> {True, False}, PlotStyle -> PointSize[0.08], PlotRange -> {{-20, 20}, {0, n + .5}}, AspectRatio -> 2,
Prolog -> Text[t, {-16, n - 1}], ImageSize -> 200],
Plot[g[t], {t, 0, 4}, PlotRange -> All,
Prolog -> {Hue[0], PointSize[0.02], Point[{t, g[t]}]},
ImageSize -> 200]}}], {t, 0, 4}]
どうでしょう。NDSolveのところが拍子抜けするくらい簡単です。