Previous ToC Next

7. 常微分方程式の初期値問題(2)

7.1. 今日の予定

今日はまず与えられた精度を達成するためにどのように刻みを調整するかとい う話をしたあと、ハミルトン系に特化した数値解法の話に移る。

7.2. 刻み幅調節と埋め込み型公式

ルンゲクッタにしても線形多段階法を使うにしても、実際に問題を解こうとい う上ではいろいろ考えないといけないことがある。もっとも重要なことは、刻 み幅 をどうやって決めるかということである。

実際問題としては、解いていくにつれて導関数の値が大きく変わるといったこ とは普通に起きる。このために、刻み幅が一定のままで計算していくというの は非常に無駄なことが多い。

無駄を減らすためには、局所誤差が小さい(解が滑らかな)ところは刻みを大 きく、逆に解が急にかわるところでは刻みを小さくしてやればよい。このよう な方法を可変刻み (variable step) 、あるいは適応刻み (adaptive step) と いう。

Figure 6: 適応的刻み幅の概念。解が急に変化する時に刻み幅が大きいまま では上手く数値解が求まらないが、変化に合わせて刻み幅を小さくすればそれらしい解が求まる。

7.2.1. RK法の場合

例えば Runge-Kutta 法の場合、 から にいくのに前に計算 した情報とかなにかを使うわけではないので、刻み幅は任意にとることができ る。問題は、どうやって刻み幅を決めるかである。

これにはいろいろな考え方があるが、通常使われるのは、局所離散化誤差を推 定して、それがある値以下になるようにするという方法である。

このための一つの考え方は、以下のようなものである。

  1. まず、ある幅 で積分する。
  2. 次に、その同じ幅を、 ずつ 2回に分けて積分する。
  3. すると、2回に分けて積分したほうがより正確な答になっていると期待できるので、二つの解の差が誤差の推定値(大きめではあるが)となっている。

誤差の推定値が求まったら、普通は次のステップを調節する。局所離散化誤差 の次数がわかっているから、それに応じてステップサイズを調整すればよい。 次数プログラムによっては、現在のステップでの誤差が予定した値よりも大き ければステップを小さくして計算し直すようになっているものもある。

実装(プログラムを作ること)が容易であることもあって、この方法はわり と広く使われている。が、かなり無駄が多い方法であるというのも確かである。 というのは、単に誤差の推定のためだけに、 50\% の無駄な計算をしているか らである。まあ、 50\% はたいしたことないという考え方もあるが、計算に何 日も掛かるというような状況なら 50\% は無視できないので、まあ、なんとか したくなるのが人情というものである。

というわけで、もうちょっとうまい方法はないかということで、以下のような ことを考えた人がいる。

一般に、 RK型の公式では、最終的な値は

の形になっている。ここで、 は全部そのまま使って、 を別の に置き換えた

で、局所離散化誤差が元の公式よりも大きいが、むやみと大きくはない(例え ば、一次次数が低い)ものがあれば、元の公式との差を誤差の推定に使うこと ができる。

この形の公式のことを埋め込み型 embedded 公式という。最初に提案した人の 名前をとって Runge-Kutta-Fehlberg 公式ということも多い。これもいろいろ 提案されているが、前回にも紹介した Dormand-Prince 公式はすべて埋め込み型であり、こ れが最近はもっとも広く使われている。精度が高いものが必要な時は、 8次の公式が使われる。これらは

  http://www.unige.ch/math/folks/hairer/software.html
から入手できる (Fotran と C がある)。

7.2.2. 線形多段階法の場合

線形多段階法の場合、 RK と違って誤差の推定は容易である。例えば、アダム ス法で予測子・修正子ペアで使うなら、予測子と修正子の差は誤差の(オーダー としては)よい推定値になっている。したがって、この部分については面倒な ことは特にない。

問題は、実際に刻みを変えるほうである。もっとも一般的な方法は、「その場 で公式を導出する」というものである。これについて以下簡単に説明する。

7.2.3. ニュートン補間とアダムス公式の一般的導出

アダムス法を使うには、要するに多項式補間の系数が出せればよい。以下、ニュー トン補間を使う方法を説明する。

補間したい関数 が、標本点 で、関数値 をとるとする。これを、以下のような形の多項式

