読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

(20) 運動方程式の級数展開による近似解法

調和振動子運動方程式は以下の通りです。

ここで、pは運動量、qは位置、mは質量、kはバネ係数です。
これは解析的に解くことができます。
初期条件を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])$