- 締切済み
次の連立微分方程式の解をRunge-Kutta法で
dR/dt=102-(R-V)*0.019 dP/dt=(R-P)*0.019 dL/dt=(P-L)*0.007 dA/dt=(L-A)*0.033 dC/dt=(A-C)*0.004 dV/dt=(C-V)*0.001 この連立微分方程式の解をRunge-Kutta法で導き出し 関数Aの時間的変動をグラフ化したいのですが C言語で作成したプログラムを実行しても想定外の結果となってしまいます (想定:t=0からt=12あたりまでAは増加し、ピークを迎えた後減少する) というのも、とある論文を基に数値シミュレーションを試みている状況なのです Runge-Kuttaの理論に関してはWikipediaやPukiwikiを参照しました プログラムの実装で参考にしたWebは http://www.geocities.jp/supermisosan/rksimultaneousequation.html http://www.330k.info/essay/Explicit-Runge-Kutta-Butcher-Tableau などの”古典的Runge-Kutta”という部分を参考にしました 実際VC++6.0で作成したプログラムも添付したいと思います Runge-Kutta法の使い方(根本)から間違っているのか? プログラムのコードが間違っているのか? ご教授いただけますでしょうか
- Nate_River
- お礼率0% (0/3)
- C・C++・C#
- 回答数3
- ありがとう数0
- みんなの回答 (3)
- 専門家の回答
みんなの回答
- akayoroshi
- ベストアンサー率50% (46/91)
Mathematicaに解かせてみました。 以下のMathematicaのスクリプトで、a(t)のグラフを描かせました。 (* By Akayoroshi on 2011-02-25 *) tmax=25 solve=NDSolve[{ r'[t]==102-(r[t]-v[t])*0.019, r[0]==0, p'[t]==(r[t]-p[t])*0.019, p[0]==0, l'[t]==(p[t]-l[t])*0.007, l[0]==0, a'[t]==(l[t]-a[t])*0.033, a[0]==0, c'[t]==(a[t]-c[t])*0.004, c[0]==0, v'[t]==(c[t]-v[t])*0.001, v[0]==0 },{r,p,l,a,c,v}, {t,0,tmax}] Plot[Evaluate[{r[t],p[t],l[t],a[t],c[t],v[l]} /. solve ], {t,0,tmax}] Plot[Evaluate[a[t] /. solve ], {t,0,tmax}] t=12のあたりにはピークは出ませんでした。 質問文中のC++での処理結果は読み取れませんが、この結果とは違っていますか。
- Tacosan
- ベストアンサー率23% (3656/15482)
念の為確認しておきたいんだけど, 「想定:t=0からt=12あたりまでAは増加し、ピークを迎えた後減少する」というのは正しいんでしょうか? ちゃんと理論的に解析できてますか? この「想定」が間違っているとしたら, プログラムについて検討しても無意味かもしれないですよ.
- Tacosan
- ベストアンサー率23% (3656/15482)
なんか下の方に模様が見えるんだけど, ひょっとしてこれが「プログラムである」ということでしょうか? そもそも初期値が分からなければ「想定外の結果」がどのような結果であるかもわからないのでコメントも難しいが, 一応「解析解」は (数値的に) 求まるはずなのでそれと比較したらどう?
関連するQ&A
- 微分方程式の問題お願いします
微分方程式の問題お願いします インダクタンスL、容量C、抵抗Rと電源を直列に接続する この電気回路を流れる電流I(t)は次の2階微分方程式を満たすとき以下の問いに答えよ L*d^2I/dt^2+R*dI/dt+I/C=dV/dt・・・(1) (1)V(t)=Eのとき(1)の一般解を求めよ (2)V(t)=sint、 L=R=C=1とする (1)の特解を,I(t)=acost+bsintとしたとき、aとbを求めよ この微分方程式を解け (3)dy/dx=-(x^3+4x^3y^3)/(y^2+3x^4y^2) (3)しかできなくて x^4/4+x^4y^3+y^3/3=Cとなりました 残りの問題お願いします
- ベストアンサー
- 数学・算数
- 微分方程式(ルンゲ=クッタ)がうまくいきません
C言語で物体の自由落下運動(空気抵抗を考慮する)を微分方程式を用いて解かなければならないのですが、どうしてもグラフの形が曲線になりません。どこが間違っているのか指摘お願いします。 空気抵抗を考慮しているので、落下速度はどんどんある一定の値に近づいていくはずなのですが、グラフは直線になってしまいます。いろいろと方法を試しましたがどれも外れました。 写真は赤が今の状態(間違っています)、黒が理想形です。 ここではルンゲ=クッタ法を使用しています。 #include <stdio.h> #include <math.h> double yd(double t,double y){ double m = 50.0; 質量 double g = 9.8; 重力加速度 double k = 0.01; 空気抵抗の定数係数 double a= g - ((k *y)/m); 運動方程式 return(a); } 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)); a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4)); return(a); } void main(){ double y = 0.0; double h = 0.001; double t = 0.0; int i; for(i = 0; i < 1000; i++){ t += h; y = runge(t,y,h); printf("%f %f\n",t,y); } }
- ベストアンサー
- C・C++・C#
- 近似計算するプログラム
近似計算法である4次のRunge-Kutta法を用い、プロット点を一斉にdos画面に表示したい場合、全体的にどのようなプログラムにすれば良いのでしょうか?最終的にはその多数のプロット点をエクセルに出力し、グラフを作りたいと思っています。 何を近似するのかは気にしていないので、色々探したところ、あるサイトから4次のRunge-Kutta法の一部分のプログラムを見つけ、以下のを使用することにしたのです。プログラムは全くの無知なので、どのくらいの量になるのか想像出来ませんが、もし短い容易なプログラムになるのであれば載せて頂けると有難いです。コンパイラはVBC++です。 反復 { v1dt = v(x[j] ) * dt v2dt = v(x[j] + v1dt/2 ) * dt v3dt = v(x[j] + v2dt/2 ) * dt v4dt = v(x[j] + v3dt ) * dt x[j+1] = x[j] + (v1dt + 2*v2dt + 2*v3dt + v4dt)/6 }
- ベストアンサー
- C・C++・C#
- 三元連立微分方程式の解き方教えてください。
dx/dt=(A-C)/A*y*z dy/dt=((C-A)x-L)/A*z dz/dt=L/C*y この連立方程式なのですが,x,y,z以外は全て定数です。 まずは,ひとつの変数のみの方程式を作ろうと思いzだけの式を作りました。 z*z'''-z''*z'+((A-C)/A)^2*z^3*z'=0 という式になりました('は微分回数です。) そして更に左辺は(z*z'')のtで微分した形に近いんじゃないかと思い解いていったのですが・・・・ここらへんからもうさっぱりでどうすればよいのかわからなくなってしまいました。もしよろしければこういう三元の連立で変数がyzみたいに積の項がある方程式のよい解き方を教えてください。パソコンで計算してみたらx,yは円の運動をしました。
- ベストアンサー
- 数学・算数
- 方程式の計算を教えることになったのですが...
y'=y+x^2 , y(0)=1 のときのy(3)の値をEuler法およびRunge-Kutta法により求める問題を教えることになったのですが、分かりやすくプログラムを教えるためにはどういう風に教えればいいのでしょうか? 皆さんの考えを教えてください。
- ベストアンサー
- その他(プログラミング・開発)
- 連立微分方程式
点P(x,y)は連立微分方程式 dx/dt=y dy/dt=-x を満たすものとする。t=0で原点以外の点から出発した点P(x,y)は、tが増加するにつれてどのようにふるまうか述べよ。図を用いてもよい。 この問題の解き方がよく分かりません。 連立微分方程式について、色々な文献を見てみたのですが、どうもいまいちです。 上の連立方程式を2つともdt=のかたちにして、dx/y=dy/-xという式にし、変数を分離して両辺を積分して・・・すると、x^2+y^2=Cという式に なりました。 円の方程式っぽいです。 でも、tは消えてしまい・・・ よく分からなくなってきました。 そもそもここまでの解き方も自分は間違っているのでしょうか?? ご意見やヒント、解答ヨロシクお願いしますm(_ _)m
- 締切済み
- 数学・算数
補足
回答ありがとうございます 以下が作成したプログラムソースです #include <stdio.h> #include <math.h> double f1(double t,double r,double v,double p,double l,double a,double c); double f2(double t,double r,double v,double p,double l,double a,double c); double f3(double t,double r,double v,double p,double l,double a,double c); double f4(double t,double r,double v,double p,double l,double a,double c); double f5(double t,double r,double v,double p,double l,double a,double c); double f6(double t,double r,double v,double p,double l,double a,double c); int main(void) { double t,r,v,p,l,a,c,dt,tmax; double k0[6],k1[6],k2[6],k3[6]; FILE *output; dt=0.01; tmax=30.0; r=v=p=l=a=c=0.0; //初期値 output=fopen("output.txt","w"); for(t=0.0;t<tmax;t+=dt) { k0[0]=dt*f1(t,r,v,p,l,a,c); k0[1]=dt*f2(t,r,v,p,l,a,c); k0[2]=dt*f3(t,r,v,p,l,a,c); k0[3]=dt*f4(t,r,v,p,l,a,c); k0[4]=dt*f5(t,r,v,p,l,a,c); k0[5]=dt*f6(t,r,v,p,l,a,c); k1[0]=dt*f1(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); k1[1]=dt*f2(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); k1[2]=dt*f3(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); k1[3]=dt*f4(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); k1[4]=dt*f5(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); k1[5]=dt*f6(t+dt/2.0,r+k0[0]/2.0,v+k0[1]/2.0,p+k0[2]/2.0,l+k0[3]/2.0,a+k0[4]/2.0,c+k0[5]/2.0); 字数制限のため再度追加入力します。