/*

代码作者：不详

代码整理者：设计天下 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;

void adams_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],*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

{

c[i]=y[i];

}

do

{

a[1]=a[0]=hh/2.0;

a[3]=a[2]=hh;

for (i=0; i

{

g[i]=y[i];

y[i]=c[i];

}

dt=lens/m; t=x;

for (j=0; j

{

(*ptr)(t,y,n,d);

for (i=0; i

{

b[i]=y[i];

e[i]=y[i];

}

for (k=0; k<=2; k++)

{

for (i=0; 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

{

y[i]=b[i]+hh*d[i]/6.0;

}

t+=dt;

}

p=0.0;

for (i=0; i

{

q=fabs(y[i]-g[i]);

if (q>p)

{

p=q;

}

}

hh/=2.0;

m<<=1;

}

while (p>=eps) ;

/*释放动态分配的内存*/

free(g);

free(b);

free(c);

free(d);

free(e);

return;

}

void adams_lis(FODE2P ap)

{

int    n,steps,i,j,m;

double lens,eps,t,t0,q;

double *b,*e,*s,*g,*d,*y,*z;

n=ap->n;

steps=ap->steps;

b=malloc(4*n*sizeof(double));

e=malloc(n*sizeof(double));

s=malloc(n*sizeof(double));

g=malloc(n*sizeof(double));

d=malloc(n*sizeof(double));

lens=ap->lens;

eps=ap->eps;

y=ap->y;

z=ap->z;

t0=t=ap->t;

for (i=0; i

{

z[i*steps]=y[i];

}

(*ap->ptr)(t,y,n,d);

for (i=0; i

{

b[i]=d[i];

}

for (i=1; i<4; i++)

{

if (i

{

t=t0+i*lens;

adams_method_lis(n,eps,t,lens,y,ap->ptr);

for (j=0; j

{

z[j*steps+i]=y[j];

}

(*ap->ptr)(t,y,n,d);

for (j=0; j

{

b[i*n+j]=d[j];

}

}

}

for (i=4; i

{

for (j=0; 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*steps+i-1]+lens*q/24.0;

b[j]=b[n+j];

b[n+j]=b[2*n+j];

b[2*n+j]=b[3*n+j];

}

t=t0+i*lens;

(*ap->ptr)(t,y,n,d);

for (m=0; m

{

b[3*n+m]=d[m];

}

for (j=0; 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*steps+i-1]+lens*q/24.0;

z[j*steps+i]=y[j];

}

(*ap->ptr)(t,y,n,d);

for (m=0; m

{

b[3*n+m]=d[m];

}

}

/*释放动态分配的内存*/

free(b);

free(e);

free(s);

free(g);

free(d);

return;

}

void adams_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];

d[2]=-2.0*y[2];

return;

}

int main()

{

double y[3] = {0.0, 1.0, 2.0};

double z[3][11]={0};

FODE2  fa = {3, 11, 0.1, 0.0001, 0.0, y, (double*)z, adams_ptr_lis};

int i,j;

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