a.outの結果がおかしい件について

質問

以下のSEIRモデルの解を求めるプログラムをコンパイルしa.outすると、
結果が
0.000000,200.000000,40.000000,30.000000,30.000000,20.000000,8.000000,12.000000,10.000000,350.000000
40.100000,275.397813,152.832528,37.778276,45.918538,0.000000,0.000000,3.615719,1.024043,516.566916
しか表示されず、またseirすべての合計は350で一定であるはずなのに516となっています。
どこがおかしいのでしょうか？
以下にプログラムを示します。
#include <stdio.h>
#include <math.h>

double alf=0.33, bet=0.0875, gam=0.2; ///alf:移行率,bet:感染率,gam:除外率
double f1(double t1,double a0,double a1,double a2,double a3){return -bet*a0*a2;}
double f2(double t1,double a0,double a1,double a2,double a3){return bet*a0*a2-alf*a1;}
double f3(double t1,double a0,double a1,double a2,double a3){return alf*a1-gam*a2;}
double f4(double t1,double a0,double a1,double a2,double a3){return gam*a2;} //箱Aの関数
double g1(double t1,double b0,double b1,double b2,double b3){return -bet*b0*b2;}
double g2(double t1,double b0,double b1,double b2,double b3){return bet*b0*b2-alf*b1;}
double g3(double t1,double b0,double b1,double b2,double b3){return alf*b1-gam*b2;}
double g4(double t1,double b0,double b1,double b2,double b3){return gam*b2;} //箱Bの関数