という形に書くことを考える。これを、 までをとったものは <$x_0 \cdots x_{p}$> までを通る最低次補間多項式になっているように構成すること にする。で、低い次数から順に作っていく。付け加える新しい項は、それまで に使った点すべてで0になるので、新しい点で標本と一致するように の 値を決めればいい。

まず、 については、 とすればいいのは明らかであろう。 次に であるが、

から

となる。同様に、 を求めると、

ただし、

である。

さて、段々式が繁雑になるが、もうちょっとの辛抱である。上の形は、

というふうにも書き直せることに注意しよう。これは、を通る補間式 の一次の係数と、 を通る補間式の係数の差みたいなものになっている。

ここで、天下りに、 階差商 (divided difference) <$D_{k,l} \quad (0\le l \le n-k)$> というものを導入する。これは以下のように定義される。

実は、このように定義された が、

を満たす、求める係数であることを示すことが出来る。

証明には、が、 から始まる 点を通る補間多項式の係数 であるということを用いる。定義により

である。 ここで、帰納 法を使って証明することにすると、

は、 を通る 次補間多項式であり、また、

は、 を通る 次補間多項式であるとしてよい。こ の2つの線形結合によって、 を通る補間多項式を作るに は、

とすればよい。ここで、は、それぞれ の最高次の係数であることに注意すると、 の最高次の係数(これは に等しい) は、で与えられることがわかる。

さて、

という多項式を考えると、これは と一致し、 さらに最高次の係数も一致しているので、結局 であることがわ かる。つまり、 が求める補間多項式であり、 がその最高 次の係数ということになる。

と、こんな風に差商を作っておくと、これから多項式の係数を機械的に計算し て積分して、、、、ということができる。さらに、差商の形でデー タをとっておくと、上の漸化式を使ってを順番に更新することができる。 これは Krogh 型の公式と呼ばれ、かなり広い範囲の問題について、{\bf もっ とも効率が良い公式}であることがわかっている。

なお、上の方法で時間刻みを変えられるなら、出発公式は不必要になることに 注意しておく。つまり、最初は1次の予測子、2次の修正子から出発して、ステッ プ毎に次数をあげていけばいい。ステップサイズは次数を考慮して決めておく。

7.3. ハミルトン系向けの方法

7.3.1. 簡単な例題

古典力学系のもっとも簡単な例といえば、一次元調和振動子、つまり、運動方 程式では

である。 これの解はもちろん単振動するものであり、それは固有値が純虚数であるから である。

さて、このような、固有値が純虚数である、いいか えれば中立安定な場合、安定性解析はちょっと厄介になる。というのは、通常 の意味で安定というのは、数値解がいつかは原点にいってしまうということ、 いいかえればエネルギーが保存しないということを意味するからである。とい うのは、安定であるとは、線型差分方程式の固有値の絶対値が全て 1 より小 さいということだからである。

つまり、安定領域のある公式では、ほとんどのもので虚軸上の一定範囲では数 値解が「安定」になってしまう。つまり、本当の解はいつまでも振動を続ける のに、数値解は減衰してしまうのである。もちろん、タイムステップを大きく とれば、安定領域から外れるのが普通であり、振動しながら 広がっていくような解になる。安定領域に虚軸を含まないような公式(例えば 前進オイラー公式)の場合では、かならず軌道が広がっていくことになる。つ まり、虚数軸が安定性限界になっていない限り、必ず軌道から一方的に外れて いってしまうのである。

安定領域とはそもそも何かという話をしておく。線形常微分方程式

を刻み で積分した時にi、これは結局 という形の一 般解を持つ。ここで、すべての固有値 の絶対値が1を超えないよう な、複素平面上での の領域のことを安定性領域という。

が複素数の場合も考えるのは、元の微分方程式が連立方程式なら固有値が 実数でない場合もありえるからである。

つまるところ、線形系では安定性限界が虚数軸である(数値解に純虚数の固有 値がある)というのが、数値解が厳密に振動的であるための必要十分条件とい うことになる。この時、数値解はもとの系の性質をある意味で良く表している といえるであろう。

実は多変数の場合や非線形の場合でも本質的な事情は変わらない。ただし、こ ちらはまだいろいろ理論的にはっきりわかっていないことがいろいろある。が、 大雑把にいうと同じような事情が成り立っている。つまり、もとが周期解であ るときに、大抵の積分法では周期解から一方的にずれていく、つまり、保存量 であるはずのエネルギー等が保存しなくなる。

