a.outの結果がおかしい件について -SEIRモデルの解を求めるプログラ-

２次元配列を関数の引数にする方法について

以下のプログラムはSEIRモデル(https://ja.wikipedia.org/wiki/SEIR%E3%83%A2%E3%83%87%E3%83%AB)を２つたて、それらを相互作用させたときのA,Bそれぞれの解をルンゲクッタ法を用いて求めています。ここで一次元配列a[i],b[j]を２次元配列a[2][4](i:領域A,B,j:SEIR)とし使用した場合、関数の定義のところなどどのように書けばよいのでしょうか？ C言語初心者でして、お見苦しい点多々あると思いますが、何卒ご教示お願いします。。。 #include <stdio.h> #include <math.h> double alf=0.05, bet=0.005, gam=0.05; ///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 = 20.0; //時間初期値 a[0] = 2000.0; a[1] = 300.0; a[2] = 400.0; a[3] = 300.0; //箱A初期値（a[0]:感受性人口、a[1]:潜伏人口、a[2]:感染人口、a[3]:隔離人口） b[0] = 1000.0; b[1] = 200.0; b[2] = 300.0; b[3] = 500.0; //箱B初期値（b[0]:感受性人口、b[1]:潜伏人口、b[2]:感染人口、b[3]:隔離人口） printf("%4.2f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f\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.05*a[0]; da[1]=-0.05*a[1]; da[2]=-0.01*a[2]; da[3]=-0.05*a[3]; //AからBへの移動量 db[0]=-0.05*b[0]; db[1]=-0.05*b[1]; db[2]=-0.02*b[2]; db[3]=-0.05*b[3]; //BからAへの移動量 if (t1 < 15.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.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f,%10.3f\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; }

コンパイルエラーの原因と対策について。

以下のプログラムは、SEIRモデルをルンゲクッタ法を使いA,Bそれぞれについて解を求めているのですが、どうしてもコンパイルエラーが出てしまいます。字数制限でエラー内容は書けませんがご希望があれば補足いたします。原因と改良すべき点をご教示お願いします。 #include <stdio.h> #include <math.h> #define alf 0.33 #define bet 0.0875 #define gam 0.2 double f1(double t1,double a[0],double a[1],double a[2],double a[3]); double f2(double t1,double a[0],double a[1],double a[2],double a[3]); double f3(double t1,double a[0],double a[1],double a[2],double a[3]); double f4(double t1,double a[0],double a[1],double a[2],double a[3]); //箱Aの関数 double g1(double t1,double b[0],double b[1],double b[2],double b[3]); double g2(double t1,double b[0],double b[1],double b[2],double b[3]); double g3(double t1,double b[0],double b[1],double b[2],double b[3]); double g4(double t1,double b[0],double b[1],double b[2],double b[3]); //箱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; } double f1(double t1,double a[0],double a[1],double a[2],double a[3]) { double r; r=-bet*a[0]*a[2]; return r; } double f2(double t1,double a[0],double a[1],double a[2],double a[3]) { double r; r=bet*a[0]*a[2]-alf*a[1]; return r; } double f3(double t1,double a[0],double a[1],double a[2],double a[3]) { double r; r=alf*a[1]-gam*a[2]; return r; } double f4(double t1,double a[0],double a[1],double a[2],double a[3]) { double r; r=gam*a[2]; return r; } double g1(double t1,double b[0],double b[1],double b[2],double b[3]) { double r; r=-bet*b[0]*b[2]; return r; } double g2(double t1,double b[0],double b[1],double b[2],double b[3]) { double r; r=bet*b[0]*b[2]-alf*b[1]; return r; } double g3(double t1,double b[0],double b[1],double b[2],double b[3]) { double r; r=alf*b[1]-gam*b[2]; return r; } double g4(double t1,double b[0],double b[1],double b[2],double b[3]) { double r; r=gam*b[2]; return r; }

配列を用いた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; }

大学のC言語の課題で

