lsoda.bak

GESPI 2.0动态系统模拟工具

            }
         }            /*   end if ( ialth == 0 )   */
         if ( ialth > 1 || l == lmax ) {
            endstoda();
            break;
         }
         yp1 = yh[lmax];
         for ( i = 1 ; i <= n ; i++ )
            yp1[i] = acor[i];
         endstoda();
         break;
      }       /*   end if ( dsm <= 1. )   */
/*
   The error test failed.  kflag keeps track of multiple failures.
   Restore tn and the yh array to their previous values, and prepare
   to try the step again.  Compute the optimum step size for this or
   one lower.  After 2 or more failures, h is forced to decrease
   by a factor of 0.2 or less.
*/
      else {
         kflag--;
         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];
            }
         rmax = 2.;
         if ( fabs( h ) <= hmin * 1.00001 ) {
            kflag = -1;
            hold = h;
            jstart = 1;
            break;
         }
         if ( kflag > -3 ) {
            rhup = 0.;
            orderswitch( &rhup, dsm, &pdh, &rh, &orderflag );
            if ( orderflag == 1 || orderflag == 0 ) {
               if ( orderflag == 0 )
                  rh = min( rh, 0.2 );
               rh = max( rh, hmin / fabs( h ) );
               scaleh( &rh, &pdh );
            }
            if ( orderflag == 2 ) {
               resetcoeff();
               rh = max( rh, hmin / fabs( h ) );
               scaleh( &rh, &pdh );
            }
            continue;
         }     /*   if ( kflag > -3 )   */
/*
   Control reaches this section if 3 or more failures have occurred.
   If 10 failures have occurred, exit with kflag = -1.
   It is assumed that the derivatives that have accumulated in the
   yh array have errors of the wrong order.  Hence the first
   derivative is recomputed, and the order is set to 1.  Then
   h is reduced by a factor of 10, and the step is retried,
   until it succeeds or h reaches hmin.
*/
         else {
            if ( kflag == -10 ) {
               kflag = -1;
               hold = h;
               jstart = 1;
               break;
            }
            else {
               rh = 0.1;
               rh = max( hmin / fabs( h ) , rh );
               h *= rh;
               yp1 = yh[1];
               for ( i = 1 ; i <= n ; i++ )
                  y[i] = yp1[i];
               (*f)( neq, tn, y, savf );
               nfe++;
               yp1 = yh[2];
               for ( i = 1 ; i <= n ; i++ )
                  yp1[i] = h * savf[i];
               ipup = miter;
               ialth = 5;
               if ( nq == 1 )
                  continue;
               nq = 1;
               l = 2;
               resetcoeff();
               continue;
            }
         }     /*   end else -- kflag <= -3 */
      }     /*   end error failure handling   */
   }      /*   end outer while   */

}           /*   end stoda   */


static void
   ewset( itol, rtol, atol, ycur )

int itol;
double *rtol, *atol, *ycur;

{
   int i;

   switch ( itol ) {
   case 1 :
      for ( i = 1 ; i <= n ; i++ )
         ewt[i] = rtol[1] * fabs( ycur[i] ) + atol[1];
      break;
   case 2 :
      for ( i = 1 ; i <= n ; i++ )
         ewt[i] = rtol[1] * fabs( ycur[i] ) + atol[i];
      break;
   case 3 :
      for ( i = 1 ; i <= n ; i++ )
         ewt[i] = rtol[i] * fabs( ycur[i] ) + atol[1];
      break;
   case 4 :
      for ( i = 1 ; i <= n ; i++ )
         ewt[i] = rtol[i] * fabs( ycur[i] ) + atol[i];
      break;
   }

}           /*   end ewset   */


void xsetf( mflag )
int mflag;

/*
   xsetf sets the flag to control the printing of messages by lsoda.
   mflag = 0 means do not print (danger.. this risks losing valuable
             information).
   mflag = 1 means print (the default)
   A call to xsetf may be made at any time and will take effect
   immediately.
*/

{

   prfl = mflag;

}           /*   end xsetf   */

void intdy( t, k, dky, iflag )
int k, *iflag;
double *dky, t;

/*
   intdy computes interpolated values of the k-th derivative of the
   dependent variable vector y, and stores it in dky.  This routine
   is called within the package with k = 0 and *t = tout, but may
   also be called by the user for any k up to the current order.
   ( See detailed instructions in the usage documentation. )

   The computed values in dky are gotten by interpolation using the
   Nordsieck history array yh.  This array corresponds uniquely to a
   vector-valued polynomial of degree nqcur or less, and dky is set
   to the k-th derivative of this polynomial at t.
   The formula for dky is

             q
   dky[i] = sum c[k][j] * ( t - tn )^(j-k) * h^(-j) * yh[j+1][i]
            j=k

   where c[k][j] = j*(j-1)*...*(j-k+1), q = nqcur, tn = tcur, h = hcur.
   The quantities nq = nqcur, l = nq+1, n = neq, tn, and h are declared
   static globally.  The above sum is done in reverse order.
   *iflag is returned negative if either k or t is out of bounds.
*/

