lsoda.c

GESPI 2.0动态系统模拟工具

	       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   */


/*
   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.
*/
void xsetf( int mflag )
{
   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++ )
	    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.;