次の連立微分方程式の解を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法の使い方（根本）から間違っているのか？

プログラムのコードが間違っているのか？

ご教授いただけますでしょうか

回答No.3

akayoroshi

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++での処理結果は読み取れませんが、この結果とは違っていますか。

回答No.2

Tacosan

念の為確認しておきたいんだけど, 「想定：t=0からt=12あたりまでAは増加し、ピークを迎えた後減少する」というのは正しいんでしょうか? ちゃんと理論的に解析できてますか?
この「想定」が間違っているとしたら, プログラムについて検討しても無意味かもしれないですよ.

回答No.1

Tacosan

なんか下の方に模様が見えるんだけど, ひょっとしてこれが「プログラムである」ということでしょうか?
そもそも初期値が分からなければ「想定外の結果」がどのような結果であるかもわからないのでコメントも難しいが, 一応「解析解」は (数値的に) 求まるはずなのでそれと比較したらどう?

補足

回答ありがとうございます
以下が作成したプログラムソースです

#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);

字数制限のため再度追加入力します。