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/rksimultane …
http://www.330k.info/essay/Explicit-Runge-Kutta- …
などの”古典的Runge-Kutta”という部分を参考にしました
実際VC++6.0で作成したプログラムも添付したいと思います
Runge-Kutta法の使い方(根本)から間違っているのか?
プログラムのコードが間違っているのか?
ご教授いただけますでしょうか
A 回答 (3件)
- 最新から表示
- 回答順に表示
No.3
- 回答日時:
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
- 回答日時:
念の為確認しておきたいんだけど, 「想定:t=0からt=12あたりまでAは増加し、ピークを迎えた後減少する」というのは正しいんでしょうか? ちゃんと理論的に解析できてますか?
この「想定」が間違っているとしたら, プログラムについて検討しても無意味かもしれないですよ.
No.1
- 回答日時:
なんか下の方に模様が見えるんだけど, ひょっとしてこれが「プログラムである」ということでしょうか?
そもそも初期値が分からなければ「想定外の結果」がどのような結果であるかもわからないのでコメントも難しいが, 一応「解析解」は (数値的に) 求まるはずなのでそれと比較したらどう?
この回答への補足
回答ありがとうございます
以下が作成したプログラムソースです
#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);
字数制限のため再度追加入力します。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 物理学 相対論的運動方程式 1 2022/07/04 06:20
- Excel(エクセル) 1つのファイルを3つのフォルダにファイル名を【明日の日付】にして、コピーをしたい 2 2022/12/21 17:43
- 数学 dx/dt=x-2y +e^t dy/dt=-3x +2y+1 初期値[1,0] [x,y] この連 3 2023/05/15 18:23
- 数学 x=r・cosθの2回微分 θ=ωtとすると? 5 2022/05/10 23:53
- 物理学 物体に一定の大きさfの力をx軸の正の向きに加える。またこの物体には抵抗係数がγの速度に比例する抵抗力 2 2023/07/06 04:01
- 物理学 問a(t)=dv/dt=−0.5v(t) 微分方程式を解いてv(t)を求めよ ∫dv/v=−∫0.5 6 2022/05/23 21:13
- 数学 単振り子とルンゲ・タック法 1 2022/07/15 00:05
- 数学 dx/dt = |y| , dy/dt = x (-∞<t<∞) をとけ 1 2022/09/17 09:56
- 数学 ラプラス変換について 3 2022/10/13 22:18
- 数学 知ってる式から応用させて 微分式を作りたいですが 例えば タンクのなかにVm3という液体がありますし 6 2023/03/09 21:45
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
PICマイコンのコピー(クローン...
-
読み込み中にアクセス違反が発...
-
めんどくさがり屋はプログラマ...
-
VBでコマンドラインから引数を...
-
ニュートン法
-
VC++コンソールアプリでウイン...
-
Notepad++の関数リスト表示でC...
-
Google カレンダーの商用利用
-
ニュートン法を使って解を求め...
-
C言語プログラムについて質問で...
-
C言語でTIFファイルを読み込む...
-
C言語でのaccess violationに...
-
プログラムの行数
-
ポケットコンピュータ(初心者)
-
Excelで4096点以上のFFTの方法
-
人工生命や人工知能の研究にお...
-
C言語で移動平均のプログラムを...
-
Excelに埋め込んだVBAのプログ...
-
寿命
-
メール送信をVB.NET化するには...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
Excelに埋め込んだVBAのプログ...
-
Notepad++の関数リスト表示でC...
-
あるプログラムのコマンドライ...
-
これってほんとにみますか?
-
Excelで4096点以上のFFTの方法
-
「Outlookが他のプログラムによ...
-
自動クエリとはどういうもので...
-
VBAでユーザーフォームが自動的...
-
VBAにてメール作成した際、一部...
-
PICマイコンのコピー(クローン...
-
テキストボックスのエンターキ...
-
読み込み中にアクセス違反が発...
-
特定のwebサイトのタイトルや記...
-
未使用の変数を一括検索する方法
-
モジュール、アプリケーション...
-
COBOLの連絡領域について
-
Google カレンダーの商用利用
-
エクセルとワードをデスクトッ...
-
ドロップダウンリストの文字を...
-
binファイルってiphone専用です...
おすすめ情報