ルンゲクッタ法を用いたバネの単振動を解析する課題をやっています。以下の部分まで作ったのですが、どうしても最後のほうがわかりません。関数fとf2のreturn やメイン関数を変更すればいいとおもうのですが、どうやっても満足な結果が得られません。どなたか詳しい方教えてもらえない？ #include<stdio.h> #include<math.h> double f(double t,double y1) { return y1; } double f2(double t2,double y2) { double w=5.0,k=0.5; return -w*w*f(t2,y2)-2.0*k*y2; } int main(void) { int count,bunkatu,npr; double t_0,T; double k1,k2,k3,k4; double y,t,dt,dt_2,dt_6; t_0 = 0.0; T=10.0; bunkatu=2000; npr=1000; y=1.0; printf("\n# [t_0,T]=[%6.4lf,%6.4lf] N=%d\n",t_0,T,bunkatu); t=t_0; dt=(T-t_0)/(double)bunkatu; dt_2=dt/2.0; dt_6=dt/6.0; for(count=0;count <=bunkatu;count++) { if(count %(bunkatu/npr)==0) printf("\n %5.2lf %10.5f",t,y); k1=f2(t,y); k2=f(t+dt_2,y+k1*dt_2); k3=f(t+dt_2,y+k2*dt_2); k4=f(t+dt,y+k3*dt); y+=(k1+2.0*k2+2.0*k3+k4)*dt_6; t+=dt; } printf("\n"); return 0; }

ルンゲクッタ法で数値解析（Ｃ言語）

コンピュータ（プログラミング）のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました（Ｃ言語を使って）。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #include <stdio.h> double f1(double t,double x,double v); double f2(double t,double x,double v); int main() { double t,x,v,dt,tmax; double k1[2],k2[2],k3[2],k4[2]; x=1.0; //位置の初期値 v=0.0; //速度の初期値 dt=0.01; //刻み幅 tmax=500.0; //繰り返し最大回数 FILE *output; output=fopen("output.dat","w"); for(t=0;t<tmax;t+=dt) { k1[0]=dt*f1(t,x,v); k1[1]=dt*f2(t,x,v); k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]); k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]); x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0; v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0; fprintf(output,"%f %f %f\n",t,x,v); } fclose(output); return 0; } double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); } ---- このソースは単振動のもので、減衰振動のときは、最後の方の ---- double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); ---- の部分が変わるのだと思うのですが、よく分かりません。減衰振動は、mx"=-kx-rx'で表されるので、この式を変形（？）したものが、入るのかな、という予想程度にしか分かりません。ソースをどのように変えれば減衰振動を解析できるのでしょうか。どなたか詳しい方、教えてください。お願いします。

ルンゲクッタのプロクラム実行できません

プログラムの実行ができません。どこを直したらよいのか教えて下さい。よろしくお願いします。 #include <stdio.h> #include <stdlib.h> #include <math.h> double *dvector(long i,long j); void free_dvector(double *a,long i); double func(double x,double y); double *rk4(double y0,double a,double b, int n,double(*f)()); int main(void) { double *y,h,a=0.0, b=1.0, y0=1.0; int i,n; printf("分割数を入力して下さい------->"); scanf("%d",&n); y=dvector(0,n); y=rk4(y0,a,b,n,func); h=(b-a)/n; for(i=0;i<=n;i++) { printf("x=%f \t y=%f \n",a+i*h,y[i]); } return 0; } double *rk4(double y0,double a, double b,int n,double (*f)()) { double k1,k2,k3,k4,h,*y,x; int i; y=dvector(0,n); h=(b-a)/n; y[0]=y0; x=a; for(i=0;i<n;i++) { k1=f(x,y[i]); k2=f(x+h/2.0,y[i]+h*k1/2.0); k3=f(x+h/2.0,y[i]+h*k2/2.0); k4=f(x+h,y[i]+h*k3); y[i+1]=y[i]+h/6.0*(k1+2.0*k2+2.0*k3+k4); x +=h; } return y; free_dvector(y,0); }

微分方程式（ルンゲ＝クッタ）がうまくいきません

