|
これだけではおもしろくないので、惑星の運動をアニメーション してみましょう。実際やってみると、今回のテーマである数値計算 の必要性を再び感じることになります。 |
| c) 多項式方程式の数値解 -- NSolve |
|
アニメーションを行うには、軌跡の式（ｒとθの関係式）ではなく、 座標を時間の関数として表現しておかなければなりません。 これは、ケプラーの第三法則と幾何学的な考察より行うことができます。 詳細はここを読んでください。 以下のような方程式が得られます。
φは位置パラメーターで時間に依存しており時刻と第一式の関係を 満たしています。実際の座標(x,y)はさらに第二、三式によりφの 関数で表されています。 φが時間の関数なので、形式的にはこれでxとyが時間の関数として求 まったことになります。しかしアニメーションを行うには、実際にｔ を陽に含む形に書き換える必要があり、それには第一式 （ケプラー方程式）を解かなければなりません。 これをSolveを使って解いてみましょう。
Mathematica に「解析解がないようです」と怒られています。 Solve は純粋な代数方程式ばかりでなく、Sin などの関数を含む方程式も解いてその解析解を求めることができます が、当然のことながら一般解が存在する場合に限られます。 ここで解きたい方程式は超越方程式なので一般解が存在しません。 仕方がないので今回の主題である数値計算を行うことにしましょう。 DSolveに対して NDSolve があったように、Solve に対しては NSolve があるのですが、これは 多項式の方程式に対してしか使えない という制限があります。そこで、Sinφを Taylor 展開してから、 NSolve を使ってみます。 また、数値計算なので時刻を具体的に 与えなければ解いてくれません。 まず、適当な時刻での位置をを求めてみます。 時刻が1周期の32分の1の場合、各パラメータを
とします。 3行目の1式目がケプラーの第三法則であり、3式目は時刻が1周期の 32分の1を表します。 (1周期をを32分割してアニメーションを行うためのものです。) また、下の2式が幾何学的考察 より求まったｘとｙの式です。
次にケプラー方程式
です。 4次の項以上を無視するには Normal を使います。
これを使ってケプラー方程式を解くと、
と計算できます。3次方程式なので解が3つ出てきています。 実数解である３つ目のφの値より位置 (x,y) が求まります。 もちろん、この値は Taylor 展開による近似誤差を含んでいます。 NSolveは、Solveで解けない、もしくは時間 がかかるような場合、もし欲しいものが数値解で充分ならば これを高速に計算することができますが、今の場合これでは 見通しが悪いので、次にFindRoot を用いた別の方法 を試してみることにします。 |
| d)Newton法を用いた方程式の数値解 -- FindRoot |
|
FindRoot はNewton法というアルゴリズムを利用して、 方程式の数値解を求める関数です。 これを使えば、解きたい方程式が超越方程式であっても、その近似解を求めることができます。 Newton法は簡単なアルゴリズムで数値計算の教科書には必ず載っています。 その原理や計算するための表式、他のアルゴリズムとの比較などを調べて、 レポートを提出してください。 ここでは、Newton法に関して以下の点を抑えておけば結構です。
Newton法は上で述べたように、アルゴリズム中で導関数を利用するため、 導関数が形式的に求まらない場合には使うことができません。 実は Mathematica の FindRootという関数はもう少し賢くて、 導関数が求まらない場合には、割線（secant）法を利用することができ ます。 その場合、初期値が2つ必要になります。 この方法についてもレポートで調べてみてください。レポートではまず、ニュートン法、割線法がそれぞれどのような方法なのかを解説してください。また、各々の場合で\( \sin(x) = x/2 \)を例として近似解が真の解に近づいていく様子も具体的に説明して下さい。レポートの採点方針は、何も知らない人が読んで理解できるかどうかに重点を置くことにします。これを読んだら初めての人でも、ニュートン法、割線法のどちらも理解できるようになる、というレポートを目指してください。 それではMathematicaに戻って、実際に計算をやってみましょう。この場合は導関数が計算できますから Newton法、 すなわち初期値は1つで充分です。
今度は数式には近似を用いていないので数値的な誤差しか含みません。 しかしFindRoot は、Solve や NSolveと違い、 すべての解を計算してはくれません。 初期値（上の例ではφ=0）から出発して、最初に到達したものを解 として表示します。念のため、別の初期値から出発してみます。
同じ結果が得られました。今のパラメータでは解が一つしかない ことはグラフの交点を見てみるとわかると思います。 こうして求まったφの値を上記の x, y に代入して、その点を軌跡のグラフにプロットしてみます。
この時刻の惑星の位置を表すために、塗りつぶされた円盤を Disk[{x, y}, r]をもちいて位置(x,y) に描きました。 ここで、Disk[{x, y}, r]は 2次元グラフィックス要素であり、半径 r で中心が (x,y) にある塗りつぶされた円盤を表します。 2次元グラフィックス要素ですので、 関数Graphics の中で用います。これを Circleに代えれば、 塗りつぶされていない円盤になります。似たようなものに Point、Line、Rectangle、Polygon、Raster などがあります。 |
| e)アニメーション |
|
では、いよいよアニメーション化してみましょう。 先ほどの場合は時刻が１周期の32分の1の時の位置パラメータφ を求めました。 これを時刻が1周期の32分の2、32分の3、...のそれぞれの位置 パラメーターφを求め(x,y)を1周期分得れば良いのですが、上 の作業を手作業で32回繰り返すのは大変です。 幸い、Mathematica には繰り返し実行を行わせる関数Do が用意されていますので、これを利用しましょう。 Do[expr[n],{n,32}]とすれば、 n を 1 から 32 まで 1 ずつ増やしながら、nに依存した数式 expr[n] を実行してくれます。 従って次に、n を使って時刻を動かしながら FindRoot を解くことを考えます。 以下のように変更します。 ここでは、 時刻が、n t = n T/32 で1周期の32分の n になるようにし、 また、FindRoot の解の探索の初期値をε=0 の時の解から 出発するようにしています。
とすれば、 32個グラフがプロットできるはずです。Mathematica の古いバージョン では、描画の速度が遅くアニメーションのようにみえました。 古い時代（といってもほんの数年前ですが）の話はやめて、 Animate を使って みましょう：
面積速度一定の法則が実感できるのではないでしょうか？ |
のSin[φ]を NDSolve できるよう Taylor 展開します。
Taylor 展開は関数 Seriesを使って行うことができます。
φ=0のまわりで展開した 3次の項までがほしければ、