もちろん、エネルギーが厳密には保存しないということが問題であるかどうか というのも実は難しい問題である。というのは、保存しないといってももちろ ん数値解の精度では保存しているわけだから、それでかまわないのではないか とも考えられるからである。

7.3.2. リープフロッグ公式

ハミルトン系のシミュレーションの例には例えば以下のようなものがある

これらはいずれも、古典的な多体問題として定式化される。プラズマの場合は 磁場も入るのでちょっと面倒だが、それを無視すれば運動方程式が

つまり、ある粒子 の加速度は他のすべての粒子からの力の合計というこ とになる。もちろん、古典力学系は他にもいろいろあるが、基本的なアプロー チは似たようなものである。

このような問題に対して、もっとも普通に使われる公式は以下の leapfrog といわれるものである。

これでは速度と位置がずれた時間でしか定義されない が、出発用公式として

を使い、さらに終了用公式として

を使うことで最初と最後を合わせることが出来る。この形は、実は

と数学的に等価である(証明してみること)また、以下のようにも書ける

さらにまた、速度を消去して, , の関係式の形で 書いてあることもあるかもしれない。なお、すべて違う名前がついていたり するが、結果は丸め誤差を別にすれば完全に同じである。時々、これらが違う ものであるかのように書いてある教科書があるので注意すること。

この公式は局所誤差が 、大域誤差がである。

さて、局所誤差という観点からはこれは決して良い公式というわけではないが、 現実にはこの公式は非常に広く使われている。

これは、別にもっとよい方法を知らないからとかではなく、実はこの方法がい くつかの意味で非常に良い方法であるからである。いくつかの意味とは、例え ばエネルギーや角運動量のような保存量が非常によく保存するということであ る。これらの保存量については多くの場合に誤差がある程度以上増えない。

この、ある程度以上誤差が増えないというのは目覚しい性質である。通常の方 法では、エネルギーの誤差は時間に比例して増えていく。従って、長時間計算 をしようとすればそれだけ正確な計算をする必要がある。ところが、エネルギー の誤差は溜っていかないのならば、かならずしも精度を上げる必要はないとも 考えられる。

もちろん、エネルギーが保存していればそれだけで計算が正しいということに はならない。が、理論的には、いくつかの重要な結果が得られている。まず、

1.leapfrog は symplectic method のもっとも簡単なものの一つである。

  1. leapfrog は symmetric method のもっとも簡単なものの一つである。

Symplectic method については、以下のようなことが知られている

  1. symplectic method は、すくなくともある種のハミルトニアンに対して使った場合に、それに近い別のハミルトン系に対する厳密解を与えることがある。
  2. 周期解を持つハミルトン系に対して使った場合に、どんな量でも誤差が最悪で時間に比例してしか増えない。
  3. 時間刻を変えると上のようなことは成り立たなくなる

これに対し、 symmetric method については以下のようなことが知られている

  1. 周期解を持つ時間対称な系に対して使った場合に、どんな量でも誤差が最悪で時間に比例してしか増えない。
  2. 時間刻を変えてもうまくいくようにすることも出来る。

以下、まず数値例でこれらの性質を確認し、それからなぜこのようなうまい性 質を持つのかについて直観的な説明を与えよう(厳密な説明/証明には力学系 の積分可能性に関するかなり深い知識が必要となる)。

7.3.3. 数値例

\subsubsection{調和振動子}

能書きを聞いていても良くわからないので、実例を見てみよう。まず、簡単な 例として調和振動子

を leap frog と 2階のルンゲクッタで解いた例を示す。初期条件は <$(x,v) = (1,0)$> で、時間刻みは である。

軌道とエネルギーを図に示す。非常に特徴的なのは、 leap frog ではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている ことである。

Figure 7: 調和振動子の数値積分。軌道。破線は2次のルンゲクッタ、実線は leapfrog の結果

Figure 8: 調和振動子の数値積分。エネルギー。破線は2次のルンゲクッタ、実線は leapfrog の結果

ルンゲクッタでは単調に変化するというのは前に説明した通り (2階のルンゲクッタでは虚軸を安定領域に含まないため)である。 さて、これに対し、 leapfrog ではエネルギーが変化していないが、これはど ういうことなのであろう?

