c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole Publ., 1990 c c Section 7.5 c c adaptive quadrature by simpson rule c c c file: adapt.f c parameter (n = 5) dimension v(6,n) logical iprt data a,b,epsi/0.0,1.0,5.0e-5/ data iprt/.true./ f(x) = sqrt(x) c print * print *,' Adaptive Quadrature' print *,' Section 7.5, Kincaid-Cheney' print * print 4,'k','v(1,n)','v(2,n)','v(3,n)','v(4,n)','v(5,n)','v(6,n)' c c initialize everything, c particularly the first column vector in the stack. c delta = b-a sigma = 0.0 h = delta/2.0 c = 0.5*(a + b) k = 1 abar = f(a) bbar = f(b) cbar = f(c) S = (abar + 4.0*cbar + bbar)*h/3.0 v(1,1) = a v(2,1) = h v(3,1) = abar v(4,1) = cbar v(5,1) = bbar v(6,1) = S if (iprt) print 5,k,(v(i,1),i=1,6) c 2 continue if ( (k .lt. 1) .or. (k .gt. n) ) go to 3 c c take the last column off the stack and process it. c h = 0.5*v(2,k) y = v(1,k) + h ybar = f(y) Star = ( v(3,k) + 4.0*ybar + v(4,k) )*h/3.0 z = v(1,k) + 3.0*h zbar = f(z) SStar = ( v(4,k) + 4.0*zbar + v(5,k) )*h/3.0 c tmp = Star + SStar - v(6,k) if ( abs(tmp) .le. 30.*epsi*h/delta ) then c if the tolerance is being met, add partial integral c and take a new vector from the bottom of the stack. c sigma = sigma + Star + SStar + tmp/15.0 if (iprt) print *,' k=',k,' sigma=',sigma k = k-1 c if (k .le. 0) then print *,' value of integral is ',sigma stop endif else if (k .ge. n) then c if the stack has reached it maximum number of columns (n), c then stop and report a failure. c print *,' method fails' stop endif c c if the tolerance is not being met, subdivide the interval c and create two new vectors in the stack, c one of which overwrites the vector just processed. if (iprt) print *,' split' vbar = v(5,k) v(2,k) = h v(5,k) = v(4,k) v(4,k) = ybar v(6,k) = Star if (iprt) print 5,k,(v(i,k),i=1,6) c k = k+1 v(1,k) = v(1,k-1) + 2.0*h v(2,k) = h v(3,k) = v(5,k-1) v(4,k) = zbar v(5,k) = vbar v(6,k) = SStar if (iprt) print 5,k,(v(i,k),i=1,6) endif c c having created two new vectors in the stack, c pick off the last one and start the proces anew. c go to 2 c 3 continue c 4 format (a2,a12,a14,a14,a14,a14,a14) 5 format (i2,1x,6(e13.6,1x)) stop end