今日は、「硬い」方程式系について、それはどういうものかと、ど うやって解くかという話をする。

さて、硬い微分方程式とは何かということだが、例えば以下のような化学反応 系を考えてみる

係数がたくさんあって面倒なので、以下 <>, <> とし、初 期条件として

なお、この計算例では、保存則

化学反応系のような非線形の方程式系では、上の例のように速度定数が大きく 違うということはそれほど珍しいことではない。上の例では であ るが、これが何桁も大きい場合というのもある。例えば、化学反応といっても 天体物理で出てくるような核化学反応を考えると、速度定数が 10桁くらい違 うものが出てくることはそれほど珍しくない。

このような、非常に時間スケールが違う方程式が連立しているような系を「硬 い」系と呼ぶ。こういった場合にどういう方法を使うべきかというのが今日考 えていきたいことである。

最初に見たように、例えば化学反応系などでタイムスケールが非常に違うと具 合の悪いことが起きるわけであるが、この具合の悪さというのは具体的にはな んだろうか？また、これはどういう意味で「硬い」のであろうか？

例によって、方程式が線形の場合について考える。非線形の場合は、局所的に 線形化することで硬さを定義することになる。非斉次の線形常微分方程式

の一般解は

さて、この硬さというのが問題かどうかということであるが、仮にある公式の 絶対安定領域が左半平面 全体を含んでいれば、硬かろうがなん だろうが問題ではないということはすぐにわかる。この、絶対安定領 域が左半平面を含むという性質を、ダールキストは A-安定性と名付けた。

A-安定性については、ダールキスト自身による次の否定的な結果がある

• 陽的ルンゲ・クッタ法は A-安定でない。
• 陽的線形多段階法は A-安定でない。
• 陰的線形多段階法でA-安定であるものの次数は 2 以下である。
• A-安定である陰的線形多段階法のうち、局所誤差がもっとも小さいものは台形公式である。

但し、陰的ルンゲ・クッタ公式については以下の素晴らしい結果が知ら れている

• 対称な 段 次公式（陰的ガウス公式）はすべて A-安定である。

陰的 RK 法はどんなものかという話を簡単にしておく。陰的 RK には陽的 RK 以上に無限にいろんなバリエーションがあり得るが、特に硬い方 程式に使われる方法はほとんど collocation method といわれるものになって いる。Collocation method の基本的な考え方は、「数値積分の公式をそのま まなんとか常微分方程式の解法に置き換えるために、数値積分を使って RK 法 の途中の点での解を求める」というものである。

といっても、これではなんだかわかりませんね。例によってまずもっとも簡単 な場合というのを考えてみる。

簡単な数値積分の公式といえば、中点公式か台形公式である。どちらも2次精 度である。台形公式はそのまま数値解法になってしまうので、以下、自明では ない例として中点公式を考える。

中点公式は、

この困難は、数値積分が近似多項式の積分として与えられたということを思い 出せば一応解決できる。中点公式の場合、数値積分ではこれはを で近似するということに相当した。同様なアイディアを微分 方程式に適用すれば、微分方程式の場合、 を定数で近似し、解で ある y を一次式で近似する、つまり

数値積分と同じように、途中にとる点の数（RKでは段数に相当）をあげれば近 似の次数が上がって精度をあげられる。普通の考え方は、点が個なら 次多項式ができる。ここでは途中の点の 座標を と 書き、構成された次多項式を と書くことにすると での「解」 について

と書くと非常に大変そうだが、係数をあらかじめ計算しておけば上の手続きが 陰的ルンゲクッタ公式として表現できる。

途中の点のとりかたであるが、普通に考えつく方法は等間隔のラグランジュ補 間で、両端を含むようなものである。これは点の数 に対して次数が で、まあ悪くはない。また、両端の点を含むので、左端は計算しなくていいこ とになりちょっと嬉しい。

が、数値積分については、うまく点をとると 個で最大 次を達成できると いう素晴らしい結果が知られている。もっとも高い次数を達成できるのは、個 の点を次ガウス・ルジャンドル多項式の零点にとるもので、これを普通陰 的ガウス公式という。これは2点で4次とか4点で8次の精度が達成できるなかな か結構な公式である。

さらに結構なことに、上に述べたA-安定になっているわけである。

さらに、ルジャンドル多項式の性質（対称性）から、陰的ガウス公式は時間反 転に対して対称であり、さらにシンプレクティックでもあるとい うことがわかっている。なお、これらのことからわかるように、陰的ガウス公 式では安定領域がちょうど左半平面全体であり、虚軸が境界になっている。

この、虚軸が境界になっているというのはまあ結構な性質ではあるが、問題に よってはもっと強力に振動を押え込みたいということもある。このような場合 には、次数が1下がるが Radau 法のような、右端の点を含む公式を使うことが ある。このへんのコードも前に紹介した Web サイトにある。

前節で見たように陰的ガウス法は A-安定であり、基本的にはこれを使えば方 程式が硬いということによる問題は起きないといっていい。が、実用上は、こ の方法にはあまり嬉しくない特性、つまり、大規模な非線形代数方程式を解か なければいけないという特性がある。

