これは解析的に解くことができます。
初期条件をp(0)=0, q(0)=1として、Maximaに解かせてみましょう。
assume(m > 0, k > 0)$ atvalue(q(t), t = 0, 1)$ atvalue(p(t), t = 0, 0)$ desolve([diff(q(t),t)=p(t)/m, diff(p(t),t)=-k^2*q(t)],[q(t),p(t)]);
使用した機能
角速度がであることがわかります。
運動方程式は解析的に解けるとは限りません。時間に関して級数展開することで、近似解を求める方法をやってみます。
pとqの級数係数それぞれが致すべき式を求めてみましょう。
/* 5.1.5m */ coeff_equations(eq, x) := block( [n : hipow(eq, x) - 1, i, result : []], /* Taylor級数の最高次は誤差があるので省く。 */ for i : 0 thru n do result : cons(coeff(lhs(eq), x, i) = coeff(rhs(eq), x, i), result), result)$ q : sum(a[i] * t ^ i, i, 1, inf) + 1$ p : sum(b[i] * t ^ i, i, 1, inf)$ deg : 15$ qexp : taylor(q, t, 0, deg)$ pexp : taylor(p, t, 0, deg)$ f : append(coeff_equations(''diff(qexp, t) = pexp / m, t), coeff_equations(''diff(pexp, t) = -m * omega ^ 2 * qexp, t));
使用した機能
続いて、各係数を代数的に解きましょう。
/* 5.1.6m */ aList : makelist(a[i], i, 1, deg)$ bList : makelist(b[i], i, 1, deg)$ solute : solve(f, append(aList, bList));
答えを級数に戻しましょう。
/* 5.1.7m */ qqq : qexp, solute[1]; ppp : pexp, solute[1];
使用した機能
視覚化してみましょう。
/* 5.1.9m */ m : 1.5$ w : 1.0$ plot2d([parametric, qqq, ppp, [t, 0, 2*%pi], [nticks, 200]], [x, -2.0, 2.0], [y, -1.5, 1.5])$
参考文献
- Crandall, Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ) p.81-p.86