lsoda.c

GESPI 2.0动态系统模拟工具

   yp1 = yh[1];
   for ( i = 1 ; i <= n ; i++ )
      y[i] = yp1[i];
   (*f)( neq, tn, y, savf );
   nfe++;
/*
   If indicated, the matrix P = I - h * el[1] * J is reevaluated and
   preprocessed before starting the corrector iteration.  ipup is set
   to 0 as an indicator that this has been done.
*/
   while ( 1 ) {
      if ( *m == 0 ) {
	 if ( ipup > 0 ) {
	    prja( neq, y, f );
	    ipup = 0;
	    rc = 1.;
	    nslp = nst;
	    crate = 0.7;
	    if ( ierpj != 0 ) {
	       corfailure( told, rh, ncf, corflag );
	       return;
	    }
	 }
	 for ( i = 1 ; i <= n ; i++ )
	    acor[i] = 0.;
      }   /*   end if ( *m == 0 )   */
      if ( miter == 0 ) {
/*
   In case of functional iteration, update y directly from
   the result of the last function evaluation.
*/
	 yp1 = yh[2];
	 for ( i = 1 ; i <= n ; i++ ) {
	    savf[i] = h * savf[i] - yp1[i];
	    y[i] = savf[i] - acor[i];
	 }
	 *del = vmnorm( n, y, ewt );
	 yp1 = yh[1];
	 for ( i = 1 ; i <= n ; i++ ) {
	    y[i] = yp1[i] + el[1] * savf[i];
	    acor[i] = savf[i];
	 }
      }      /*   end functional iteration   */
/*
   In the case of the chord method, compute the corrector error,
   and solve the linear system with that as right-hand side and
   P as coefficient matrix.
*/
      else {
	 yp1 = yh[2];
	 for ( i = 1 ; i <= n ; i++ )
	    y[i] = h * savf[i] - ( yp1[i] + acor[i] );
	 solsy( y );
	 *del = vmnorm( n, y, ewt );
	 yp1 = yh[1];
	 for ( i = 1 ; i <= n ; i++ ) {
	    acor[i] += y[i];
	    y[i] = yp1[i] + el[1] * acor[i];
	 }
      }   /*   end chord method   */
/*
   Test for convergence.  If *m > 0, an estimate of the convergence
   rate constant is stored in crate, and this is used in the test.

   We first check for a change of iterates that is the size of
   roundoff error.  If this occurs, the iteration has converged, and a
   new rate estimate is not formed.
   In all other cases, force at least two iterations to estimate a
   local Lipschitz constant estimate for Adams method.
   On convergence, form pdest = local maximum Lipschitz constant
   estimate.  pdlast is the most recent nonzero estimate.
*/
      if ( *del <= 100. * pnorm * ETA )
	 break;
      if ( *m != 0 || meth != 1 ) {
	 if ( *m != 0 ) {
	    rm = 1024.0;
	    if ( *del <= ( 1024. * *delp ) )
	       rm = *del / *delp;
	    rate = max( rate, rm );
	    crate = max( 0.2 * crate, rm );
	 }
	 dcon = *del * min( 1., 1.5 * crate ) / ( tesco[nq][2] * conit );
	 if ( dcon <= 1. ) {
	    pdest = max( pdest, rate / fabs( h * el[1] ) );
	    if ( pdest != 0. )
	       pdlast = pdest;
	    break;
	 }
      }
/*
   The corrector iteration failed to converge.
   If miter != 0 and the Jacobian is out of date, prja is called for
   the next try.   Otherwise the yh array is retracted to its values
   before prediction, and h is reduced, if possible.  If h cannot be
   reduced or mxncf failures have occured, exit with corflag = 2.
*/
      (*m)++;
      if ( *m == maxcor || ( *m >= 2 && *del > 2. * *delp ) ) {
	 if ( miter == 0 || jcur == 1 ) {
	    corfailure( told, rh, ncf, corflag );
	    return;
	 }
	 ipup = miter;
/*
   Restart corrector if Jacobian is recomputed.
*/
	 *m = 0;
	 rate = 0.;
	 *del = 0.;
	 yp1 = yh[1];
	 for ( i = 1 ; i <= n ; i++ )
	    y[i] = yp1[i];
	 (*f)( neq, tn, y, savf );
	 nfe++;
      }
/*
   Iterate corrector.
*/
      else {
	 *delp = *del;
	 (*f)( neq, tn, y, savf );
	 nfe++;
      }
   }   /*   end while   */
}       /*   end correction   */


static void
   corfailure( told, rh, ncf, corflag )

int *ncf, *corflag;
double *told, *rh;

{
   int j, i1, i;

   *ncf++;
   rmax = 2.;
   tn = *told;
   for ( j = nq ; j >= 1 ; j-- )
      for ( i1 = j ; i1 <= nq ; i1++ ) {
	 yp1 = yh[i1];
	 yp2 = yh[i1+1];
	 for ( i = 1 ; i <= n ; i++ )
	    yp1[i] -= yp2[i];
      }
   if ( fabs( h ) <= hmin * 1.00001 || *ncf == mxncf ) {
      *corflag = 2;
      return;
   }
   *corflag = 1;
   *rh = 0.25;
   ipup = miter;

}              /*   end corfailure   */


