配列を用いたC言語プログラミングについて

質問

以下のルンゲクッタ法を用いたプログラムに配列などを使いさらに短くしたいのですが
どのような方法が有りますか？

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

double f1(double t1,double w,double x,double y,double z);
double f2(double t1,double w,double x,double y,double z);
double f3(double t1,double w,double x,double y,double z);
double f4(double t1,double w,double x,double y,double z); //箱Aの関数

double g1(double t1,double a,double b,double c,double d);
double g2(double t1,double a,double b,double c,double d);
double g3(double t1,double a,double b,double c,double d);
double g4(double t1,double a,double b,double c,double d); //箱Bの関数

int main(void)
{
double t1,w,x,y,z,a,b,c,d,dt,t1max,t2max,lam,gam,lat,dw,dx,dy,dz,da,db,dc,dd ;

double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[4] ; ///宣言

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

w = 200.0;
x = 40.0;
y = 30.0;
z = 30.0; ///箱A初期値（w：感受性人口、x:潜伏人口、y:感染人口、z:隔離人口）

a = 20.0;
b = 8.0;
c = 12.0;
d = 10.0; ///箱B初期値（a:感受性人口、b:潜伏人口、c:感染人口,d:隔離人口）

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

k1[0]=dt*f1(t1,w,x,y,z);
k1[1]=dt*f2(t1,w,x,y,z);
k1[2]=dt*f3(t1,w,x,y,z);
k1[3]=dt*f4(t1,w,x,y,z);
k2[0]=dt*f1(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0);
k2[1]=dt*f2(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0);
k2[2]=dt*f3(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0);
k2[3]=dt*f4(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0);
k3[0]=dt*f1(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0);
k3[1]=dt*f2(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0);
k3[2]=dt*f3(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0);
k3[3]=dt*f4(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0);
k4[0]=dt*f1(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]);
k4[1]=dt*f2(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]);
k4[2]=dt*f3(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]);
k4[3]=dt*f4(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]); ///箱Aルンゲクッタ

l1[0]=dt*g1(t1,a,b,c,d);
l1[1]=dt*g2(t1,a,b,c,d);
l1[2]=dt*g3(t1,a,b,c,d);
l1[3]=dt*g4(t1,a,b,c,d);
l2[0]=dt*g1(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0);
l2[1]=dt*g2(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0);
l2[2]=dt*g3(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0);
l2[3]=dt*g4(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0);
l3[0]=dt*g1(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0);
l3[1]=dt*g2(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0);
l3[2]=dt*g3(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0);
l3[3]=dt*g4(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0);
l4[0]=dt*g1(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]);
l4[1]=dt*g2(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]);
l4[2]=dt*g3(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]);
l4[3]=dt*g4(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]); ///箱Bルンゲクッタ

w=w+((k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0);
x=x+((k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0);
y=y+((k1[2]+2.0*k2[2]+2.0*k3[2]+k4[2])/6.0);
z=z+((k1[3]+2.0*k2[3]+2.0*k3[3]+k4[3])/6.0);

a=a+((l1[0]+2.0*l2[0]+2.0*l3[0]+l4[0])/6.0);
b=b+((l1[1]+2.0*l2[1]+2.0*l3[1]+l4[1])/6.0);
c=c+((l1[2]+2.0*l2[2]+2.0*l3[2]+l4[2])/6.0);
d=d+((l1[3]+2.0*l2[3]+2.0*l3[3]+l4[3])/6.0);
}
return 0;
}

回答No.3

trapezium · ベストアンサー

#2ですが補足しときますね。コンパイルは通りますが動作するかまでは検証してません。

typedef double (*func_t)(double, double, double, double, double);
typedef double box_t[5][4];
void box(double t1, double dt, func_t f[], box_t k)
{
int i;
#define arg2(n,i) (k[0][i] + k[n][i] / 2.0)
#define arg(i) (k[0][i] + k[3][i])
for (i = 0; i < 4; i++)
k[1][i]=dt*f[i](t1,k[0][0],k[0][1],k[0][2],k[0][3]);
for (i = 0; i < 4; i++)
k[2][i]=dt*f[i](t1+dt/2.0,arg2(1,0),arg2(1,1),arg2(1,2),arg2(1,3));
for (i = 0; i < 4; i++)
k[3][i]=dt*f[i](t1+dt/2.0,arg2(2,0),arg2(2,1),arg2(2,2),arg2(2,3));
for (i = 0; i < 4; i++)
k[4][i]=dt*f[i](t1+dt,arg(0),arg(1),arg(2),arg(3));
for (i = 0; i < 4; i++)
k[0][i] += ((k[1][i]+2.0*k[2][i]+2.0*k[3][i]+k[4][i])/6.0);
}

