ch7.h

将C语言的常用程序集移植到VC开发环境的源代码

                for (j=0; j<=n-1; j++)
                  { qq=c[j]-9.0*(b[j]-2.0*a[j]/9.0)/8.0;
                    qq=d[j]-8.0*qq;
                    y[j]=y[j]+hh*qq/6.0;
                  }
                t=t+dt;
              }
            p=0.0;
            for (j=0; j<=n-1; j++)
              { qq=fabs(y[j]-v[j]);
                if (qq>p) p=qq;
              }
            hh=hh/2.0; nn=nn+nn;
          }
        for (j=0; j<=n-1; j++) z[j*k+i]=y[j];
      }
    free(a); free(b); free(c); free(d); free(u); free(v);
    return;
}
/////////////////////////////////////////////////////////////
static void rkt1(double t,double h,int n,double y[],double b[],double d[],double e[],
				 void (*gpbs1f)(double t,double y[],int n,double d[]))
{ 
	//extern void gpbs1f();
    int i,k;
    double a[4],tt;
    a[0]=h/2.0; a[1]=a[0]; a[2]=h; a[3]=h;
    gpbs1f(t,y,n,d);
    for (i=0; i<=n-1; i++) { b[i]=y[i]; e[i]=y[i];}
    for (k=0; k<=2; k++)
      { for (i=0; i<=n-1; i++)
          { y[i]=e[i]+a[k]*d[i];
            b[i]=b[i]+a[k+1]*d[i]/3.0;
          }
        tt=t+a[k];
        gpbs1f(tt,y,n,d);
      }
    for (i=0; i<=n-1; i++)
      y[i]=b[i]+h*d[i]/6.0;
    return;
}

