#include "stdafx.h"
#include <stdlib.h>
#include <math.h>
#include "DifferentialEguation.h"
#include "MatrixAlgo.h"
#include "LinearEquation.h"
//全区间积分的特雷纳方法
void gtnr2(double t, double h, int n, double y[], int k, double z[])
{
int i,j;
double a,s,aa,bb,dd,g,dy,dy1,*d,*p,*w,*q,*r;
d=(double *)malloc(n*sizeof(double));
p=(double *)malloc(n*sizeof(double));
w=(double *)malloc(4*n*sizeof(double));
q=(double *)malloc(4*n*sizeof(double));
r=(double *)malloc(4*n*sizeof(double));
a=t;
for (j=0; j<=n-1; j++)
z[j*k]=y[j];
for (i=1; i<=k-1; i++)
{
t=a+(i-1)*h;
for(j=0; j<=n-1; j++)
w[j]=y[j];
gtnr2f(t,y,n,d);
for(j=0; j<=n-1; j++)
{
q[j]=d[j];
y[j]=w[j]+h*d[j]/2.0;
w[n+j]=y[j];
}
s=t+h/2.0;
gtnr2f(s,y,n,d);
for (j=0; j<=n-1; j++)
{
q[n+j]=d[j];
y[j]=w[j]+h*d[j]/2.0;
w[n+n+j]=y[j];
}
gtnr2f(s,y,n,d);
for (j=0; j<=n-1; j++)
q[n+n+j]=d[j];
for (j=0; j<=n-1; j++)
{
aa=q[n+n+j]-q[n+j];
bb=w[n+n+j]-w[n+j];
if (-aa*bb*h>0.0)
{
p[j]=-aa/bb;
dd=-p[j]*h;
r[j]=exp(dd);
r[n+j]=(r[j]-1.0)/dd;
r[n+n+j]=(r[n+j]-1.0)/dd;
r[3*n+j]=(r[n+n+j]-1.0)/dd;
}
else
p[j]=0.0;
if (p[j]<=0.0)
g=q[n+n+j];
else
{
g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
}
w[3*n+j]=w[j]+g*h;
y[j]=w[3*n+j];
}
s=t+h;
gtnr2f(s,y,n,d);
for(j=0; j<=n-1; j++)
q[3*n+j]=d[j];
for(j=0; j<=n-1; j++)
{
if (p[j]<=0.0)
{
dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
dy=(dy+q[n+n+n+j])*h/6.0;
}
else
{
dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]+p[j]*w[n+j]);
dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
dy=dy*r[n+n+j]+q[j]*r[n+j];
dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
dy=(dy+4.0*dy1*r[n+n+n+j])*h;
}
y[j]=w[j]+dy;
z[j*k+i]=y[j];
}
}
free(d);
free(p);
free(w);
free(q);
free(r);
return;
}
//全区间积分的连分式法
void gpbs2(double t, double h, int n, double y[], double eps, int k, double z[])
{
int i,j,kk,l,m,nn,it;
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);
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);
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;
}
//全区间积分的双边法
void ggjfq(double t, double h, int n, double y[], double eps, int k, double z[])
{
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);
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;
}
//积分一步的变步长欧拉方法
void gelr2(double t, double h, double y[], int n, double eps)
{
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 gtnr1(double t, double h, int n, double y[])
{
int j;
double s,aa,bb,dd,g,dy,dy1,*d,*p,*w,*q,*r;
w=(double *)malloc(4*n*sizeof(double));
q=(double *)malloc(4*n*sizeof(double));
r=(double *)malloc(4*n*sizeof(double));
d=(double *)malloc(n*sizeof(double));
p=(double *)malloc(n*sizeof(double));
for (j=0; j<=n-1; j++)
w[j]=y[j];
gtnr1f(t,y,n,d);
for (j=0; j<=n-1; j++)
{
q[j]=d[j];
y[j]=w[j]+h*d[j]/2.0;
w[n+j]=y[j];
}
s=t+h/2.0;
gtnr1f(s,y,n,d);
for (j=0; j<=n-1; j++)
{
q[n+j]=d[j];
y[j]=w[j]+h*d[j]/2.0;
w[n+n+j]=y[j];
}
gtnr1f(s,y,n,d);
for (j=0; j<=n-1; j++)
q[n+n+j]=d[j];
for (j=0; j<=n-1; j++)
{
aa=q[n+n+j]-q[n+j];
bb=w[n+n+j]-w[n+j];
if (-aa*bb*h>0.0)
{
p[j]=-aa/bb;
dd=-p[j]*h;
r[j]=exp(dd);
r[n+j]=(r[j]-1.0)/dd;
r[n+n+j]=(r[n+j]-1.0)/dd;
r[3*n+j]=(r[n+n+j]-1.0)/dd;
}
else
p[j]=0.0;
if (p[j]<=0.0)
g=q[n+n+j];
else
{
g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
}
w[3*n+j]=w[j]+g*h;
y[j]=w[3*n+j];
}
s=t+h;
gtnr1f(s,y,n,d);
for (j=0; j<=n-1; j++)
q[3*n+j]=d[j];
for (j=0; j<=n-1; j++)
{
if (p[j]<=0.0)
{
dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
dy=(dy+q[n+n+n+j])*h/6.0;
}
else
{
dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]+p[j]*w[n+j]);
dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
dy=dy*r[n+n+j]+q[j]*r[n+j];
dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
dy=(dy+4.0*dy1*r[n+n+n+j])*h;
}
y[j]=w[j]+dy;
}
free(d);
free(p);
free(w);
free(q);
free(r);
return;
}
//积分一步的连分式法
void gpbs1(double t, double h, int n, double y[], double eps)
{
int i,j,k,m,nn,it;
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);
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);
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;
}
//全区问积分的变步长基尔方法
void ggil2(double t, double h, double y[], int n, double eps, int k, double z[])
{
int i,j,m,kk,ii;
double aa,hh,x,p,dt,r,s,t0,qq,*g,*q,*d,*u,*v;
double a[4]={0.5,0.29289321881,1.7071067812,0.166666667};
double b[4]={2.0,1.0,1.0,2.0};
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 gelr1(double t, double y[], int n, double h, int k, double z[])
{
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 grkt2(double t, double h, double y[], int n, double eps)
{
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 gmrsn(double t, double h, int n, double y[], double eps, int k, double z[])
{
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);
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;
}
//积分一步的变步长基尔方法
void ggil1(double t, double h, double y[], int n, double eps, double q[])
{
int i,j,k,m,ii;
double x,p,hh,r,s,t0,dt,qq,*d,*u,*v,*g;
double a[4]={0.5,0.29289321881,1.7071067812,0.166666667};
double b[4]={2.0,1.0,1.0,2.0};
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 gdfte(double a, double b, double ya, double yb, int n, double y[])
{
int i,k,nn,m1;
double z[4],h,x,*g,*d;
g=(double *)malloc(6*n*sizeof(double));
d=(double *)malloc(2*n*sizeof(double));
h=(b-a)/(n-1.0);
nn=2*n-1;
g[0]=1.0;
g[1]=0.0;
y[0]=ya;
y[n-1]=yb;
g[3*n-3]=1.0;
g[3*n-4]=0.0;
for (i=2; i<=n-1; i++)
{
x=a+(i-1)*h;
gdftef(x,z);
k=3*(i-1)-1;
g[k]=z[0]-h*z[1]/2.0;
g[k+1]=h*h*z[2]-2.0*z[0];
g[k+2]=z[0]+h*z[1]/2.0;
y[i-1]=h*h*z[3];
}
m1=3*n-2;
atrde(g,n,m1,y);
h=h/2.0;
g[0]=1.0;
g[1]=0.0;
d[0]=ya;
d[nn-1]=yb;
g[3*nn-3]=1.0;
g[3*nn-4]=0.0;
for (i=2; i<=nn-1; i++)
{
x=a+(i-1)*h;
gdftef(x,z);
k=3*(i-1)-1;
g[k]=z[0]-h*z[1]/2.0;
g[k+1]=h*h*z[2]-2.0*z[0];
g[k+2]=z[0]+h*z[1]/2.0;
d[i-1]=h*h*z[3];
}
m1=3*nn-2;
atrde(g,nn,m1,d);
for (i=2; i<=n-1; i++)
{
k=i+i-1;
y[i-1]=(4.0*d[k-1]-y[i-1])/3.0;
}
free(g);
free(d);
return;
}
//全区间积分的定步长维梯方法
void gwity(double t, double y[], int n, double h, int k, double z[])
{
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[])
{
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 ghamg(double t, double h, int n, double y[], double eps, int k, double z[])
{
int i,j,m;
double a,q,*b,*d,*u,*v,*w,*g;
b=(double *)malloc(4*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));
g=(double *)malloc(n*sizeof(double));
a=t;
for (i=0; i<=n-1; i++)
z[i*k]=y[i];
ghamgf(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;
rkt5(t,h,y,n,eps);
for (m=0; m<=n-1; m++)
z[m*k+i]=y[m];
ghamgf(t,y,n,d);
for (m=0; m<=n-1; m++)
b[i*n+m]=d[m];
}
for (i=0; i<=n-1; i++)
u[i]=0.0;
for (i=4; i<=k-1; i++)
{
for (j=0; j<=n-1; j++)
{
q=2.0*b[3*n+j]-b[n+n+j]+2.0*b[n+j];
y[j]=z[j*k+i-4]+4.0*h*q/3.0;
}
for (j=0; j<=n-1; j++)
y[j]=y[j]+112.0*u[j]/121.0;
t=a+i*h;
ghamgf(t,y,n,d);
for (j=0; j<=n-1; j++)
{
q=9.0*z[j*k+i-1]-z[j*k+i-3];
q=(q+3.0*h*(d[j]+2.0*b[3*n+j]-b[n+n+j]))/8.0;
u[j]=q-y[j];
z[j*k+i]=q-9.0*u[j]/121.0;
y[j]=z[j*k+i];
b[n+j]=b[n+n+j];
b[n+n+j]=b[n+n+n+j];
}
ghamgf(t,y,n,d);
for (m=0; m<=n-1; m++)
b[3*n+m]=d[m];
}
free(b);
free(d);
free(u);
free(v);
free(w);
free(g);
return;
}
//积分刚性方程组的吉尔方法
int ggear(double a, double b, double hmin, double hmax, double h, double eps, int n, double y0[], int k, double t[], double z[])
{
int kf,jt,nn,nq,i,m,irt1,irt,j,nqd,idb;
int iw,j1,j2,nt,nqw,l,*iis,*jjs;
double aa[7],hw,hd,rm,t0,td,r,dd,pr1,pr2,pr3,rr;
double enq1,enq2,enq3,eup,e,edwn,bnd,r1;
double pp[7][3]={ {2.0,3.0,1.0},{4.5,6.0,1.0},
{7.333,9.167,0.5},{10.42,12.5,0.1667},
{13.7,15.98,0.04133},{17.15,1.0,0.008267},
{1.0,1.0,1.0}};
double *d,*p,*s,*s02,*ym,*er,*yy,*y;
d=(double *)malloc(n*sizeof(double));
p=(double *)malloc(n*n*sizeof(double));
s=(double *)malloc(10*n*sizeof(double));
s02=(double *)malloc(n*sizeof(double));
ym=(double *)malloc(n*sizeof(double));
er=(double *)malloc(n*sizeof(double));
yy=(double *)malloc(n*sizeof(double));
y=(double *)malloc(8*n*sizeof(double));
iis=(int *)malloc(n*sizeof(int));
jjs=(int *)malloc(n*sizeof(int));
aa[1]=-1.0;
jt=0;
nn=0;
nq=1;
t0=a;
for (i=0; i<=8*n-1; i++)
y[i]=0.0;
for (i=0; i<=n-1; i++)
{
y[i*8]=y0[i];
yy[i]=y[i*8];
}
ggearf(t0,yy,n,d);
for (i=0; i<=n-1; i++)
y[i*8+1]=h*d[i];
hw=h;
m=2;
for (i=0; i<=n-1; i++)
ym[i]=1.0;
l20:
irt=1;
kf=1;
nn=nn+1;
t[nn-1]=t0;
for (i=0; i<=n-1; i++)
z[i*k+nn-1]=y[i*8];
if ((t0>=b)||(nn==k))
{
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
return(kf);
}
for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
s[i*10+j]=y[i*8+j];
hd=hw;
if (h!=hd)
{
rm=h/hd;
irt1=0;
rr=fabs(hmin/hd);
if (rm<rr)
rm=rr;
rr=fabs(hmax/hd);
if (rm>rr)
rm=rr;
r=1.0;
irt1=irt1+1;
for (j=1; j<=m-1; j++)
{
r=r*rm;
for (i=0; i<=n-1; i++)
y[i*8+j]=s[i*10+j]*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y[i*8]=s[i*10];
idb=m;
}
nqd=nq;
td=t0;
rm=1.0;
if (jt>0)
goto l80;
l60:
switch (nq){
case 1:
aa[0]=-1.0;
break;
case 2:
aa[0]=-2.0/3.0;
aa[2]=-1.0/3.0;
break;
case 3:
aa[0]=-6.0/11.0;
aa[2]=aa[0];
aa[3]=-1.0/11.0;
break;
case 4:
aa[0]=-0.48;
aa[2]=-0.7;
aa[3]=-0.2;
aa[4]=-0.02;
break;
case 5:
aa[0]=-120.0/274.0;
aa[2]=-225.0/274.0;
aa[3]=-85.0/274.0;
aa[4]=-15.0/274.0;
aa[5]=-1.0/274.0;
break;
case 6:
aa[0]=-720.0/1764.0;
aa[2]=-1624.0/1764.0;
aa[3]=-735.0/1764.0;
aa[4]=-175.0/1764.0;
aa[5]=-21.0/1764.0;
aa[6]=-1.0/1764.0;
break;
default:
{
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
return(-2);
}
}
m=nq+1;
idb=m;
enq2=0.5/(nq+1.0);
enq3=0.5/(nq+2.0);
enq1=0.5/(nq+0.0);
eup=pp[nq-1][1]*eps;
eup=eup*eup;
e=pp[nq-1][0]*eps;
e=e*e;
edwn=pp[nq-1][2]*eps;
edwn=edwn*edwn;
if (edwn==0.0)
{
for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
y[i*8+j]=s[i*10+j];
h=hd;
nq=nqd;
jt=nq;
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
return(-4);
}
bnd=eps*enq3/(n+0.0);
iw=1;
if (irt==2)
{
r1=1.0;
for (j=1; j<=m-1; j++)
{
r1=r1*r;
for (i=0; i<=n-1; i++)
y[i*8+j]=y[i*8+j]*r1;
}
idb=m;
for (i=0; i<=n-1; i++)
if (ym[i]<fabs(y[i*8]))
ym[i]=fabs(y[i*8]);
jt=nq;
goto l20;
}
l80:
t0=t0+h;
for (j=2; j<=m; j++)
for (j1=j; j1<=m; j1++)
{
j2=m-j1+j-1;
for (i=0; i<=n-1; i++)
y[i*8+j2-1]=y[i*8+j2-1]+y[i*8+j2];
}
for (i=0; i<=n-1; i++)
er[i]=0.0;
j1=1;
nt=1;
for (l=0; l<=2; l++)
{
if ((j1!=0)&&(nt!=0))
{
for (i=0; i<=n-1; i++)
yy[i]=y[i*8];
ggearf(t0,yy,n,d);
if (iw>=1)
{
for (i=0; i<=n-1; i++)
yy[i]=y[i*8];
ggears(t0,yy,n,p);
r=aa[0]*h;
for (i=0; i<=n-1; i++)
for (j=0; j<=n-1; j++)
p[i*n+j]=p[i*n+j]*r;
for (i=0; i<=n-1; i++)
p[i*n+i]=1.0+p[i*n+i];
iw=-1;
jjs[0]=brinv(p,n);
j1=jjs[0];
}
if (jjs[0]!=0)
{
for (i=0; i<=n-1; i++)
s02[i]=y[i*8+1]-d[i]*h;
for (i=0; i<=n-1; i++)
{
dd=0.0;
for (j=0; j<=n-1; j++)
dd=dd+s02[j]*p[i*n+j];
s[i*10+8]=dd;
}
nt=n;
for (i=0; i<=n-1; i++)
{
y[i*8]=y[i*8]+aa[0]*s[i*10+8];
y[i*8+1]=y[i*8+1]-s[i*10+8];
er[i]=er[i]+s[i*10+8];
if (fabs(s[i*10+8])<=(bnd*ym[i]))
nt=nt-1;
}
}
}
}
if (nt>0)
{
t0=td;
if ((h>(hmin*1.00001))||(iw>=0))
{
if (iw!=0)
rm=0.25*rm;
iw=1;
irt1=2;
rr=fabs(hmin/hd);
if (rm<rr)
rm=rr;
rr=fabs(hmax/hd);
if (rm>rr)
rm=rr;
r=1.0;
for (j=1; j<=m-1; j++)
{
r=r*rm;
for (i=0; i<=n-1; i++)
y[i*8+j]=s[i*10+j]*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y[i*8]=s[i*10];
idb=m;
goto l80;
}
for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
y[i*8+j]=s[i*10+j];
h=hd;
nq=nqd;
jt=nq;
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
return(-3);
}
dd=0.0;
for (i=0; i<=n-1; i++)
dd=dd+(er[i]/ym[i])*(er[i]/ym[i]);
iw=0;
if (dd<=e)
{
if (m>=3)
for (j=2; j<=m-1; j++)
for (i=0; i<=n-1; i++)
y[i*8+j]=y[i*8+j]+aa[j]*er[i];
kf=1;
hw=h;
if (idb>1)
{
idb=idb-1;
if (idb<=1)
for (i=0; i<=n-1; i++)
s[i*10+9]=er[i];
for (i=0; i<=n-1; i++)
if (ym[i]<fabs(y[i*8]))
ym[i]=fabs(y[i*8]);
jt=nq;
goto l20;
}
}
if (dd>e)
{
kf=kf-2;
if (h<=(hmin*1.00001))
{
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
hw=h;
jt=nq;
return(-1);
}
t0=td;
if (kf<=-5)
{
if (nq==1)
{
for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
y[i*8+j]=s[i*10+j];
h=hd;
nq=nqd;
jt=nq;
free(d);
free(p);
free(s);
free(s02);
free(ym);
free(er);
free(yy);
free(iis);
free(jjs);
free(y);
return(-4);
}
for (i=0; i<=n-1; i++)
yy[i]=y[i*8];
ggearf(t0,yy,n,d);
r=h/hd;
for (i=0; i<=n-1; i++)
{
y[i*8]=s[i*10];
s[i*10+1]=hd*d[i];
y[i*8+1]=s[i*10+1]*r;
}
nq=1;
kf=1;
goto l60;
}
}
pr2=log(dd/e);
pr2=enq2*pr2;
pr2=exp(pr2);
pr2=1.2*pr2;
pr3=1.0e+20;
if (nq<7)
if (kf>-1)
{
dd=0.0;
for (i=0; i<=n-1; i++)
{
pr3=(er[i]-s[i*10+9])/ym[i];
dd=dd+pr3*pr3;
}
pr3=log(dd/eup);
pr3=enq3*pr3;
pr3=exp(pr3);
pr3=1.4*pr3;
}
pr1=1.0e+20;
if (nq>1)
{
dd=0.0;
for (i=0; i<=n-1; i++)
{
pr1=y[i*8+m-1]/ym[i];
dd=dd+pr1*pr1;
}
pr1=log(dd/edwn);
pr1=enq1*pr1;
pr1=exp(pr1);
pr1=1.3*pr1;
}
if (pr2<=pr3)
{
if (pr2>pr1)
{
r=1.0e+04;
if (pr1>1.0e-04)
r=1.0/pr1;
nqw=nq-1;
}
else
{
nqw=nq;
r=1.0e+04;
if (pr2>1.0e-04)
r=1.0/pr2;
}
}
else
{
if (pr3<pr1)
{
r=1.0e+04;
if (pr3>1.0e-04)
r=1.0/pr3;
nqw=nq+1;
}
else
{
r=1.0e+04;
if (pr1>1.0e-04)
r=1.0/pr1;
nqw=nq-1;
}
}
idb=10;
if (kf==1)
if (r<1.1)
{
for (i=0; i<=n-1; i++)
if (ym[i]<fabs(y[i*8]))
ym[i]=fabs(y[i*8]);
jt=nq;
goto l20;
}
if (nqw>nq)
for (i=0; i<=n-1; i++)
y[i*8+nqw]=er[i]*aa[m-1]/(m+0.0);
m=nqw+1;
if (kf==1)
{
irt=2;
rr=hmax/fabs(h);
if (r>rr)
r=rr;
h=h*r;
hw=h;
if (nq==nqw)
{
r1=1.0;
for (j=1; j<=m-1; j++)
{
r1=r1*r;
for (i=0; i<=n-1; i++)
y[i*8+j]=y[i*8+j]*r1;
}
idb=m;
for (i=0; i<=n-1; i++)
if (ym[i]<fabs(y[i*8]))
ym[i]=fabs(y[i*8]);
jt=nq;
goto l20;
}
nq=nqw;
goto l60;
}
rm=rm*r;
irt1=3;
rr=fabs(hmin/hd);
if (rm<rr)
rm=rr;
rr=fabs(hmax/hd);
if (rm>rr)
rm=rr;
r=1.0;
for (j=1; j<=m-1; j++)
{
r=r*rm;
for (i=0; i<=n-1; i++)
y[i*8+j]=s[i*10+j]*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y[i*8]=s[i*10];
idb=m;
if (nqw==nq)
goto l80;
nq=nqw;
goto l60;
}
//全区间积分的阿当姆斯预报-校正法
void gadms(double t, double h, int n, double y[], double eps, int k, double z[])
{
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);
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;
}
void rkt1(double t, double h, int n, double y[], double b[], double d[], double e[])
{
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 rkt2(double t, double h, int n, double y[], double b[], double d[], double e[])
{
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 rkt3(double t, double h, double y[], int n, double eps)
{
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 rkt4(double t, double h, double y[], int n, double eps)
{
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 rkt5(double t, double h, double y[], int n, double eps)
{
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++)
{
ghamgf(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];
ghamgf(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 gelr1f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void gelr2f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void gwityf(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void grkt1f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void grkt2f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
return;
}
void ggil1f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void ggil2f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void gmrsnf(double t, double y[], int n, double d[])
{
double q;
t = t;
n = n;
q = 60.0 * (0.06 + t * (t - 0.6));
d[0] = q * y[1];
d[1] = -q * y[0];
return;
}
void gpbs1f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
return;
}
void gpbs2f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = -y[1];
d[1] = y[0];
return;
}
void ggjfqf(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = -y[1];
d[1] = y[0];
return;
}
void gadmsf(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = -y[2];
return;
}
void ghamgf(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = y[1];
d[1] = -y[0];
d[2] = y[2];
return;
}
void gtnr1f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = -21.0 * y[0] + 19.0 * y[1] - 20.0 * y[2];
d[1] = 19.0 * y[0] - 21.0 * y[1] + 20.0 * y[2];
d[2] = 40.0 * y[0] - 40.0 * y[1] - 40.0 * y[2];
return;
}
void gtnr2f(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = -21.0 * y[0] + 19.0 * y[1] - 20.0 * y[2];
d[1] = 19.0 * y[0] - 21.0 * y[1] + 20.0 * y[2];
d[2] = 40.0 * y[0] - 40.0 * y[1] - 40.0 * y[2];
return;
}
void ggearf(double t, double y[], int n, double d[])
{
t = t;
n = n;
d[0] = -21.0 * y[0] + 19.0 * y[1] - 20.0 * y[2];
d[1] = 19.0 * y[0] - 21.0 * y[1] + 20.0 * y[2];
d[2] = 40.0 * y[0] - 40.0 * y[1] - 40.0 * y[2];
return;
}
void ggears(double t, double y[], int n, double p[3][3])
{
t = t;
n = n;
y[0] = y[0];
p[0][0] = -21.0;
p[0][1] = -19.0;
p[0][2] = -20.0;
p[1][0] = 19.0;
p[1][1] = -21.0;
p[1][2] = 20.0;
p[2][0] = 40.0;
p[2][1] = -40.0;
p[2][2] = -40.0;
return;
}
void gdftef(double x, double z[4])
{
z[0] = -1.0;
z[1] = 0.0;
z[2] = 2.0/(x * x);
z[3] = 1.0/x;
return;
}
----根据《C语言算法程序集》整理
扫一扫
7万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
