この方程式、stiffなんでしょうか・・・

質問

僕は以下のような連立微分方程式を数値計算で解こうと試みています。
現在微分方程式の解を級数展開と変数変換とルンゲクッタ法で求めようとしています。
しかし、どのように初期値をいじっても必ず解が発散してしまいます。
色々と自分で調べてみる過程で微分方程式を数値計算で解くときにstiff問題というものが生じることがあると知りました。
解の挙動が急に変わる場合や境界条件に大きく依存する場合などで途中からルンゲクッタ法などで数値計算できなくなり、解が発散してしまうという条件がぴったり自分の場合と一致しました。
以下の連立微分方程式はstiffなのでしょうか？

(x-1)^4*f ''(x)+(2x-1)(x-1)^3/x*f '(x)-{ ((x-1)/x+A(x))^2+f^2-1 }f=0,
(x-1)^4*A''(x)+(2x-1)(x-1)^3/x*A'(x)-(x-1)^2/x^2A(x)-f^2(A+(x-1)/x)=0

stiffだとルンゲクッタ法は使えなくなってしまうんですよね・・・・？

回答No.4

FT56F001 · ベストアンサー

>　x=0周りでf(x)=f0+f1x+f2*x^2+f3*x^3+f4*x^4などと級数解を仮定しておき、
x=0近傍で作った漸近解が，x→1まで解けるとは思えないので，
解が途中でツギハギになるのは仕方ないと思います。
また，漸近解として問題になる/調べなければならないのは，
x→0で発散する解，すなわち
f(x)=f0+f_(-1)/x+f_(-2)/x^2+f_(-3)/x^3+・・・の方でしょう。

元の微分方程式，
f''(r)+1/r*f'(r)-{ (1/r-A(r))^2+f^2-1 }f=0,
A''(r)+1/r*A'(r)-1/r^2A(r)-f^2(A-1/r)=0
は，ベッセルの微分方程式
u''(r)+(1/r)*u'(r)+(1-ν^2/r^2)*u(r)=0

に似ているので，ベッセル関数の性質から何か分かるかもしれませんね。
第１種J(r)はr=0で有限，第２種N(r)はr→0で発散の性質は似ていそうです。

とにかく，数値計算で特異点近くまで解こうとしてもムリでしょう。
まずは特異点周りの漸近的性質をじっくり調べる，
次に特異点を除去する方法があれば考える，といった段取りでしょうか。

お礼

回答ありがとうございました。
x=0まわりの級数解を初期条件にしてルンゲクッタで解いていく方法でやっていたんですが、
x=1までその級数解を使うわけではなくあくまで初期条件としてなのでちゃんと解が出ると思っていたんですがそうではないんですね？

分数式を考察するのもしてみたんですけど項数が足りなかったのかもしれませんね。

ベッセル関数に似たような挙動を確かに示しています。
同じではないですが、初期条件次第ではかなり似た形になっています。
以下のような定常解にしたいのですがベッセル関数の特性だとこのようにはならなさそうですね・・・
http://www.google.co.jp/imgres?q=nielsen+olesen+vortex&hl=ja&sa=X&rlz=1T4ADSA_jaJP386JP387&biw=1147&bih=614&tbm=isch&prmd=imvns&tbnid=FTIp5T87sN2yMM:&imgrefurl=http://iopscience.iop.org/1475-7516/2005/02/003/fulltext/&docid=-Jt68YqLwMTJJM&imgurl=http://ej.iop.org/images/1475-7516/2005/02/003/Full/9383603.jpg&w=547&h=424&ei=tG5DT8u0Mu-NmQWZmYXaBw&zoom=1&iact=hc&vpx=849&vpy=146&dur=6145&hovh=198&hovw=255&tx=127&ty=146&sig=101522157301032231749&page=1&tbnh=136&tbnw=175&start=0&ndsp=18&ved=0CFQQrQMwBQ

回答No.3

FT56F001

あと，手計算でx=0近傍およびx=1近傍の漸近的な解を求める，という方法があります。

