• 締切済み

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

二階の常微分方程式をプログラムをつかって数値的に解きたいのですが 初期条件をいれると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など小さい値をいれてみても途中の計算で発散してしまいます．

みんなの回答

noname#108554
noname#108554
回答No.1

この問題では、y’=zとして、zの解析解を求めることができます。 zの解析解はx=1で発散しますので、 >ちなみにy'(0)=0.001など小さい値をいれてみても途中の計算で発散してしまいます． これは当然ではないかと。なお、 >初期条件をいれるとk1の分母が０になってしまうばあい， http://sci.tea-nifty.com/blog/2008/08/excel_f62a.html のように、原点付近でなんとかしてベキ展開して 計算をスタートさせるのがいい方法だと思います。

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

関連するQ&A

  • ルンゲクッタ法の二階微分方程式（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)ｄｘとなってしまうんですがどうなんでしょう？ おそらくどこかで勘違いしてると思うんです。 長い質問になってしまいましたがどうかご教授のほどよろしくお願いします。

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

    陽的な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 と仮定します． よろしくお願いいたします！！

  • 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
  • いろいろなルンゲクッタがありますが違いはあるのでしょうか？

    ルンゲクッタをインターネットで検索すると様々なルンゲクッタの式が出てきます． 例えば，y'=f(x,y)に対して (1) k1=f(x,y) k2=f(x+dx/2,y+dx/2*k1) k3=f(x+dx/2,y+dx/2*k2) k4=f(x+dx,y+dx*k3) y(n+1)=y(n)+dx/6*(k1+2k2+2k3+k4) (2) k1=f(x,y) k2=f(x+dx/2,y+k1/2) k3=f(x+dx/2,y+k2/2) k4=f(x+dx,y+k3) y(n+1)=y(n)+dx/6*(k1+2k2+2k3+k4) など、以上にあげたもののほかにいくつか出てきます． この(1)と(2)の違いはk1～k4を求める際に k2=f(x+dx/2,y+dx/2*k1) k2=f(x+dx/2,y+k1/2) と,f(x,y)のyの方にdxがかけられている場合とかけられていない場合があります． この違いは何なんでしょうか？詳しい方がいらっしゃいましたらお願いします． ちなみに簡単な微分方程式を上記の両方で解いてみたのですが同じ結果になりました。

  • ルンゲクッタのプロクラム実行できません

    プログラムの実行ができません。 どこを直したらよいのか教えて下さい。 よろしくお願いします。 #include <stdio.h> #include <stdlib.h> #include <math.h> double *dvector(long i,long j); void free_dvector(double *a,long i); double func(double x,double y); double *rk4(double y0,double a,double b, int n,double(*f)()); int main(void) { double *y,h,a=0.0, b=1.0, y0=1.0; int i,n; printf("分割数を入力して下さい------->"); scanf("%d",&n); y=dvector(0,n); y=rk4(y0,a,b,n,func); h=(b-a)/n; for(i=0;i<=n;i++) { printf("x=%f \t y=%f \n",a+i*h,y[i]); } return 0; } double *rk4(double y0,double a, double b,int n,double (*f)()) { double k1,k2,k3,k4,h,*y,x; int i; y=dvector(0,n); h=(b-a)/n; y[0]=y0; x=a; for(i=0;i<n;i++) { k1=f(x,y[i]); k2=f(x+h/2.0,y[i]+h*k1/2.0); k3=f(x+h/2.0,y[i]+h*k2/2.0); k4=f(x+h,y[i]+h*k3); y[i+1]=y[i]+h/6.0*(k1+2.0*k2+2.0*k3+k4); x +=h; } return y; free_dvector(y,0); }

  • ルンゲクッタ法（RK4）で斜方投射の問題を解く

    下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法（RK4）で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k：空気抵抗、m:物体の質量 g：重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2=h/2.0; k1=yd(t,y); k2=yd(t+h2,y+(h2*k1)); k3=yd(t+h2,y+(h2*k2)); k4=yd(t+h,y+(h*k3)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }

  • 数値解析法

    この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; } }

  • ルンゲクッタ法で数値解析（Ｃ言語）

    コンピュータ（プログラミング）のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。 単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。 書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました（Ｃ言語を使って）。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #include <stdio.h> double f1(double t,double x,double v); double f2(double t,double x,double v); int main() { double t,x,v,dt,tmax; double k1[2],k2[2],k3[2],k4[2]; x=1.0; //位置の初期値 v=0.0; //速度の初期値 dt=0.01; //刻み幅 tmax=500.0; //繰り返し最大回数 FILE *output; output=fopen("output.dat","w"); for(t=0;t<tmax;t+=dt) { k1[0]=dt*f1(t,x,v); k1[1]=dt*f2(t,x,v); k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]); k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]); x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0; v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0; fprintf(output,"%f %f %f\n",t,x,v); } fclose(output); return 0; } double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); } ---- このソースは単振動のもので、減衰振動のときは、最後の方の ---- double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); ---- の部分が変わるのだと思うのですが、よく分かりません。 減衰振動は、mx"=-kx-rx'で表されるので、この式を変形（？）したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。

  • mathematicaでリストの格納

    mathematicaでTable関数で作成したリスト {{x1, y1, z1, f(x1,y1,z1)}, {x2, y2, z2, f(x2,y2,z2)}, ... , {xn, yn, zn, f(xn,yn,zn)}} 中のx1～xnまでの各成分とy1～ynまで(、z1～znまで、 f(x1,y1,z1)～f(xn, yn, zn)までの各成分)をそれぞれ配列に格納するにはどうすればいいのでしょうか？(C言語のようにループ文で配列に格納することはできないのでしょうか？) もしくは、行列中で列の成分を取り出すことはできますか？ どなたか解法を示していただければ幸いです。

  • 制約条件の下でラグランジェの未定乗数法を用いる計算について

    条件 g(x,y) = x + y - B = 0 の下で、f(x,y) = xy の最大値を求める問題を解いている途中でわからなくなってしまったので質問させていただきます。 この問題をラグランジェの未定乗数法を用いて解いていくと L = f(x,y) + λg(x,y) = xy + λ(x + y -B)とおいて L1 = y + λ = 0 - (1) L2 = x + λ = 0 - (2) L3 = x + y - B =0 - (3) (1),(2),(3)を展開して 　x = -B/2 　y = -B/2 が求まりした。 上記の値はあくまで候補であり、最大値となっているかどうかを調べる必要があるので、ヘッセの行列式を用いて計算しようとしました。 しかし、それぞれの行と列にどのような値を代入して、どのような計算をしていけばいいのかわからないので質問させていただきました。 よろしくお願いします。

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

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

質問する