前々回前回と質点系の振動の話を書きました。
詳細は以前の記事を参照してください。

前回の「多質点振動の弾性振動」の場合は解こうとする対象がマトリクスなので、fortranのプログラミングに慣れた人でも作成までに3日程度かかることを紹介しました。

今回はこの問題をMathematicaで解こうとするとどういうアプローチになるのかを紹介します。
1.png


上の図の振動モデルは最終的に下記の式として整理することができます。
1_siki.png
この式を導く流れは前回の記事を参照してください。

さて、上記の式は2階の常微分方程式を解くという問題に帰着します。
しかもその対象はマトリクス。

そう、そんなときこそMathematica!
Mathematicaならマトリクスを対象とした2階の常微分方程式を解くことなんて一発です。

実際、ヒューリンクスブログの執筆者で、主にMathematicaについて書いている守谷くんに「振動方程式って2階の常微分方程式でマトリクスはこうなるんだよ」って伝えたら3時間ほどでプログラミングしてくれました。

そのサンプルが下記です。

先ず、入力する波形を定義します。
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のところが拍子抜けするくらい簡単です。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/435
内容に対しての関連性がみられないものは削除する場合があります

コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報