int main(void) {
double t1,dt,t1max;
double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[4],a[4],b[4],da[4],db[4] ;
double (*pf[4])(double,double,double,double,double)={f1,f2,f3,f4};
double (*pg[4])(double,double,double,double,double)={g1,g2,g3,g4};
int i; //宣言

t1 = 0.0;
dt = 0.1;
t1max = 40.0; //時間初期値

a[0] = 200.0;
a[1] = 40.0;
a[2] = 30.0;
a[3] = 30.0; //箱A初期値（a[0]:感受性人口、a[1]:潜伏人口、a[2]:感染人口、a[3]:隔離人口）

b[0] = 20.0;
b[1] = 8.0;
b[2] = 12.0;
b[3] = 10.0; //箱B初期値（b[0]:感受性人口、b[1]:潜伏人口、b[2]:感染人口、b[3]:隔離人口）

printf("%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n",t1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);

for(t1=0.0;t1<=t1max;t1+=dt) {

for (i=0;i<4;i++) k1[i]=dt*pf[i](t1,a[0],a[1],a[2],a[3]);
for (i=0;i<4;i++) k2[i]=dt*pf[i](t1+dt/2.0,a[0]+k1[0]/2.0,a[1]+k1[1]/2.0,a[2]+k1[2]/2.0,a[3]+k1[3]/2.0);
for (i=0;i<4;i++) k3[i]=dt*pf[i](t1+dt/2.0,a[0]+k2[0]/2.0,a[1]+k2[1]/2.0,a[2]+k2[2]/2.0,a[3]+k2[3]/2.0);
for (i=0;i<4;i++) k4[i]=dt*pf[i](t1+dt,a[0]+k3[0],a[1]+k3[1],a[2]+k3[2],a[3]+k3[3]);
for (i=0;i<4;i++) l1[i]=dt*pg[i](t1,b[0],b[1],b[2],b[3]);
for (i=0;i<4;i++) l2[i]=dt*pg[i](t1+dt/2.0,b[0]+l1[0]/2.0,b[1]+l1[1]/2.0,b[2]+l1[2]/2.0,b[3]+l1[3]/2.0);
for (i=0;i<4;i++) l3[i]=dt*pg[i](t1+dt/2.0,b[0]+l2[0]/2.0,b[1]+l2[1]/2.0,b[2]+l2[2]/2.0,b[3]+l2[3]/2.0);
for (i=0;i<4;i++) l4[i]=dt*pg[i](t1+dt,b[0]+l3[0],b[1]+l3[1],b[2]+l3[2],b[3]+l3[3]);}
for (i=0;i<=3;i++) a[i]=a[i]+((k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6.0);
for (i=0;i<=3;i++) b[i]=b[i]+((l1[i]+2.0*l2[i]+2.0*l3[i]+l4[i])/6.0); //4次のルンゲクッタ

da[0]=-0.8*a[0];
da[1]=-0.8*a[1];
da[2]=-0.2*a[2];
da[3]=-0.5*a[3]; //AからBへの移動量

db[0]=-b[0];
db[1]=-b[1];
db[2]=-0.7*b[2];
db[3]=-0.9*b[3]; //BからAへの移動量

if (t1 < 21.0) { //移動開始
for (i=0;i<=3;i++) {
a[i]=a[i]+da[i];
b[i]=b[i]-db[i]; ////AからBへ移動後のAにおけるSEIRの数,AからBへ移動後のBにおけるSEIRの数
} }
else {
for (i=0;i<=3;i++) {
a[i]=a[i]-da[i]; ////BからAへの移動後のAにおけるSEIRの数
b[i]=b[i]+db[i]; ////BからAへの移動後のBにおけるSEIRの数
}
printf("%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n",t1+0.1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);
}
return 0;
}

回答No.3

f272 · ベストアンサー

やりたいことが，よくわからないが...
AからBへ移動とかBからAへの移動というのは，一体いつ移動するの？21日までは各タイムステップでAからBへ移動して，その後は各タイムステップでBからAへ移動するのだろうか？
もしそうなら
#include <stdio.h>
#include <math.h>

double alf=0.33, bet=0.0875, gam=0.2; ///alf:移行率,bet:感染率,gam:除外率
double f1(double t1,double a0,double a1,double a2,double a3){return -bet*a0*a2;}
double f2(double t1,double a0,double a1,double a2,double a3){return bet*a0*a2-alf*a1;}
double f3(double t1,double a0,double a1,double a2,double a3){return alf*a1-gam*a2;}
double f4(double t1,double a0,double a1,double a2,double a3){return gam*a2;} //箱Aの関数
double g1(double t1,double b0,double b1,double b2,double b3){return -bet*b0*b2;}
double g2(double t1,double b0,double b1,double b2,double b3){return bet*b0*b2-alf*b1;}
double g3(double t1,double b0,double b1,double b2,double b3){return alf*b1-gam*b2;}
double g4(double t1,double b0,double b1,double b2,double b3){return gam*b2;} //箱Bの関数

int main(void) {
double t1,dt,t1max;
double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[4],a[4],b[4],da[4],db[4] ;
double (*pf[4])(double,double,double,double,double)={f1,f2,f3,f4};
double (*pg[4])(double,double,double,double,double)={g1,g2,g3,g4};
int i; //宣言

t1 = 0.0;
dt = 0.1;
t1max = 40.0; //時間初期値

a[0] = 200.0;
a[1] = 40.0;
a[2] = 30.0;
a[3] = 30.0; //箱A初期値（a[0]:感受性人口、a[1]:潜伏人口、a[2]:感染人口、a[3]:隔離人口）

b[0] = 20.0;
b[1] = 8.0;
b[2] = 12.0;
b[3] = 10.0; //箱B初期値（b[0]:感受性人口、b[1]:潜伏人口、b[2]:感染人口、b[3]:隔離人口）

printf("%4.2f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n",t1+0.1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);

for(t1=0.0;t1<=t1max;t1+=dt) {

for (i=0;i<4;i++) k1[i]=dt*pf[i](t1,a[0],a[1],a[2],a[3]);
for (i=0;i<4;i++) k2[i]=dt*pf[i](t1+dt/2.0,a[0]+k1[0]/2.0,a[1]+k1[1]/2.0,a[2]+k1[2]/2.0,a[3]+k1[3]/2.0);
for (i=0;i<4;i++) k3[i]=dt*pf[i](t1+dt/2.0,a[0]+k2[0]/2.0,a[1]+k2[1]/2.0,a[2]+k2[2]/2.0,a[3]+k2[3]/2.0);
for (i=0;i<4;i++) k4[i]=dt*pf[i](t1+dt,a[0]+k3[0],a[1]+k3[1],a[2]+k3[2],a[3]+k3[3]);
for (i=0;i<4;i++) l1[i]=dt*pg[i](t1,b[0],b[1],b[2],b[3]);
for (i=0;i<4;i++) l2[i]=dt*pg[i](t1+dt/2.0,b[0]+l1[0]/2.0,b[1]+l1[1]/2.0,b[2]+l1[2]/2.0,b[3]+l1[3]/2.0);
for (i=0;i<4;i++) l3[i]=dt*pg[i](t1+dt/2.0,b[0]+l2[0]/2.0,b[1]+l2[1]/2.0,b[2]+l2[2]/2.0,b[3]+l2[3]/2.0);
for (i=0;i<4;i++) l4[i]=dt*pg[i](t1+dt,b[0]+l3[0],b[1]+l3[1],b[2]+l3[2],b[3]+l3[3]);
for (i=0;i<=3;i++) a[i]=a[i]+((k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6.0);
for (i=0;i<=3;i++) b[i]=b[i]+((l1[i]+2.0*l2[i]+2.0*l3[i]+l4[i])/6.0); //4次のルンゲクッタ

da[0]=-0.8*a[0];
da[1]=-0.8*a[1];
da[2]=-0.2*a[2];
da[3]=-0.5*a[3]; //AからBへの移動量

db[0]=-b[0];
db[1]=-b[1];
db[2]=-0.7*b[2];
db[3]=-0.9*b[3]; //BからAへの移動量

if (t1 < 21.0) { //移動開始
for (i=0;i<=3;i++) {
a[i]=a[i]+da[i];
b[i]=b[i]-da[i]; ////AからBへ移動後のAにおけるSEIRの数,AからBへ移動後のBにおけるSEIRの数
}
} else {
for (i=0;i<=3;i++) {
a[i]=a[i]-db[i]; ////BからAへの移動後のAにおけるSEIRの数
b[i]=b[i]+db[i]; ////BからAへの移動後のBにおけるSEIRの数
}
}
printf("%4.2f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n",t1+0.1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);
}
return 0;
}
こんな感じだろう。

お礼

いつもお世話になっております。エラーなく正常に動きました！
Aを例えば駅、Bを学校とし、ある時間（放課後）を過ぎると逆にBからAへ行くような
シュミレーションをしてみたのですが、いかがでしょうか？
また質問させていただいた際には是非ご回答いただけると幸いです。
ありがとうございました。

回答No.2

chie65535

＞しか表示されず、

プログラムの波括弧の対応がおかしいので「初期値」と「最終結果」のみしか表示されません。

貴方は「t1が0から40.1まで0.1刻みで401回分、表示される筈」でプログラムを書いたつもりですが、そうなってはいません。

＞またseirすべての合計は350で一定であるはずなのに516となっています。

感染が進むうち「出生」と「死亡」が発生するので、最終結果が「人口増加」で終わっただけの話です。

このプログラムでは「各種パラメータの設定値」により「全員死滅して終わり」になったり、「出生よりも感染による死亡が多くなり人口減少して終わり」になったり、「出生と感染による死亡が釣り合って人口が変わらずに終わり」になったり、「感染による死亡よりも出生が多くなり人口増加して終わり」になったりします。

＞どこがおかしいのでしょうか？

「波括弧の位置」がおかしいです。

以下のプログラムと見比べて下さい（「波括弧」を追加しただけです）

#include <stdio.h>
#include <math.h>

double alf=0.33, bet=0.0875, gam=0.2; ///alf:移行率,bet:感染率,gam:除外率
double f1(double t1,double a0,double a1,double a2,double a3){return -bet*a0*a2;}
double f2(double t1,double a0,double a1,double a2,double a3){return bet*a0*a2-alf*a1;}
double f3(double t1,double a0,double a1,double a2,double a3){return alf*a1-gam*a2;}
double f4(double t1,double a0,double a1,double a2,double a3){return gam*a2;} //箱Aの関数
double g1(double t1,double b0,double b1,double b2,double b3){return -bet*b0*b2;}
double g2(double t1,double b0,double b1,double b2,double b3){return bet*b0*b2-alf*b1;}
double g3(double t1,double b0,double b1,double b2,double b3){return alf*b1-gam*b2;}
double g4(double t1,double b0,double b1,double b2,double b3){return gam*b2;} //箱Bの関数

int main(void) {
　double t1,dt,t1max;
　double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[4],a[4],b[4],da[4],db[4] ;
　double (*pf[4])(double,double,double,double,double)={f1,f2,f3,f4};
　double (*pg[4])(double,double,double,double,double)={g1,g2,g3,g4};
　int i; //宣言

　t1 = 0.0;
　dt = 0.1;
　t1max = 40.0; //時間初期値

　a[0] = 200.0;
　a[1] = 40.0;
　a[2] = 30.0;
　a[3] = 30.0; //箱A初期値（a[0]:感受性人口、a[1]:潜伏人口、a[2]:感染人口、a[3]:隔離人口）

　b[0] = 20.0;
　b[1] = 8.0;
　b[2] = 12.0;
　b[3] = 10.0; //箱B初期値（b[0]:感受性人口、b[1]:潜伏人口、b[2]:感染人口、b[3]:隔離人口）

　printf("%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n",t1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);

　for(t1=0.0;t1<=t1max;t1+=dt) {

　　for (i=0;i<4;i++) k1[i]=dt*pf[i](t1,a[0],a[1],a[2],a[3]);
　　for (i=0;i<4;i++) k2[i]=dt*pf[i](t1+dt/2.0,a[0]+k1[0]/2.0,a[1]+k1[1]/2.0,a[2]+k1[2]/2.0,a[3]+k1[3]/2.0);
　　for (i=0;i<4;i++) k3[i]=dt*pf[i](t1+dt/2.0,a[0]+k2[0]/2.0,a[1]+k2[1]/2.0,a[2]+k2[2]/2.0,a[3]+k2[3]/2.0);
　　for (i=0;i<4;i++) k4[i]=dt*pf[i](t1+dt,a[0]+k3[0],a[1]+k3[1],a[2]+k3[2],a[3]+k3[3]);
　　for (i=0;i<4;i++) l1[i]=dt*pg[i](t1,b[0],b[1],b[2],b[3]);
　　for (i=0;i<4;i++) l2[i]=dt*pg[i](t1+dt/2.0,b[0]+l1[0]/2.0,b[1]+l1[1]/2.0,b[2]+l1[2]/2.0,b[3]+l1[3]/2.0);
　　for (i=0;i<4;i++) l3[i]=dt*pg[i](t1+dt/2.0,b[0]+l2[0]/2.0,b[1]+l2[1]/2.0,b[2]+l2[2]/2.0,b[3]+l2[3]/2.0);
　　for (i=0;i<4;i++) l4[i]=dt*pg[i](t1+dt,b[0]+l3[0],b[1]+l3[1],b[2]+l3[2],b[3]+l3[3]);}
　　for (i=0;i<=3;i++) a[i]=a[i]+((k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6.0);
　　for (i=0;i<=3;i++) b[i]=b[i]+((l1[i]+2.0*l2[i]+2.0*l3[i]+l4[i])/6.0); //4次のルンゲクッタ

　　da[0]=-0.8*a[0];
　　da[1]=-0.8*a[1];
　　da[2]=-0.2*a[2];
　　da[3]=-0.5*a[3]; //AからBへの移動量

　　db[0]=-b[0];
　　db[1]=-b[1];
　　db[2]=-0.7*b[2];
　　db[3]=-0.9*b[3]; //BからAへの移動量

　　if (t1 < 21.0) { //移動開始
　　　for (i=0;i<=3;i++) {
　　　　a[i]=a[i]+da[i];
　　　　b[i]=b[i]-db[i]; ////AからBへ移動後のAにおけるSEIRの数,AからBへ移動後のBにおけるSEIRの数
　　　}
　　} else {
　　　for (i=0;i<=3;i++) {
　　　　a[i]=a[i]-da[i]; ////BからAへの移動後のAにおけるSEIRの数
　　　　b[i]=b[i]+db[i]; ////BからAへの移動後のBにおけるSEIRの数
　　　}
　　}
　　printf("%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n",t1+0.1,a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3],a[0]+a[1]+a[2]+a[3]+b[0]+b[1]+b[2]+b[3]);
　}
　return 0;
}

お礼

printfがelse内に入ってましたね。
ありがとうございました。

回答No.1

maiko0333

間違いっちゅうかdouble を使う限り、小数点以下の数字は近似値にしかなりませんよ。