static void
   solsy( y )

double *y;

/*
   This routine manages the solution of the linear system arising from
   a chord iteration.  It is called if miter != 0.
   If miter is 2, it calls dgesl to accomplish this.
   If miter is 5, it calls dgbsl.

   y = the right-hand side vector on input, and the solution vector
       on output.
*/


{
   iersl = 0;
   if ( miter != 2 ) {
      if ( prfl )
	 printf( "\nsolsy -- miter != 2" );
      return;
   }

   if ( miter == 2 )
      dgesl( wm, n, ipvt, y, 0 );
   return;

}          /*   end solsy   */


static void methodswitch( double dsm, double pnorm, double *pdh, double *rh )
{
   int lm1, lm1p1, lm2, lm2p1, nqm1, nqm2;
   double rh1, rh2, rh1it, exm2, dm2, exm1, dm1, alpha, exsm;

/*
   We are current using an Adams method.  Consider switching to bdf.
   If the current order is greater than 5, assume the problem is
   not stiff, and skip this section.
   If the Lipschitz constant and error estimate are not polluted
   by roundoff, perform the usual test.
   Otherwise, switch to the bdf methods if the last step was
   restricted to insure stability ( irflag = 1 ), and stay with Adams
   method if not.  When switching to bdf with polluted error estimates,
   in the absence of other information, double the step size.

   When the estimates are ok, we make the usual test by computing
   the step size we could have (ideally) used on this step,
   with the current (Adams) method, and also that for the bdf.
   If nq > mxords, we consider changing to order mxords on switching.
   Compare the two step sizes to decide whether to switch.
   The step size advantage must be at least ratio = 5 to switch.
*/
   if ( meth == 1 ) {
      if ( nq > 5 )
	 return;
      if ( dsm <= ( 100. * pnorm * ETA ) || pdest == 0. ) {
	 if ( irflag == 0 )
	    return;
	 rh2 = 2.;
	 nqm2 = min( nq, mxords );
      }
      else {
	 exsm = 1. / ( double ) l;
	 rh1 = 1. / ( 1.2 * pow( dsm, exsm ) + 0.0000012 );
	 rh1it = 2. * rh1;
	 *pdh = pdlast * fabs( h );
	 if ( ( *pdh * rh1 ) > 0.00001 )
	    rh1it = sm1[nq] / *pdh;
	 rh1 = min( rh1, rh1it );
	 if ( nq > mxords ) {
	    nqm2 = mxords;
	    lm2 = mxords + 1;
	    exm2 = 1. / ( double ) lm2;
	    lm2p1 = lm2 + 1;
	    dm2 = vmnorm( n, yh[lm2p1], ewt ) / cm2[mxords];
	    rh2 = 1. / ( 1.2 * pow( dm2, exm2 ) + 0.0000012 );
	 }
	 else {
	    dm2 = dsm * ( cm1[nq] / cm2[nq] );
	    rh2 = 1. / ( 1.2 * pow( dm2, exsm ) + 0.0000012 );
	    nqm2 = nq;
	 }
	 if ( rh2 < ratio * rh1 )
	    return;
      }
/*
   The method switch test passed.  Reset relevant quantities for bdf.
*/
      *rh = rh2;
      icount = 20;
      meth = 2;
      miter = jtyp;
      pdlast = 0.;
      nq = nqm2;
      l = nq + 1;
      return;
   }   /*   end if ( meth == 1 )   */
/*
   We are currently using a bdf method, considering switching to Adams.
   Compute the step size we could have (ideally) used on this step,
   with the current (bdf) method, and also that for the Adams.
   If nq > mxordn, we consider changing to order mxordn on switching.
   Compare the two step sizes to decide whether to switch.
   The step size advantage must be at least 5/ratio = 1 to switch.
   If the step size for Adams would be so small as to cause
   roundoff pollution, we stay with bdf.
*/
   exsm = 1. / ( double ) l;
   if ( mxordn < nq ) {
      nqm1 = mxordn;
      lm1 = mxordn + 1;
      exm1 = 1. / ( double ) lm1;
      lm1p1 = lm1 + 1;
      dm1 = vmnorm( n, yh[lm1p1], ewt ) / cm1[mxordn];
      rh1 = 1. / ( 1.2 * pow( dm1, exm1 ) + 0.0000012 );
   }
   else {
      dm1 = dsm * ( cm2[nq] / cm1[nq] );
      rh1 = 1. / ( 1.2 * pow( dm1, exsm ) + 0.0000012 );
      nqm1 = nq;
      exm1 = exsm;
   }
   rh1it = 2. * rh1;
   *pdh = pdnorm * fabs( h );
   if ( ( *pdh * rh1 ) > 0.00001 )
      rh1it = sm1[nqm1] / *pdh;
   rh1 = min( rh1, rh1it );
   rh2 = 1. / ( 1.2 * pow( dsm, exsm ) + 0.0000012 );
   if ( ( rh1 * ratio ) < ( 5. * rh2 ) )
      return;
   alpha = max( 0.001, rh1 );
   dm1 *= pow( alpha, exm1 );
   if ( dm1 <= 1000. * ETA * pnorm )
      return;
/*
   The switch test passed.  Reset relevant quantities for Adams.
*/
   *rh = rh1;
   icount = 20;
   meth = 1;
   miter = 0;
   pdlast = 0.;
   nq = nqm1;
   l = nq + 1;

}     /*   end methodswitch   */