C言語で物体の自由落下運動（空気抵抗を考慮する）を微分方程式を用いて解かなければならないのですが、どうしてもグラフの形が曲線になりません。どこが間違っているのか指摘お願いします。空気抵抗を考慮しているので、落下速度はどんどんある一定の値に近づいていくはずなのですが、グラフは直線になってしまいます。いろいろと方法を試しましたがどれも外れました。写真は赤が今の状態（間違っています）、黒が理想形です。ここではルンゲ＝クッタ法を使用しています。 #include <stdio.h> #include <math.h> double yd(double t,double y){ double m = 50.0; 質量 double g = 9.8; 重力加速度 double k = 0.01; 空気抵抗の定数係数 double a= g - ((k *y)/m); 運動方程式 return(a); } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2 = h / 2.0; k1 = yd(t,y); k2 = yd(t+h2,y+(h2*k1)); k3 = yd(t+h2,y+(h2*k2)); k4 = yd(t+h,y+(h*k3)); a = y + (h/6.0) * ((k1 + (2.0 * k2) + (2.0 * k3) + k4)); return(a); } void main(){ double y = 0.0; double h = 0.001; double t = 0.0; int i; for(i = 0; i < 1000; i++){ t += h; y = runge(t,y,h); printf("%f %f\n",t,y); } }

Javaプログラムのフローチャートについて

下に記述したのは4次のルンゲクッタ法のJavaによるプログラム例です。このプログラムの簡単なフローチャートを作成してプレゼンしたいのですが、作成した経験がなく、本等をみてもいまいちわかりません。どなたかご教授いただきたく、お願い申し上げます。 // 2階常微分方程式に対する Runge-Kutta 法による積分メソッド public void Runge_Kutta(double x[], double y[]){ double xi, vi, k1, k2, k3, k4, l1, l2, l3, l4; for(int i = 0; i < iMax-1; i++){ xi = x[i]; vi=y[i]; k1 = dt*vi; l1 = dt*fun(vi, xi); k2 = dt*(vi+0.5*l1); l2 = dt*fun(vi+0.5*l1, xi+0.5*k1); k3 = dt*(vi+0.5*l2); l3 = dt*fun(vi+0.5*l2, xi+0.5*k2); k4 = dt*(vi+l3); l4 = dt*fun(vi+l3, xi+k3); x[i+1] = xi+(k1+2.0*(k2+k3)+k4)/6.0; y[i+1] = vi+(l1+2.0*(l2+l3)+l4)/6.0; if(y[i] < ymin){ Max=i; break;} if(Math.abs(y[i+1]-y[i]) < diff && Math.abs(x[i+1]-x[i])< diff){Max=i; break;} } }

セグメントエラー

#include <stdio.h> #define S 10.0 #define L 10.0 void solve(int m,double a[][10],double f[]); main() { double t,me[10][10],fe[10]; int i,x; t=(S*L)/420; me[0][0]=t*156;me[0][1]=t*22*L;me[0][2]=t*54;me[0][3]=t*(-13)*L; me[1][0]=me[0][1];me[1][1]=t*4*L*L;me[1][2]=t*13*L;me[1][3]=t*(-3)*L*L; me[2][0]=me[0][2];me[2][1]-me[1][2];me[2][2]=t*156;me[2][3]=t*(-22)*L; me[3][0]=me[0][3];me[3][1]=me[1][3];me[2][3]=me[3][2];me[3][3]=t*4*L*L; for(i=0;i<4;i++){ fe[i]=1.0; } x=4; solve(x,me,fe); } void solve(int m,double a[][10],double f[]) { /* m=number of unknowns */ int i,j,k; double aa; /* forward elimination */ for(i=0;i<m-1;i++){ for(j=i+1;j<m;j++){ aa=a[j][i]/a[i][i]; f[j]-=aa*f[i]; for(k=i+1;k<m;j++){ a[j][k]-=aa*a[i][k]; } } } /* backward substitution */ f[m-1]/=a[m-1][m-1]; for(i=m-2;i>=0;i--){ for(j=i+1;j<m;j++){ f[i]-=a[i][j]*f[j]; } f[i]/=a[i][i]; } } 以上のようなプログラムについて、コンパイルはできるのですが、実行時に「セグメントエラー」と出てしまいます。たぶん配列の引渡しに原因があるようなのですが、どなたか詳しい原因（対処法）をご教授願いたいと思います。

お初に質問させていただきます。最近、Excelのマクロを始め、下のよう

