全区间积分的双边法(常微分方程组的求解)

/*
代码作者：不详
代码整理者：设计天下 MySDN网站 算法天下工作室

功能：全区间积分的双边法(常微分方程组的求解)
*/

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
/*全区间积分的定步长欧拉方法*/
/*全区间积分的维梯方法*/
/*全区间积分的定步长龙格-库塔方法*/
typedef struct _fode {
 int        n;      /*微分方程组中方程个数，也是未知函数的个数*/
 int        steps;  /*积分步数（包括起始点这一步）*/
 double     lens;   /*积分的步长*/
 double     t;      /*对微分方程进行积分的起始点 */
 double    *y;      /*存放n个未知函数在起始点t处的函数值*/
 double    *z;      /*返回steps个积分点（包括起始点）上的未知函数值*/
 void     (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名（由用户自编）*/
} FODE, *FODEP;

/*全区间积分的变步长默森方法*/
/*全区间积分的双边法*/
/*全区间积分的阿当姆斯预报校正法*/
/*全区间积分的哈明方法*/
typedef struct _fode2 {
 int        n;      /*微分方程组中方程个数，也是未知函数的个数*/
 int        steps;  /*积分步数（包括起始点这一步）*/
 double     lens;   /*积分的步长*/
 double     eps;    /*控制精度要求*/
 double     t;      /*对微分方程进行积分的起始点 */
 double    *y;      /*存放n个未知函数在起始点t处的函数值*/
 double    *z;      /*返回steps个积分点（包括起始点）上的未知函数值*/
 void     (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名（由用户自编）*/
} FODE2, *FODE2P;

/*积分一步的特雷纳方法*/
typedef struct _tlode {
 int        n;      /*微分方程组中方程个数，也是未知函数的个数*/
 double     lens;   /*积分的步长*/
 double     t;      /*对微分方程进行积分的起始点 */
 double    *y;      /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
 void     (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名（由用户自编）*/
} TODE, *TODEP;

/*积分一步的变步长欧拉方法*/
/*积分一步的变步长龙格-库塔方法*/
/*积分一步的连分式法*/
typedef struct _eode {
 int        n;      /*微分方程组中方程个数，也是未知函数的个数*/
 double     lens;   /*积分的步长*/
 double     eps;    /*积分的精度要求*/
 double     t;      /*对微分方程进行积分的起始点 */
 double    *y;      /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
 void     (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名（由用户自编）*/
} EODE, *EODEP;

/*积分一步的变步长基尔方法*/
typedef struct _eode2 {
 int        n;      /*微分方程组中方程个数，也是未知函数的个数*/
 double     lens;   /*积分的步长*/
 double     eps;    /*积分的精度要求*/
 double     t;      /*对微分方程进行积分的起始点 */
 double    *y;      /*存放n个未知函数在起始点t处的函数值。返回t+lens点处的n个未知函数值*/
 double    *q;      /*在主函数第一次调用本函数时，应赋值0，
        以后每调用一次函数（即每积分一步），将由本函数的返回值以便循环使用。*/
 void     (*ptr)(); /*指向计算微分方程组中各方程右端函数值的函数名（由用户自编）*/
} EODE2, *EODE2P;

/*二阶微分方程边值问题的数值解法*/
typedef struct _dode {
 int        n;      /*求解区间[a,b]的等分点数(包括左端点a与右端点b)*/
 double     a;      /*求解区间的左端点*/
 double     b;      /*求解区间的右端点求*/
 double     ya;     /*未知函数在求解区间左端点处的函数值y(a)*/
 double     yb;     /*未知函数在求解区间右端点处的函数值y(b)*/
 double    *y;      /*返回n个等距离散点上的未知函数值*/
 void     (*ptr)(); /*指向计算二阶微分方程中函数值的函数名(由用户自编)*/
} DODE, *DODEP;

static void bilateral_method_lis(int n,double eps,double t,double lens,double y[],void (*ptr)())
{
 int    m,i,j,k;
 double hh,p,dt,x,tt,q,a[4];
 double *g,*b,*c,*d,*e;

 g=malloc(n*sizeof(double));
 b=malloc(n*sizeof(double));
 c=malloc(n*sizeof(double));
 d=malloc(n*sizeof(double));
 e=malloc(n*sizeof(double));
 hh=lens;
 m=1;
 x=t;
 for (i=0; i<n; i++)
 {
  c[i]=y[i];
 }

 do
 {
  a[1]=a[0]=hh/2.0;
  a[3]=a[2]=hh;
  for (i=0; i<n; i++)
  {
   g[i]=y[i];
   y[i]=c[i];
  }
  dt=lens/m; t=x;
  for (j=0; j<m; j++)
  {
   (*ptr)(t,y,n,d);
   for (i=0; i<n; i++)
   {
    b[i]=y[i]; e[i]=y[i];
   }
   for (k=0; k<3; k++)
   {
    for (i=0; i<n; i++)
    {
     y[i]=e[i]+a[k]*d[i];
     b[i]+=a[k+1]*d[i]/3.0;
    }
    tt=t+a[k];
    (*ptr)(tt,y,n,d);
   }
   for (i=0; i<n; i++)
   {
    y[i]=b[i]+hh*d[i]/6.0;
   }
   t+=dt;
  }
  p=0.0;
  for (i=0; i<n; i++)
  {
   q=fabs(y[i]-g[i]);
   if (q>p)
   {
    p=q;
   }
  }
  hh=hh/2.0;
  m<<=1;
 }
 while (p>=eps);
 /*释放动态分配的内存*/
 free(g); free(b); free(c); free(d); free(e);
 return;
}

void bilateral_lis(FODE2P ap)
{
 int    n,steps,i,j;
 double lens,eps,t,t0,qq;
 double *d,*p,*u,*v,*w,*y,*z;

 n=ap->n;
 steps=ap->steps;
 d=malloc(n*sizeof(double));
 p=malloc(n*sizeof(double));
 u=malloc(n*sizeof(double));
 v=malloc(n*sizeof(double));
 w=malloc(n*sizeof(double));
 lens=ap->lens;
 eps=ap->eps;
 y=ap->y;
 z=ap->z;
 for (i=0; i<n; i++)
 {
  p[i]=0.0;
  z[i*steps]=y[i];
 }
 t0=t=ap->t;
 (*ap->ptr)(t,y,n,d);
 for (j=0; j<n; j++)
 {
  u[j]=d[j];
 }
 bilateral_method_lis(n,eps,t,lens,y,ap->ptr);
 t+=lens;
 (*ap->ptr)(t,y,n,d);
 for (j=0; j<n; j++)
 {
  z[j*steps+1]=y[j];
  v[j]=d[j];
 }
 for (j=0; j<n; j++)
 {
  p[j]=-4.0*z[j*steps+1]+5.0*z[j*steps]+2.0*lens*(2.0*v[j]+u[j]);
  y[j]=p[j];
 }
 t+=lens;
 (*ap->ptr)(t,y,n,d);
 for (j=0; j<n; j++)
 {
  qq=2.0*lens*(d[j]-2.0*v[j]-2.0*u[j])/3.0
   +4.0*z[j*steps+1]-3.0*z[j*steps];
  z[j*steps+2]=(p[j]+qq)/2.0;
  y[j]=z[j*steps+2];
 }
 for (i=3; i<steps; i++)
 {
  t=t0+(i-1)*lens;
  (*ap->ptr)(t,y,n,d);
  for (j=0; j<=n-1; j++)
  {
   u[j]=v[j];
   v[j]=d[j];
  }
  for (j=0; j<n; j++)
  {
   qq=-4.0*z[j*steps+i-1]+5.0*z[j*steps+i-2];
   p[j]=qq+2.0*lens*(2.0*v[j]+u[j]);
   y[j]=p[j];
  }
  t+=lens;
  (*ap->ptr)(t,y,n,d);
  for (j=0; j<n; j++)
  {
   qq=2.0*lens*(d[j]-2.0*v[j]-2.0*u[j])/3.0
    +4.0*z[j*steps+i-1]-3.0*z[j*steps+i-2];
   y[j]=(p[j]+qq)/2.0;
   z[j*steps+i]=y[j];
  }
 }
 /*释放动态分配的内存*/
 free(d);
 free(p);
 free(u);
 free(v);
 free(w);
 return;
}

void bilateral_ptr_lis(double t,double y[],int n,double d[])
{
 t=t;
 n=n;
 d[0]=-2.0*y[1];
 d[1]=2.0*y[0];
 return;
}

int main()
{
 double y[2]={0.0, 1.0};
 double z[2][11]={0};
 FODE2  fa={2, 11, 0.1, 0.00001, 0.1, y, (double*)z, bilateral_ptr_lis};
 int    i;

 bilateral_lis(&fa);
 printf("/n");
 for (i=0; i<fa.steps; i++)
 {
  printf("t=%7.3f  y(0)=%e  y(1)=%e/n", i*fa.lens,z[0][i],z[1][i]);
 }
 printf("/n");
 return 0;
}