vsrk4.c

数值方法求解常微分方程的C实现[TEXT为主程序,其余方法(即C文件)放入Include文件夹即可]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
extern double f( double x, double y );
extern double Y_initial;
extern int N;
void vsrk4(double x_min, double x_max, double y[])
{
double rk4( double x, double y, double hh );
int i,k,j;
double xx,h,yn,ynh,delta;
double x[51];
static double hmin=0.0005; 
static double hmax=0.1;
static double eps=1.e-8;
h=hmax;
xx=x_min;
i=0;
x[0]=xx;
y[0]=Y_initial;
while( xx < x_max )
  {
again:  if( xx+h > x_max ) h=x_max-xx;
  yn=rk4(xx,y[i],h);
  ynh=rk4(xx,y[i],h/2.0);
  ynh=rk4(xx+h/2.0,ynh,h/2.0);
  delta=fabs((yn-ynh)/15.0);
  if( ( delta < eps )||( h == hmin ) )
    {
    xx=xx+h;
    i=i+1;
    x[i]=xx;
    y[i]=ynh;
    if( i >=50 ) goto endd;
    if( delta < 0.05*eps )
      {
      h=2.0*h;
      if( h > hmax ) h=hmax;
      }
    }
  else
    {
    h=h/2.0;
    if( h < hmin ) h=hmin;
    goto again;
    }
  }
endd:k=i;
printf("     xn           yn\n");
for(j=0;j<=k;j++)
   printf("%12.8f %12.8f\n",x[j],y[j]);
getch();
}