• 締切済み

常微分方程式の数値解法： 陰解法のルンゲクッタ法の公式について

陽的なRunge-Kutta法は, y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4) ただし，k1=hf(xi,yi), k2=hf(xi+h/2, yi+k1/2) k3=hf(xi+h/2, yi+k2/2) k4=hf(xi+h,yi+k3) と表すことは x=[x0,xn]の範囲でルンゲ・クッタ法により数値的にy(x)を解きたいのですが， 解こうとしている問題の初期条件がx=xnの時，y=0となっており，陰的(xnからx0に向かって)に解かなければならないのだろう，と考えています． 上記の公式でhの代わりに-hを入れて，プログラムを走らせても求めたい結果に大きな差異が生じてしまい困っています． そこで，陰的なRunge-kutta法の公式には，陽的な解き方と比較してどのような修正をすればよいか，教えてください． ちなみに解きたい微分方程式は， d2y/dx2 = x と仮定します． よろしくお願いいたします！！

みんなの回答

回答No.1

u = -x として、 u0,u1,...,un = -xn,...,-x0 とすれば、陽的解法で問題ないと思いますが。

aminofish
質問者

お礼

そうなんですか！ 陰的に解く必要はないのですね． もう一度教えていただいた通りに，計算し直してみます． 素早い回答くださいましてありがとうございます．

  1. カテゴリ
  2. 学問・教育
  3. 数学・算数

