请各位学长帮忙修改个程序！
下面这个程序是用龙格库塔法求解微分方程的，哪位学长能帮忙把它改成用数组形式的，我不太会，想学习一下，我编的这个太繁杂可，谢谢！
/*龙格库塔法求解微分方程*/
#include "stdio.h"
#include "math.h"
#define A 0 /*定义初值和步长*/
#define B 0.002
#define C 1
#define D 1
#define E 1
double f(double,double,double,double); /* 定义函数*/
double g(double,double,double,double);
double e(double,double,double,double);
main()
{
double x0=A, h=B,
y01=C, u01=D, v01=E,
j01,j02,j03,j04,k01,k02,k03,k04,l01,l02,l03,l04;
int i,n=10000;
FILE *cfPtr;
if((cfPtr=fopen("luolunzi.dat","w"))==NULL) /*数据要写入的文件*/
printf("File could not be opened\n");
else
{
for (i=1;i<=n;i++)
{ /*套用公式开始计算*/
x0=x0+h;
j01=f(x0,y01,u01,v01);
k01=g(x0,y01,u01,v01);
l01=e(x0,y01,u01,v01);
j02=f(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);
k02=g(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);
l02=e(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);
j03=f(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
k03=g(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
l03=e(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
j04=f(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
k04=g(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
l04=e(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
y01=y01+h/6.0*(j01+2.0*j02+2.0*j03+j04);
u01=u01+h/6.0*(k01+2.0*k02+2.0*k03+k04);
v01=v01+h/6.0*(l01+2.0*l02+2.0*l03+l04);
printf("%15.7f\t%15.7f\t%15.7f\t\n",x0,y01,u01,v01);
fprintf(cfPtr,"%15.7f\t%15.7f\t%15.7f\t\n",x0,y01,u01,v01);
}
fclose(cfPtr);
}
}
double f(double s1,double q1,double u1,double r1) /*定义的微分方程*/
{
double z1;
z1=-10*q1+(10*u1);
return z1;
}
double g(double s2,double q2,double u2,double r2)
{
double z2;
z2=28*q2-u2-q2*r2;
return z2;
}
double e(double s3,double q3,double u3,double r3)
{
double z3;
z3=(-8/3)*r3+q3*u3;
return z3;
}
[解决办法]
楼主的意思是把每次计算出来的x0,y01,u01,v01用数组的形式保存下来,是吗?
[解决办法]
整理修改如下：
#include <stdio.h> #include <math.h> #define A 0 /*定义初值和步长*/ #define B 0.002 #define C 1 #define D 1 #define E 1 #define N 100double f1(double s1,double q1,double u1,double r1) //f { return (-10*q1 + 10*u1); } double f2(double s2,double q2,double u2,double r2) //g{ return (28*q2 - u2 - q2*r2);} double f3(double s3,double q3,double u3,double r3) //e{ return ((-8/3.0)*r3 + q3*u3);} int main(){ double result[3] = {C, D, E}; //y,u,v double tmp[3][4]; //j,k,i //初始化 double x = A, h = B; double (* p[3])(double, double, double, double) = {f1, f2, f3}; int i, j; FILE *cfPtr; if((cfPtr = fopen("luolunzi.dat","w")) == NULL) printf("File could not be opened\n"); else { for(i = 1; i <= N; i++) { x += h; for(j = 0; j < 3; j++) { tmp[j][0] = p[j](x, result[0], result[1], result[2]); } for(j = 0; j < 3; j++) { tmp[j][1] = p[j](x + h/2.0, result[0]+h/2.0*tmp[0][0], result[1]+h/2.0*tmp[1][0], result[2]+h/2.0*tmp[2][0]); } for(j = 0; j < 3; j++) { tmp[j][2] = p[j](x + h/2.0, result[0]+h/2.0*tmp[0][1], result[1]+h/2.0*tmp[1][1], result[2]+h/2.0*tmp[2][1]); } for(j = 0; j < 3; j++) { tmp[j][3] = p[j](x + h, result[0]+h*tmp[0][1], result[1]+h*tmp[1][1], result[2]+h*tmp[2][1]); } printf("%15.7f\t", x); for(j = 0; j < 3; j++) { result[j] += h/6.0*(tmp[j][0]+2.0*tmp[j][1]+2.0*tmp[j][2]+tmp[j][3]); printf("%15.7f\t", result[j]); } fprintf(cfPtr,"%15.7f\t%15.7f\t%15.7f\t%15.7f\n", x, result[0], result[1], result[2]); } fclose(cfPtr); } return 0;}
[解决办法]
太麻烦了，先这样了，不合适的地方，lz自己改改吧。
/*龙格库塔法求解微分方程*/ #include "stdio.h" #include "math.h" #define A 0 /*定义初值和步长*/ #define B 0.002 #define C 1 #define D 1 #define E 1 #define ARRAY_NUM(a) (sizeof(a) / sizeof((a)[0]))typedef double (*FUNC)(double,double,double,double);struct R_K_PARA { double x; double y; double u; double v;};typedef struct R_K_PARA R_K_PARA_t;double f(double,double,double,double); /* 定义函数*/ double g(double,double,double,double); double e(double,double,double,double); FUNC r_k_func[] = { f, g, e,};#define R_K_NUM ARRAY_NUM(r_k_func)R_K_PARA_t r_k(FUNC* func, int fn, R_K_PARA_t pa, double h);int main(void) { double h = B; int i, n = 10000; R_K_PARA_t pa = {A, C, D, E}, ret; FILE *cfPtr; if ((cfPtr = fopen("luolunzi.dat", "w")) == NULL) /*数据要写入的文件*/ printf("File could not be opened\n"); else { for (i = 1; i <= n; i++) { /*套用公式开始计算*/ pa.x += h; pa = r_k(r_k_func, R_K_NUM, pa, h); printf("%15.7f\t%15.7f\t%15.7f\t%15.7f\t\n", pa.x, pa.y, pa.u, pa.v); fprintf("%15.7f\t%15.7f\t%15.7f\t%15.7f\t\n", ret.x, ret.y, ret.u, ret.v); } fclose(cfPtr); } return 0;}R_K_PARA_t r_k(FUNC* func, int fn, R_K_PARA_t pa, double h){ int i, j; double rk[] = {0, 1 / 2.0, 1 / 2.0, 1}; double jkl[ARRAY_NUM(rk)][R_K_NUM]; for (i = 0; i < ARRAY_NUM(rk); i++) { for (j = 0; j < R_K_NUM; j++) { if (i == 0) { jkl[i][j] = func[j](pa.x, pa.y, pa.u, pa.v); } else { jkl[i][j] = func[j](pa.x + h * rk[i], pa.y + h * rk[i] * jkl[i - 1][0], pa.u + h * rk[i] * jkl[i - 1][1], pa.v + h * rk[i] * jkl[i - 1][2]); } } } pa.y += (jkl[0][0] + jkl[1][0] * 2.0 + jkl[2][0] * 2.0 + jkl[3][0]) * h / 6.0; pa.u += (jkl[0][1] + jkl[1][1] * 2.0 + jkl[2][1] * 2.0 + jkl[3][1]) * h / 6.0; pa.v += (jkl[0][2] + jkl[1][2] * 2.0 + jkl[2][2] * 2.0 + jkl[3][2]) * h / 6.0; return pa;}double f(double s1,double q1,double u1,double r1) /*定义的微分方程*/ { double z1; z1=-10*q1+(10*u1); return z1; } double g(double s2,double q2,double u2,double r2) { double z2; z2=28*q2-u2-q2*r2; return z2; } double e(double s3,double q3,double u3,double r3) { double z3; z3=(-8/3)*r3+q3*u3; return z3; }