adams_pece_impl.h

TSTOOL应用软件,内有说明文件,可以适当修改,应用比较方便.

	i__ = *k - j;
/* L125: */
	v[i__] -= alpha[j + 1] * v[i__ + 1];
    }

/*   update v(*) and set w(*) */

L130:
    limit1 = kp1 - *ns;
    temp5 = alpha[*ns];
    i__1 = limit1;
    for (iq = 1; iq <= i__1; ++iq) {
	v[iq] -= temp5 * v[iq + 1];
/* L135: */
	w[iq] = v[iq];
    }
    g[nsp1] = w[1];

/*   compute the g(*) in the work vector w(*) */

L140:
    nsp2 = *ns + 2;
    if (kp1 < nsp2) {
	goto L199;
    }
    i__1 = kp1;
    for (i__ = nsp2; i__ <= i__1; ++i__) {
	limit2 = kp2 - i__;
	temp6 = alpha[i__ - 1];
	i__2 = limit2;
	for (iq = 1; iq <= i__2; ++iq) {
/* L145: */
	    w[iq] -= temp6 * w[iq + 1];
	}
/* L150: */
	g[i__] = w[1];
    }
L199:
/*       ***     end block 1     *** */

/*       ***     begin block 2     *** */
/*   predict a solution p(*), evaluate derivatives using predicted */
/*   solution, estimate local error at order k and errors at orders k, */
/*   k-1, k-2 as if constant step size were used. */
/*                   *** */

/*   change phi to phi star */

    if (*k < nsp1) {
	goto L215;
    }
    i__1 = *k;
    for (i__ = nsp1; i__ <= i__1; ++i__) {
	temp1 = beta[i__];
	i__2 = neqn;
	for (l = 1; l <= i__2; ++l) {
/* L205: */
	    phi[l + i__ * phi_dim1] = temp1 * phi[l + i__ * phi_dim1];
	}
/* L210: */
    }

/*   predict solution and differences */

L215:
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
	phi[l + kp2 * phi_dim1] = phi[l + kp1 * phi_dim1];
	phi[l + kp1 * phi_dim1] = 0.;
/* L220: */
	p[l] = 0.;
    }
    i__1 = *k;
    for (j = 1; j <= i__1; ++j) {
	i__ = kp1 - j;
	ip1 = i__ + 1;
	temp2 = g[i__];
	i__2 = neqn;
	for (l = 1; l <= i__2; ++l) {
	    p[l] += temp2 * phi[l + i__ * phi_dim1];
/* L225: */
	    phi[l + i__ * phi_dim1] += phi[l + ip1 * phi_dim1];
	}
/* L230: */
    }
    if (*nornd) {
	goto L240;
    }
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
	tau = *h__ * p[l] - phi[l + phi_dim1 * 15];
	p[l] = y[l] + tau;
/* L235: */
	phi[l + (phi_dim1 << 4)] = p[l] - y[l] - tau;
    }
    goto L250;
L240:
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
/* L245: */
	p[l] = y[l] + *h__ * p[l];
    }
L250:
    xold = *x;
    *x += *h__;
    absh = abs(*h__);
    f(x, &p[1], &yp[1]);

/*   estimate errors at orders k,k-1,k-2 */

    erkm2 = 0.;
    erkm1 = 0.;
    erk = 0.;
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
	temp3 = 1. / wt[l];
	temp4 = yp[l] - phi[l + phi_dim1];
	if (km2 < 0) {
	    goto L265;
	} else if (km2 == 0) {
	    goto L260;
	} else {
	    goto L255;
	}
L255:
/* Computing 2nd power */
	d__1 = (phi[l + km1 * phi_dim1] + temp4) * temp3;
	erkm2 += d__1 * d__1;
L260:
/* Computing 2nd power */
	d__1 = (phi[l + *k * phi_dim1] + temp4) * temp3;
	erkm1 += d__1 * d__1;
L265:
/* Computing 2nd power */
	d__1 = temp4 * temp3;
	erk += d__1 * d__1;
    }
    if (km2 < 0) {
	goto L280;
    } else if (km2 == 0) {
	goto L275;
    } else {
	goto L270;
    }
L270:
    erkm2 = absh * sig[km1] * gstr[km2 - 1] * sqrt(erkm2);
L275:
    erkm1 = absh * sig[*k] * gstr[km1 - 1] * sqrt(erkm1);
L280:
    temp5 = absh * sqrt(erk);
    err = temp5 * (g[*k] - g[kp1]);
    erk = temp5 * sig[kp1] * gstr[*k - 1];
    knew = *k;

/*   test if order should be lowered */

    if (km2 < 0) {
	goto L299;
    } else if (km2 == 0) {
	goto L290;
    } else {
	goto L285;
    }
L285:
    if (max(erkm1,erkm2) <= erk) {
	knew = km1;
    }
    goto L299;