もちろん、線形多段階法のところで考えたような逐次（直接）代入法で答が収束すれば、 これはたいした問題ではない。しかし、実はこれはうまくいかない。つまり、 陰的ルンゲ・クッタに対する逐次代入が収束する条件は リプシッツ連続条件

それではいったいどうすればいいかということだが、要するに逐次代入なんて いう手抜きな方法を使わないで、もうちょっとまっとうな方法を使えばいい。通常使われるのはニュートン法による反復、つまり、方程式を局所的に線形化し、 その線形方程式の解を新しい近似解にする方法である。

原理を後退オイラー法の場合について説明しよう。もちろん、高次の陰的ルン ゲクッタやあとに述べる陰的線形多段階法でも理屈は同じである。後退オイラー 法の場合、解くべき方程式は

すぐにわかるように、ニュートン法は方程式が線形であれば硬さに無関係に一 回で収束する。つまり、逐次代入法とは違ってA-安定性を保てる。

方程式が非線形の場合にうまくいくか、また、どれくらいの反復回 数でうまくいくかというのはなかなか難しい問題である。良く知られているよ うに、ニュートン法は二次収束する、つまり、初期推定が十分に良ければ誤差 が一回反復するとまえの誤差の2乗程度まで小さくなるという性質（これは、2 次の項の大きさからすぐに証明できる）があるので、うまい初期値が得られれ ば反復回数は少ない。

が、うまい初期値をとる一般的な方法はない。例えば陽的ルンゲクッタや陽的 な線形多段階法で初期値を作ることも考えられるが、これは系が硬い時には悪 い推定値を出す可能性もあるからである。比較的うまくいくのは、なにも考え ないで とする、また陰的ガウス公式であれば途中の 値についてもそうする方法である。

なお、もともとのニュートン法ではヤコビアンを反復毎に計算し直すが、これ らの方法を実装する時には、「収束しているうちはヤコビアンをそのまま使う」 という方法で計算量を減らすのが普通である。

ここまでで見たように、ガウス法などの陰的ルンゲクッタは硬い方程式にたい して非常によい性質をもつが、一つ大きな問題がある。それは、ニュートン法 で解く方程式の大きさが、元の方程式の本数 と公式の段数 の積 に比例するということである。

ニュートン法で方程式を解くときに、 LU分解（普通のガウスの消去法のこと と思っていい）を使うとすれば、計算量が行列の次元の3乗で増える。つまり、 公式の段数の3乗ということになる。これでは、せっかく高い次数の公式が作 れるといっても実用性は低い。

ニュートン法ででてくる線形化した方程式自体を、LU分解ではなく適当な反復 法、たとえば Gauss-Seidel 法、 SOR、共役勾配(CG)法といったもので解くこ とも考えられるが、これらを使うとまたA-安定性が失われたりすることになる。

これらの問題をあるていど回避するため、非常にさまざまな公式群が研究され ている。以下にその例をあげる

前に述べたように、計算量と精度の関係を見た時に線形多段階法はルンゲクッ タに比べてかなり良いのが普通である。特に陰的公式の時には、解くべき代数 方程式が次数が増えても大きくならない。従って、線形多段階法で比較的よい性 質をもつものはないかということを考える。

ダールキストの結果によって、A-安定では2次以上にならないということがわ かってしまっているので、A-安定の条件のほうをちょっと緩めてみる。緩め方 にはいろいろあるが、ギア (C. W. Gear) の提案は、図\ref{fig:sstable}の ように原点近くで虚部が大きいところでは不安定でもいいということにしよう というものである。

係数の表は適当な参考書をみること。これは、現在でも多くのライブラリやパッ ケージソフトで使われている。

実際のプログラムの例は例えば http://www.netlib.org/ode/cvode.tar.gz といったものを見てみること。

これは無数にある。名前だけをあげると、ローゼンブロック法、 DIRK, MIRK, SIRK, SDIRK といったものがあり、現在盛んに研究されている。

まあ、実用上は、完全陰的なガウスや Radau 公式では問題があるという場合 でなければ、わざわざこれらの方法を使う意味はあまりないように思う。

常微分方程式の話はこれでおしまいにして、最適化の話。なお、12/4 は休講 にします。

陰的中点公式の安定性領域を求めよ。

van der Pol 方程式

の数値解を、 の時に適当な陰的公式と陽的公式 （後退オイラーと前進オイラーでもOK）で求めてみよ。

2段4次の陰的ガウス公式について

1. ルジャンドル多項式の 0 点をとるということから公式を導け。

2. 実際に線形系に適用して、4次精度であることと A-安定であることを確かめよ。

3. 2変数以上の場合に適用出来るプログラムを作り、適当な非線形方程式を解いてみよ。例えば van der Pol 方程式で極度に非線形性の強い場合について、4次の古典的ルンゲクッタと解の精度、信頼性等を比較してみよ。

   Previous

   ToC

   Next