void gpbs1(double t,double h,int n,double y[],double eps,
		void (*gpbs1f)(double t,double y[],int n,double d[]))
{ 
	int i,j,k,m,nn,it;
    //void rkt1();
    double x,hh,dd,q,p,g[10],*b,*d,*u,*v,*w,*e;
    b=(double*)malloc(10*n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    w=(double*)malloc(n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    for (j=0; j<=n-1; j++) v[j]=y[j];
    x=t; nn=1; hh=h; g[0]=hh;
    rkt1(x,hh,n,y,w,d,e,gpbs1f);
    for (j=0; j<=n-1; j++)
      { b[j]=y[j]; u[j]=y[j];}
    k=1; it=1;
    while (it==1)
      { nn=nn+nn; hh=hh/2.0; it=0;
        g[k]=hh;
        for (j=0; j<=n-1; j++) y[j]=v[j];
        t=x;
        for (j=0; j<=nn-1; j++)
          { rkt1(t,hh,n,y,w,d,e,gpbs1f);
            t=t+hh;
          }
        for (j=0; j<=n-1; j++)
          { dd=y[j]; m=0;
            for (i=0; i<=k-1; i++)
              if (m==0)
                { q=dd-b[i*n+j];
                  if (fabs(q)+1.0==1.0) m=1;
                  else dd=(g[k]-g[i])/q;
                }
            b[k*n+j]=dd;
            if (m!=0) b[k*n+j]=1.0e+35;
          }
        for (j=0; j<=n-1; j++)
          { dd=0.0;
            for (i=k-1; i>=0; i--)
              dd=-g[i]/(b[(i+1)*n+j]+dd);
            y[j]=dd+b[j];
          }
        p=0.0;
        for (j=0; j<=n-1; j++)
          { q=fabs(y[j]-u[j]);
            if (q>p) p=q;
          }
        if ((p>=eps)&&(k<7))
          { for (j=0; j<=n-1; j++) u[j]=y[j];
            k=k+1; it=1;
          }
      }
    free(b); free(d); free(u); free(v); free(w); free(e);
    return;
}
/////////////////////////////////////////////////////////////
static void rkt2(double t,double h,int n,double y[],double b[],double d[],double e[],
				 void (*gpbs2f)(double t,double y[],int n,double d[]))
{ 
	//extern void gpbs2f();
    int i,k;
    double a[4],tt;
    a[0]=h/2.0; a[1]=a[0]; a[2]=h; a[3]=h;
    gpbs2f(t,y,n,d);
    for (i=0; i<=n-1; i++) { b[i]=y[i]; e[i]=y[i];}
    for (k=0; k<=2; k++)
      { for (i=0; i<=n-1; i++)
          { y[i]=e[i]+a[k]*d[i];
            b[i]=b[i]+a[k+1]*d[i]/3.0;
          }
        tt=t+a[k];
        gpbs2f(tt,y,n,d);
      }
    for (i=0; i<=n-1; i++)
      y[i]=b[i]+h*d[i]/6.0;
    return;
}
void gpbs2(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*gpbs2f)(double t,double y[],int n,double d[]))
{ 
	int i,j,kk,l,m,nn,it;
    //void rkt2();
    double x,hh,dd,tt,q,p,g[10],*b,*d,*u,*v,*w,*e;
    b=(double*)malloc(10*n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    w=(double*)malloc(n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    for (i=1; i<=k-1; i++)
    { 
		for (j=0; j<=n-1; j++) v[j]=y[j];
        x=t+(i-1)*h; nn=1; hh=h; g[0]=hh;
        rkt2(x,hh,n,y,w,d,e,gpbs2f);
        for (j=0; j<=n-1; j++)
          { b[j]=y[j]; u[j]=y[j];}
        kk=1; it=1;
        while (it==1)
          { nn=nn+nn; hh=hh/2.0; it=0;
            g[kk]=hh;
            for (j=0; j<=n-1; j++) y[j]=v[j];
            tt=x;
            for (j=0; j<=nn-1; j++)
              { rkt2(t,hh,n,y,w,d,e,gpbs2f);
                tt=tt+hh;
              }
            for (j=0; j<=n-1; j++)
              { dd=y[j]; l=0;
                for (m=0; m<=kk-1; m++)
                  if (l==0)
                    { q=dd-b[m*n+j];
                      if (fabs(q)+1.0==1.0) l=1;
                      else dd=(g[kk]-g[m])/q;
                    }
                b[kk*n+j]=dd;
                if (l!=0) b[kk*n+j]=1.0e+35;
              }
            for (j=0; j<=n-1; j++)
              { dd=0.0;
                for (m=kk-1; m>=0; m--)
                  dd=-g[m]/(b[(m+1)*n+j]+dd);
                y[j]=dd+b[j];
              }
            p=0.0;
            for (j=0; j<=n-1; j++)
              { q=fabs(y[j]-u[j]);
                if (q>p) p=q;
              }
            if ((p>=eps)&&(kk<7))
              { for (j=0; j<=n-1; j++) u[j]=y[j];
                kk=kk+1; it=1;
              }
          }
        for (j=0; j<=n-1; j++)
          z[j*k+i]=y[j];
	}
    free(b); free(d); free(u); free(v); free(w); free(e);
    return;
}
/////////////////////////////////////////////////////////////
static void rkt3(double t,double h,double y[],int n,double eps,
				 void (*ggjfqf)(double t,double y[],int n,double d[]))
{
	//extern void ggjfqf();
    int m,i,j,k;
    double hh,p,dt,x,tt,q,a[4],*g,*b,*c,*d,*e;
    g=(double*)malloc(n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    c=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    hh=h; m=1; p=1.0+eps; x=t;
    for (i=0; i<=n-1; i++) c[i]=y[i];
    while (p>=eps)
      { a[0]=hh/2.0; a[1]=a[0]; a[2]=hh; a[3]=hh;
        for (i=0; i<=n-1; i++)
          { g[i]=y[i]; y[i]=c[i];}
        dt=h/m; t=x;
        for (j=0; j<=m-1; j++)
          { ggjfqf(t,y,n,d);
            for (i=0; i<=n-1; i++) 
              { b[i]=y[i]; e[i]=y[i];}
            for (k=0; k<=2; k++)
              { for (i=0; i<=n-1; i++)
                  { y[i]=e[i]+a[k]*d[i];
                    b[i]=b[i]+a[k+1]*d[i]/3.0;
                  }
                tt=t+a[k];
                ggjfqf(tt,y,n,d);
              }
            for (i=0; i<=n-1; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
          }
        p=0.0;
        for (i=0; i<=n-1; i++)
          { q=fabs(y[i]-g[i]);
            if (q>p) p=q;
          }
        hh=hh/2.0; m=m+m;
      }
    free(g); free(b); free(c); free(d); free(e);
    return;
}
void ggjfq(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*ggjfqf)(double t,double y[],int n,double d[]))
{ 
	//extern void ggjfqf();
    //void rkt3();
    int i,j;
    double a,qq,*d,*p,*u,*v,*w;
    d=(double*)malloc(n*sizeof(double));
    p=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    w=(double*)malloc(n*sizeof(double));
    for (i=0; i<=n-1; i++) {p[i]=0.0; z[i*k]=y[i];}
    a=t;
    ggjfqf(t,y,n,d);
    for (j=0; j<=n-1; j++) u[j]=d[j];
    rkt3(t,h,y,n,eps,ggjfqf);
    t=a+h;
    ggjfqf(t,y,n,d);
    for (j=0; j<=n-1; j++)
      { z[j*k+1]=y[j]; v[j]=d[j];}
    for (j=0; j<=n-1; j++)
      { p[j]=-4.0*z[j*k+1]+5.0*z[j*k]+2.0*h*(2.0*v[j]+u[j]);
        y[j]=p[j];
      }
    t=a+2.0*h;
    ggjfqf(t,y,n,d);
    for (j=0; j<=n-1; j++)
      { qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
        qq=qq+4.0*z[j*k+1]-3.0*z[j*k];
        z[j*k+2]=(p[j]+qq)/2.0;
        y[j]=z[j*k+2];
      }
    for (i=3; i<=k-1; i++)
      { t=a+(i-1)*h;
        ggjfqf(t,y,n,d);
        for (j=0; j<=n-1; j++)
          { u[j]=v[j]; v[j]=d[j];}
        for (j=0; j<=n-1; j++)
          { qq=-4.0*z[j*k+i-1]+5.0*z[j*k+i-2];
            p[j]=qq+2.0*h*(2.0*v[j]+u[j]);
            y[j]=p[j];
          }
        t=t+h;
        ggjfqf(t,y,n,d);
        for (j=0; j<=n-1; j++)
          { qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
            qq=qq+4.0*z[j*k+i-1]-3.0*z[j*k+i-2];
            y[j]=(p[j]+qq)/2.0;
            z[j*k+i]=y[j];
          }
      }
    free(d); free(p); free(u); free(v); free(w);
    return;
}
/////////////////////////////////////////////////////////////
static void rkt4(double t,double h,double y[],int n,double eps,
				 void (*gadmsf)(double t,double y[],int n,double d[]))
{ 
	//extern void gadmsf();
    int m,i,j,k;
    double hh,p,dt,x,tt,q,a[4],*g,*b,*c,*d,*e;
    g=(double*)malloc(n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    c=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    hh=h; m=1; p=1.0+eps; x=t;
    for (i=0; i<=n-1; i++) c[i]=y[i];
    while (p>=eps)
      { a[0]=hh/2.0; a[1]=a[0]; a[2]=hh; a[3]=hh;
        for (i=0; i<=n-1; i++)
          { g[i]=y[i]; y[i]=c[i];}
        dt=h/m; t=x;
        for (j=0; j<=m-1; j++)
          { gadmsf(t,y,n,d);
            for (i=0; i<=n-1; i++) 
              { b[i]=y[i]; e[i]=y[i];}
            for (k=0; k<=2; k++)
              { for (i=0; i<=n-1; i++)
                  { y[i]=e[i]+a[k]*d[i];
                    b[i]=b[i]+a[k+1]*d[i]/3.0;
                  }
                tt=t+a[k];
                gadmsf(tt,y,n,d);
              }
            for (i=0; i<=n-1; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
          }
        p=0.0;
        for (i=0; i<=n-1; i++)
          { q=fabs(y[i]-g[i]);
            if (q>p) p=q;
          }
        hh=hh/2.0; m=m+m;
      }
    free(g); free(b); free(c); free(d); free(e);
    return;
}
void gadms(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*gadmsf)(double t,double y[],int n,double d[]))
{ 
	//extern void gadmsf();
    //void rkt4();
    int i,j,m;
    double a,q,*b,*e,*s,*g,*d;
    b=(double*)malloc(4*n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    s=(double*)malloc(n*sizeof(double));
    g=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    a=t;
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    gadmsf(t,y,n,d);
    for (i=0; i<=n-1; i++) b[i]=d[i];
    for (i=1; i<=3; i++)
      if (i<=k-1)
        { t=a+i*h;
          rkt4(t,h,y,n,eps,gadmsf);
          for (j=0; j<=n-1; j++) z[j*k+i]=y[j];
          gadmsf(t,y,n,d);
          for (j=0; j<=n-1; j++) b[i*n+j]=d[j];
        }
    for (i=4; i<=k-1; i++)
      { for (j=0; j<=n-1; j++)
          { q=55.0*b[3*n+j]-59.0*b[2*n+j];
            q=q+37.0*b[n+j]-9.0*b[j];
            y[j]=z[j*k+i-1]+h*q/24.0;
            b[j]=b[n+j];
            b[n+j]=b[n+n+j];
            b[n+n+j]=b[n+n+n+j];
          }
        t=a+i*h;
        gadmsf(t,y,n,d);
        for (m=0; m<=n-1; m++) b[n+n+n+m]=d[m];
        for (j=0; j<=n-1; j++)
          { q=9.0*b[3*n+j]+19.0*b[n+n+j]-5.0*b[n+j]+b[j];
            y[j]=z[j*k+i-1]+h*q/24.0;
            z[j*k+i]=y[j];
          }
        gadmsf(t,y,n,d);
        for (m=0; m<=n-1; m++) b[3*n+m]=d[m];
      }
    free(b); free(e); free(s); free(g); free(d);
    return;
}
/////////////////////////////////////////////////////////////
static void rkt5(double t,double h,double y[],int n,double eps,
				 void (*ghamgf)(double t,double y[],int n,double d[]))
{ 
	//extern void ghamgf();