配列を用いた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;
}
補足
お知恵を貸してくださりありがとうございます。 K[i_]w->F[i,w]でなくて K[i_]->F[i,w]/wを用いるのは面白いなと思いました。 ただ、説明不足だったと思いますが、 (K[1]K[1]+k[2]) w not equal K[1]w K[1] w + K[2] w であり、 (K[1]K[1]+k[2]) w = K[1]K[1]w +K[2]w =F[1,K[1]w] + F[2,w] =F[1,F[1,w]] + F[2,w] というように汎関数的に適用するようなルールになっていて 難しいのです。 線形汎関数ですので、 F[i,a*p(w)]=a*F[i,p(w)] などとしたいのですが、いちいちa*の後のパターンp(w)を指定して あげて変形しています。もっといい方法がないのかなあと思っています。