微分方程式の数値解法

2階の微分方程式で、
例えば、初期条件が、y(0)=0、y(1)=0のような
「globalな」ものだった場合って、どのように数値的にとけばよいのでしょう？
普通、数値解法って、y(0)=0、y'(1)=0のような
局所的に折れ線をつなげていくじゃないですか・・・？
質問の場合には、その解法だと恐ろしく面倒ですよね。

No.1 ベストアンサー 回答者： chukanshi 回答日時： 2002/03/08 19:24

このような場合には、

y(0)=0,y'(0)=k
と置いて、kをいろいろな値にとって
微分方程式を数値的に解きます。
そして、y(1)=0となる条件の元で、
kの値を決めてやります。そうして
求まったy(x)が解です。

プログラムでかくと、

定義
eps=10E-8
y10=-10 :y'(0)の初期値
step=1.E-2 :y'(0)を変化させる１ステップの大きさ

y00=0 :初期条件
y01=0 :初期条件

本体

While(k<10)

サブルーチン
（y00=0,y10=kを初期条件として微分方程式を解く。
解いた結果、y(1)をy0endとして返す。）

if(|y0end-y01|<eps)

Then End（プログラム終了）
else(k=k+step)

End while

みたいな感じで解きます。

アルゴリズム的には、「shooting method」といいます。
この言葉は、日本語の教科書には余り出ていませんが、
英語の数値計算の辞書のように厚い本には、
「shooting method」と出ています。

たとえば、こんな本には出ています。

Numerical Recipes in Fortran
William H. Press (編集), Saul A. Teukolsky (著), William T. Vetterling (著)
ハードカバー - 963 p / 2nd 版 (1992/01/15)
Cambridge University Press ; ISBN: 052143064X

読んでみてください。

この回答へのお礼

ありがとうございます。
やっぱりそういう方法しかないですよね。
なんか超絶技巧な技があるのかと思いましたが。

お礼日時：2002/03/09 12:03

No.2 回答者： tgb 回答日時： 2002/03/09 21:34

　線形問題にしか適用できませんが、学生の頃、境界値問題

の解法の１つとして何かの本に載っていたのを見たことがあ
ります。この方法は両端の境界条件でｙやｙ'を組み合わせて
指定する場合にも応用できます。

・ｙ（０）＝０となる解ｙａを１つ求めます。
・ｙ（１）＝０となる解ｙｂを１つ求めます。
・ｙ＝ａ＊ｙａ＋ｂ＊ｙｂとして
　境界条件を満たすようにａ、ｂを決定します。
（実質的にはそれぞれ反対側の境界条件を満たすよう係数を
　決定している）

　境界条件の値を０とする解を選ぶことにより、係数を決める
時片側だけに寄与するようにしたのがミソですが、質問のよう
な両端で０と言う条件に対してはうまく行かないようです。そ
れ以外ならすっきりした解法だと思いますがいかがでしょうか？

この回答へのお礼

なるほど！確かに線形なら可能ですね。
ありがとうございます。

お礼日時：2002/03/10 06:37

No.3 回答者： chukanshi 回答日時： 2002/03/10 00:23

>やっぱりそういう方法しかないですよね。

>なんか超絶技巧な技があるのかと思いましたが。

あとで、自分で回答を読んで、わかりにくかったかな、と思いましたが、
わかっていただけてよかったです。
というよりも、誰でも真っ先に思いつく方法ですよね。
英語で、「shooting method」とちゃんと名前がついていること自体
始めは驚きでした。

でも、周囲の数値計算の専門の仲間に聞いても、これ以外の方法は
みんな知らないって言うので、私は論文（理論物理です）の
微分方程式の数値計算をこの方法でやって、英文専門誌に投稿して受理
されましたから、認められた方法だと今では思っています。

tgbさんのように、なるほど、両側から解いていって、間で継ぎ合わせる
という方法もありますよね。
収束の速さや、解の精度を挙げるには、いろいろと検討しなければ
ならないかもしれません。

私もいまだに、超絶技巧があるのでは？と調べたり、研究していますが、
新しい方法は、見つかりません。

このページに画期的な方法が現れるのを期待していますし、
ibm_111さんが、新しい方法を手に入れられたら、こちらが
教えていただきたいほどです。

No.4 回答者： motsuan 回答日時： 2002/03/11 17:36

基本的に境界値問題のことをおっしゃっているのだとおもいますが

線形であれば 線型方程式を解く方法もあると思います。
Σa_n(x) （∂ /∂ｘ）^n f(x) = Ｃ
というようなものであれば、結局、各場所の値 f(x_j) (j=1,2,...,N)を未知数として、
方程式から境界条件を除いて適当な数の連立方程式がえられますよね（適切な問題であれば）。
たとえば、端の点の値が与えられているのであれば、１階の微分方程式であれば完全にきまって、
2個の境界条件と微分方程式を差分式になおしたN-2個の連立方程式が得られて
結局、線形連立方程式に帰着できることがわかります。

でも、線形じゃないばあいはやっぱり「いろいろやってみる」方式のようで、
Mathematica などでは（マニュアルによると）
非剛性アダムズ(Adams)法と剛性ギア(Gear)法とを交互に切り替えるのだそうです。
http://www.wolfram.com/products/calculationcente …
ちなみに
線型境界値問題の場合はゲルファンド・ロクチェフスキー(Gel'fand-Lokutsiyevskii)追跡法
だそうです。
確認できなかったのですが、方程式には硬い軟らかいのがあるそうで
http://grape.c.u-tokyo.ac.jp/~makino/kougi/syste …
局所的にそれを切り替えて収束をはやくしているようです。
非線形数値計画法とかの（準）ニュートン法の微分方程式版みたいなものなのでしょうか？

あと、Ｌ（f(x) ）＝０を微分方程式として
発展方程式　(∂/∂t ) f(x) = Ｌ（f(x) ）
＋「境界では境界を満たすようにする」
を考えて、（良く収束するとして）、
発展方程式の時間を進めて
定常解　(∂/∂t ) f(x) ＝０を解とするという考え方もありますよね。