lsoda.bak

GESPI 2.0动态系统模拟工具

            wm[i][j] = ( acor[i] - savf[i] ) * fac;
         y[j] = yj;
      }
      nfe += n;
/*
   Compute norm of Jacobian.
*/
      pdnorm = fnorm( n, wm, ewt ) / fabs( hl0 );
/*
   Add identity matrix.
*/
      for ( i = 1 ; i <= n ; i++ )
         wm[i][i] += 1.;
/*
   Do LU decomposition on P.
*/
      dgefa( wm, n, ipvt, &ier );
      if ( ier != 0 )
         ierpj = 1;
      return;
   }

}      /*   end prja   */


static double
   vmnorm( n, v, w )

int n;
double *v, *w;

/*
   This function routine computes the weighted max-norm
   of the vector of length n contained in the array v, with weights
   contained in the array w of length n.

   vmnorm = max( i = 1, ..., n ) fabs( v[i] ) * w[i].
*/

{
   int i;
   double vm;

   vm = 0.;
   for ( i = 1 ; i <= n ; i++ )
      vm = max( vm, fabs( v[i] ) * w[i] );
   return vm;

}                  /*   end vmnorm   */


static double
   fnorm( n, a, w )

int n;
double **a, *w;

/*
   This subroutine computes the norm of a full n by n matrix,
   stored in the array a, that is consistent with the weighted max-norm
   on vectors, with weights stored in the array w.

      fnorm = max(i=1,...,n) ( w[i] * sum(j=1,...,n) fabs( a[i][j] ) / w[j] )
*/

{
   int i, j;
   double an, sum, *ap1;

   an = 0.;
   for ( i = 1; i <= n ; i++ ) {
      sum = 0.;
      ap1 = a[i];
      for ( j = 1 ; j <= n ; j++ )
         sum += fabs( ap1[j] ) / w[j];
      an = max( an, sum * w[i] );
   }
   return an;

}     /*   end fnorm   */


static double
   bnorm()

{
}   /*   end bnorm   */

static void correction( int neq, double *y, void (*f)(), int *corflag, 
                        double pnorm, double *del, double *delp,
                        double *told, int *ncf, double *rh, int *m )
/*
   *corflag = 0 : corrector converged,
              1 : step size to be reduced, redo prediction,
              2 : corrector cannot converge, failure flag.
*/
{
   int i;
   double rm, rate, dcon;

/*
   Up to maxcor corrector iterations are taken.  A convergence test is
   made on the r.m.s. norm of each correction, weighted by the error
   weight vector ewt.  The sum of the corrections is accumulated in the
   vector acor[i].  The yh array is not altered in the corrector loop.
*/

   *m = 0;
   *corflag = 0;
   rate = 0.;
   *del = 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   */


st