物理学の具体的な計算にMaximaを使って、物理学もMaximaも同時に学んでしまいましょう。
今回はMaximaを使って微分方程式を解いてみたいと思います。
一般に、未知関数とその導関数を含んだ方程式を微分方程式といいます。
一変数関数の場合を常微分方程式、多変数関数の場合を偏微分方程式と呼びます。
今回は力学で現れる具体的な常微分方程式をいくつか解いてみます。
Maximaのごく基本的な使い方については以下の記事を参照してください:
不定積分
まず微分方程式を解くにあたって関数の積分を復習しておきましょう。
integrate(exp(-a*x), x);
などとします(ただし積分定数は省かれてしまいます)。
少しの間、この不定積分を出力する関数を I と定義してこれを用います:
I(f):=integrate(f, x)$
例えば
assume(notequal(n,-1))$ I(x^n);
が求まります。また
makelist( I(cos(x)^n), n, 1, 5);
として の積分が求められます。あるいは
I(1/x);
として
がわかります。
さてウォーミングアップは軽く済ませて、早速、微分方程式を解いていきましょう。
微分方程式
1階線型常微分方程式
簡単な例として次の1階線型常微分方程式
を考えましょう。ここで は の関数、 は定数とします。
この微分方程式は変数分離法によって解くことができます:
より
両辺の指数関数をとって
と解くことができます。
integrate を使って微分方程式を解く
さて、これをMaximaで解かせるには、不定積分で使った integrate という組み込み関数を使って
integrate( 1/y(x)*diff(y(x),x) = -a, x );
として
がわかります(不定積分の場合と違って積分定数も含めて結果を出力します)。
ちなみに
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 でしょうか?)
裳華房の『解析学概論』(矢野・石原)は、常微分方程式(だけでなく、ベクトル解析、複素解析、フーリエ解析)に関してよくまとめられた解説と例題、さまざまな練習問題が載っており、効率よく応用数学が学べる理工系学部生におすすめのテキストです:
テキストに載っている問題を、手計算で解けるようになることはもとより、Maximaを利用して解いてみるのも、Maximaを有効に活用するための練習として有用です。
自由落下
で記述されます。
ここで は鉛直方向の座標(上を正の方向にとりました)、
は質点の質量、 は重力加速度です。
この常微分方程式を解くには integrate を2回使って
kill(all)$ integrate( diff(x(t), t, 2) = -g, t); integrate( %, t);
とすることで
がわかります。ここで は初速度、 は初期位置に相当する積分定数です。
さて、解を求める別な方法として desolve を使って
kill(all)$ desolve( diff(x(t), t, 2) = -g, x(t) );
によっても解を求めることができます。さらに、2階微分方程式を解く
斜方投射
水平右向きを 軸の正方向、鉛直上向きを 軸の正方向にとります。運動方程式は
です。いま
とします。この状況設定を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$
とします。 成分、 成分の運動方程式を積分すると(ここでは ode2 を使っています)
GSolx: ode2(EoMx, x, t); GSoly: ode2(EoMy, y, t);
により一般解が
と求まります。
あとは初期条件から4つの積分定数を決定すれば良いのですが、
その処理も ic2 という組み込み関数によってMaximaに実装されています:
ic2 を使って一般解に初期条件を課す
いま初期条件として、時刻 に原点 から初速度
で質点が飛び出したとします。この初期条件を課した解が 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));
によって
と求められます。
運動の軌跡
さて、解から時刻 を消去して運動の軌跡を調べてみます:
Traj:subst(rhs(solve(Solx, t)[1]), t, Soly);
より運動の軌跡は についての2次関数
となり、放物線を描くことがわかります。
ここで、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));
より
となって、放物線の頂点が
で与えられることがわかります。また、質点が再び地面に着く位置は
solve(rhs(Traj)=0, x);
により
であることがわかります。よって飛距離を最大とする初速度の投射角度は となります。
この角度 で初速 を与えたとき、質点の最高到達点は
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);
により
と算出されます。ただし重力加速度 を としました。
以上の運動を 秒ごとにプロットしたのが以下の図です:
この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次元調和振動子
水平右向きを 軸の正方向にとります。運動方程式は
です。この状況設定を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)));
として、一般解が
とわかります。ここで は積分定数、 は角振動数で で与えられます。
atvalue を使って一般解に初期条件を課す
いま初期条件として、時刻 に位置 から初速度 で運動を始めたとします。
この初期条件を課した解を desolve を使って求めるには、組み込み atvalue を使います:
atvalue(x(t), t=0, x0)$ atvalue('diff(x(t), t), t=0, 0)$ Solx: expand(desolve(EoMx, x(t)));
と入力することにより、解が
と与えられます。単振動の運動の様子は次のようなものです:
一方は
をプロットした点、他方の点は
です。単振動が回転運動の射影であることがわかります。
Lorentz力による運動
電場 、磁場 の中を運動する荷電粒子(電荷: )は
という運動方程式にしたがいます。ここで は位置ベクトルで
です。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]])$
と入力して、外積を定義しておきます。外場として時間変動しない電磁場
を考えたとき、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 を使っています。
いま時刻 での初期条件を
とすると解が
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 を用いました。
以上の初期条件のもとで、位置ベクトル、速度ベクトルが
と求まります。詳細はランダウ・リフシッツ『場の古典論』セクション22が参考になります。