int main(void)
{
box_t k, l; ///宣言
func_t f[] = {f1, f2, f3, f4};
func_t g[] = {g1, g2, g3, g4};
double t1 = 0.0;
double const dt = 0.3;
double const t1max = 40.0; //時間初期値

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

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

for(t1=0.0;t1<=t1max;t1+=dt) {
box(t1, dt, f, k);
box(t1, dt, g, l);
}
return 0;
}

お礼

補足有り難うございます。
２次元配列を使えばかなり短くできそうですね！
C初心者なものでこれから２次元配列について勉強し、活用していこうかと思います。
大変参考になりました。ありがとうございました。

回答No.2

trapezium

A, B のセットは同じ演算パラメータですよね？
それならワンセット分で関数化したらいいでしょう。

現状だと変数が無駄にばらばらなので配列にまとめて、関数の f1～4(), g1～4() なども
typedef double (*func_t)(double, double, double, double, double);
func_t f[] = {f1, f2, f3, f4};
func_t g[] = {g1, g2, g3, g4};
としておけば

> k1[0]=dt*f1(t1,w,x,y,z);
> k1[1]=dt*f2(t1,w,x,y,z);
> k1[2]=dt*f3(t1,w,x,y,z);
> k1[3]=dt*f4(t1,w,x,y,z);
を例えば
double k[5][4];
int i;
for (i = 0; i < 3; i++)
k[1][i]=dt*f[i](t1, k[0][0], k[0][1], k[0][2], k[0][3]);
以下 k2～4, wxyz 部分も同様にできそうな印象です。k[0][0] が w に相当
あとはマクロを適当に活用するのも可読性が良くなります。

お礼

２次元配列ですか！
これは参考になりそうです。
ありがとうございました。

回答No.1

f272

こんな感じ？
#include <stdio.h>
#include <math.h>

double alf=0.33, bet=0.0875, gam=0.2;

double f1(double t1,double w,double x,double y,double z){return -bet*w*y;}
double f2(double t1,double w,double x,double y,double z){return bet*w*y-alf*x;}
double f3(double t1,double w,double x,double y,double z){return alf*x-gam*y;}
double f4(double t1,double w,double x,double y,double z){return gam*y;} //箱Aの関数
double g1(double t1,double a,double b,double c,double d){return -bet*a*c;}
double g2(double t1,double a,double b,double c,double d){return bet*a*c-alf*b;}
double g3(double t1,double a,double b,double c,double d){return alf*b-gam*c;}
double g4(double t1,double a,double b,double c,double d){return gam*c;} //箱Bの関数

int main(void)
{
double t1,w,x,y,z,a,b,c,d,dt,t1max;

double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[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.03;//0.3;
t1max = 40.0; //時間初期値

w = 200.0;
x = 40.0;
y = 30.0;
z = 30.0; ///箱A初期値（w：感受性人口、x:潜伏人口、y:感染人口、z:隔離人口）

a = 20.0;
b = 8.0;
c = 12.0;
d = 10.0; ///箱B初期値（a:感受性人口、b:潜伏人口、c:感染人口,d:隔離人口）

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

for (i=0;i<4;i++) k1[i]=dt*pf[i](t1,w,x,y,z);
for (i=0;i<4;i++) k2[i]=dt*pf[i](t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0);
for (i=0;i<4;i++) k3[i]=dt*pf[i](t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0);
for (i=0;i<4;i++) k4[i]=dt*pf[i](t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]);

for (i=0;i<4;i++) l1[i]=dt*pg[i](t1,a,b,c,d);
for (i=0;i<4;i++) l2[i]=dt*pg[i](t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0);
for (i=0;i<4;i++) l3[i]=dt*pg[i](t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0);
for (i=0;i<4;i++) l4[i]=dt*pg[i](t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]);

w=w+((k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0);
x=x+((k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0);
y=y+((k1[2]+2.0*k2[2]+2.0*k3[2]+k4[2])/6.0);
z=z+((k1[3]+2.0*k2[3]+2.0*k3[3]+k4[3])/6.0);

a=a+((l1[0]+2.0*l2[0]+2.0*l3[0]+l4[0])/6.0);
b=b+((l1[1]+2.0*l2[1]+2.0*l3[1]+l4[1])/6.0);
c=c+((l1[2]+2.0*l2[2]+2.0*l3[2]+l4[2])/6.0);
d=d+((l1[3]+2.0*l2[3]+2.0*l3[3]+l4[3])/6.0);

//printf("%f %f %f %f %f %f\n", t1,w,x,y,z,w+x+y+z);
//printf("%f %f %f %f %f %f\n", t1,a,b,c,d,a+b+c+d);
}
return 0;
}

お礼

SEIRの微分方程式をそこに書けるとは知りませんでした！
実は単純に２つの箱それぞれの解を求めているのではなく、この後相互作用を書くつもりでしたが字数的に厳しかったので省きました。
ありがとうございました。