/*
   This routine returns from stoda to lsoda.  Hence freevectors() is
   not executed.
*/

static void endstoda( void )
{
   double r;
   int i;

   r = 1. / tesco[nqu][2];
   for ( i = 1 ; i <= n ; i++ )
      acor[i] *= r;
   hold = h;
   jstart = 1;

}      /*   end endstoda   */


static void
   orderswitch( rhup, dsm, pdh, rh, orderflag )

int *orderflag;
double *rhup, dsm, *pdh, *rh;


/*
   Regardless of the success or failure of the step, factors
   rhdn, rhsm, and rhup are computed, by which h could be multiplied
   at order nq - 1, order nq, or order nq + 1, respectively.
   In the case of a failure, rhup = 0. to avoid an order increase.
   The largest of these is determined and the new order chosen
   accordingly.  If the order is to be increased, we compute one
   additional scaled derivative.

   orderflag = 0  : no change in h or nq,
	       1  : change in h but not nq,
	       2  : change in both h and nq.
*/

{
   int newq, i;
   double exsm, rhdn, rhsm, ddn, exdn, r;

   *orderflag = 0;

   exsm = 1. / ( double ) l;
   rhsm = 1. / ( 1.2 * pow( dsm, exsm ) + 0.0000012 );

   rhdn = 0.;
   if ( nq != 1 ) {
      ddn = vmnorm( n, yh[l], ewt ) / tesco[nq][1];
      exdn = 1. / ( double ) nq;
      rhdn = 1. / ( 1.3 * pow( ddn, exdn ) + 0.0000013 );
   }
/*
   If meth = 1, limit rh accordinfg to the stability region also.
*/
   if ( meth == 1 ) {
      *pdh = max( fabs( h ) * pdlast, 0.000001 );
      if ( l < lmax )
	 *rhup = min( *rhup, sm1[l] / *pdh );
      rhsm = min( rhsm, sm1[nq] / *pdh );
      if ( nq > 1 )
	 rhdn = min( rhdn, sm1[nq-1] / *pdh );
      pdest = 0.;
   }
   if ( rhsm >= *rhup ) {
      if ( rhsm >= rhdn ) {
	 newq = nq;
	 *rh = rhsm;
      }
      else {
	 newq = nq - 1;
	 *rh = rhdn;
	 if ( kflag < 0 && *rh > 1. )
	    *rh = 1.;
      }
   }
   else {
      if ( *rhup <= rhdn ) {
	 newq = nq - 1;
	 *rh = rhdn;
	 if ( kflag < 0 && *rh > 1. )
	    *rh = 1.;
      }
      else {
	 *rh = *rhup;
	 if ( *rh >= 1.1 ) {
	    r = el[l] / ( double ) l;
	    nq = l;
	    l = nq + 1;
	    yp1 = yh[l];
	    for ( i = 1 ; i <= n ; i++ )
	       yp1[i] = acor[i] * r;
	    *orderflag = 2;
	    return;
	 }
	 else {
	    ialth = 3;
	    return;
	 }
      }
   }
/*
   If meth = 1 and h is restricted by stability, bypass 10 percent test.
*/
   if ( meth == 1 ) {
      if ( ( *rh * *pdh * 1.00001 ) < sm1[newq] )
	 if ( kflag == 0 && *rh < 1.1 ) {
	    ialth = 3;
	    return;
	 }
   }
   else {
      if ( kflag == 0 && *rh < 1.1 ) {
	 ialth = 3;
	 return;
      }
   }
   if ( kflag <= -2 )
      *rh = min( *rh, 0.2 );
/*
   If there is a change of order, reset nq, l, and the coefficients.
   In any case h is reset according to rh and the yh array is rescaled.
   Then exit or redo the step.
*/
   if ( newq == nq ) {
      *orderflag = 1;
      return;
   }
   nq = newq;
   l = nq + 1;
   *orderflag = 2;

}      /*   end orderswitch   */


static void resetcoeff( void )
/*
   The el vector and related constants are reset
   whenever the order nq is changed, or at the start of the problem.
*/
{
 int i;
 double *ep1;

 ep1 = elco[nq];
 for ( i=1; i<=l; i++ )
  el[i] = ep1[i];
 rc = rc * el[1] / el0;
 el0 = el[1];
 conit = 0.5 / ( double ) ( nq + 2 );
}     /*   end resetcoeff   */