ok.c

本方程为二阶边值问题的数值算法之打靶法计算

/*********************************************************************
 *本方程为二阶边值问题的数值算法之打靶法计算
 
 *程序在Dev-C++ 4.9.8.0下编译成功 
 
 *本程序在windows xp pro 下运行成功
 **********************************************************************/
 
 
#include  "stdio.h"

#define  M  100

void init_step(void);
void ans_display(void);
void compute(int flag);

  double a=1;
  double h=0.1;
  double lambda=0.5,tau,T,N,J;
  double xj[400],u0j[400],u1j[400];
  int j,n;
  
int main()
{
  int a,b,N;
  double e,h,Tk;
  double z1,z2,bata_arfa,y1[400],y2[400],x[100];
  int k=1,i;
  double k11,k12,k21,k22,k31,k32,k41,k42,k011,k012,k021,k022,k031,k032,k041,k042;
  b=0;a=-1,e=0.05;
  bata_arfa=8;
  N=100;
  h=(double)(b-a)/N;
  /* printf("h= %lf\n",h);*/
  Tk=bata_arfa/(b-a);
      /*计算前半个方程的算法*/
      y1[0]=0;y2[0]=Tk;z1=0;z2=1;
      for(i=0;i<N+1;i++)
      {
         x[i]=a+i*h;
         k11=h*y2[i];
         k12=h*(-42*x[i]*x[i]*x[i]*x[i]*x[i]*x[i]+1);
         k21=h*(y2[i]+0.5*k12);
         k22=h*(-42*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)+1);
         k31=h*(y2[i]+0.5*k22);
         k32=h*(-42*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)+1);
         k41=h*(y2[i]+k32);
         k42=h*(-42*(x[i]+h)*(x[i]+h)*(x[i]+h)*(x[i]+h)*(x[i]+h)+1);
         y1[i+1]=y1[i]+(k11+2*k21+2*k31+k41)/6;
         y2[i+1]=y2[i]+(k12+2*k22+2*k32+k42)/6;

         k011=h*z2;
         k012=0;
         k021=h*(z2+0.5*k012);
         k022=0;
         k031=h*(z2+0.5*k022);
         k032=0;
         k041=h*(z2+k032);
         k042=0;
         z1=z1+(k011+2*k021+2*k031+k041)/6;
         z2=z2+(k012+2*k022+2*k032+k042)/6;
     }
           for(i=0;i<N+1;i++)
           printf("%f,\t %f\n",x[i],y1[i]);

      b=0;a=1;
      bata_arfa=-9;
      h=(double)(b-a)/N;
      Tk=bata_arfa/(b-a);

	  /*计算后半个方程的算法*/
      y1[0]=0;y2[0]=Tk;z1=0;z2=1;
      for(i=0;i<N;i++)
      {
         x[i]=a+i*h;
         k11=h*y2[i];
         k12=h*(42*x[i]*x[i]*x[i]*x[i]*x[i]*x[i]+2);
         k21=h*(y2[i]+0.5*k12);
         k22=h*(42*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)+2);
         k31=h*(y2[i]+0.5*k22);
         k32=h*(42*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)*(x[i]+h/2)+2);
         k41=h*(y2[i]+k32);
         k42=h*(42*(x[i]+h)*(x[i]+h)*(x[i]+h)*(x[i]+h)*(x[i]+h)+2);
         y1[i+1]=y1[i]+(k11+2*k21+2*k31+k41)/6;
         y2[i+1]=y2[i]+(k12+2*k22+2*k32+k42)/6;

         k011=h*z2;
         k012=0;
         k021=h*(z2+0.5*k012);
         k022=0;
         k031=h*(z2+0.5*k022);
         k032=0;
         k041=h*(z2+k032);
         k042=0;
         z1=z1+(k011+2*k021+2*k031+k041)/6;
         z2=z2+(k012+2*k022+2*k032+k042)/6;
        }
           for(i=N-1;i>=0;i--)
           printf("%f,\t %f\n",x[i],y1[i]);

  getchar();
}