実は、この調和振動子の場合には、 leapfrog 公式は以下の量

を保存するということを確かめることが出来る。つまり、 で与えら れる位相平面の上で考えると、 leapfrog 公式の解は上の式で与えられる楕円 の上にのっているのである。このために、エネルギーの誤差がある値よりも大 きくなり得ないことになる。

7.3.4. 非線形振動

では、非線形振動ではどうだろう?簡単な例として

を leap frog と 2階のルンゲクッタで解いた例を示す。初期条件は <$(x,v) = (1,0)$> で、時間刻みは である。

軌道とエネルギーを図に示す。調和振動の場合と同様に、 leap frog ではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている。

Figure 9: 非線形振動の数値積分。軌道。破線は2次のルンゲクッタ、実線は leapfrog の結果

Figure 10: 非線形振動の数値積分。エネルギー。破線は2次のルンゲクッタ、実線は leapfrog の結果

この場合も保存量があるので、頑張れば求まるかもしれない。

7.3.5. シンプレクティック公式

さて、 leapfrog 公式は、上で見たようにハミルトン系に対してエネルギー誤 差が有界に留まるという大きな特長がある。が、2次精度でしかない。もっと 高次の方法はないのだろうか、またあるとすればどういう原理で作れるのだろ うか。

高次の方法を構成する一つのアプローチが、シンプレクティック公式といわれ るものである。これはなにかというと、積分公式がシンプレクティック写像に なるように作るということである。

シンプレクティック写像とは、ようする に正準変換のことである。正準変換とはなにかというのは解析力学を思い出し てみて欲しいが、直観的には、力学系を不変に保つ、つまり変換まえの座標系で 求めた軌道を変換したものと、変換後の座標系での力学系の軌道が厳密に同じ になるようなものである。

ハミルトン力学系の解そのもの(ある時刻 での座標から、 での座 標に移す変換)もシンプレクティックである。まあ、だから、シンプレクティッ クになっているような積分公式は、そうでないものに比べてなんとなく力学系 の性質にあっているような気はする。

で、いいたかったことは何かっていうと、上の leapfrog 公式はこのシンプレ クティック性を満たしているということだった。詳しくは、「数理科学」 1995年6月号にのった吉田春夫さんによる解説記事でもみてもらうことにして、ここで は高次の公式にはどんなものがあるかという話をしておく。

7.3.6. 陽解法

陽解法の組み立て方はいろいろある。一つは、 RK 系の公式で、係数をシン プレクティック性を満たすように決めるということである。これはここ15年で 無数に論文がでた。4次、あるいは6次の公式としては、吉田や鈴木による作用 素分解に基づく公式が良く知られている。

これらの方法の原理は、要するに上の leap frog をタイムステップを変えて いくつか組合わせるというものである。うまくタイムステップを組み合わせる と誤差の高次の項を消すことができるわけである。3段4次の公式、7段6次の公 式等が吉田によって導かれている。

なお、実は陽解法はハミルトニアンが の形の場合にしか使え ないが、大抵の問題はこう書けることはいうまでもないであろう。

また、 RK 系の公式を力任せに構成する試みもあり、4次から8次までの公式 が作られている。計算精度という観点からはこちらのほうが leapfrog を組合 わせるものよりもいいものもある。

7.3.7. シンプレクティック公式の意味

さて、シンプレクティックであるということと、「良い」ということの間には ちゃんとした理論的な関係があるのであろうか?一応あるということになって いる。つまり、あるハミルトニアン で表される系をシンプレクティック な 次の公式を使って積分した解は、別のハミルトニアン で与えられ るシステムの厳密解になっていて、 の間に

という関係がある(を求める数列が収束すれば)ということがわかってい る。

なお、例えば上の調和振動子の場合には実際に数列が収束し、 が求まっ ている。が、これは極めて例外的な場合で、一般の場合には収束するかどうか も明らかではないようである。

収束するかどうか明らかではないのでは、使っていいことがあるという保証は ないではないかと思うかもしれない。理論的にはその通りなのであるが、実際 にはいろいろな問題に適用されて、従来の方法よりも高い精度が得られるとい うことが確認されている。

7.3.8. シンプレクティック公式の問題点と対応

