Gear算法程序

zmf_0218

请问各位大侠，有没有Gear算法的程序啊，或者告诉我一声有介绍这个程序的书籍，谢谢。

风花雪月

吉尔法求解一阶微分方程组c程序，仅供参考

10GEAR.C

1.   #include "math.h"
   
2.   #include "stdlib.h"
   
3.   #include "4rinv.c"
   
4.   int gear(a,b,hmin,hmax,h,eps,n,y0,k,t,z,ss,f)
   
5.   void  (*f)(),(*ss)();
   
6.   int n,k;
   
7.   double a,b,hmin,hmax,h,eps,y0[],t[],z[];
   
8.   { int kf,jt,nn,nq,i,m,irt1,irt,j,nqd,idb;
   
9.     int iw,j1,j2,nt,nqw,l,*iis,*jjs;
   
10.     double aa[7],hw,hd,rm,t0,td,r,dd,pr1,pr2,pr3,rr;
    
11.     double enq1,enq2,enq3,eup,e,edwn,bnd,r1;
    
12.     static double pp[7][3]={ {2.0,3.0,1.0},{4.5,6.0,1.0},
    
13.            {7.333,9.167,0.5},{10.42,12.5,0.1667},
    
14.            {13.7,15.98,0.04133},{17.15,1.0,0.008267},
    
15.            {1.0,1.0,1.0}};
    
16.     double *d,*p,*s,*s02,*ym,*er,*yy,*y;
    
17.     d=malloc(n*sizeof(double));
    
18.     p=malloc(n*n*sizeof(double));
    
19.     s=malloc(10*n*sizeof(double));
    
20.     s02=malloc(n*sizeof(double));
    
21.     ym=malloc(n*sizeof(double));
    
22.     er=malloc(n*sizeof(double));
    
23.     yy=malloc(n*sizeof(double));
    
24.     y=malloc(8*n*sizeof(double));
    
25.     iis=malloc(n*sizeof(int));
    
26.     jjs=malloc(n*sizeof(int));
    
27.     aa[1]=-1.0; jt=0; nn=0; nq=1; t0=a;
    
28.     for (i=0; i<=8*n-1; i++) y[i]=0.0;
    
29.     for (i=0; i<=n-1; i++)
    
30.        { y[i*8]=y0[i]; yy[i]=y[i*8];}
    
31.     (*f)(t0,yy,n,d);
    
32.     for (i=0; i<=n-1; i++) y[i*8+1]=h*d[i];
    
33.     hw=h; m=2;
    
34.     for (i=0; i<=n-1; i++) ym[i]=1.0;
    
35.     l20:
    
36.     irt=1; kf=1; nn=nn+1;
    
37.     t[nn-1]=t0;
    
38.     for (i=0; i<=n-1; i++) z[i*k+nn-1]=y[i*8];
    
39.     if ((t0>=b)||(nn==k))
    
40.       { free(d); free(p); free(s); free(s02);
    
41.         free(ym); free(er); free(yy); free(iis);
    
42.         free(jjs); free(y); return(kf);}
    
43.     for (i=0; i<=n-1; i++)
    
44.       for (j=0; j<=m-1; j++) s[i*10+j]=y[i*8+j];
    
45.     hd=hw;
    
46.     if (h!=hd)
    
47.       { rm=h/hd; irt1=0;
    
48.         rr=fabs(hmin/hd);
    
49.         if (rm<rr) rm=rr;
    
50.         rr=fabs(hmax/hd);
    
51.         if (rm>rr) rm=rr;
    
52.     r=1.0; irt1=irt1+1;
    
53.         for (j=1; j<=m-1; j++)
    
54.           { r=r*rm;
    
55.             for (i=0; i<=n-1; i++)
    
56.               y[i*8+j]=s[i*10+j]*r;
    
57.           }
    
58.         h=hd*rm;
    
59.         for (i=0; i<=n-1; i++)
    
60.           y[i*8]=s[i*10];
    
61.         idb=m;
    
62.       }
    
63.     nqd=nq; td=t0; rm=1.0;
    
64.     if (jt>0) goto l80;
    
65.     l60:
    
66.     switch (nq)
    
67.       { case 1: aa[0]=-1.0; break;
    
68.         case 2: aa[0]=-2.0/3.0; aa[2]=-1.0/3.0; break;
    
69.         case 3: aa[0]=-6.0/11.0; aa[2]=aa[0];
    
70.                 aa[3]=-1.0/11.0; break;
    
71.         case 4: aa[0]=-0.48; aa[2]=-0.7; aa[3]=-0.2;
    
72.                 aa[4]=-0.02; break;
    
73.         case 5: aa[0]=-120.0/274.0; aa[2]=-225.0/274.0;
    
74.                 aa[3]=-85.0/274.0; aa[4]=-15.0/274.0;
    
75.                 aa[5]=-1.0/274.0; break;
    
76.         case 6: aa[0]=-720.0/1764.0; aa[2]=-1624.0/1764.0;
    
77.                 aa[3]=-735.0/1764.0; aa[4]=-175.0/1764.0;
    
78.                 aa[5]=-21.0/1764.0; aa[6]=-1.0/1764.0;
    
79.                 break;
    
80.         default: { free(d); free(p); free(s); free(s02);
    
81.                    free(ym); free(er); free(yy);
    
82.                    free(iis); free(jjs); free(y); return(-2);
    
83.                  }
    
84.       }
    
85.     m=nq+1; idb=m;
    
86.     enq2=0.5/(nq+1.0); enq3=0.5/(nq+2.0);
    
87.     enq1=0.5/(nq+0.0);
    
88.     eup=pp[nq-1][1]*eps; eup=eup*eup;
    
89.     e=pp[nq-1][0]*eps; e=e*e;
    
90.     edwn=pp[nq-1][2]*eps; edwn=edwn*edwn;
    
91.     if (edwn==0.0)
    
92.       { for (i=0; i<=n-1; i++)
    
93.           for (j=0; j<=m-1; j++)
    
94.             y[i*8+j]=s[i*10+j];
    
95.         h=hd; nq=nqd; jt=nq;
    
96.         free(d); free(p); free(s); free(s02);
    
97.         free(ym); free(er); free(yy); free(iis);
    
98.         free(jjs); free(y); return(-4);
    
99.       }
    
100.     bnd=eps*enq3/(n+0.0);
     
101.     iw=1;
     
102.     if (irt==2)
     
103.       { r1=1.0;
     
104.         for (j=1; j<=m-1; j++)
     
105.           { r1=r1*r;
     
106.             for (i=0; i<=n-1; i++)
     
107.               y[i*8+j]=y[i*8+j]*r1;
     
108.           }
     
109.         idb=m;
     
110.         for (i=0; i<=n-1; i++)
     
111.           if (ym[i]<fabs(y[i*8]))
     
112.             ym[i]=fabs(y[i*8]);
     
113.         jt=nq;
     
114.         goto l20;
     
115.       }
     
116.     l80:
     
117.     t0=t0+h;
     
118.     for (j=2; j<=m; j++)
     
119.       for (j1=j; j1<=m; j1++)
     
120.         { j2=m-j1+j-1;
     
121.           for (i=0; i<=n-1; i++)
     
122.             y[i*8+j2-1]=y[i*8+j2-1]+y[i*8+j2];
     
123.         }
     
124.     for (i=0; i<=n-1; i++) er[i]=0.0;
     
125.     j1=1; nt=1;
     
126.     for (l=0; l<=2; l++)
     
127.       { if ((j1!=0)&&(nt!=0))
     
128.           { for (i=0; i<=n-1; i++) yy[i]=y[i*8];
     
129.             (*f)(t0,yy,n,d);
     
130.             if (iw>=1)
     
131.               { for (i=0; i<=n-1; i++) yy[i]=y[i*8];
     
132.                 (*ss)(t0,yy,n,p);
     
133.                 r=aa[0]*h;
     
134.                 for (i=0; i<=n-1; i++)
     
135.                   for (j=0; j<=n-1; j++)
     
136.                     p[i*n+j]=p[i*n+j]*r;
     
137.                 for (i=0; i<=n-1; i++)
     
138.                   p[i*n+i]=1.0+p[i*n+i];
     
139.                 iw=-1;
     
140.                 jjs[0]=rinv(p,n);
     
141.                 j1=jjs[0];
     
142.               }
     
143.             if (jjs[0]!=0)
     
144.               { for (i=0; i<=n-1; i++)
     
145.                   s02[i]=y[i*8+1]-d[i]*h;
     
146.                 for (i=0; i<=n-1; i++)
     
147.                   { dd=0.0;
     
148.                     for (j=0; j<=n-1; j++)
     
149.                       dd=dd+s02[j]*p[i*n+j];
     
150.                     s[i*10+8]=dd;
     
151.                   }
     
152.                 nt=n;
     
153.                 for (i=0; i<=n-1; i++)
     
154.                   { y[i*8]=y[i*8]+aa[0]*s[i*10+8];
     
155.                     y[i*8+1]=y[i*8+1]-s[i*10+8];
     
156.                     er[i]=er[i]+s[i*10+8];
     
157.                     if (fabs(s[i*10+8])<=(bnd*ym[i]))
     
158.                       nt=nt-1;
     
159.                   }
     
160.               }
     
161.           }
     
162.       }
     
163.     if (nt>0)
     
164.       { t0=td;
     
165.         if ((h>(hmin*1.00001))||(iw>=0))
     
166.           { if (iw!=0) rm=0.25*rm;
     
167.             iw=1; irt1=2;
     
168.             rr=fabs(hmin/hd);
     
169.             if (rm<rr) rm=rr;
     
170.             rr=fabs(hmax/hd);
     
171.             if (rm>rr) rm=rr;
     
172.             r=1.0;
     
173.             for (j=1; j<=m-1; j++)
     
174.               { r=r*rm;
     
175.                 for (i=0; i<=n-1; i++)
     
176.                   y[i*8+j]=s[i*10+j]*r;
     
177.               }
     
178.             h=hd*rm;
     
179.             for (i=0; i<=n-1; i++)
     
180.               y[i*8]=s[i*10];
     
181.             idb=m;
     
182.             goto l80;
     
183.           }
     
184.         for (i=0; i<=n-1; i++)
     
185.           for (j=0; j<=m-1; j++)
     
186.             y[i*8+j]=s[i*10+j];
     
187.         h=hd; nq=nqd; jt=nq;
     
188.         free(d); free(p); free(s); free(s02);
     
189.         free(ym); free(er); free(yy);
     
190.         free(iis); free(jjs); free(y); return(-3);
     
191.       }
     
192.     dd=0.0;
     
193.     for (i=0; i<=n-1; i++)
     
194.       dd=dd+(er[i]/ym[i])*(er[i]/ym[i]);
     
195.     iw=0;
     
196.     if (dd<=e)
     
197.       { if (m>=3)
     
198.           for (j=2; j<=m-1; j++)
     
199.             for (i=0; i<=n-1; i++)
     
200.               y[i*8+j]=y[i*8+j]+aa[j]*er[i];
     
201.         kf=1; hw=h;
     
202.         if (idb>1)
     
203.           { idb=idb-1;
     
204.             if (idb<=1)
     
205.               for (i=0; i<=n-1; i++)
     
206.                 s[i*10+9]=er[i];
     
207.             for (i=0; i<=n-1; i++)
     
208.           if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
     
209.             jt=nq;
     
210.             goto l20;
     
211.           }
     
212.       }
     
213.     if (dd>e)
     
214.       { kf=kf-2;
     
215.         if (h<=(hmin*1.00001))
     
216.           { free(d); free(p); free(s); free(s02);
     
217.             free(ym); free(er); free(yy);
     
218.             free(iis); free(jjs); free(y);
     
219.             hw=h; jt=nq; return(-1);
     
220.           }
     
221.         t0=td;
     
222.         if (kf<=-5)
     
223.           { if (nq==1)
     
224.               { for (i=0; i<=n-1; i++)
     
225.                   for (j=0; j<=m-1; j++)
     
226.                     y[i*8+j]=s[i*10+j];
     
227.                 h=hd; nq=nqd; jt=nq;
     
228.         free(d); free(p); free(s); free(s02);
     
229.         free(ym); free(er); free(yy);
     
230.         free(iis); free(jjs); free(y); return(-4);
     
231.               }
     
232.             for (i=0; i<=n-1; i++) yy[i]=y[i*8];
     
233.             (*f)(t0,yy,n,d);
     
234.             r=h/hd;
     
235.             for (i=0; i<=n-1; i++)
     
236.               { y[i*8]=s[i*10];
     
237.                 s[i*10+1]=hd*d[i];
     
238.                 y[i*8+1]=s[i*10+1]*r;
     
239.               }
     
240.             nq=1; kf=1; goto l60;
     
241.           }
     
242.       }
     
243.     pr2=log(dd/e); pr2=enq2*pr2; pr2=exp(pr2);
     
244.     pr2=1.2*pr2;
     
245.     pr3=1.0e+20;
     
246.     if (nq<7)
     
247.       if (kf>-1)
     
248.         { dd=0.0;
     
249.           for (i=0; i<=n-1; i++)
     
250.             { pr3=(er[i]-s[i*10+9])/ym[i];
     
251.               dd=dd+pr3*pr3;
     
252.             }
     
253.           pr3=log(dd/eup); pr3=enq3*pr3;
     
254.           pr3=exp(pr3); pr3=1.4*pr3;
     
255.         }
     
256.     pr1=1.0e+20;
     
257.     if (nq>1)
     
258.       { dd=0.0;
     
259.         for (i=0; i<=n-1; i++)
     
260.           { pr1=y[i*8+m-1]/ym[i];
     
261.             dd=dd+pr1*pr1;
     
262.           }
     
263.         pr1=log(dd/edwn); pr1=enq1*pr1;
     
264.         pr1=exp(pr1); pr1=1.3*pr1;
     
265.       }
     
266.     if (pr2<=pr3)
     
267.       { if (pr2>pr1)
     
268.           { r=1.0e+04;
     
269.             if (pr1>1.0e-04) r=1.0/pr1;
     
270.             nqw=nq-1;
     
271.           }
     
272.         else
     
273.           { nqw=nq; r=1.0e+04;
     
274.             if (pr2>1.0e-04) r=1.0/pr2;
     
275.           }
     
276.       }
     
277.     else
     
278.       { if (pr3<pr1)
     
279.           { r=1.0e+04;
     
280.             if (pr3>1.0e-04) r=1.0/pr3;
     
281.             nqw=nq+1;
     
282.           }
     
283.         else
     
284.           { r=1.0e+04;
     
285.             if (pr1>1.0e-04) r=1.0/pr1;
     
286.             nqw=nq-1;
     
287.           }
     
288.       }
     
289.     idb=10;
     
290.     if (kf==1)
     
291.       if (r<1.1)
     
292.         { for (i=0; i<=n-1; i++)
     
293.             if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
     
294.           jt=nq; goto l20;
     
295.         }
     
296.     if (nqw>nq)
     
297.       for (i=0; i<=n-1; i++)
     
298.         y[i*8+nqw]=er[i]*aa[m-1]/(m+0.0);
     
299.     m=nqw+1;
     
300.     if (kf==1)
     
301.       { irt=2; rr=hmax/fabs(h);
     
302.         if (r>rr) r=rr;
     
303.         h=h*r; hw=h;
     
304.         if (nq==nqw)
     
305.           { r1=1.0;
     
306.             for (j=1; j<=m-1; j++)
     
307.               { r1=r1*r;
     
308.                 for (i=0; i<=n-1; i++)
     
309.                   y[i*8+j]=y[i*8+j]*r1;
     
310.               }
     
311.             idb=m;
     
312.             for (i=0; i<=n-1; i++)
     
313.               if (ym[i]<fabs(y[i*8])) ym[i]=fabs(y[i*8]);
     
314.             jt=nq; goto l20;
     
315.           }
     
316.         nq=nqw;
     
317.         goto l60;
     
318.       }
     
319.     rm=rm*r; irt1=3;
     
320.     rr=fabs(hmin/hd);
     
321.     if (rm<rr) rm=rr;
     
322.     rr=fabs(hmax/hd);
     
323.     if (rm>rr) rm=rr;
     
324.     r=1.0;
     
325.     for (j=1; j<=m-1; j++)
     
326.       { r=r*rm;
     
327.         for (i=0; i<=n-1; i++)
     
328.           y[i*8+j]=s[i*10+j]*r;
     
329.       }
     
330.     h=hd*rm;
     
331.     for (i=0; i<=n-1; i++)
     
332.       y[i*8]=s[i*10];
     
333.     idb=m;
     
334.     if (nqw==nq) goto l80;
     
335.     nq=nqw; goto l60;
     
336.   }