関連するQ&A

  • 常微分方程式の解き方について

    常微分方程式の解き方について質問です。 Fickの第2法則の式を4次のルンゲ・クッタ法で解きたいと考えています。 4次のルンゲ・クッタ法は、 k1=f(xi,yi) k2=f(xi+h/2,yi+h/2*k1) k3=f(xi+h/2,yi+h/2*k2) k4=f(xi+h,yi+h*k3) として、 yi+1=yi+h/6*(k1+2k2+2*k3+k4)で解いています。 そこで、Fickの法則（dC/dt=D*d^2C/dx^2）を、以下の式ように変換しました。 dC/dt=D*（Cj+1+Cj）/dx^2　←これ(dC/dt)をf(xi,yi)としています。 f(xi,yi)のyiはCとして、hはdtとして解いてみました（xiはないものとして解いています） エクセルでは、 　 t0 ， t0+⊿t ， ・・・ x1: Cj,n ， Cj,n+1 ， ・・・ x2: Cj+1,n ，Cj+1,n+1， ・　　 ・ のようなセルを組み、t0+⊿tにCn+1=Cn+h/6*(k1+2k2+2*k3+k4)の式を入れています。 上記の方法で解いてみましたが、うまくいきませんでした。 説明が分かりにくく申し訳ないですが、上記のやり方で間違っているところ、アドバスなど ありましたら回答お願いします。

  • ルンゲクッタ法の二階微分方程式（Fortran)

    数値計算の演習問題で以下の二階微分方程式をルンゲ・クッタで解けという問題があります。 -y"+x^2・y=e・ｙ（eは定数、”・”は単なる掛け算） y(0)=1, y'(0)=0, 0<=x<=2までを計算せよ。 これは ｙ’=z・・・(1) 　 z'=(x^2-e)y・・・(2) この２つの連立方程式を解けばよいところまではわかります。 まず(2)を解くときにルンゲ・クッタの場合 (k1+２k2+２k3+k4)/6の項（←公式の右辺第二項）のｋ（1～4）を求めなければいけません。 質問はｋの求め方です。 本にはy'=f(x,y,z) , z'=g(x,y,z)とおけば (2)の場合だと例えばk1は k1=g(xn,yn,zn)ｄｘで計算する。と書いてあります しかしz'=(x^2-e)y（←ｚが入ってない）　なので、計算すると k1=g(xn,yn)ｄｘとなってしまうんですがどうなんでしょう？ おそらくどこかで勘違いしてると思うんです。 長い質問になってしまいましたがどうかご教授のほどよろしくお願いします。

  • 数値解析法

    このHeun法のプログラムをRunge-Kutta法にするにはどうしたらいいですか？ #include <stdio.h> #include <math.h> double f1(double y) { return y; } double f2(double y) { return -4*y; } int main(){ double a=0; double b; int m=10; int n; double h; double x,y; int k; double e; double f; double k1,k2; printf("Heun法計算例：y=e^x, y=1/e^4x\n\n"); // y = e^x b = 1; for(n=100;n<=10000;n*=100){ h = (b-a)/n; printf("y' = y: h(=dx) = %.1e (y=e^x)\n",h); x = a; y = 1; for(k=0;k<=n;k++) { x = k*h; if(k%(n/m)==0) { f = exp(x); e = fabs(y-f); printf("x=%.2f, y=%f, e^x=%f er=%.0e\n",x,y,f,e); } // Heun's method k1 = h*f1(y); k2 = h*f1(y+k1); y += (k1+k2)/2; } }

  • 常微分方程式の問題で。

    常微分方程式の問題で オイラー法、ルンゲ・クッタ法、ルンゲ・クッタ・ギル法、ミルン法、アダムス法より２つの方法を用いて解け。 y''+y'/x+y=0 初期条件y'/x=-1/2 0<=x<=4 h=0.4 という問題なのですが、基本的に解き方がわかりません。 二階微分なので連立ルンゲ・クッタで解くのかな？っと思ったのですが 結局それからどうすればいいのかわかりません。 ２つの方法とありますが、これはどの方法でも解けるものなのでしょうか？ 本当はＣ言語のプログラムを組まなければいけないのですが、 この理論がわからないのでおねがいします。

  • 偏微分方程式の数値解析

    お忙しいところすみません． 現在，FDTD法による数値解析を行っています． FDTD法は空間及び時間微分に対して中心差分を行っていますがこれを4次のルンゲクッタ法で行うことは出来るのでしょうか？ある参考書に以下のような「連立1階微分方程式」を4次のルンゲ-クッタ法で解くというのがありました． dy/dt=f(t,y,z) dz/dt=g(t,y,z) これを4次のルンゲ-クッタ法で解くと y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4) z(n+1)=z(n)+1/6*(j1+2*j2+2*j3+j4) ここで k1=Δt*f(t(n),y(n),z(n)) j1=Δt*g(t(n),y(n),z(n)) k2=Δt*f(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) j2=Δt*g(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) k3=Δt*f(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) j3=Δt*g(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) k4=Δt*f(t(n)+Δt,y(n)+k3,z(n)+j3) j4=Δt*g(t(n)+Δt,y(n)+k3,z(n)+j3) これをマクスウェルの方程式に対応させ ε∂E(t,x,y,z)/∂t=f(t,H)=∇×H μ∂H/∂t=g(t,E)=-∇×E として出来ないこともない気がするのですが． なぜこの方法を試みようとした理由は現在，通常のFDTD定式化により数値解析を行っていますが ・プログラムは正しい ・始めの部分は正常に計算出来ている ・途中から凹凸が出始め，それが成長して発散している という事態になり，数値的不安定の可能性が強いと思ったからです．オイラー法やルンゲ-クッタ法ではこのような異常振動は進み幅を小さくすれば止まる(ということが数学的に証明されている)ようです． 何卒回答の方，よろしくお願い致します．

  • ４次のルンゲクッタ公式

    　４次のルンゲクッタ法の公式 　　　yi+1＝ｙi＋1/6(k1+２k2+２k3+k4) 　 の導出の仕方がわかりません。 本やウェブで探してみましたが、省略されてしまっています。 　もし、導出方法がわかる方、詳しく載っているウェブページなどをご存知の方がいらっしゃいいましたら、教えていただけないでしょうか。 　よろしくお願いいたします。

  • ルンゲクッタ法の解き方（初期条件）

    二階の常微分方程式をプログラムをつかって数値的に解きたいのですが 初期条件をいれるとk1の分母が０になってしまうばあい，どのように計算すれば良いかわかりません．どのように通常解くのかおしえていただけませんか？ 具体例として y''=((y')^2+1)/sin(arctan(y')) 初期条件 y(0)=C1(定数） y'(0)=0 k1=0.5*h*f(xn,yn,yn') k2=0.5*h*f(xn+0.5h,yn+K,yn'+k1) (ここでK=0.5h(yn'+0.5k1)) k3=0.5*h*f(xn+0.5h,yn+K,yn'+k2) k4=0.5*h*f(xn+h,yn+L,yn'+2k3)　　(ここでL=h(yn'+ｋ３) ｋ１の値を計算すると分母が０になるので計算できません． どうぞ宜しくお願いします. ちなみにy'(0)=0.001など小さい値をいれてみても途中の計算で発散してしまいます．

  • オイラー法、ルンゲクッタ法について。

    オイラー法、ルンゲクッタ法について。 この２つについて分からない事があるので質問します。 まず、オイラーについてですが、yi+1=yi+hf(x,y)という式がテイラー展開によって求まると言われましたが、テイラー展開の2次以降の項は微少量として無視できるのは分かります。でもそもそもテイラー展開ってひとつ先の値を今の値から求まるみたいな展開でしたっけ？？というのが一つ目の質問です。 ２つ目は、オイラーの式の中のf(x,y)についてです。簡単なバネ・マス・ダンパ系を考えた時、運動方程式はm・d2x/dt2+c・dx/dt+kx=0となると思いますが、この場合のf(x,y)はどうやって求めるのでしょうか。 ３つ目はルンゲクッタそもそもについてです。 ルンゲクッタとはK1K2K3K4という係数(？)に1221という重みをかけるとyi+1が求まるそうですが、この理由がどんなサイトや本を見ても納得出来ません。 何か分かりやすい本やサイトがあれば教えて頂けないでしょうか。 以上３つの質問、回答よろしくお願いします。

  • Rubyのルンゲクッタ法がうまくいきません

    Rubyでルンゲクッタ法でy=exp(x**2)と比較しようという問題で、下のようなプログラムを組むとyの値がかなり大きくなってしまいます。どこが間違っているのでしょうか。教えてください。 #! ruby -Ks #ルンゲックッタ def df(a,b) z=2*a*b return z end n=10 x=Array.new(n) y=Array.new(n) x[0]=0.0 y[0]=1.0 h=1.0/n #任意で変更 fw=File.open("output2.txt","w") for i in 0..n k1=df(x[i],y[i])*h k2=df(x[i]+h/2,y[i]+k1/2)*h k3=df(x[i]+h/2,y[i]+k2/2)*h k4=df(x[i]+h,y[i]+k3) k=(k1+2*k2+2*k3+k4)/6 x[i+1]=x[i]+h y[i+1]=y[i]+k fw.print x[i]," ",y[i]," ","\n" end fw.close()

    • ベストアンサー
    • Ruby
  • ルンゲクッタ法について

    初期値問題 d^2y/dx^2=dy/dx-(dy/dx)^2 (0<x<5) , y(0)=1+log(3/2) , y'(0)=1/3 を１階連立系に直し、コンピュータを用いて古典的ルンゲ・クッタ公式によって数値解を求め、真の解と比較せよ。但し、ステップ幅hは1/4 , 1/8 ,...のように二通り以上とることとし、数値解の誤差の振る舞いを吟味せよ。ルンゲ・クッタ公式のプログラムも添付し出来れば解の様子をx‐y平面に図示すること。 という問題があるのですが、全然手がつきません。誰か数学の偉い方、ご指導お願いします。お礼もいたします。

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する