手始めにx=0の近傍と考えて，x以外の因子ではx=0と近似して定数に置くと，与微分方程式は
f ''(x)+1/x*f'(x)-{ (-1/x+A(x))^2+f^2-1 }f=0,
A''(x)+1/x*A'(x)-1/x^2A(x)-f^2(A-1/x)=0
と簡単になります。

x→0のとき，1/xのオーダで発散すると仮定して，A=A0/x，f=F0/xと置き，整理してみます。
f'=-F0/x^2，f''=2*F0/x^3，A'=-A0/x^2，A''=2*A0/x^3を代入すると，

F0-{(-1+A0)^2+F0^2}F0=0
F0^2*(A0-1)=0

F0=0,A0=1，すなわち
f(x)=0，A(x)=1/xがx→0における漸近解のひとつであることがわかります。
(無限大に発散する解がある以上，数値計算でx→0へ向かって解くのは，相当注意が必要です)

このようにして，特異点まわりの漸近解を手計算で求めて，数値解と継ぎ合わせていけば，
全貌が見えてくるかもしれません。

お礼

回答ありがとうございます。
実は、最初はx=0周りでf(x)=f0+f1x+f2*x^2+f3*x^3+f4*x^4などと級数解を仮定しておき、
この仮定した式を元の方程式に代入してmathematicaで級数解の係数を求め、
それを用いて、x=0まわりのf(x),f'(x)を求めてx=0.0001などを代入して初期条件にして解いています。
これも一種の漸近解となっているはずなのですが、それでもやはりx=1近傍で発散してしまいます・・・

回答No.2

FT56F001

>実はこの方程式の定義域は「0:1」なんです。
元の問題によりますが，微分方程式がおとなしい点，例えばx=1/2で初期値を仮定して，
プラス方向に1-0まで，マイナス方向に0+0まで解いて行く。
だましだまし解くため，0付近，1付近では，刻み幅をとても細かく選び直す，というテクはあります。

ただし，x=0，x=1の付近は，数学的に難しくなっていて，当然，数値的にも解けなくなる，
という感触があります。定式化や変数変換のやり方から見直すべきかもしれません。

お礼

回答ありがとうございます。
実は変数変換前の式は以下のような感じです。
この場合rは[0,∞]の定義域です。
変数変換はx=r/(r+1)で行いました。
f ''(r)+1/r*f'(r)-{ (1/r-A(r))^2+f^2-1 }f=0,
A''(r)+1/r*A'(r)-1/r^2A(r)-f^2(A-1/r)=0

r=tanxとかで変数変換してみたりもしたんですが、やはりx=1近傍で発散してしまいます・・・

回答No.1

FT56F001

あいまいな記憶+パッと見の感じ=見当はずれ
かもしれませんので，参考程度に

stiffとは，
「微分方程式の解が，数学的には存在して一意だけれど，数値的に解こうとすると条件が厳しくて難しい」といった意味合いだったと思います。

質問者さんの微分方程式，x=1やx=0で分母が0になり，f''やA''が求められなくなっていませんか。これ，微分方程式の特異点になってないかな?

数学的に解けない(解が存在しない，一意性が崩れる)とすれば，
「数値計算が難しいから，ルンゲクッタ法じゃダメで高度なテクの数値解法が必要」
といった小手先の話ではなくて，
「解がないのだから求めようとしても無理」という原理的な話になっているかも???

お礼

回答ありがとうございました。
実はこの方程式の定義域は「0:1」なんです。
元々の方程式は[0,∞]なんですけど、計算しにくいので自分で変数変換しました。
なのでx=0,1をぎりぎり通らないように数値計算してやればいいのかなって思ったのですが・・・
stiffな微分方程式に対してルンゲクッタを用いたときの挙動とそっくりなんです。
一応、数値計算で求められたグラフも存在するんですが、
そのグラフだと急激に立ちあがりを起こして、そのあとf,Ａの傾きが0になって定常状態になります。
そういう解はstiffになると調べたらありましたのでそうなのかなと・・・
ルンゲクッタだとx=1に近づくにつれて急激な発散を起こしますし・・・