全区间积分的定步长欧拉方法(常微分方程组的求解)

/*
代码作者：不详
代码整理者：设计天下 MySDN网站 算法天下工作室 LS
代号: _MySDN_SFTX_LS
网址：http://www.my-sdn.net

功能：全区间积分的定步长欧拉方法(常微分方程组的求解)

*/

#include "stdio.h"
#include "stdlib.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;

void euler_step_lis(FODEP ap)
{
 int    n,steps,i,j;
 double x,lens,t,*y,*z,*d;

 n=ap->n;
 steps=ap->steps;
 lens=ap->lens;
 t=ap->t;
 y=ap->y;
 z=ap->z;
 d=malloc(n*sizeof(double));
 for (i=0; i<n; i++)
 {
  z[i*steps]=y[i];
 }
 for (j=1; j<steps; j++)
 {
  x=t+(j-1)*lens;
  (*ap->ptr)(y,d);
  for (i=0; i<n; i++)
  {
   y[i]=z[i*steps+j-1]+lens*d[i];
  }
  x=t+j*lens;
  (*ap->ptr)(y,d);
  for (i=0; i<n; i++)
  {
   d[i]=z[i*steps+j-1]+lens*d[i];
  }
  for (i=0; i<n; i++)  
  {
   y[i]=0.5*(y[i]+d[i]);
   z[i*steps+j]=y[i];
  }
 }
 free(d);/*释放动态分配的内存*/
}
 
void step_ptr_lis(double y[],double d[])
{
 d[0]=2.0*y[1];
 d[1]=-2.0*y[0];
 d[2]=-2.0*y[2];
}
 
int main()
{
 int i,j;
 double y[3] = {-1.5, 0.0, 1.5};
 double z[3][11] = {0};
 FODE   fa = {3, 11, 0.01, 0.0, y, (double*)z, step_ptr_lis};

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