integrate2.c

二重和一次积分程序

#include <math.h>
#include "nrutil.h"
#define EPS 1.0e-6
#define JMAX 22
#define JMAXP (JMAX+1)
#define K 5
/* Here EPS is the fractional accuracy desired, as determined by the
 * extrapolation error estimate; JMAX limits the total number of steps;
 * K is the number of points used in the extrapolation. */

/* integrate.c -- integration routines from Numerical Recipes, */

/* Date: Tue May 09 09:37:14 PDT 2000 */

/* NOTES:
 *  o converted from float to double 
 *  o modified to handle functions of (exactly) 2 variables
 */

#define FUNC(x,y) ((*func)((x),(y)))

/* This routine computes the nth stage of refinement of an extended
 * trapezoidal rule. func is input as a pointer to the function to be
 * integrated between limits a and b, also input. When called with n=1,
 * the routine returns the crudest estimate of R b a f(x)dx. Subsequent
 * calls with n=2,3,...  (in that sequential order) will improve the
 * accuracy by adding 2 n-2 additional interior points. */
double 
trapzd(double (*func)(double,double), double a, double b, int n, double delta)
{
   double x,tnm,sum,del;
   static double s;
   int it,j;

   if (n == 1) {
      return (s=0.5*(b-a)*(FUNC(a,delta)+FUNC(b,delta)));
   } else {
      for (it=1,j=1;j<n-1;j++) it <<= 1;
      tnm=it;
      del=(b-a)/tnm; /* This is the spacing of the points to be added. */
         x=a+0.5*del;
      for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x,delta);
      s=0.5*(s+(b-a)*sum/tnm); /* This replaces s by its refined value. */
         return s;
   }
}

/* Returns the integral of the function func from a to b. The parameters
 * EPS can be set to the desired fractional accuracy and JMAX so that 2
 * to the power JMAX-1 is the maximum allowed number of steps.
 * Integration is performed by Simpson's rule. */
double
qsimp(double (*func)(double,double), double a, double b, double delta)
{
   double trapzd(double (*func)(double,double), double a, double b, int n, 
         double delta);
   int j;
   double s,st,ost,os;

   ost = os = -1.0e30;
   for (j=1;j<=JMAX;j++) {
      st=trapzd(func,a,b,j,delta);
      s=(4.0*st-ost)/3.0; /* Compare equation (4.2.4), above. */
      if (j > 5) /* Avoid spurious early convergence. */
         if (fabs(s-os) < EPS*fabs(os) ||
               (s == 0.0 && os == 0.0)) return s;
      os=s;
      ost=st;
   }
   nrerror("Too many steps in routine qsimp.");
}

/* Given arrays xa[1..n] and ya[1..n], and given a value x, this routine
 * returns a value y, and an error estimate dy. If P(x) is the polynomial
 * of degree N-1 such that P(xai) =yai ; i= 1, ..., n, then the
 * returned value y = P(x). */
void 
polint(double xa[], double ya[], int n, double x, double *y, double *dy)
{
   int i,m,ns=1;
   double den,dif,dift,ho,hp,w;
   double *c,*d;

   dif=fabs(x-xa[1]);
   c=dvector(1,n);
   d=dvector(1,n);
   for (i=1;i<=n;i++) { /* Find the index ns of the closest table entry, */
      if ( (dift=fabs(x-xa[i])) < dif) {
         ns=i;
         dif=dift;
      }
      c[i]=ya[i]; /* and initialize the tableau of c's and d's. */
      d[i]=ya[i];
   }
   *y=ya[ns--]; /* This is the initial approximation to y. */
   for (m=1;m<n;m++) { /* For each column of the tableau, */
      for (i=1;i<=n-m;i++) { /* loop over current c's and d's and update */
         ho=xa[i]-x;         /* them */
         hp=xa[i+m]-x;
         w=c[i+1]-d[i];
         if ( (den=ho-hp) == 0.0) nrerror("Error in routine polint");
         /* This error can occur only if two input xa's are (to
          * within roundoff) identical. */
         den=w/den;
         d[i]=hp*den; /* Here the c's and d's are updated. */
         c[i]=ho*den;
      }
      *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
      /* After each column in the tableau is completed, we decide which
       * correction, c or d, we want to add to our accumulating value of
       * y, i.e., which path to take through the tableau---forking up or
       * down. We do this in such a way as to take the most \straight
       * line" route through the tableau to its apex, updating ns
       * accordingly to keep track of where we are. This route keeps the
       * partial approximations centered (insofar as possible) on the
       * target x.  The last dy added is thus the error indication. */
   }
   free_dvector(d,1,n);
   free_dvector(c,1,n);
}

/* Returns the integral of the function func from a to b. Integration is
 * performed by Romberg's method of order 2K, where, e.g., K=2 is
 * Simpson's rule. */
double 
qromb(double (*func)(double,double), double a, double b, double delta)
{
   void polint(double xa[], double ya[], int n, double x, double *y, double *dy);
   double trapzd(double (*func)(double,double), double a, double b, int n, 
         double delta);
   double ss,dss;
   /* These store the successive trapezoidal approximations and their
    * relative stepsizes.  */
   double s[JMAXP],h[JMAXP+1]; 
   int j;

   h[1]=1.0;
   for (j=1;j<=JMAX;j++) {
      s[j]=trapzd(func,a,b,j,delta);
      if (j >= K) {
         polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);
         if (fabs(dss) <= EPS*fabs(ss)) return ss;
      }
      h[j+1]=0.25*h[j];
      /* This is a key step: The factor is 0.25 even though the stepsize
       * is decreased by only 0.5. This makes the extrapolation a
       * polynomial in h^2 as allowed by equation (4.2.1), not just a
       * polynomial in h. */
   }
   nrerror("Too many steps in routine qromb");
}