お初に質問させていただきます。最近、Excelのマクロを始め、下のようなユーザー定義関数を作成し無事動作を確認しました。そして、他の方法で実測した任意の時刻における値との差の二乗の総和が最小になるようにK(1),K(2)の値を最適化したいと思い、最小二乗の総和をワークシート上に計算し、Excelのツールのソルバーを用いて最適化を試みました。しかし、「最適解は見つかりました」と出るのですが、値に全く変化はなく最適化できていません。ユーザー定義関数を含むとソルバーは使えないのですか？もし何か解決法があれば教えてください。私が使っているExcelはExcel 2003です。 Public Function RK1(成分, 時刻かける20) As Double Application.Volatile Dim Koc As Double, Kco As Double, Eo As Integer, Ec As Integer, Y0(2) As Double Dim Y(2) As Double, K1(4) As Double, K2(4) As Double, K3(4) As Double Dim K4(4) As Double, T As Double Dim K(2) As Double, DelY(2) As Double Dim Tt(1500) As Double, TV(2, 1500) As Double For I = 1 To 2 Y0(I) = Cells(3 + I, 2).Value Const nn As Integer = 20 dx = 1# / nn T = 0# K(1) = Cells(2, 2).Value K(2) = Cells(3, 2).Value Next I For J = 1 To 300 Y(1) = Y0(1) Y(2) = Y0(2) K1(1) = F1(T, Y(1), Y(2)) * dx K1(2) = F2(T, Y(1), Y(2)) * dx ' Y(1) = Y0(1) + 0.5 * K1(1) Y(2) = Y0(2) + 0.5 * K1(2) K2(1) = F1(T, Y(1), Y(2)) * dx K2(2) = F2(T, Y(1), Y(2)) * dx ' Y(1) = Y0(1) + 0.5 * K2(1) Y(2) = Y0(2) + 0.5 * K2(2) K3(1) = F1(T, Y(1), Y(2)) * dx K3(2) = F2(T, Y(1), Y(2)) * dx ' Y(1) = Y0(1) + K3(1) Y(2) = Y0(2) + K3(2) K4(1) = F1(T, Y(1), Y(2)) * dx K4(2) = F2(T, Y(1), Y(2)) * dx ' DelY(1) = (K1(1) + 2# * K2(1) + 2# * K3(1) + K4(1)) / 6# DelY(2) = (K1(2) + 2# * K2(2) + 2# * K3(2) + K4(2)) / 6# T = T + dx Y0(1) = Y0(1) + DelY(1) Y0(2) = Y0(2) + DelY(2) Tt(J) = T TV(1, J) = Y0(1) TV(2, J) = Y0(2) Next J RK1 = TV(成分, 時刻かける20) End Function Function F1(T, a, b) Dim K(2) As Double Eo = Cells(7, 2).Value Ec = Cells(8, 2).Value K(1) = Cells(2, 2).Value K(2) = Cells(3, 2).Value F1 = (-K(1) * a + K(2) * b) * (1 - 10 ^ (-(Eo * a + Ec * b))) / (Eo * a + Ec * b) End Function Function F2(T, a, b) Dim K(2) As Double Eo = Cells(7, 2).Value Ec = Cells(8, 2).Value K(1) = Cells(2, 2).Value K(2) = Cells(3, 2).Value F2 = (K(1) * a - K(2) * b) * (1 - 10 ^ (-(Eo * a + Ec * b))) / (Eo * a + Ec * b) End Function

Xcodeでのc言語のプログクムについて。

Xcodeでのc言語のプログクムについて。ルンゲークッタ法のプログラムを書いているのですがビルドエラーになってしまい、原因が分かりません。　ビルド結果は、 "_func",referenced from: _main in main.o "_dvector",referenced from: _main in main.o _rk4 in main.o symbol(s) not found collect2: ld returned 1 exit status 　プログラムは、 #include <stdio.h> #include <stdlib.h> double *dvector(long i, long j); /* ベクトル量領域の確保*/ void free_dvector(double *a, long i); /*　領域の解放*/ double func(double x, double y); /* 関数の定義 */ /* ルンゲークッタ法 */ double *rk4( double y0, double a, double b, int n, double (*f)()); int main (void) { double *y, h, a=0.0, b=1.0, y0=1.0 ; int i, n; printf("分割式を入力してください-->"); scanf("%d",&n); y = dvector( 0, n ); /*領域の確保*/ y = rk4( y0, a, b, n, func ); /*ルンゲクッタ法*/ /*　結果の表示 */ h = (b-a)/n ; /* 刻み幅　*/ for (i = 0 ; i <= n ; i++) { printf("x=%f \t y=%f \n", a+i*h, y[i] ); } return 0; } /* ルンゲクッタ法　*/ double *rk4( double y0, double a, double b, int n, double (*f)() ) { double k1, k2, k3, k4, h, *y, x; int i; y = dvector( 0, n ); /* y[0,1,....,n]の確保　*/ h = (b-a)/n ; /* 初期値の設定　*/ y[0] = y0; x = a; /* ルンゲクッタ法　*/ for ( i = 0 ; i < n ; i++) { k1 = f(x,y[i]); k2 = f(x+h/2.0, y[i]+h*k1/2.0); k3 = f(x+h/2.0, y[i]+h*k2/2.0); k4 = f(x+h, y[i]+h*k3); y[i+1] = y[i] + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 ); x += y[i]; } return y; free_dvector(y, 0); /* 領域の解放　*/ }

