/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH7.h"  常微分方程(組)的求解*****/
#ifndef CH7_H_
#define CH7_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

//#include "ch2.h"
#include "ch1.h"
//*******************************************************************
//用戶自編函數，在具體應用中必須予以確定 (給函數指針(全局變量)賦值)，函數定義參見文件末
//*******************************************************************
void gelr1(double t,double y[],int n,double h,int k,double z[],
		void (*gelr1f)(double t,double y[],int n,double d[]));//全區間積分的定步長歐拉方法
void gelr2(double t,double h,double y[],int n,double eps,
		void (*gelr2f)(double t,double y[],int n,double d[]));//積分一步的變步長歐拉方法
void gwity(double t,double y[],int n,double h,int k,double z[],
		void (*gwityf)(double t,double y[],int n,double d[]));//全區間積分的定步長維梯方法
void grkt1(double t,double y[],int n,double h,int k,double z[],
		void (*grkt1f)(double t,double y[],int n,double d[]));//全區間積分的定步長龍格-庫塔法
void grkt2(double t,double h,double y[],int n,double eps,
		void (*grkt2f)(double t,double y[],int n,double d[]));//積分一步的變步長龍格-庫塔法
void ggil1(double t,double h,double y[],int n,double eps,double q[],
		void (*ggil1f)(double t,double y[],int n,double d[]));//積分一步的變步長基爾方法
void ggil2(double t,double h,double y[],int n,double eps,int k,double z[],
		void (*ggil2f)(double t,double y[],int n,double d[]));//全區間積分的變步長基爾方法
void gmrsn(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*gmrsnf)(double t,double y[],int n,double d[]));//全區間積分的變步長默森方法
void gpbs1(double t,double h,int n,double y[],double eps,
		void (*gpbs1f)(double t,double y[],int n,double d[]));//積分一步的連分式法
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[]));//全區間積分的連分式法
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[]));//全區間積分的雙邊法
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[]));//全區間積分的阿當姆斯預報校正法
void ghamg(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*ghamgf)(double t,double y[],int n,double d[]));//全區間積分的哈明方法
void gtnr1(double t,double h,int n,double y[],
		void (*gtnr1f)(double t,double y[],int n,double d[]));//積分一步的特雷納方法
void gtnr2(double t,double h,int n,double y[],int k,double z[],
		void (*gtnr2f)(double t,double y[],int n,double d[]));//全區間積分的特雷納方法
		
int ggear(double a,double b,double hmin,double hmax,double h,double eps,
		int n,double y0[],int k,double t[],double z[],
		void (*ggearf)(double t,double y[],int n,double d[]),
		void (*ggears)(double t,double y[],int n,double p[]));//積分剛性方程組的吉爾方法
		/* int n; double t,y[],p[3][3]*/
void gdfte(double a,double b,double ya,double yb,int n,double y[],
		void (*gdftef)(double x,double z[]));//二階微分方程邊值問題的數值數值解法