L290:
    if (erkm1 <= erk * .5) {
	knew = km1;
    }

/*   test if step successful */

L299:
    if (err <= *eps) {
	goto L400;
    }
/*       ***     end block 2     *** */

/*       ***     begin block 3     *** */
/*   the step is unsuccessful.  restore  x, phi(*,*), psi(*) . */
/*   if third consecutive failure, set order to one.  if step fails more */
/*   than three times, consider an optimal step size.  double error */
/*   tolerance and return if estimated step size is too small for machine */
/*   precision. */
/*                   *** */

/*   restore x, phi(*,*) and psi(*) */

    *phase1 = FALSE_;
    *x = xold;
    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
	temp1 = 1. / beta[i__];
	ip1 = i__ + 1;
	i__2 = neqn;
	for (l = 1; l <= i__2; ++l) {
/* L305: */
	    phi[l + i__ * phi_dim1] = temp1 * (phi[l + i__ * phi_dim1] - phi[
		    l + ip1 * phi_dim1]);
	}
/* L310: */
    }
    if (*k < 2) {
	goto L320;
    }
    i__1 = *k;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* L315: */
	psi[i__ - 1] = psi[i__] - *h__;
    }

/*   on third failure, set order to one.  thereafter, use optimal step */
/*   size */

L320:
    ++ifail;
    temp2 = .5;
    if ((i__1 = ifail - 3) < 0) {
	goto L335;
    } else if (i__1 == 0) {
	goto L330;
    } else {
	goto L325;
    }
L325:
    if (p5eps < erk * .25) {
	temp2 = sqrt(p5eps / erk);
    }
L330:
    knew = 1;
L335:
    *h__ = temp2 * *h__;
    *k = knew;
    if (abs(*h__) >= fouru * abs(*x)) {
	goto L340;
    }
    *crash = TRUE_;
    d__1 = fouru * abs(*x);
    *h__ = d_sign(&d__1, h__);
    *eps += *eps;
    return 0;
L340:
    goto L100;
/*       ***     end block 3     *** */

/*       ***     begin block 4     *** */
/*   the step is successful.  correct the predicted solution, evaluate */
/*   the derivatives using the corrected solution and update the */
/*   differences.  determine best order and step size for next step. */
/*                   *** */
L400:
    *kold = *k;
    *hold = *h__;

/*   correct and evaluate */

    temp1 = *h__ * g[kp1];
    if (*nornd) {
	goto L410;
    }
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
	rho = temp1 * (yp[l] - phi[l + phi_dim1]) - phi[l + (phi_dim1 << 4)];
	y[l] = p[l] + rho;
/* L405: */
	phi[l + phi_dim1 * 15] = y[l] - p[l] - rho;
    }
    goto L420;
L410:
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
/* L415: */
	y[l] = p[l] + temp1 * (yp[l] - phi[l + phi_dim1]);
    }
L420:
    f(x, &y[1], &yp[1]);

/*   update differences for next step */

    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
	phi[l + kp1 * phi_dim1] = yp[l] - phi[l + phi_dim1];
/* L425: */
	phi[l + kp2 * phi_dim1] = phi[l + kp1 * phi_dim1] - phi[l + kp2 * 
		phi_dim1];
    }
    i__1 = *k;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = neqn;
	for (l = 1; l <= i__2; ++l) {
/* L430: */
	    phi[l + i__ * phi_dim1] += phi[l + kp1 * phi_dim1];
	}
/* L435: */
    }

/*   estimate error at order k+1 unless: */
/*     in first phase when always raise order, */
/*     already decided to lower order, */
/*     step size not constant so estimate unreliable */

    erkp1 = 0.;
    if (knew == km1 || *k == 12) {
	*phase1 = FALSE_;
    }
    if (*phase1) {
	goto L450;
    }
    if (knew == km1) {
	goto L455;
    }
    if (kp1 > *ns) {
	goto L460;
    }
    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
/* L440: */
/* Computing 2nd power */
	d__1 = phi[l + kp2 * phi_dim1] / wt[l];
	erkp1 += d__1 * d__1;
    }
    erkp1 = absh * gstr[kp1 - 1] * sqrt(erkp1);

/*   using estimated error at order k+1, determine appropriate order */
/*   for next step */

    if (*k > 1) {
	goto L445;
    }
    if (erkp1 >= erk * .5) {
	goto L460;
    }
    goto L450;
L445:
    if (erkm1 <= min(erk,erkp1)) {
	goto L455;
    }
    if (erkp1 >= erk || *k == 12) {
	goto L460;
    }

/*   here erkp1 .lt. erk .lt. dmax1(erkm1,erkm2) else order would have */
/*   been lowered in block 2.  thus order is to be raised */

/*   raise order */

L450:
    *k = kp1;
    erk = erkp1;
    goto L460;

