/* Adams-Bashforth-Moulton integration formulas. Reference: Shampine, L. F. and M. K. Gordon, _Computer Solution of Ordinary Differential Equations_, W. H. Freeman, 1975. Program by Steve Moshier. */ # include # include # include "ode_moshier.h" /* Divided differences */ #define N 19 /* Predictor coefficients */ double precof[N] = { 1.0, 1.0 / 2.0, 5.0 / 12.0, 3.0 / 8.0, 251.0 / 720.0, 95.0 / 288.0, 19087.0 / 60480.0, 5257.0 / 17280.0, 1070017.0 / 3628800.0, 25713.0 / 89600.0, 26842253.0 / 95800320.0, 4777223.0 / 17418240.0, 703604254357.0 / 2615348736000.0, 106364763817.0 / 402361344000.0, 1166309819657.0 / 4483454976000.0, 2.5221445e7 / 9.8402304e7, 8.092989203533249e15 / 3.201186852864e16, 8.5455477715379e13 / 3.4237292544e14, 1.2600467236042756559e19 / 5.109094217170944e19, }; /* Corrector coefficients */ double corcof[N] = { 1.0, -1.0 / 2.0, -1.0 / 12.0, -1.0 / 24.0, -19.0 / 720.0, -3.0 / 160.0, -863.0 / 60480.0, -275.0 / 24192.0, -33953.0 / 3628800.0, -8183.0 / 1036800.0, -3250433.0 / 479001600.0, -4671.0 / 788480.0, -13695779093.0 / 2615348736000.0, -2224234463.0 / 475517952000.0, -132282840127.0 / 31384184832000.0, -2639651053.0 / 689762304000.0, 1.11956703448001e14 / 3.201186852864e16, 5.0188465e7 / 1.5613165568e10, 2.334028946344463e15 / 7.86014494949376e17, }; /******************************************************************************/ void divdif ( double vec[], int k, double *diffn ) /******************************************************************************/ /* Compute zeroth through kth backward differences of the data in the input array */ /* double vec[], input array of k+1 data items; int k; double *diffn, output array of ith differences; */ { double diftbl[N]; double *p, *q; double y; int i, o; /* Copy the given data (zeroth difference) into temp array */ p = diftbl; q = vec; for( i=0; i<=k; i++ ) *p++ = *q++; /* On the first outer loop, k-1 first differences are calculated. * These overwrite the original data in the temp array. */ o = k; for( o=k; o>0; o-- ) { p = diftbl; q = p; for( i=0; i c D f . * n+1 n - i n * i=0 * * This subroutine returns the summation term, not multiplied by h. * The coefficients c_i of the integration formula are either * precof[] or corcof[], given above. */ /******************************************************************************/ double intpol( double *diffn, double *coeffs, int k ) /******************************************************************************/ /* double *diffn; array of backward differences double *coeffs; coefficients of integration formula int k; differences used go up to order k-1 */ { double s; int i; s = 0.0; coeffs += k; diffn += k; for( i=0; i= nequat in adstep() below. If neq > nequat, unused space * is left for extra variables that are not actually integrated. */ /******************************************************************************/ void adstart ( double h, double yn[], double work[], int neq, int ord, double t ) /******************************************************************************/ { double *p; int j; hstep = h; ccor = hstep * precof[ord]; jstep = 0; dsiz = ord + 2; asiz = neq * dsiz; order = ord; ordp1 = ord + 1; p = work; dv = p; p += asiz; dvp = p; p += dsiz; vp = p; p += neq; yp = p; p += neq; vn = p; p += neq; delta = p; p += neq; sdelta = p; p += neq; y0 = p; func ( t, yn, vn ); p = dv; for( j=0; j e ) e = e0; pdv += dsiz; } done: time += hstep; *t = time; return( e ); }