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

このQ&Aに関連する最新のQ&A

A 回答 (4件)

このような場合には、


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)を変化させる1ステップの大きさ

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

読んでみてください。
    • good
    • 0
この回答へのお礼

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

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

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


線形であれば 線型方程式を解く方法もあると思います。
Σa_n(x) (∂ /∂x)^n f(x) = C
というようなものであれば、結局、各場所の値 f(x_j) (j=1,2,...,N)を未知数として、
方程式から境界条件を除いて適当な数の連立方程式がえられますよね(適切な問題であれば)。
たとえば、端の点の値が与えられているのであれば、1階の微分方程式であれば完全にきまって、
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 …
局所的にそれを切り替えて収束をはやくしているようです。
非線形数値計画法とかの(準)ニュートン法の微分方程式版みたいなものなのでしょうか?

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

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


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

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

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

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

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

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

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


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

・y(0)=0となる解yaを1つ求めます。
・y(1)=0となる解ybを1つ求めます。
・y=a*ya+b*ybとして
 境界条件を満たすようにa、bを決定します。
(実質的にはそれぞれ反対側の境界条件を満たすよう係数を
 決定している)

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

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

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

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Qx1=(1,1,1),x2=(1,1,-1),x3=(1,-1,-1)をC^3の基底,{y1,y2,y3}がその双対基底でx=(0,1,0)の時,y1(x),y

[問] ベクトルx1=(1,1,1),x2=(1,1,-1),x3=(1,-1,-1)をC^3の基底とする。
{y1,y2,y3}がその双対基底でx=(0,1,0)の時、
y1(x),y2(x),y3(x)を求めよ。

という問題の解き方をお教え下さい。

双対基底とは
{f;fはF線形空間VからFへの線形写像}
という集合(これをV*と置く)において、
V(dimV=nとする)の一組基底を{v1,v2,…,vn}とすると
fi(vj)=δij(:クロネッカーのデルタ)で定めるV*の部分集合
{f1,f2,…,fn}はV*の基底となる。これを{v1,v2,…,vn}の双対基底と呼ぶ。

まず、
C^3の次元は6(C^3の基底は(1,0,0),(0,1,0),(0,0,1),(i,0,0),(0,i,0),(0,0,i))
だと思うので上記のx1,x2,x3は基底として不足してると思うのです(もう3ベクトル必要?)。

うーん、どのようにしたらいいのでしょうか?

Aベストアンサー

>C^3の次元は6(

これが間違え.
「x1=(1,1,1),x2=(1,1,-1),x3=(1,-1,-1)をC^3の基底」
といってるんだから,係数体はRではなく,C.

あとは定義にしたがって,
dualな基底を書き下せばいいだけ.
y1(x1)=1,y1(x2)=y1(x3)=0であって
v=ax1+bx2+cx2と表わせるわけだし,
v=(v1,v2,v3)とすれば,a,b,cはv1,v2,v3で表現できる
#単なる基底変換の問題.

Qx+y+z=0,2x^2+2y^2-z^2=0のとき,x=yであることを証明せよ。

クリックありがとうございます(∩´∀`)∩

 ★x+y+z=0,2x^2+2y^2-z^2=0のとき,x=yであることを証明せよ。

この問題について説明をお願いします。

Aベストアンサー

おおざっぱな説明になりますが、左の式を
z=-x-y
として、それを右の式のzに代入します。
それを展開してまとめると
x^2-2xy+y^2=0
という式になります。
あとはこれを因数分解すれば
(x-y)^2=0
となるので、x=yという答えがでます。
与えられた条件がほかになければこれでいいはずです。

Q線形です (1)を x+3y-2z=0 x-2y+4z=0 x^2+y^2+z^2=1をもちいて 答

線形です
(1)を
x+3y-2z=0
x-2y+4z=0
x^2+y^2+z^2=1をもちいて
答えが+-の答えになりました
(2)では外せきが8,-6,-5となり
おおきさの5ルート5で割ると
+-の答えにはなりませんでした
どちらが正しいのでしょうか?

Aベストアンサー

外積からでてきた単位べクトルは、外積の定義から、ベクトルa、bに垂直ですよね。
だからそれと正反対のベクトルも、ベクトルa、bに垂直な単位ベクトルだから、これも答えに入れれば
よいのです。つまり外積から出した単位ベクトルの各成分に(-1)をかけた成分のベクトルも答えに
なります。そしてこうして出した2つのベクトルは、先に内積で出した2つのベクトルと一致します。

Qx^2-y^2+x+3y-2=0 ⇔(x+y-1)(x-y+2)=0にする方法

教えてください!!いま二次曲線を学んでるのですけど、x^2-y^2+x+3y-2=0 
⇔(x+y-1)(x-y+2)=0にする方法を
教えてください!!
なぜかというと、私は
(x+1/2)^2-(y-3/2)=4としかできません!

Aベストアンサー

(⇒)
xの次数でそろえます。

x^2-y^2+x+3y-2=0

x^2+x-(y^2-3y+2)=0

()の中を因数分解します。

x^2+x-(y-1)(y-2)=0

全体を因数分解します。

(x+y-1)(x-y+2)=0

※yの次数でそろえてもできます。
※因数分解の仕方は教科書がわかりやすいと思うので、
教科書を参照してください。


(←)
一つずつ掛け合わせて展開していきましょう。

(x+y-1)(x-y+2)=0

x^2-y^2+x+3y-2=0

Qy=0.4431e^(-0.708x) x=0 の時 y=0.4431 y=0.4431/2 こ

y=0.4431e^(-0.708x)

x=0 の時 y=0.4431

y=0.4431/2
ここの式からわかりません!

y/0.4431=e^(-0.708x)
-0.708x=ln(y/0.4431)
y=0.4431/2を代入すると
-0.708x=ln(0.4431/(2×0.4431))
-0.708x=ln(1/2)
-0.708x=-ln2
0.708x=ln2
x=ln2/0.708
=0.693147181/0.708
=0.979=58.7秒

Aベストアンサー

y=0.4431/2 ← 初期値の半分です

y=0.4431e^(-0.708x)
y/0.4431=e^(-0.708x)
これは0.4431を左辺に移しただけ

y=e^x の逆関数 x=ln(y)
を使ってxを求める、ここから高校あたりの数学です。

-0.708x=ln(y/0.4431)

逆関数にした後、 y=0.4431/2を代入すると
-0.708x=ln(0.4431/(2×0.4431))
対数関数の中身の分数を単純に通分して
-0.708x=ln(1/2)

対数関数の公式で
ln(1/b)=ln1-lnb となる
ln1は底がeですが対数は、底>0の時、真数が1の時、0となるので
ln1=0
-0.708x=ln(1/2) ← ln(1/b)=ln1-lnb
-0.708x=ln1-ln2 ← ln1=0
-0.708x=-ln2
0.708x=ln2
x=ln2/0.708 ← ここで関数電卓でも対数表でも使ってln2の値を求める
=0.693147181/0.708
=0.979=58.7秒


人気Q&Aランキング

おすすめ情報