Gear算法程序
请问各位大侠，有没有Gear算法的程序啊，或者告诉我一声有介绍这个程序的书籍，谢谢。 吉尔法求解一阶微分方程组c程序，仅供参考10GEAR.C
#include "math.h"
#include "stdlib.h"
#include "4rinv.c"
int gear(a,b,hmin,hmax,h,eps,n,y0,k,t,z,ss,f)
void(*f)(),(*ss)();
int n,k;
double a,b,hmin,hmax,h,eps,y0[],t[],z[];
{ int kf,jt,nn,nq,i,m,irt1,irt,j,nqd,idb;
int iw,j1,j2,nt,nqw,l,*iis,*jjs;
double aa,hw,hd,rm,t0,td,r,dd,pr1,pr2,pr3,rr;
double enq1,enq2,enq3,eup,e,edwn,bnd,r1;
static double pp={ {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=malloc(n*sizeof(double));
p=malloc(n*n*sizeof(double));
s=malloc(10*n*sizeof(double));
s02=malloc(n*sizeof(double));
ym=malloc(n*sizeof(double));
er=malloc(n*sizeof(double));
yy=malloc(n*sizeof(double));
y=malloc(8*n*sizeof(double));
iis=malloc(n*sizeof(int));
jjs=malloc(n*sizeof(int));
aa=-1.0; jt=0; nn=0; nq=1; t0=a;
for (i=0; i<=8*n-1; i++) y=0.0;
for (i=0; i<=n-1; i++)
{ y=y0; yy=y;}
(*f)(t0,yy,n,d);
for (i=0; i<=n-1; i++) y=h*d;
hw=h; m=2;
for (i=0; i<=n-1; i++) ym=1.0;
l20:
irt=1; kf=1; nn=nn+1;
t=t0;
for (i=0; i<=n-1; i++) z=y;
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=y;
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=s*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y=s;
idb=m;
}
nqd=nq; td=t0; rm=1.0;
if (jt>0) goto l80;
l60:
switch (nq)
{ case 1: aa=-1.0; break;
case 2: aa=-2.0/3.0; aa=-1.0/3.0; break;
case 3: aa=-6.0/11.0; aa=aa;
aa=-1.0/11.0; break;
case 4: aa=-0.48; aa=-0.7; aa=-0.2;
aa=-0.02; break;
case 5: aa=-120.0/274.0; aa=-225.0/274.0;
aa=-85.0/274.0; aa=-15.0/274.0;
aa=-1.0/274.0; break;
case 6: aa=-720.0/1764.0; aa=-1624.0/1764.0;
aa=-735.0/1764.0; aa=-175.0/1764.0;
aa=-21.0/1764.0; aa=-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*eps; eup=eup*eup;
e=pp*eps; e=e*e;
edwn=pp*eps; edwn=edwn*edwn;
if (edwn==0.0)
{ for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
y=s;
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=y*r1;
}
idb=m;
for (i=0; i<=n-1; i++)
if (ym<fabs(y))
ym=fabs(y);
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=y+y;
}
for (i=0; i<=n-1; i++) er=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=y;
(*f)(t0,yy,n,d);
if (iw>=1)
{ for (i=0; i<=n-1; i++) yy=y;
(*ss)(t0,yy,n,p);
r=aa*h;
for (i=0; i<=n-1; i++)
for (j=0; j<=n-1; j++)
p=p*r;
for (i=0; i<=n-1; i++)
p=1.0+p;
iw=-1;
jjs=rinv(p,n);
j1=jjs;
}
if (jjs!=0)
{ for (i=0; i<=n-1; i++)
s02=y-d*h;
for (i=0; i<=n-1; i++)
{ dd=0.0;
for (j=0; j<=n-1; j++)
dd=dd+s02*p;
s=dd;
}
nt=n;
for (i=0; i<=n-1; i++)
{ y=y+aa*s;
y=y-s;
er=er+s;
if (fabs(s)<=(bnd*ym))
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=s*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y=s;
idb=m;
goto l80;
}
for (i=0; i<=n-1; i++)
for (j=0; j<=m-1; j++)
y=s;
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/ym)*(er/ym);
iw=0;
if (dd<=e)
{ if (m>=3)
for (j=2; j<=m-1; j++)
for (i=0; i<=n-1; i++)
y=y+aa*er;
kf=1; hw=h;
if (idb>1)
{ idb=idb-1;
if (idb<=1)
for (i=0; i<=n-1; i++)
s=er;
for (i=0; i<=n-1; i++)
if (ym<fabs(y)) ym=fabs(y);
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=s;
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=y;
(*f)(t0,yy,n,d);
r=h/hd;
for (i=0; i<=n-1; i++)
{ y=s;
s=hd*d;
y=s*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-s)/ym;
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/ym;
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<fabs(y)) ym=fabs(y);
jt=nq; goto l20;
}
if (nqw>nq)
for (i=0; i<=n-1; i++)
y=er*aa/(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=y*r1;
}
idb=m;
for (i=0; i<=n-1; i++)
if (ym<fabs(y)) ym=fabs(y);
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=s*r;
}
h=hd*rm;
for (i=0; i<=n-1; i++)
y=s;
idb=m;
if (nqw==nq) goto l80;
nq=nqw; goto l60;
}
4rinv.c
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
int rinv(a,n)
int n;
double a[];
{ int *is,*js,i,j,k,l,u,v;
double d,p;
is=malloc(n*sizeof(int));
js=malloc(n*sizeof(int));
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs(a);
if (p>d) { d=p; is=i; js=j;}
}
if (d+1.0==1.0)
{ free(is); free(js); printf("err**not inv\n");
return(0);
}
if (is!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is*n+j;
p=a; a=a; a=p;
}
if (js!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js;
p=a; a=a; a=p;
}
l=k*n+k;
a=1.0/a;
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j; a=a*a;}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
a=a-a*a;
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; a=-a*a;}
}
for (k=n-1; k>=0; k--)
{ if (js!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js*n+j;
p=a; a=a; a=p;
}
if (is!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is;
p=a; a=a; a=p;
}
}
free(is); free(js);
return(1);
}
回复 2楼 的帖子
谢谢，仔细研究一番！
页:
[1]