各位高手帮我看看这个问题?
我略懂C,现在要计算一个微分方程的问题,以往都是用matlab的,但是速度太慢.下面是程序,书上的.但是为了要求精度,我需要取很多点,也就是步长要取的很小(就是下面红色的语句).
但是步长取小了后出现这样的错误 "错误 rkuttas.c 6: 数组太小" 请问怎么修改子函数???
子函数
/* rkuttas.c: Fourth-order Runge-kutta method for equations. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
extern void f( double x, double yy[M], double ff[M] );
void rkuttas( double x0, double x1, double y[M][N] )
{
int i,j;
double x,k1[M],k2[M],k3[M],k4[M],yy[M],ff[M],h;
h=(x1-x0)/(N-1.0);
for(i=0;i<N-1;i++)
{
x=x0+i*h;
for(j=0;j<M;j++) yy[j]=y[j][i];
f(x,yy,ff);
for(j=0;j<M;j++) k1[j]=h*ff[j];
for(j=0;j<M;j++) yy[j]=y[j][i]+0.5*k1[j];
f(x+0.5*h,yy,ff);
for(j=0;j<M;j++) k2[j]=h*ff[j];
for(j=0;j<M;j++) yy[j]=y[j][i]+0.5*k2[j];
f(x+0.5*h,yy,ff);
for(j=0;j<M;j++) k3[j]=h*ff[j];
for(j=0;j<M;j++) yy[j]=y[j][i]+k3[j];
f(x+h,yy,ff);
for(j=0;j<M;j++) k4[j]=h*ff[j];
for(j=0;j<M;j++)
y[j][i+1]=y[j][i]+(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/6.0;
}
}
主函数
/* rkuttasm.c: Fourth-order Runge-kutta method for equations. */
#define N 22220
#define M 3
#include "rkuttas.c"
void f( double x, double yy[M], double ff[M] );
extern void rkuttas( double x0, double x1, double y[M][N] );
void main( )
{
int i,j;
double x0,x1,h,x,y0[M],y[M][N];
x0=0.0; x1=1;
y0[0]=-1.0; y[0][0]=y0[0];
y0[1]=0.0; y[1][0]=y0[1];
y0[2]=1.0; y[2][0]=y0[2];
system("cls");
printf("Fourth-order Runge-Kutta method for equations:\n");
rkuttas(x0,x1,y);
h=(x1-x0)/(N-1.0);
printf(" xn ");
for(j=0;j<M;j++) printf("y[%d] ",j+1);
printf("\n");
for(i=0;i<N;i++)
{
x=x0+i*h;
printf("%10.6f",x);
for(j=0;j<M;j++) printf("%12.8f",y[j][i]);
printf("\n");
}
getch();
}
/********* function f( x, yy[M], ff[M] ) ********/
void f( double x, double yy[M], double ff[M] )
{
ff[0]=yy[1];
ff[1]=-yy[0];
ff[2]=-yy[2];
}