{
   int i, ic, j, jj, jp1;
   double c, r, s, tp;

   *iflag = 0;
   if ( k < 0 || k > nq ) {
      if ( prfl )
         printf( "\nintdy -- k = %d illegal", k );
      *iflag = -1;
      return;
   }
   tp = tn - hu - 100. * ETA * ( tn + hu );
   if ( ( t - tp ) * ( t - tn ) > 0. ) {
      if ( prfl ) {
         printf( "\nintdy -- t = %g illegal", t );
         printf( "\n         t not in interval tcur - hu to tcur" );
      }
      *iflag = -2;
      return;
   }

   s = ( t - tn ) / h;
   ic = 1;
   for ( jj = l - k ; jj <= nq ; jj++ )
      ic *= jj;
   c = ( double ) ic;
   yp1 = yh[l];
   for ( i = 1 ; i <= n ; i++ )
      dky[i] = c * yp1[i];
   for ( j = nq - 1 ; j >= k ; j-- ) {
      jp1 = j + 1;
      ic = 1;
      for ( jj = jp1 - k ; jj <= j ; jj++ )
         ic *= jj;
      c = ( double ) ic;
      yp1 = yh[jp1];
      for ( i = 1 ; i <= n ; i++ )
         dky[i] = c * yp1[i] + s * dky[i];
   }
   if ( k == 0 )
      return;
   r = pow( h, ( double ) ( -k ) );
   for ( i = 1 ; i <= n ; i++ )
      dky[i] *= r;

}      /*   end intdy   */


static void
   cfode( meth )

int meth;

{
   int i, nq, nqm1, nqp1;
   double agamq, fnq, fnqm1, pc[13], pint, ragq,
          rqfac, rq1fac, tsign, xpin;
/*
   cfode is called by the integrator routine to set coefficients
   needed there.  The coefficients for the current method, as
   given by the value of meth, are set for all orders and saved.
   The maximum order assumed here is 12 if meth = 1 and 5 if meth = 2.
   ( A smaller value of the maximum order is also allowed. )
   cfode is called once at the beginning of the problem, and
   is not called again unless and until meth is changed.

   The elco array contains the basic method coefficients.
   The coefficients el[i], 1 < i < nq+1, for the method of
   order nq are stored in elco[nq][i].  They are given by a generating
   polynomial, i.e.,

      l(x) = el[1] + el[2]*x + ... + el[nq+1]*x^nq.

   For the implicit Adams method, l(x) is given by

      dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1),   l(-1) = 0.

   For the bdf methods, l(x) is given by

      l(x) = (x+1)*(x+2)*...*(x+nq)/k,

   where   k = factorial(nq)*(1+1/2+...+1/nq).

   The tesco array contains test constants used for the
   local error test and the selection of step size and/or order.
   At order nq, tesco[nq][k] is used for the selection of step
   size at order nq-1 if k = 1, at order nq if k = 2, and at order
   nq+1 if k = 3.
*/
   if ( meth == 1 ) {
      elco[1][1] = 1.;
      elco[1][2] = 1.;
      tesco[1][1] = 0.;
      tesco[1][2] = 2.;
      tesco[2][1] = 1.;
      tesco[12][3] = 0.;
      pc[1] = 1.;
      rqfac = 1.;
      for ( nq = 2 ; nq <= 12 ; nq++ ) {
/*
   The pc array will contain the coefficients of the polynomial

      p(x) = (x+1)*(x+2)*...*(x+nq-1).

   Initially, p(x) = 1.
*/
         rq1fac = rqfac;
         rqfac = rqfac / ( double ) nq;
         nqm1 = nq - 1;
         fnqm1 = ( double ) nqm1;
         nqp1 = nq + 1;
/*
   Form coefficients of p(x)*(x+nq-1).
*/
         pc[nq] = 0.;
         for ( i = nq ; i >= 2 ; i-- )
            pc[i] = pc[i-1] + fnqm1 * pc[i];
         pc[1] = fnqm1 * pc[1];
/*
   Compute integral, -1 to 0, of p(x) and x*p(x).
*/
         pint = pc[1];
         xpin = pc[1] / 2.;
         tsign = 1.;
         for ( i = 2 ; i <= nq ; i++ ) {
            tsign = -tsign;
            pint += tsign * pc[i] / ( double ) i;
            xpin += tsign * pc[i] / ( double ) ( i + 1 );
         }
/*
   Store coefficients in elco and tesco.
*/
         elco[nq][1] = pint * rq1fac;
         elco[nq][2] = 1.;
         for ( i = 2 ; i <= nq ; i++ )
            elco[nq][i+1] = rq1fac * pc[i] / ( double ) i;
         agamq = rqfac * xpin;
         ragq = 1. / agamq;
         tesco[nq][2] = ragq;
         if ( nq < 12 )
            tesco[nqp1][1] = ragq * rqfac / ( double ) nqp1;
         tesco[nqm1][3] = ragq;
      }      /*   end for   */
      return;
   }      /*   end if ( meth == 1 )   */

/*
   meth = 2.
*/
   pc[1] = 1.;
   rq1fac = 1.;
/*
   The pc array will contain the coefficients of the polynomial

      p(x) = (x+1)*(x+2)*...*(x+nq).

   Initially, p(x) = 1.
*/
   for ( nq = 1 ; nq <= 5 ; nq++ ) {
      fnq = ( double ) nq;
      nqp1 = nq + 1;
/*
   Form coefficients of p(x)*(x+nq).
*/
      pc[nqp1] = 0.;
      for ( i = nq + 1 ; i >= 2 ; i-- )
         pc[i] = pc[i-1] + fnq * pc[i];
      pc[1] *= fnq;
/*
   Store coefficients in elco and tesco.
*/
      for ( i = 1 ; i <= nqp1 ; i++ )
         elco[nq][i] = pc[i] / pc[2];
      elco[nq][2] = 1.;
      tesco[nq][1] = rq1fac;
      tesco[nq][2] = ( ( double ) nqp1 ) / elco[nq][1];
      tesco[nq][3] = ( ( double ) ( nq + 2 ) ) / elco[nq][1];
      rq1fac /= fnq;
   }
   return;

}       /*   end cfode   */