//*******************************************************************
void gelr1(double t,double y[],int n,double h,int k,double z[],
		void (*gelr1f)(double t,double y[],int n,double d[]))
{ 
	//extern void gelr1f();
    int i,j;
    double x,*d;
    d=(double*)malloc(n*sizeof(double));
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    for (j=1; j<=k-1; j++)
      { x=t+(j-1)*h;
        gelr1f(x,y,n,d);
        for (i=0; i<=n-1; i++)
          y[i]=z[i*k+j-1]+h*d[i];
        x=t+j*h;
        gelr1f(x,y,n,d);
        for (i=0; i<=n-1; i++)
          d[i]=z[i*k+j-1]+h*d[i];
        for (i=0; i<=n-1; i++)
          { y[i]=(y[i]+d[i])/2.0;
            z[i*k+j]=y[i];
          }
      }
    free(d); return;
}
/////////////////////////////////////////////////////////////
void gelr2(double t,double h,double y[],int n,double eps,
		void (*gelr2f)(double t,double y[],int n,double d[]))
{ 
	//extern void gelr2f();
    int i,j,m;
    double hh,p,x,q,*a,*b,*c,*d;
    a=(double*)malloc(n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    c=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    hh=h; m=1; p=1.0+eps;
    for (i=0; i<=n-1; i++) a[i]=y[i];
    while (p>=eps)
      { for (i=0; i<=n-1; i++)
          { b[i]=y[i]; y[i]=a[i];}
        for (j=0; j<=m-1; j++)
          { for (i=0; i<=n-1; i++) c[i]=y[i];
            x=t+j*hh;
            gelr2f(x,y,n,d);
            for (i=0; i<=n-1; i++)
              y[i]=c[i]+hh*d[i];
            x=t+(j+1)*hh;
            gelr2f(x,y,n,d);
            for (i=0; i<=n-1; i++)
              d[i]=c[i]+hh*d[i];
            for (i=0; i<=n-1; i++)
              y[i]=(y[i]+d[i])/2.0;
          }
        p=0.0;
        for (i=0; i<=n-1; i++)
          { q=fabs(y[i]-b[i]);
            if (q>p) p=q;
          }
        hh=hh/2.0; m=m+m;
      }
    free(a); free(b); free(c); free(d);
    return;
}
/////////////////////////////////////////////////////////////
void gwity(double t,double y[],int n,double h,int k,double z[],
		void (*gwityf)(double t,double y[],int n,double d[]))
{ 
	//extern void gwityf();
    int i,j;
    double x,*a,*d;
    a=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    gwityf(t,y,n,d);
    for (j=1; j<=k-1; j++)
      { for (i=0; i<=n-1; i++)
          a[i]=z[i*k+j-1]+h*d[i]/2.0;
        x=t+(j-0.5)*h;
        gwityf(x,a,n,y);
        for (i=0; i<=n-1; i++)
          { d[i]=2.0*y[i]-d[i];
            z[i*k+j]=z[i*k+j-1]+h*y[i];
          }
      }
    free(a); free(d);
    return;
}
/////////////////////////////////////////////////////////////
void grkt1(double t,double y[],int n,double h,int k,double z[],
		void (*grkt1f)(double t,double y[],int n,double d[]))
{ 
	//extern void grkt1f();
    int i,j,l;
    double a[4],tt,*b,*d;
    b=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    a[0]=h/2.0; a[1]=a[0];
    a[2]=h; a[3]=h;
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    for (l=1; l<=k-1; l++)
      { grkt1f(t,y,n,d);
        for (i=0; i<=n-1; i++) b[i]=y[i];
        for (j=0; j<=2; j++)
          { for (i=0; i<=n-1; i++)
              { y[i]=z[i*k+l-1]+a[j]*d[i];
                b[i]=b[i]+a[j+1]*d[i]/3.0;
              }
            tt=t+a[j];
            grkt1f(tt,y,n,d);
          }
        for (i=0; i<=n-1; i++)
          y[i]=b[i]+h*d[i]/6.0;
        for (i=0; i<=n-1; i++)
          z[i*k+l]=y[i];
        t=t+h;
      }
    free(b); free(d);
    return;
}
/////////////////////////////////////////////////////////////
void grkt2(double t,double h,double y[],int n,double eps,
		void (*grkt2f)(double t,double y[],int n,double d[]))
{ 
	//extern void grkt2f();
    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++)
          { grkt2f(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];
                grkt2f(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 ggil1(double t,double h,double y[],int n,double eps,double q[],
		void (*ggil1f)(double t,double y[],int n,double d[]))
{
	//extern void ggil1f();
    int i,j,k,m,ii;
    double x,p,hh,r,s,t0,dt,qq,*d,*u,*v,*g;
    static double a[4]={0.5,0.29289321881,
                        1.7071067812,0.166666667};
    static double b[4]={2.0,1.0,1.0,2.0};
    static double c[4],e[4]={0.5,0.5,1.0,1.0};
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    g=(double*)malloc(n*sizeof(double));
    for (i=0; i<=2; i++) c[i]=a[i];
    c[3]=0.5;
    x=t; p=1.0+eps; hh=h; m=1;
    for (j=0; j<=n-1; j++) u[j]=y[j];
    while (p>=eps)
      { for (j=0; j<=n-1; j++)
          { v[j]=y[j]; y[j]=u[j]; g[j]=q[j];}
        dt=h/m; t=x;
        for (k=0; k<=m-1; k++)
          { ggil1f(t,y,n,d);
            for (ii=0; ii<=3; ii++)
              { for (j=0; j<=n-1; j++)
                  d[j]=d[j]*hh;
                for (j=0; j<=n-1; j++)
                  { r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
                    y[j]=y[j]+r;
                    s=g[j]+3.0*r;
                    g[j]=s-c[ii]*d[j];
                  }
                t0=t+e[ii]*hh;
                ggil1f(t0,y,n,d);
              }
            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; m=m+m;
      }
    for (j=0; j<=n-1; j++) q[j]=g[j];
    free(g); free(d); free(u); free(v);
    return;
}
/////////////////////////////////////////////////////////////
void ggil2(double t,double h,double y[],int n,double eps,int k,double z[],
		void (*ggil2f)(double t,double y[],int n,double d[]))
{ 
	//extern void ggil2f();
    int i,j,m,kk,ii;
    double aa,hh,x,p,dt,r,s,t0,qq,*g,*q,*d,*u,*v;
    static double a[4]={0.5,0.29289321881,
                        1.7071067812,0.166666667};
    static double b[4]={2.0,1.0,1.0,2.0};
    static double c[4],e[4]={0.5,0.5,1.0,1.0};
    q=(double*)malloc(n*sizeof(double));
    g=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    for (i=0; i<=2; i++) c[i]=a[i];
    c[3]=0.5;
    aa=t;
    for (i=0; i<=n-1; i++) 
      { z[i*k]=y[i]; q[i]=0.0;}
    for (i=2; i<=k; i++)
      { x=aa+(i-2)*h; m=1; hh=h;
        p=1.0+eps;
        for (j=0; j<=n-1; j++) u[j]=y[j];
        while (p>=eps)
          { for (j=0; j<=n-1; j++)
              { v[j]=y[j]; y[j]=u[j]; g[j]=q[j];}
            dt=h/m; t=x;
            for (kk=0; kk<=m-1; kk++)
              { ggil2f(t,y,n,d);
                for (ii=0; ii<=3; ii++)
                  { for (j=0; j<=n-1; j++)
                      d[j]=d[j]*hh;
                    for (j=0; j<=n-1; j++)
                      { r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
                        y[j]=y[j]+r;
                        s=g[j]+3.0*r;
                        g[j]=s-c[ii]*d[j];
                      }
                    t0=t+e[ii]*hh;
                    ggil2f(t0,y,n,d);
                  }
                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; m=m+m;
          }
        for (j=0; j<=n-1; j++)
          { q[j]=g[j]; z[j*k+i-1]=y[j];}
      }
    free(q); free(g); free(d); free(u); free(v);
    return;
}
/////////////////////////////////////////////////////////////
void gmrsn(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*gmrsnf)(double t,double y[],int n,double d[]))
{ 
	//extern void gmrsnf();
    int i,j,m,nn;
    double aa,bb,x,hh,p,dt,t0,qq,*a,*b,*c,*d,*u,*v;
    a=(double*)malloc(n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    c=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    aa=t;
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    for (i=1; i<=k-1; i++)
      { x=aa+(i-1)*h; nn=1; hh=h;
        for (j=0; j<=n-1; j++) u[j]=y[j];
        p=1.0+eps;
        while (p>=eps)
          { for (j=0; j<=n-1; j++)
              { v[j]=y[j]; y[j]=u[j];}
            dt=h/nn; t=x;
            for (m=0; m<=nn-1; m++)
              { gmrsnf(t,y,n,d);
                for (j=0; j<=n-1; j++)
                  { a[j]=d[j]; y[j]=y[j]+hh*d[j]/3.0;}
                t0=t+hh/3.0;
                gmrsnf(t0,y,n,d);
                for (j=0; j<=n-1; j++)
                  { b[j]=d[j]; y[j]=y[j]+hh*(d[j]-a[j])/6.0;}
                gmrsnf(t0,y,n,d);
                for (j=0; j<=n-1; j++)
                  { b[j]=d[j];
                    bb=(d[j]-4.0*(b[j]+a[j]/4.0)/9.0)/8.0;
                    y[j]=y[j]+3.0*hh*bb;
                  }
                t0=t+hh/2.0;
                gmrsnf(t0,y,n,d);
                for (j=0; j<=n-1; j++)
                  { c[j]=d[j];
                    qq=d[j]-15.0*(b[j]-a[j]/5.0)/16.0;
                    y[j]=y[j]+2.0*hh*qq;
                  }
                t0=t+hh;
                gmrsnf(t0,y,n,d);