09. 二階微分はどうよ: Newton運動方程式
Photo by Kamil Pietrzak on Unsplash
時間方向に二階微分が!
これまでは
\[ \frac{du}{dt} = F(u,t) \]
というタイプの、時間方向の微分が一階しかない時間発展問題の常微分方程式を扱ってきた.
しかしもちろんのことながら、われわれが取り組みたい問題はこういうものばかりではなく、例えば
\[ \frac{d^2 u}{dt^2} = F(u,\frac{du}{dt}, t) \]
といったような、時間方向の微分が二階の問題や、三階以上の問題も考えうる.
こういう場合にどうしたらよいか、それを学んでみよう.
いわゆる Newton 運動方程式
時間方向に二階微分が入る時間発展問題の常微分方程式といえば、高校時代に学ぶ Newton の運動方程式がまずは挙げられよう. そこで、ここではそれを扱ってみよう.
具体的には、バネでぶら下げられている重りの挙動を例として考えてみよう. バネ定数を $k$, 重りの重さを $m$, 重力定数を $g$ とすれば、重りの高さ位置 $h(t)$ $(t: $時間$)$ についての Newton の運動方程式は
\[ m \frac{d^2 h}{dt^2} = -kh -mg \]
となる. むろん、解をきちんと考えるためにはさらに初期値等の情報が必要で、例えば
\[ \left\{ \begin{array}{ccl} h(0) & = & h_{ini}, \cr \displaystyle \frac{dh}{dt}(0) & = & v_{ini} \end{array} \right. \]
などという情報を付与しないといけない.
さて、これを例に、しばらく考察を進めてみよう.
この問題は実は手計算で厳密解が求まる. 書いてしまうと, \[ \left\{ \begin{array}{ccl} h(t) & = & C \cos(\alpha t) + S \sin(\alpha t) + \beta , \cr \alpha & := & \sqrt{k/m}, \cr \beta & := & - mg/k \end{array} \right. \label{eq:strict-solution-form} \] という感じだ. ただし,定数 $C$, $S$ は初期条件で定まる量だね.
だから,上のように $h_{ini}$, $v_{ini}$ という初期条件を与える場合はさらに具体的に, \[ \left\{ \begin{array}{ccl} h(t) & = & C \cos(\alpha t) + S \sin(\alpha t) + \beta , \cr C & = & h_{ini} - \beta , \cr S & = & v_{ini}/\alpha , \cr \alpha & := & \sqrt{k/m} , \cr \beta & := & - mg/k , \cr \end{array} \right. \] と求まってしまう. まあ今回は「この厳密解を知らない体(ていと読む)で数値計算する」話になっているのでそのつもりで以下を読もう.
(計算その前に) 保存量
計算について考える前に、この系の保存量を考えてみよう. これは高校の物理で学んだ通りで、全エネルギー = 運動エネルギー + 位置エネルギー が保存されるはずだ.
そして、運動エネルギーは
\[ \frac{1}{2} m \left( \frac{dh}{dt} \right)^2 \]
と表される.
位置エネルギーについては,位置エネルギー = "力の積分 on 逆らって移動した距離"" だったので
\[ \int_0^h (kx + mg) dx = \frac{1}{2} kh^2 + mgh \]
となることは覚えているだろうか? まあそこは思い出してもらうとして、結論としてはこの計算によって全エネルギーは
\[ E(h(t)) = \frac{1}{2} m \left( \frac{dh}{dt} \right)^2 + \frac{1}{2} kh^2 + mgh \label{eq:energy-original} \]
となるわけなので、これがこの系の保存量である…はずだな.
高階微分問題への近似方法のアプローチ
さて,この問題をどう数値計算していくか,今から見ていこう. いくつかアプローチが考えられるので,先に分類を見せておこう.
(分類) | 問題そのものを変形しない | 問題そのものを変形する |
---|---|---|
離散化を工夫しない | アプローチ 1 | アプローチ 3 |
離散化を工夫する | アプローチ 2 | アプローチ 4 |
それぞれ利点,欠点があるのでそのあたりに気をつけつつ,読んでいこう.
アプローチ 1: そのまま離散化
高階微分が含まれる時間発展問題へのアプローチはいくつか考えられるが、 いちばんストレートなのは
「方程式はそのままに、離散近似も素直に」
という手法だろう.
先のバネの問題に対しては,方程式はいじらずに,そして,離散化は 近似解を $h_n \cong h(n\Delta t)$, $n = 0,1, \cdots$ として
\[ \left\{ \begin{array}{ccl} \displaystyle m \left( \frac{h_{n+1} - 2h_n + h_{n-1} }{\Delta t^2} \right) & = & -k h_n - mg, \,\, n = 2, 3, \cdots , \cr \mbox{初期条件: }\, \, h_0 & \cong & h(0), \cr \mbox{初期条件: }\, \, h_1 & \cong & h(\Delta t) \end{array} \right. \]
という近似式を提案する、などだ. プログラムしやすいように少し変形すると次の式になる.
\[ \left\{ \begin{array}{ccl} h_{n+1} & = & \displaystyle \left( 2 - \frac{k\Delta t^2}{m} \right) h_n - h_{n-1} - g\Delta t^2 , \,\, n = 2, 3, \cdots , \cr \mbox{初期条件: }\, \, h_0 & \cong & h(0), \cr \mbox{初期条件: }\, \, h_1 & \cong & h(\Delta t) \end{array} \right. \]
この場合,初期条件の1つ目の具体化はあまり問題ないが,2つ目については $h(\Delta t)$ の値がわかっていないので困ってしまう. そこで,1ステップだけ他の方法で近似値を求めるなどのアイディアが必要になる. たとえば Euler 法的な近似を使うならば,まず下記のように考えて,
\[ \left\{ \begin{array}{ccl} h_0 & = & h_{ini} \, , \cr \mbox{Euler 法的に $v$ を近似: }\,\, v_{ini} & = & \displaystyle \frac{h_1 - h_0}{\Delta t} \end{array} \right. \]
これを整理した次の式を採用することになるだろう. これなら $h_{ini}$, $v_{ini}$ を使って $h_0$, $h_1$ が作れる.
\[ \left\{ \begin{array}{ccl} h_0 & = & h_{ini} \, , \cr h_1 & = & h_0 + v_{ini}\,\Delta t \end{array} \right. \]
これはあくまで「例」であって、良い離散化かどうかは全く考慮していない…
実際にやってみよう: アプローチ1
では早速,実際に動かしてみよう.
各種定数や初期値等をまず定義して…
|
|
上の初期値 h_ini
, v_ini
の値から,
エネルギーの定義式 eq.$(\ref{eq:energy-original})$ を用いて
すでに系のエネルギー値が計算できることに注意しておこう.
この場合,エネルギー値は
\[
E = 0.5
\]
となる.確かめておこう.
さて話を戻して,次に,
$h_n$ を h_curr
,
$h_{n-1}$ を h_old
と書いて,新しい値 $h_{n+1}$ を返す関数を定義しよう.
|
|
あとは初期条件の処理も合わせて,プログラムを全体を書いてしまおう.
|
|
近似解をプロットしよう.
手順はいつものとおり,
|
|
としてから,わかりやすいように時間軸方向もデータを丁寧に用意して、
|
|
これまたいつものように次の命令で良い.
|
|
ふむ,悪くない数値解のように見える. 重りが振動している様子が出力されている様子だよね.
ここでこっそり,「知らない振りをしている厳密解」と比べてみよう.まずプログラムは次のような感じ.
|
|
数値解と厳密解を同時にプロットしてみよう.
|
|
よほどひどい数値計算でない限り,まあ重なるよね. というわけで,「誤差」をプロットしよう.
|
|
ふむ,近似解は結構あっているように見えたが,桁数で言えば 2桁ぐらいしか正しくないということのようだ.
上の例では,Julia の
broadcast
という機能を使っている.「ピリオド」がついている箇所がそうだね.
よくわからない人は,次のようにしても同じことができるのでそちらで理解しよう.
|
|
さて,この数値計算でエネルギーはどうなっただろうか. 確かめてみよう.
まず,エネルギーの計算式を作りたい…のだが, エネルギーの式 eq.$(\ref{eq:energy-original})$ を見るとわかるように,その定義には $dh/dt$ が含まれていて,$h(t)$ の近似値である $h_0, h_1, h_2, \cdots$ しか計算していない上の数値計算だと「どうやって $dh/dt$ の近似値を求めよう?」となってしまう.
まあしかたないので,「なるべく適切に近似」するように自分で決めるしかないね.
例えば,この数値計算法は(初期値の処理を除いて)時間方向に対称なので,エネルギーも時間方向に対称な近似をしたほうが適切だろうと考えて,一つの例として \[ \left( \left. \frac{dh}{dt} \right|_n \right)^2 \cong \displaystyle \frac{1}{2} \left\{ \left( \frac{h_{n+1} - h_n}{\Delta t} \right)^2 + \left( \frac{h_{n} - h_{n-1}}{\Delta t} \right)^2 \right\} \]
とし,次の式を使ってエネルギー値の近似とする方法があるね.ここではそうしておこう. \[ E_n = \displaystyle \frac{m}{4} \left\{ \left( \frac{h_{n+1} - h_n}{\Delta t} \right)^2 + \left( \frac{h_{n} - h_{n-1}}{\Delta t} \right)^2 \right\} + \frac{1}{2}kh_n^2 + mgh_n \]
これは $h_{n+1}$, $h_n$, $h_{n-1}$ の 3つの値から $t = n\Delta t$ でのエネルギー値 $E_n$ が求まるという形になっている.
そこで,たとえばこれを h_new
, h_curr
, h_old
の 3つの値をもらって E_n
相当を返す関数として作っておこう.
|
|
では,これを使ってエネルギー値を計算しながらもう一度計算しよう.
|
|
さてエネルギーはどうかな?
まず数字を見てみよう(ちなみに,h_n
よりもデータが 2つ少ないよね).
|
|
499-element Array{Float64,1}:
1.235025
1.2350745025000003
1.23024653235025
1.2207337255597261
1.206915643119367
⋮
1.0859048467179022
1.110261220668761
1.13472017191485
1.158305788301444
1.1800770037347217
あれ? 今の状況だと $E=0.5$ のはず.これはまずいな.
プロットでも確認しよう.
|
|
値もおかしいし,上下に 10% ぐらいの振動もある. この様子だと,少なくともこのエネルギーの近似値はあまり信用できないなあ.
う~ん,状況をもう少し見るために,Δt
の値を小さくして(誤差が小さくなると思われる)再計算して様子を見よう.
|
|
まずは近似解 h_n
の確認だ.
|
|
さっきのグラフと比較しても,近似解は悪くなさそうだということがわかる.
一応誤差も見ておこう.
|
|
う~ん,誤差は 2桁小さくなっているな.Δt
を2桁小さくしたら誤差も同じように2桁小さくなった,ということで,これは
この計算法は精度でいえば一次精度, つまり
\[
\mbox{error} \cong O(\Delta)
\]
ということかなあ.
計算法そのものは時間方向に対称なので二次精度
($\mbox{error} \cong O(\Delta^2)$)
なんだけど,初期値の作り方が Euler的で一次精度になってしまっているのが原因なんだろうな.
あまり良くないねえ.
さて,ではエネルギーはどうだろうか.
近似エネルギー値 E_n
値をまずはみてみよう.
|
|
50000-element Array{Float64,1}:
0.5049240149002502
0.5049240197529657
0.5049240245581272
0.5049240293157156
0.5049240340257121
⋮
0.5049210620561304
0.5049210784579437
0.5049210948237981
0.5049211111535099
0.5049211274472656
おや? 値がだいぶ $E = 0.5$ に近くなったな. 振動も小さくなったようだ. 初期値を引いた差をプロットしてみよう.
|
|
4桁ぐらい保たれるようになったな.まあなんとか及第点か.
全体としては,う~ん,そうだなあ,このアプローチは数値解そのものは近似は悪くなさそうだが,エネルギー値など「系の性質」をみるような量を計算してみると状況は良くないね.
この場合は Δt
をかなり小さくとらないとダメだねえ.
なので,このアプローチ 1 については以下のような感じだ.
- 数値計算の為の式は作りやすい.
- プログラムも簡単だ.
- 初期値の処理に困るかも…, エネルギー値の近似に困るかも… 等々,理論との整合性を保つのに困るシーンがありうる.
- 数値計算の結果は悪くないが,精度は低い.エネルギー値などはあまり良くない.エネルギー値を調べるなら
Δt
を相当小さくしないとダメ.
というわけで,このアプローチ 1 は「ひどく悪いわけではないが,オススメしにくいアプローチ」というところかな.
アプローチ 2: 式はそのままに,離散化に工夫を加える
アプローチ 1 の改善版として考えられるアプローチは、「方程式はそのままに、離散化は工夫して」という手法だろう. これは悪くないアプローチだが、説明が少し必要だ.そこで詳しくは次回以降にするとしよう.
アプローチ 3: 式を変形! 離散化は素直に
次に考えられるアプローチは、「従属変数を新たに必要なだけ導入して方程式を一階常微分方程式に変形し、その常微分方程式に対する汎用解法を適用する」という手法だ. これはこうした問題の王道なので,以下の内容をよく理解しておこう.
話は簡単で、従属変数を新たに一つ導入すると全体の微分階数が 1階下げられるので(方程式は連立になる), これを必要なだけ繰り返せば良い. 具体的に上の例で示してみよう.
まず、従属変数として、
\[ v(t) := \frac{dh}{dt} \]
を新たに導入する.すると元の方程式そのものは
\[ m\frac{d v}{dt} = -kh -mg \]
と書き換えることができて微分階数を一つ下げられる. よって全体としては,連立ではあるが一階常微分方程式である次の形
\[ \left\{ \begin{array}{rcl} \displaystyle \frac{d h}{dt} & = & v, \cr \displaystyle m \frac{d v}{dt} & = & -kh -mg \end{array} \right. \]
に帰着できるというわけだ.
そして,この変形後の問題はこれまで見てきた一階微分方程式そのものであるので, Euler 法でも時間対称離散化でも Runge-Kutta 法でも,よく知られた汎用解法を適用できる という利点がある.これは大きいね.
この「変形後」の全エネルギーは,$h$ だけではなく $v$ も使った
\[ E(h,v) = \frac{1}{2} m v^2 + \frac{1}{2} kh^2 + mgh \]
と表現した方が良さそうだ.というのも,この式の中に微分が登場しないので扱いやすそうだからね.
実際にやってみよう: Euler 法
Euler 法だと、近似値を $h_n \cong h(n\Delta t)$, $v_n \cong v(n\Delta t)$ として, まあ普通は近似式は次のようになるかな.
\[ \left\{ \begin{array}{rcl} \displaystyle \frac{h_{n+1} - h_n}{\Delta t} & = & v_n, \cr \displaystyle m \left( \frac{v_{n+1} - v_n}{\Delta t} \right) & = & -k h_n -mg . \end{array} \right. \]
プログラムに向いた形に書き換えると次のような感じだ.
\[ \left\{ \begin{array}{rcl} h_{n+1} & = & h_n + \Delta t \, v_n, \cr v_{n+1} & = & \displaystyle v_n - \Delta t \, (\frac{k}{m} h_n + g), \end{array} \right. \]
あとはプログラミングだ.
Euler 法の 1step を関数としていつものように定義しよう.
|
|
エネルギーの計算式も定義しておいて、
|
|
早速動かしてみようではないか!
|
|
近似解をプロットしてみよう.
|
|
う~ん、Euler法の欠点がもろに出ているね. 保存系なのに振幅が大きくなっている… 悪い意味で,厳密解と比較するまでもないな,こりゃ.
まあ「良くない」だろうけど念の為にエネルギーもプロットして確認しておこう.
|
|
う~ん、エネルギーがすごい勢いで「増えている」.こりゃだめだ.
念のために相図も描いておこう.
|
|
本来は保存系であることが明白なのに、数値解は一周してもとに戻ってこない.やはりまずいねえ…
さて,アプローチ 1 の例と同様に,この Euler法も Δt
を十分小さくするとまあ悪くない近似解を出すかもしれない,という期待はできる.
確かめよう.
まず,計算は以下のようになるね.
|
|
まず近似解のプロットだ.
|
|
数値解そのものはさっきよりはだいぶいいね. 誤差も見てみよう.
|
|
誤差はひどいね.こんなに Δt
を小さくしたのに,h_n
の値は2桁ぐらいしかあってないということだね.
エネルギーも見てみよう.
|
|
さっきよりはマシだが,$E = 0.5$ から純粋にエネルギーが増大してしまっている. 相図で様子を確かめよう.
|
|
じりじりと軌道が大きくなっているようだな… 戻ってくるだけ,アプローチ1 の方がマシかな.
というわけで,いまのところ「アプローチ3 + Euler法」は以下のような感じかな.
- 理論との整合性はいいね.
- 数値計算のために式は簡単に作れる.
- プログラムも簡単だ.
- (上で言及してないが)実行速度も速いはず.
- ただ,数値解の誤差はでかいし,本質的に「悪い」.
- エネルギーなどの「系の性質」に関わる量の近似計算も悪いね.
Δt
を小さくすると誤差は小さくなるが,挙動が「戻ってこない」ので保存系の計算には向いてないね.
というわけで,この「アプローチ3 + Euler法」は正直おすすめできないかな. アプローチ 1 より悪いんじゃないかな.
Runge-Kutta 法なら?
というわけで、このアプローチの「汎用解法」の部分を Runge-Kutta 法に切り替えてみよう. つまり、もとの微分方程式を、左辺 = 一階微分 の形式に書き直した
\[ \left\{ \begin{array}{rcl} \displaystyle \frac{d h}{dt} & = & v, \cr \displaystyle \frac{d v}{dt} & = & \displaystyle -\frac{k}{m} h -g \end{array} \right. \]
に対して Runge-Kutta 法を適用するわけだ.
あとはいつもどおりなので、さくさくやってみよう.
まず、微分方程式の右辺を定義する.
|
|
そしてそれを元に Runge-Kutta 法を適用したものを作ってしまう.
|
|
エネルギー計算式も呼び出しやすい形式にしておこう.
|
|
あとは実際に計算するだけだ.
|
|
プロットしやすいように結果を切り出して…
|
|
プロットする.
|
|
おお、これは最初から悪くない.悪くないぞ!
誤差もみてみよう.
|
|
だんだん大きくなっているのは気になるが,誤差の値は小さいね.h
の値は 5桁ぐらい正しい,ということなので,うん,悪くないね.
エネルギーもプロットしよう.
|
|
一方的にさがる傾向はあるが,3桁ぐらいはなんとか保っているか.Δt
を小さくしたわけではないので,まあ悪くないというところかな.
もう少し良くてもいいかなあとも思うが.
相図も描いてみよう.
|
|
ふむ、目で見る限り図の上では一周して戻ってきていると言ってよさそうだ.悪くないね.
というわけで「アプローチ3 + Runge-Kutta法」は,
- 理論との整合性はいいね.
- 数値計算のための式作成は簡単だ.
- プログラムも簡単だ.
- (言及してないが)実行速度は 「Euler法に比べ4倍遅いだけ」で,結構速いだろう.
- 数値解の誤差はけっこう小さいぞ.
- エネルギーなどの系の性質に関する量の近似計算も悪くない.
- 保存系だと,数値計算の状況もほぼきちんと「戻ってくる」.
という感じなので,この「アプローチ3 + Runge-Kutta法」は間違いなく「オススメの方法」だね.
アプローチ 4: 式を変形!離散化も工夫するぞ!
最後に次に考えられるアプローチは、
- 従属変数を新たに必要なだけ導入して方程式を一階常微分方程式に変形
- その常微分方程式に対する 特別な 解法を設計、適用
という手法で、これもまあ王道の一つといえよう.
この場合の 特別 というのはもちろん文脈依存なので、以下に「エネルギー保存性を再現するという意味で特別な」例を示すのでそれを見て雰囲気をつかもう.
まずは、先に示したように、 $v = dh/dt$ という従属変数を一つ新たに導入したときの全エネルギーは
\[ E(h,v) = \frac{1}{2} m v^2 + \frac{1}{2} k h^2 + mgh \]
となるはずだ.
そして、離散化版のエネルギーについては、ある時点 $t = n\Delta t$ での数値解 $h_n, v_n$ に対してやはりそっくりな式
\[ E_d(h_n, v_n) = \frac{1}{2} m v_n^2 + \frac{1}{2} k h_n^2 + mg h_n \]
であると定義しよう.
そして、天下りであるが、例えば、次のような数値計算式で近似値 $h_{n+1}, v_{n+1}$ を計算すると、$E_d$ は変化しない. つまり、$E_d(h_{n+1}, v_{n+1}) = E_d(h_n, v_n)$ となる.
\[ \left\{ \begin{array}{rcl} \displaystyle \frac{h_{n+1} - h_n}{\Delta t} & = & \displaystyle \frac{v_{n+1} + v_n}{2}, \cr \displaystyle \frac{v_{n+1} - v_n}{\Delta t} & = & \displaystyle -\frac{k}{m} \left( \frac{h_{n+1} + h_n}{2} \right) - g \end{array} \right. \]
なぜこの数値計算方法でうまくいくかは,レポート課題になっているよ.
こいつは幸いにも線形な式なので、キレイに変形できてしまう. やってみると、まず、簡単のために $\tau = \Delta/2$ とおいてから,
\[ \left\{ \begin{array}{rcl} \displaystyle h_{n+1} - \tau v_{n+1} & = & \displaystyle h_n + \tau v , \cr \displaystyle \frac{k}{m}\tau h_{n+1} + v_{n+1} & = & \displaystyle -\frac{k}{m}\tau h_n + v_n -2\tau g \end{array} \right. \]
となり、これは線形代数の記号では
\[ \left( \begin{array}{cc} 1 & -\tau \cr \displaystyle \frac{k}{m}\tau & 1 \end{array} \right) % \left( \begin{array}{c} h_{n+1} \cr v_{n+1} \end{array} \right) = \left( \begin{array}{cc} 1 & \tau \cr \displaystyle -\frac{k}{m}\tau & 1 \end{array} \right) \left( \begin{array}{c} h_n \cr v_n \end{array} \right) - \left( \begin{array}{c} 0 \cr 2\tau g \end{array} \right) \]
となる. よって、定数部分をわかりやすく
\[ A = \left(\begin{array}{cc} 1 & -\tau \cr \displaystyle \frac{k}{m}\tau & 1 \end{array}\right), \]
\[ B = \left(\begin{array}{cc} 1 & \tau \cr \displaystyle -\frac{k}{m}\tau & 1 \end{array}\right), \]
\[ \boldsymbol{b} = \left(\begin{array}{c} 0 \cr 2\tau g \end{array}\right) \]
とおくと、この数値計算式は
\[ A \, \left( \begin{array}{c} h_{n+1} \cr v_{n+1} \end{array} \right) = B \, \left( \begin{array}{c} h_n \cr v_n \end{array}\right) - \boldsymbol{b} \]
という、連立一次方程式で書けることになる. 計算上,右辺は既知のベクトルになるので,この連立一次方程式を解けば $(h_{n+1}, v_{n+1})$ が手に入るというシンプルな形であることがわかる.
さてあとはプログラムだ.
まず必要な記号の定義をしよう. 行列やベクトルは、下記のように、手書きの数式と同じように書いておくと紛れがないのでなるべくそうしよう.
|
|
すると、この近似計算式の 1step 分はいきなり次のように書けてしまう.
|
|
逆行列 $A^{-1}$ を使っていないことに注意しよう.連立一次方程式を解いて答えを求めるために逆行列を計算するのは愚の骨頂というもので、理由がない限りやってはいけない.スカラーの場合だって、$24 / 3 = 24 * 0.333333.. = 8$ なんて計算手順だと計算量も多いし不要な誤差も入ってくるのがわかるだろう?
さて、あとはいつものように計算するだけだ.
|
|
プロット用にデータを切り出して、
|
|
まずは近似解そのものをプロット.
|
|
うむ.まあうまくいっているようだな.
誤差もみてみよう.
|
|
う~ん,誤差についてはあまり良くないね.h
としては2桁ぐらいしか合ってないね.
ではエネルギーはどうだろう.
|
|
GKS: Possible loss of precision in routine SET_WINDOW
値の変化が乏しすぎてプロットに警告が出ている. というわけで初期値からのズレをプロットさせよう.
|
|
おお! エネルギー値は 11桁ほど保たれているぞ.ふむ、これは良い.
そして少し面白い.というのも, 数値解としては厳密解からそれなりにずれているのに,エネルギー値はほぼ正しく保存している,ということだからねえ.
一応、相図もプロットしておこう.明らかに,目で見る範囲ではきれいに戻ってくると思うぞ.
|
|
というわけで、このアプローチ4の「謎の近似式」だと不思議な良さがあるね.こんな感じだ.
- 理論との整合性はいいね.
- 数値計算のための式作成はどうなっている? → レポートで!
- プログラムはまあ簡単な方かな.
- (言及してないが)実行速度は少し遅そうだ.非線形問題なら非線形連立方程式を解く必要があるからね.
- 数値解の誤差はそんなに小さいわけではないね.
- エネルギーなどの系の性質に関する量の近似計算はびっくりするぐらい良い.
- 保存系だと,数値計算の状況もきちんと「戻ってくる」.
という感じなので,この「アプローチ4」もまあ「オススメ」して良さそうかな.
問題はどうやってこの謎の近似式を作成するか,だね.それはまずは今回のレポートで取り組んでみよう.
レポート
下記要領でレポートを出してみよう.
- e-mail にて,
- 題名を 2020-numerical-analysis-report-09 として,
- 教官宛(アドレスは web の "TOP" を見よう)に,
- 自分の学籍番号と名前を必ず書き込んで,
- 内容はテキストでも良いし,pdf などの電子ファイルを添付しても良いので,
下記の内容を実行して,結果や解析,感想等をレポートとして提出しよう.
- ここで扱った問題のエネルギー $E(h)$ が保存されることを手計算で確認しよう.
具体的には,$E(h)$ eq.$(\ref{eq:energy-original})$ の時間微分$\displaystyle \frac{dE}{dt}$ を手で計算して、ゼロになることを確認すればいい. なお,この際に厳密解の関数形 eq.$(\ref{eq:strict-solution-form})$ などは知らないという前提で計算すること. 逆に言えば「このエネルギー保存性は厳密解の関数形と関係ない」のだ. - アプローチ 4の解法が、なぜこんなによくエネルギー $E(h,v)$ を保つのか、自分なりに説明してみよう.
ヒント: 前回のレポートによく似た問題があったね.そこのヒントを参考にしよう. - 「アプローチ 3 + Runge-Kutta法」と「アプローチ 4」では
Δt
を 0.001 に小さくした計算を行っていない.
そこで, この 2つの解法でΔt = 0.001
として計算を行い,Δt = 0.1
のときに比べて誤差がどれくらい小さくなるかみてみよう.