さて、式(\ref{eq:shadow}) をみるとわかるように、シンプレクティック公式 に付随するハミルトニアン には時間刻み が入っている。従って、 をふらふら変えると も変わって、結果的に求まった数値解は一つの 力学系の軌道ではないなんか変なものになってしまう。ということは、可変時 間刻みにするとシンプレクティック公式はうまく働かないのではないかという ことが想像される。

じつはその通りで、例えばシンプレクティックで埋め込み型のルンゲクッタと いうものを作って、実際に時間刻みを変えてみた人がいる。その結果、普通の ルンゲクッタよりも良くならないということが発見された (1992年頃)。問題 によってはこれは致命的な欠陥となるので、さまざまな対応法が精力的に研究 されているのが現状である。が、あまり一般的にうまくいく方法というのは見 つかっていないようである。つまり、 万能な公式というのはシンプレクティック公式に関する限り知られていない。

もう一つのシンプレクティック公式の問題点は、「写像」であるということか ら一段法、具体的にはルンゲクッタ型の公式であるので、局所誤差に対する計 算量という観点からは線形多段階法に比べて必ず悪いということである。

これらを解決しようという研究はいろいろあるが、一つの方向が以下に述べる対称型の公式である。

7.3.9. 対称型公式とは

さて、まず対称型公式とはなにかということだが、まず時間反転に対する対称性というものを定義しておこう。時間反 転に対する対称性とは、常微分方程式

に対する一段法

が、

を満たすこと、つまり、直観的には、ある微分方程式系があって、それを1ス テップ分数値解を求めたとする。で、そこから逆に戻ると厳密に元の値に戻る ということである。(ここでは丸め誤差はないものとする)

で、一段法の場合には、この意味で対称なものを対称型公式ということにする。

具体例で考えてみよう。例えば前進オイラー法は対称型ではない。これは、いっ た先で導関数を計算すれば、一般にもとのところでの導関数とは違うから当然 であろう。前進オイラー法が対称ではないのだから、その逆写像である後退オ イラー法も対称型ではないことになる。

では、対称型にはどういうものがあるかということだが、明らかに対称型であ るものとしては、台形公式

がある。これが上の対称性を満たしていることはいうまでもない。

さて、なぜこの型の公式を考えるかということであるが、 実験的には以下のようなことが知られている。

前のほうの性質はシンプレクティック公式と同じであるが、後の方の性質はシ ンプレクティック公式よりもある意味でよいものである。

ここではとりあえず対称型公式にはどんなものがあ るかという話をしておこう。

さて、一階の方程式に対する一段法という制限をつけた場合、つまり、ルンゲ クッタ法の場合、対称な公式はかならず陰的公式になる。逆に、陰的公式であ れば対称型のものを作るのは容易である。つまり、ルンゲクッタ法を与える行 列 とかベクトル , が対称になっていればいい。

一階の系に対する一段法では、陰的 RK のほかにうまい方法があるわけではな い。それでは、

以下、順番に考えていくことにする。

7.3.10. ハミルトン系用の陽的対称型RK公式

実は、すでに述べたように、 leapfrog 公式は対称型でありしかも陽公 式である。で、 leapfrog の組合わせで作られる公式も、実は すべて対称型になっている。

7.3.11. 対称型線形多段法

線形多段階法でも、対称型というものを考えることはできる。但し、これら の公式は、今のところ

ということが知られており、実用になっていない。これは 国立天文台の福島さんのグループにより精力的に研究が進められている。

7.3.12. エルミート型公式

一般に対称型というわけではないが、4次の対称型公式を導く特別な方法とし てエルミート型(あるいはエルミート・オブレヒホフ型)と呼ばれる公式のク ラスについてここで触れておこう。

エルミート型公式は、線形多段階法のアダムス型公式の一つの一般化である。 どのような意味において一般化であるかというと、アダムス型と同じように補 間多項式を積分することで新しい時刻での値を求めるところは同じである。違 うところは、アダムス型公式では補間多項式を作るのに関数の値を使う(ラグ ランジュ・ニュートン補間)が、エルミート型公式では、関数の導関数の値も 使うのである(エルミート補間)。

もっとも単純な発想としては、テイラー法というものがある。つまり、微分方 程式

というものがあったとしたら、

という関係式を使って の全微分を求めることが出来、同様に高次の導関数 もどんどん計算していくことが出来る(微分が式で書ければ)。これらを使っ て直接テイラー級数を評価して近似解を求めるのである。

