pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Classical Music

Maximaで微分方程式を解く(力学の問題を題材にして)

物理学の具体的な計算にMaximaを使って、物理学もMaximaも同時に学んでしまいましょう。

今回はMaximaを使って微分方程式を解いてみたいと思います。

一般に、未知関数とその導関数を含んだ方程式を微分方程式といいます。

一変数関数の場合を常微分方程式、多変数関数の場合を偏微分方程式と呼びます。

今回は力学で現れる具体的な常微分方程式をいくつか解いてみます。

Maximaのごく基本的な使い方については以下の記事を参照してください:

pianofisica.hatenablog.com


不定積分

まず微分方程式を解くにあたって関数の積分を復習しておきましょう。

関数  y=f(x)不定積分

 \displaystyle{F(x)=\int dx \, f(x) +C}

 C積分定数)はMaximaでは

integrate(exp(-a*x), x);

などとします(ただし積分定数は省かれてしまいます)。

少しの間、この不定積分を出力する関数を I と定義してこれを用います:

I(f):=integrate(f, x)$

例えば

assume(notequal(n,-1))$
I(x^n);

として  x^n不定積分

 \displaystyle{\int dx \, x^n =\frac{x^{n+1}}{n+1}+C \qquad (n\neq-1)}

が求まります。また

makelist( I(cos(x)^n), n, 1, 5);

として  f(x)=\cos^n x \ (n=1,2,3,\dots)積分が求められます。あるいは

I(1/x);

として

 \displaystyle{\int dx \, \frac1x =\log x+C }

がわかります。

さてウォーミングアップは軽く済ませて、早速、微分方程式を解いていきましょう。



微分方程式

1階線型常微分方程式

簡単な例として次の1階線型常微分方程式

 \displaystyle{\frac{dy}{dx}=-ay}

を考えましょう。ここで  y x の関数、 a は定数とします。

この微分方程式は変数分離法によって解くことができます:

 \displaystyle{\int dx\, \frac{dy}{dx}\,\frac{1}{y}=-a\int dx}

より

 \displaystyle{\log y=-ax + c \quad (c:\text{積分定数})}

両辺の指数関数をとって

 \displaystyle{y=C\exp\left(-ax\right)  \quad (C:\text{積分定数})}

と解くことができます。

integrate を使って微分方程式を解く

さて、これをMaximaで解かせるには、不定積分で使った integrate という組み込み関数を使って

integrate( 1/y(x)*diff(y(x),x) = -a, x );

として

 \displaystyle{\log y=-ax + c \quad (c:\text{積分定数})}

がわかります(不定積分の場合と違って積分定数も含めて結果を出力します)。

ちなみに

integrate( diff(y(x),x) = -a*y(x), x );

とした場合には答えを出力せず、うまくいかないようです。

desolve を使って微分方程式を解く

さて、微分方程式を解く実装がMaximaには用意されています:すると前者のような書き方

desolve( diff(y(x),x)/y(x) = -a, y(x) );

ではうまく計算できず、後者のような書き方

desolve( diff(y(x),x) = -a*y(x), y(x) );

はうまく処理できます(このあたりは trial and error でしょうか?)

自由落下

1次元自由落下は2階常微分方程式運動方程式

 \displaystyle{m\frac{d^2x}{dt^2}=-mg}

で記述されます。

ここで  x \ [{\rm m}] は鉛直方向の座標(上を正の方向にとりました)、

 m \ [{\rm kg}] は質点の質量、 g \ [{\rm m/s^2}] は重力加速度です。

この常微分方程式を解くには integrate を2回使って

kill(all)$
integrate( diff(x(t), t, 2) = -g, t);
integrate( %, t);

とすることで

 \displaystyle{x(t)=-\frac12 g t^2+C_1 t+C_2}

がわかります。ここで  C_1 は初速度、 C_2 は初期位置に相当する積分定数です。

さて、解を求める別な方法として desolve を使って

kill(all)$
desolve( diff(x(t), t, 2) = -g, x(t) );

によっても解を求めることができます。さらに、2階微分方程式を解く

ode2 を使って微分方程式を解く
ode2( diff(x(t), t, 2) = -g, x(t), t );

によっても解くことができます。

以上見たように、Maxima微分方程式を解くにあたっては、主に

integrate、desolve、ode2 を使うという3つの方法が考えられます。

求めたい関数や変数の書式、うまく処理できる場合が異なっていることに注意が必要です。

斜方投射

水平右向きを  x 軸の正方向、鉛直上向きを  y 軸の正方向にとります。運動方程式

 \displaystyle{m\frac{d^2x}{dt^2}=F_x \qquad m\frac{d^2y}{dt^2}=F_y}

です。いま

 \displaystyle{F_x=0 \qquad F_y=-mg}

とします。この状況設定をMaximaに入力するには

kill(all)$ assume(m>0,g>0)$ declare([m,g],constant)$ 
depends([x,y],t)$
Fx:0$    Fy:-m*g$    
EoMx: m*diff(x,t,2)=Fx$
EoMy: m*diff(y,t,2)=Fy$