/*   lower order */

L455:
    *k = km1;
    erk = erkm1;

/*   with new order determine appropriate step size for next step */

L460:
    hnew = *h__ + *h__;
    if (*phase1) {
	goto L465;
    }
    if (p5eps >= erk * two[*k]) {
	goto L465;
    }
    hnew = *h__;
    if (p5eps >= erk) {
	goto L465;
    }
    temp2 = (double) (*k + 1);
    d__1 = p5eps / erk;
    d__2 = 1. / temp2;
    r__ = pow_dd(&d__1, &d__2);
/* Computing MAX */
    d__1 = .5, d__2 = min(.9,r__);
    hnew = absh * max(d__1,d__2);
/* Computing MAX */
    d__2 = hnew, d__3 = fouru * abs(*x);
    d__1 = max(d__2,d__3);
    hnew = d_sign(&d__1, h__);
L465:
    *h__ = hnew;
    return 0;
/*       ***     end block 4     *** */
} /* step_ */



template<class ExplicitODE>
int Adams<ExplicitODE>::intrp(const double x, const double *y, const double xout, double *yout, double *ypout, const long neqn, const long kold, 
	double *phi, double *psi)
{
    /* Initialized data */

    static double g[13] = { 1. };
    static double rho[13] = { 1. };

    /* System generated locals */
    long phi_dim1, phi_offset, i__1, i__2;

    /* Local variables */
    double term, temp1, temp2, temp3;
    long i__, j, l;
    double gamma, w[13];
    long limit1;
    double psijm1, hi;
    long ki, jm1;
    double eta;
    long kip1;


/*   the methods in subroutine  step  approximate the solution near  x */
/*   by a polynomial.  subroutine  intrp  approximates the solution at */
/*   xout  by evaluating the polynomial there.  information defining this */
/*   polynomial is passed from  step  so  intrp  cannot be used alone. */

/*   this code is completely explained and documented in the text, */
/*   computer solution of ordinary differential equations:  the initial */
/*   value problem  by l. f. shampine and m. k. gordon. */

/*   input to intrp -- */

/*   all floating point variables are double precision */
/*   the user provides storage in the calling program for the arrays in */
/*   the call list */
/*   and defines */
/*      xout -- point at which solution is desired. */
/*   the remaining parameters are defined in  step  and passed to  intrp */
/*   from that subroutine */

/*   output from  intrp -- */

/*      yout(*) -- solution at  xout */
/*      ypout(*) -- derivative of solution at  xout */
/*   the remaining parameters are returned unaltered from their input */
/*   values.  integration with  step  may be continued. */

    /* Parameter adjustments */
    phi_dim1 = neqn;
    phi_offset = 1 + phi_dim1 * 1;
    phi -= phi_offset;

    --ypout;
    --yout;
    --y;
    --psi;

    /* Function Body */

    hi = xout - x;
    ki = kold + 1;
    kip1 = ki + 1;

/*   initialize w(*) for computing g(*) */

    i__1 = ki;
    for (i__ = 1; i__ <= i__1; ++i__) {
		temp1 = (double) i__;
		w[i__ - 1] = 1. / temp1;
    }
    term = 0.;

/*   compute g(*) */

    i__1 = ki;
	
    for (j = 2; j <= i__1; ++j) {
		jm1 = j - 1;
		psijm1 = psi[jm1];
		gamma = (hi + term) / psijm1;
		eta = hi / psijm1;
		limit1 = kip1 - j;
		i__2 = limit1;
		for (i__ = 1; i__ <= i__2; ++i__) {
	    	w[i__ - 1] = gamma * w[i__ - 1] - eta * w[i__];
		}
		g[j - 1] = w[0];
		rho[j - 1] = gamma * rho[jm1 - 1];
		term = psijm1;
    }

/*   interpolate */

    i__1 = neqn;
    for (l = 1; l <= i__1; ++l) {
		ypout[l] = 0.;
		yout[l] = 0.;
    }
	
    i__1 = ki;
    
	for (j = 1; j <= i__1; ++j) {
		i__ = kip1 - j;
		temp2 = g[i__ - 1];
		temp3 = rho[i__ - 1];
		i__2 = neqn;
		for (l = 1; l <= i__2; ++l) {
	    	yout[l] += temp2 * phi[l + i__ * phi_dim1];
	    	ypout[l] += temp3 * phi[l + i__ * phi_dim1];
		}
    }

    i__1 = neqn;

    for (l = 1; l <= i__1; ++l) {
		yout[l] = y[l] + hi * yout[l];
    }
    return 0;
} /* intrp_ */

double NumericalAlgorithm::machine_precision()
{
	double halfu = 0.5;

	while(((double)(1.0 + halfu)) > 1.0) {
		halfu *= 0.5;
	}

	return 2 * halfu;
}	

#endif