高次の導関数が簡単に求められる場合(例えば線形な系とか)には、この方法 はかなり強力である。線形多段階法やルンゲクッタに比べて誤差項の係数がずっ と小さく、タイムステップが大きくとれるからである。

が、多くの問題において、高階導関数の直接計算は現実的ではない。これは、 計算量が指数関数的に爆発するのが普通だからである。とはいえ、の一階 導関数や二階導関数くらいならば指数関数的といってもたかがしれている。し かしながら、これではテイラー法では2次や3次の公式しか作れない。実用的な 公式としては、(leap frog のような特別に良い性質をもっているとかいうの でなければ)もうちょっと高次のものが必要である。

そこで考えられるのが、高階導関数を使ってさらにルンゲ・クッタ型とか線形 多段階法のような公式を作ることである。例によってもっとも簡単な場合とい うのを考えると、それは のところで の一階導関数 の値を使うもの(高次導関数を使わないものであれば台形公式に相当)の一般 化ということになろう。

台形公式は、2次の陰的アダムス型公式と考えることができる。つまり、 を線 形補間してそれを積分しているわけである。これに対し、 と <$f' = df/dx$> を使う補間というのはどういうものかということを考えてみると、こ れは 3次多項式を構成できることがわかる。具体的には、導関数の近似多項式 を

としたとき、

という形で係数を求めることが出来る。さらに、上の近似多項式を積分してや れば、結局

という形の公式が得られる。これは台形公式と同じく陰的公式である。対称型 であることは明らかであろう。

この公式はプログラムが簡単であるわりには精度が高く、また対称型であるというよ い性質をもつこともあり、良く使われるようになりつつあるようである。収束 させるための反復には普通の直接代入が使える。

7.3.13. 対称型一段公式における時間刻みの変更

対称型一段公式には、対称性を保ったままで時間刻みを変更できるというシン プレクティック公式では(普通の意味では)なかった良い性質がある。以下で は、どのようにしてそれが可能であるかを簡単に述べておく。

今、任意の対称型一段法があって、その時間刻みを与える関数として が与えられているとする。ここで はなんでもいいが、これも一段法 である、つまり前のステップの値とかを必要としないものであるとする。

一般に、 によって刻み幅を決めて求めていった数値解は、時間反転 に対する対称性を満たしていない。つまり、1ステップ積分してから、逆に戻 すと、タイムステップが変わってしまうためにもとのところに戻らない。この ことのために、例えば周期軌道の場合に誤差が周期的にならないということが 起こるわけである。

刻み幅を変えつつ時間反転に対する対称性を保つ一つの方法は(これ以外の方 法もあるかもしれないが、今のところ知られていない)、 自体に対 称性を持たせることである。つまり、一段法を

という形に書いた時、

がなりたつことを保証するように を決めるのである。具体的には、上の 対称性が成り立っていないような時間刻みの式 があった時に、

によって対称化した時間刻みを作ればいいことになる。この場合時間刻み自体 が陰的に決まることになるが、とにかく対称性という要求は満たされるのであ る。

7.4. 練習

プログラムが必要なものはプログラムを提出 すること。

7.4.1. 練習1

2次元ケプラー問題

について、古典ルンゲクッタと leap frog の両方で、適当な初期条件(解が 楕円になるものをとること)から出発して 10, 1000, 周期後の解の厳 密解からの誤差が刻み幅のどのような関数になるか調べよ。ここでは刻み幅は 固定とする。

7.4.2. 練習2

4次のシンプレクティック公式は、 が刻み幅 の leap frog を表す作用素であるとして、以下のように書ける。

つまり、リープフロッグでまずちょっと行き過 ぎて、次に逆に戻って、最後に最初と同じだけ進んで次の時刻に行くというも のである。これをプログラムし、 1 と同様な解析を行なえ。

7.4.3. 練習3

4次のエルミート公式についても同様な解析を行なえ。説明されている ものは陰的公式なので、予測子に何を使えるかは各自検討すること。

7.4.4. 練習4

2.1 で述べた古典的 RK 法で誤差の調整をする方法をプログラムし、離 心率の大きなケプラー問題について刻み一定の場合にくらべてどの程度得にな るか調べてみよ。
Previous ToC Next