static void
   scaleh( rh, pdh )

double *rh, *pdh;

{
   double r;
   int j, i;
/*
   If h is being changed, the h ratio rh is checked against
   rmax, hmin, and hmxi, and the yh array is rescaled.  ialth is set to
   l = nq + 1 to prevent a change of h for that many steps, unless
   forced by a convergence or error test failure.
*/
   *rh = min( *rh, rmax );
   *rh = *rh / max( 1., fabs( h ) * hmxi * *rh );
/*
   If meth = 1, also restrict the new step size by the stability region.
   If this reduces h, set irflag to 1 so that if there are roundoff
   problems later, we can assume that is the cause of the trouble.
*/
   if ( meth == 1 ) {
      irflag = 0;
      *pdh = max( fabs( h ) * pdlast, 0.000001 );
      if ( ( *rh * *pdh * 1.00001 ) >= sm1[nq] ) {
         *rh = sm1[nq] / *pdh;
         irflag = 1;
      }
   }
   r = 1.;
   for ( j = 2 ; j <= l ; j++ ) {
      r *= *rh;
      yp1 = yh[j];
      for ( i = 1 ; i <= n ; i++ )
         yp1[i] *= r;
   }
   h *= *rh;
   rc *= *rh;
   ialth = l;

}     /*   end scaleh   */


static void prja( int neq, double *y, void (*f)() )
{
   int i, i1, i2, ier, ii, j, j1, jj, lenp,
       mba, mband, meb1, meband, ml3, np1;
   double con, fac, hl0, r, r0, yi, yj, yjj;
/*
   prja is called by stoda to compute and process the matrix
   P = I - h * el[1] * J, where J is an approximation to the Jacobian.
   Here J is computed by finite differencing.
   J, scaled by -h * el[1], is stored in wm.  Then the norm of J ( the
   matrix norm consistent with the weighted max-norm on vectors given
   by vmnorm ) is computed, and J is overwritten by P.  P is then
   subjected to LU decomposition in preparation for later solution
   of linear systems with p as coefficient matrix.  This is done
   by dgefa if miter = 2, and by dgbfa if miter = 5.
*/
   nje++;
   ierpj = 0;
   jcur = 1;
   hl0 = h * el0;
/*
   If miter = 2, make n calls to f to approximate J.
*/
   if ( miter != 2 ) {
      if ( prfl )
         printf( "\nprja -- miter != 2" );
      return;
   }

   if ( miter == 2 ) {
      fac = vmnorm( n, savf, ewt );
      r0 = 1000. * fabs( h ) * ETA * ( ( double ) n ) * fac;
      if ( r0 == 0. )
         r0 = 1.;
      for ( j = 1 ; j <= n ; j++ ) {
         yj = y[j];
         r = max( sqrteta * fabs( yj ), r0 / ewt[j] );
         y[j] += r;
         fac = -hl0 / r;
         (*f)( neq, tn, y, acor );
         for ( i = 1 ; i <= n ; i++ )