ルンゲクッタ法（RK4）で斜方投射の問題を解く

下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法（RK4）で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k：空気抵抗、m:物体の質量 g：重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2=h/2.0; k1=yd(t,y); k2=yd(t+h2,y+(h2*k1)); k3=yd(t+h2,y+(h2*k2)); k4=yd(t+h,y+(h*k3)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }

数値解析法

このHeun法のプログラムをRunge-Kutta法にするにはどうしたらいいですか？ #include <stdio.h> #include <math.h> double f1(double y) { return y; } double f2(double y) { return -4*y; } int main(){ double a=0; double b; int m=10; int n; double h; double x,y; int k; double e; double f; double k1,k2; printf("Heun法計算例：y=e^x, y=1/e^4x\n\n"); // y = e^x b = 1; for(n=100;n<=10000;n*=100){ h = (b-a)/n; printf("y' = y: h(=dx) = %.1e (y=e^x)\n",h); x = a; y = 1; for(k=0;k<=n;k++) { x = k*h; if(k%(n/m)==0) { f = exp(x); e = fabs(y-f); printf("x=%.2f, y=%f, e^x=%f er=%.0e\n",x,y,f,e); } // Heun's method k1 = h*f1(y); k2 = h*f1(y+k1); y += (k1+k2)/2; } }

微分方程式のプログラムについて

次の問題を解こうとしましたが、プログラムが動かなかったので実行できませんでした。今、非常に困っています。分かりにくいプログラムですみません。問　温度Taの空気中に置かれている物体の温度Tは、次の常微分方程式に従う。　　　dT/dt=-k(T-Ta) t=0でT=100[℃]、Ta=20[℃]、k=1×10^-2[s^-1]として、物体の温度が60℃にまで冷却される時間を求めよ。厳密解とも比較すること。自分が作ったプログラムは #include <stdio.h> #include <math.h> double k = 1e-2, tmpa=20.0; double f(double t, double tmp) { return -k*(tmp-tmpa); } int main(void) { double t0 = 0.0, tmp0 = 100, t1, tmp1, ttmp1; double k1, k2, k3, k4, h; printf("刻み幅="); scanf("%lf", h); t1 = t0+h; ttmp1=tmpa+(tmp0-tmpa)+exp(-k*(t0-t0)); printf("%8.4lf %8.4lf %8.4lf\n", t0, tmp0, ttmp1); while(t0 <= 100.0){ k1 = f(t0, tmp0); k2 = f(t0+0.5*h, tmp0+0.5*h*k1); k3 = f(t0+0.5*h, tmp0+0.5*h*k2); k4 = f(t0+h, tmp0+h*k3); tmp1 = tmp0+h*(k1+2.0*k2+2.0*k3+k4)/6.0; printf("%8.4lf %8.4lf %8.4lf\n", t1, tmp1, ttmp1); if(tmp1 <= 60.0) break; t0 = t1; tmp0 = tmp1; t1 = t0+h; } return 0; } あと厳密解はT=e^-kt(T0-Ta)+Ta求められるそうです。変な間違いが多いと思いますが、このプログラムのどこが間違っているか教えてください。厳密解だけでもいいです。お願いします。

C++言語のプログラムをfortranに変換