とします。 x 成分、 y 成分の運動方程式積分すると(ここでは ode2 を使っています)

GSolx: ode2(EoMx, x, t);
GSoly: ode2(EoMy, y, t);

により一般解が

 \displaystyle{x(t)=B_1 t+B_2\qquad y(t)=-\frac12 g t^2+C_1 t+C_2}

と求まります。

あとは初期条件から4つの積分定数を決定すれば良いのですが、

その処理も ic2 という組み込み関数によってMaximaに実装されています:

ic2 を使って一般解に初期条件を課す

いま初期条件として、時刻  t=0 に原点  (x=0, \ y=0) から初速度

 \displaystyle{v_x(t=0)=v\cos\theta, \qquad v_y(t=0)=v\sin\theta}

で質点が飛び出したとします。この初期条件を課した解が ic2 を使って

Solx: ic2(GSolx, t=0, x=0, 'diff(x, t)=v*cos(theta));
Soly: ic2(GSoly, t=0, y=0, 'diff(y, t)=v*sin(theta));

によって

 \displaystyle{x(t)=vt\cos\theta, \qquad y(t)=-\frac12gt^2+vt\sin\theta}

と求められます。

運動の軌跡

さて、解から時刻  t を消去して運動の軌跡を調べてみます:

Traj:subst(rhs(solve(Solx, t)[1]), t, Soly);

より運動の軌跡は  x についての2次関数

 \displaystyle{y={{x\sin \theta}\over{\cos \theta}}-{{gx^2}\over{2v^2\cos^2\theta}}}

となり、放物線を描くことがわかります。

ここで、2次式を平方完成する関数をMaximaで作って

a(f):=ratcoeff(f, x, 2)$ b(f):=ratcoeff(f, x, 1)$ c(f):=ratcoeff(f, x, 0)$
complsq(f):=a(f)*(x+b(f)/(2*a(f)))^2+c(f)-b(f)^2/(4*a(f))$

運動の軌跡の式を平方完成すると

complsq(rhs(Traj));

より

 \displaystyle{y={{v^2\sin^2\theta}\over{2g}}-{{g}\over{2v^2\cos^2\theta}}\left(x-{{v^2\sin\theta\cos\theta}\over{g}}\right)^2}

となって、放物線の頂点が

 \displaystyle{(x_{\rm max},y_{\rm max})=\left({{v^2\sin\theta\cos\theta}\over{g}},{{v^2\sin^2\theta}\over{2g}}\right)}

で与えられることがわかります。また、質点が再び地面に着く位置は

solve(rhs(Traj)=0, x);

により

 \displaystyle{x_{\rm farthest}={{2v^2\sin\theta\cos\theta}\over{g}}={{v^2\sin2\theta}\over{g}}}

であることがわかります。よって飛距離を最大とする初速度の投射角度は  \displaystyle{\theta_{\rm max}=\frac{\pi}{4}} となります。

この角度  \displaystyle{\theta_{\rm max}} で初速  80\ [{\rm km/h}]=22.2\ [{\rm m/s}] を与えたとき、質点の最高到達点は

v: 80/3.6$
theta: %pi/4$
g: 9.8$
float(v^2/g*sin(theta)*cos(theta));
float(v^2/(2*g)*sin(theta)^2);

により

 \displaystyle{(x_{\rm max},y_{\rm max})=\left(25.2\ [ {\rm m} ],12.6\ [ {\rm m} ]\right)}

と算出されます。ただし重力加速度  g 9.8\ [ {\rm m/s^2} ] としました。

以上の運動を  0.05 秒ごとにプロットしたのが以下の図です:
f:id:pianofisica:20190119121825g:plain

このGIFアニメの作成方法は以下の通りです:まずMaxima

for i:0 step 1 thru 70 do(
     b:concat("/保存先ディレクトリのパス/", i, ".eps"),
     t : float(0.05*i),
     position_x: [v*t*cos(theta)],
     position_y: [-1/2*g*t^2+v*t*sin(theta)],
	plot2d ([discrete,position_x,position_y],[style,points],[x, -2, 53], [y,-2,15],
	            [plot_format, gnuplot], [gnuplot_term, ps], [gnuplot_out_file, b]));

として複数のEPSファイルを作成します。

これらの1つ1つの画像が各時刻におけるスナップショットになっています。

コンピュータにImageMagickをインストールしたうえで、

カレントディレクトリをEPSファイルの保存してあるディレクトリに移動し、ターミナルから

$ convert -layers optimize -loop 0 -delay 10 *.eps anim.gif

としてGIFアニメを作成しました。

GIFアニメの作成過程で入力した投射角度や初速を変更すれば、

それに応じて運動の軌跡が変わっていく様子を観察することができます。



1次元調和振動子

水平右向きを  x 軸の正方向にとります。運動方程式

 \displaystyle{m\frac{d^2x}{dt^2}=F_x \qquad F_x=-kx}

です。この状況設定をMaximaに入力して解を求めるには

kill(all)$ assume(m>0,k>0)$ declare([m,k],constant)$ 
Fx:-k*x(t)$ 
EoMx: m*diff(x(t),t,2)=Fx$
GSolx: expand(desolve(EoMx, x(t)));

として、一般解が

 \displaystyle{x(t)=A\sin(\omega t)+B\cos(\omega t)}

とわかります。ここで  A, \ B積分定数 \displaystyle{\omega} は角振動数で  \displaystyle{\omega=\sqrt{\frac{k}{m}}} で与えられます。

atvalue を使って一般解に初期条件を課す

いま初期条件として、時刻  t=0 に位置  x=x_0 から初速度  v_x=0 で運動を始めたとします。

この初期条件を課した解を desolve を使って求めるには、組み込み atvalue を使います:

atvalue(x(t), t=0, x0)$
atvalue('diff(x(t), t), t=0, 0)$
Solx: expand(desolve(EoMx, x(t)));

と入力することにより、解が

 \displaystyle{x(t)=x_0\cos(\omega t)}

と与えられます。単振動の運動の様子は次のようなものです:
f:id:pianofisica:20190119133035g:plain

一方は

 x=\cos t \quad y=0

をプロットした点、他方の点は

 x=\cos t \quad y=\sin t

です。単振動が回転運動の射影であることがわかります。

Lorentz力による運動

電場  \vec{E}、磁場  \vec{B} の中を運動する荷電粒子(電荷 q )は

 \displaystyle{\frac{d^2\vec{r}}{dt^2}=q\left(\vec{E}+\frac{d\vec{r}}{dt}\times \vec{B}\right)}

という運動方程式にしたがいます。ここで  \vec{r} は位置ベクトルで

 \displaystyle{\vec{r}=\left(\begin{array}{c}  x \\ y \\ z \end{array} \right)}

です。Lorentz力はベクトルの外積を含むので、あらかじめMaxima

kill(all)$ load(eigen)$ 
exprod(a,b):= columnvector([
a[2][1]*b[3][1]-a[3][1]*b[2][1],
a[3][1]*b[1][1]-a[1][1]*b[3][1],
a[1][1]*b[2][1]-a[2][1]*b[1][1]])$

と入力して、外積を定義しておきます。外場として時間変動しない電磁場

 \displaystyle{\vec{E}=\left(\begin{array}{c}  0 \\ E_y \\ E_z \end{array} \right) \qquad \vec{B}=\left(\begin{array}{c}  0 \\ 0 \\ B_z \end{array} \right)}

を考えたとき、Lorentz力は

r:columnvector([x(t),y(t),z(t)])$
declare([Ey,Ez,Bz,q,m],constant)$ assume(Ey>0,Ez>0,Bz>0,q>0,m>0)$
E:columnvector([0,Ey,Ez])$
B:columnvector([0,0,Bz])$
F:q*(E+exprod(diff(r,t),B));

で求められます。運動方程式の一般解が

EoM[i]:= m*diff(r[i][1],t,2)=F[i][1]$
GSol:expand(desolve( [EoM[1], EoM[2], EoM[3]], [x(t), y(t), z(t)]))$
i:1$ while i<=3 do(print(GSol[i]), i:i+1);
i:1$ while i<=3 do(print(diff(GSol[i],t)), i:i+1);

によって出力されます。ここで微分方程式が連立式であることから desolve を使っています。

いま時刻  t=0 での初期条件を

 \displaystyle{\vec{r}(t=0)=\left(\begin{array}{c}  0 \\ 0 \\ 0 \end{array} \right) \qquad \frac{d\vec{r}}{dt}(t=0)=\left(\begin{array}{c}  0 \\ 0 \\ 0 \end{array} \right)}

とすると解が

atvalue(x(t), t=0, 0)$ atvalue('diff(x(t), t), t=0, 0)$
atvalue(y(t), t=0, 0)$ atvalue('diff(y(t), t), t=0, 0)$
atvalue(z(t), t=0, 0)$ atvalue('diff(z(t), t), t=0, 0)$
Sol:expand(desolve( [EoM[1], EoM[2], EoM[3]], [x(t), y(t), z(t)]))$
i:1$ while i<=3 do(print(Sol[i]), i:i+1);
i:1$ while i<=3 do(print(diff(Sol[i],t)), i:i+1);

で求められます。

ここで desolve を使ったことに応じて初期条件の入力に atvalue を用いました。

以上の初期条件のもとで、位置ベクトル、速度ベクトルが

 \displaystyle{\vec{r}(t)=\left( \frac{E_y}{B_z}\left(t-\frac{m}{qB_z}\sin\left(\frac{qB_z}{m}t\right)\right), \ \frac{mE_y}{qB_z^2}\left(1-\cos\left(\frac{qB_z}{m}t\right)\right), \ \frac{qE_z}{2m}t^2  \right) \\ \vec{v}(t)=\left( \frac{E_y}{B_z}\left(1-\cos\left(\frac{qB_z}{m}t\right)\right), \ \frac{E_y}{B_z}\sin\left(\frac{qB_z}{m}t\right), \ \frac{qE_z}{m}t  \right)}

と求まります。詳細はランダウ・リフシッツ『場の古典論』セクション22が参考になります。

キーワード微分方程式運動方程式Maxima、自由落下、Lorentz力