复制代码

4rinv.c

1.   #include "stdlib.h"
   
2.   #include "math.h"
   
3.   #include "stdio.h"
   
4.   int rinv(a,n)
   
5.   int n;
   
6.   double a[];
   
7.   { int *is,*js,i,j,k,l,u,v;
   
8.     double d,p;
   
9.     is=malloc(n*sizeof(int));
   
10.     js=malloc(n*sizeof(int));
    
11.     for (k=0; k<=n-1; k++)
    
12.       { d=0.0;
    
13.         for (i=k; i<=n-1; i++)
    
14.         for (j=k; j<=n-1; j++)
    
15.           { l=i*n+j; p=fabs(a[l]);
    
16.             if (p>d) { d=p; is[k]=i; js[k]=j;}
    
17.           }
    
18.         if (d+1.0==1.0)
    
19.           { free(is); free(js); printf("err**not inv\n");
    
20.             return(0);
    
21.           }
    
22.         if (is[k]!=k)
    
23.           for (j=0; j<=n-1; j++)
    
24.             { u=k*n+j; v=is[k]*n+j;
    
25.               p=a[u]; a[u]=a[v]; a[v]=p;
    
26.             }
    
27.         if (js[k]!=k)
    
28.           for (i=0; i<=n-1; i++)
    
29.             { u=i*n+k; v=i*n+js[k];
    
30.               p=a[u]; a[u]=a[v]; a[v]=p;
    
31.             }
    
32.         l=k*n+k;
    
33.         a[l]=1.0/a[l];
    
34.         for (j=0; j<=n-1; j++)
    
35.           if (j!=k)
    
36.             { u=k*n+j; a[u]=a[u]*a[l];}
    
37.         for (i=0; i<=n-1; i++)
    
38.           if (i!=k)
    
39.             for (j=0; j<=n-1; j++)
    
40.               if (j!=k)
    
41.                 { u=i*n+j;
    
42.                   a[u]=a[u]-a[i*n+k]*a[k*n+j];
    
43.                 }
    
44.         for (i=0; i<=n-1; i++)
    
45.           if (i!=k)
    
46.             { u=i*n+k; a[u]=-a[u]*a[l];}
    
47.       }
    
48.     for (k=n-1; k>=0; k--)
    
49.       { if (js[k]!=k)
    
50.           for (j=0; j<=n-1; j++)
    
51.             { u=k*n+j; v=js[k]*n+j;
    
52.               p=a[u]; a[u]=a[v]; a[v]=p;
    
53.             }
    
54.         if (is[k]!=k)
    
55.           for (i=0; i<=n-1; i++)
    
56.             { u=i*n+k; v=i*n+is[k];
    
57.               p=a[u]; a[u]=a[v]; a[v]=p;
    
58.             }
    
59.       }
    
60.     free(is); free(js);
    
61.     return(1);
    
62.   }

复制代码

zmf_0218

谢谢，仔细研究一番！