const int N = 100; const double q = 10.0, dt = 0.00001, Dm = 30.0, t0 = 2.0, K = 1.0, pi = 3.141592, f = 3.0; double C[N], dC[N]; double dx = q/N; for (int i = 0; i < N; ++i) C[i] = 0;　// 初期条件 for (double t = 0; t < t0; t += dt) { 　C[0] = 0;　// 境界条件1 　C[N - 1] = K*sin(2*pi*f*t);　// 境界条件2 　for (int i = 1; i < N - 1; ++i) dC[i] = (Dm*(C[i + 1] - 2*C[i] + C[i - 1])/(dx*dx))*dt; 　for (int i = 1; i < N - 1; ++i) C[i] += dC[i]; } for (int i = 0; i < N; ++i) cout << C[i] << endl; // t = t0 このプログラムをfortranに変換できる方いますか？

ストリータ・ヘルプス式

ストリータ・ヘルプス式で、ｄD/dt=K1L-K2DからｄL/ｄｔ＝-K1L (t=t0 の時、　L=L0) を用いて、 D= {K1*L0*(e^-k1t-e~k2t)/ K2-K1}+Do*e^-k2t (t=0 の時、　D=Do）を導く方法を教えてください。色々調べたのですが、途中過程が省略されていて、わからないので、できるだけ詳しくお願いします…。

-1.#IND00と出てしまうのですが・・・

vc++(2010)でガウスの消去法を使って連立方程式を解く、というプログラムを組みました。正直中身自体をきちんと理解していないので、組みながら理解しようと思って組んだのですが、結果に-1.#IND00と出てしまいました。これは何なのでしょうか？下記がそのプログラムです。 #include<stdio.h> #include<math.h> #define N 3 int main(void){ double a[N][N+1] = {{2,3,4}, {3,-2,5}, {5,4,-7}}; double b[1][N+1] = {20,14,-8}; double x[N]; int i = 0; int j = 0; int k = 0; int l = 0; int pivot = 0; double p = 0; double q = 0; double m = 0; for(i = 0; i < N; i++){ x[i] = 0; } for(i = 0; i < N; i++){ pivot = i; for(l = i; l < N; l++){ if(fabs(a[l][i]) > m){ m = fabs(a[l][i]); pivot = l; } } if(pivot != i){ for(j = 0; j < N; j++){ b[0][j] = a[i][j]; a[i][j] = a[pivot][j]; a[pivot][j] = b[0][j]; } } } for(k = 0; k < N; k++){ a[k][j] = a[k][j] / p; a[k][k] = 1; for(j = k; j < N; j++){ a[k][j] = a[k][j] / p; } for(i = k+1; i < N; i++){ q = a[i][k]; for(j = k+1; j < N; j++){ a[i][j] = a[i][j] - q*a[k][j]; } a[i][k] = 0; } } for(i = N-1; i >= 0; i--){ x[i] = a[i][N]; for(j = N-1; j > i; j--){ x[i] = x[i] - a[i][j] * x[j]; } } // for(i = 0; i < N; i++){ // for(j = 0; j < N+1; j++){ // printf("%.1f", a[i][j]); // } // printf("\n"); // } // printf("解\n"); // for(i = 0; i < N; i++){ // printf("%f\n", x[i]); // } return 0; } 最後のコメントにしてある行は解を表す時と行列を表す時で使い分けているので、実際はどちらかを外して使用しています。

配列のエラーに関して

java言語を用いて，Householder変換を用いた固有値の数値計算に挑戦してみました．しかし，次のようなエラーが発生し上手くいきません．どなたかこの問題を解決するためにお力をかしていただけないでしょうか． ----------エラー内容-------------------------------------------------------------------------------- Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: Index 0 out of bounds for length 0 at Out.Mhouse(House.java:90) at House.main(House.java:10) ---------------------------------------------------------------------------------------------------- //Householder変換 public class House{ public static void main(String[] args){ double[][] A = new double[3][3]; int n = A.length; Out out = new Out(); for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ if(j < n-1){ System.out.print(out.Mhouse(A)[i][j] + " "); }else if (j == n-1) System.out.println(out.Mhouse(A)[i][j]); }; }; }; }; class Out{ double[][] outpro(double[] x){ int n; n = x.length; double[][] A = new double[n][n]; for(int i = 0;i < n;i++ ){ for(int j = 0;j < n;j++){ A[i][j] = x[i] * x[j]; } } return A; }; double[][] Msca(double a,double[][] A){ int n = A.length; for(int i = 0;i < n; i++){ for(int j = 0;j < n;j++){ A[i][j] = a * A[i][j]; } } return A; }; double selfpro(double[] x){ double a = 0; int n = x.length; for(int i = 0;i < n; i++){ a = a + x[i] * x[i]; }; return a; }; double[] minus(double[] x, double[] y){ int n = x.length; double[] z = new double[n]; for(int i = 0;i < n;i++){ z[i] = x[i] - y[i]; }; return z; }; double[][] house_1(double[] x){ int n = x.length; double[][] A = new double[n][n]; for(int i=0;i < n;i++){ for(int j = 0;j < n;j++){ if(i == j){ A[i][j] = 1 - Msca(2/selfpro(x),outpro(x))[i][j]; }else{ A[i][j] = - Msca(2/selfpro(x),outpro(x))[i][j]; }; }; }; return A; }; double[][] house_2(double[] x){ double[][] z = new double[1][1]; z[1][1] = 1 - 2; return z; }; double[][] Mhouse(double[][] A){ int n = A.length; double[][] H = new double[n][n]; for(int i = 0;i < n;i++){ double[] x = new double[n-i]; double[] y = new double[n-i]; double[][][] L = new double[i][n-i][n-i]; for(int j = 0;j < n-i;j++){ x[j] = A[i][i+j]; if(j == 0){ y[j] = 1; }else{ y[j] = 0; }; x[j] = y[j] - x[j]; }; if(i < n-1){ L[i] = house_1(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }else if(i == n-1){ L[i] = house_2(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }; }; double[][] B = new double[n][n]; for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ for(int k = 0;k < n;k++){ B[i][j] = H[i][k] * A[k][j]; }; }; }; return A; }; };

二次元配列の引数渡し

二次元配列を関数の引数として渡し、その配列を戻り値として呼び出し元の変数に返したいです。具体的にはルンゲクッタ法で微分する関数の引数を行列で扱いたいです。 for (i=0;i<N;i++){ for (j=0;j<N;j++){ y[i][j]=rand(); } } 　　K1=h*ff(t,y); K2=h*ff(t+h/2,y+K1/2); K3=h*ff(t+h/2,y+K2/2); K4=h*ff(t+h,y+K3); Ka=(K1+2*K2+2*K3+K4)/6; double ff(int t,double y[][]){ ここで return y[][]; のような形で二次元配列を戻り値として変数に返したいのです。 } Ｃ言語初心者なのでいまいちよくわかりません、宜しくお願いいたします。

c言語でDFTのプログラムを作成したのですが

c言語でDFTのプログラムを作成しました。以下にソースを載せます。 #include<stdio.h> #include<math.h> #include<stdlib.h> #include<time.h> #define PI 3.141592653589793 #define N 64 //データ数 DFT(double result[]){ int i,k; double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部 double a,b; for(k=0;k<N;++k){ a=b=0; for(i=0;i<N;++i){ a=a+result[i]*cos(2*PI*i*k/N);　　　　　 b=b+(-1.0)*result[i]*sin(2*PI*i*k/N); } A[k]=a/N; B[k]=b/N; } for(i=0;i<N;++i){ printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示 T=T+(0.1/N); } } main(){ int i; double T=0; double Sampdata; double result[N]; for(i=0;i<N;++i){ Sampdata=5*sin(20*PI*T);　　　　　　//0~0.1秒間をN個にサンプリング result[i]=Sampdata; //サンプリングデータを代入 T=T+(0.1/N); } clock_t start,end; //処理時間計測開始 start=clock(); DFT(result); end=clock(); printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC); //処理時間表示 } 元信号には5sin(20πt)の値を入れています。この信号は周期は0.1secです。これでフーリエ変換を行うとデータ数N/2を中心に対称なデータが出てくるのですが、処理が終わるのが早い気がするんです。例えば2^15個のデータで実行しても2分もかからずに処理が終わってしまいます。一応、対称性が出てるとはいえ、終わるのが早すぎる気がするのですが、おかしい所があれば教えていただけると嬉しいです。よろしくお願いします。