Section 43-04.html

Unit 8: Numerical Methods

Chapter 43: Numeric Integration

Section 43.4: adaptive quadrature

Copyright

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(student):
read(`C:/Program Files/Maple 6.01/pvac.txt`):

Warning, the protected names norm and trace have been redefined and unprotected

>

Motivation

Suppose Simpson's rule requires points to evaluate the integral

to an accuracy of . It may well be that in one portion of the interval the function varies more rapidly than in another, so that the number of function evaluations is actually being determined by only a small part of . If so, it might be more efficient to integrate over each subsection of , each time using just enough points appropriate for that subsection. Perhaps the total number of function evaluations can thereby be reduced.

By way of illustration, consider the function

> f := cosh(sqrt(1+x+2*x^2));

>

and the definite integral

> Lambda := Int(f,x=-2..3);

>

which evaluates to

> lambda := evalf(Lambda);

>

Using Simpson's rule to approximate with an error no greater than takes 51 function-evaluations, as seen from the following computation.

> for k from 24 to 25 do
evalf(simpson(f,x=-2..3,2*k)-lambda);
od;

>

With , the error is still too large, but with , so the total number of function-evaluations is 51, the error drops under the required bound.

Suppose we split the interval exactly in half, forming the two subintervals and , then integrate separately over the two subintervals, permitting an error of at most on each subinterval. On each subinterval we have the integrals

> Lambda[1] := Int(f,x=-2..1/2);
Lambda[2] := Int(f,x=1/2..3);

>

whose values are, respectively,

> lambda[1] := evalf(Lambda[1]);
lambda[2] := evalf(Lambda[2]);

>

To determine the number of function-evaluations needed by Simpson's rule on each subinterval, perform the following calculations. To approximate , obtain

> for k from 8 to 9 do
Left_error[2*k] := evalf(simpson(f,x=-2..1/2,2*k)-lambda[1]);
od;

>

so that 19 function-evaluations are needed to make the error smaller than .

To approximate , obtain

> for k from 14 to 15 do
Right_error[2*k] := evalf(simpson(f,x=1/2..3,2*k)-lambda[2]);
od;

>

so that 31 function-evaluations are needed to make the error smaller than . Hence, a total of function-evaluations are needed with the interval-splitting technique. (Since the two subintervals are contiguous, the function need be computed at the midpoint of just once, not twice.) This is a slight savings over the original tally of , but our intent was merely to provide a motivation for the interval-splitting that constitutes the adaptive quadrature strategy.

>

Error-Test Mechanism

A key ingredient of an adaptive scheme is the error-test mechanism. In Section 4.1, in the context of numeric solutions of differential equations, we derived the error-test mechanism

valid if the error in is . Thus, the error in the more accurate approximation (i.e., the one made with stepsize ), is measured by the difference between the more accurate -approximation, and the less-accurate -approximation, all divided by . This very calculation was reviewed in Section 42.2 in the context of Richardson extrapolation.

Since Simpson's composite rule is , the factor in the denominator of the error-test estimate would be .

We will show how the interval-splitting strategy is coupled with this error-test expression to produce an adaptive, or self-directed, numeric integration scheme.

>

Naive Maple Implementation of Adaptive Quadrature

We next sketch a naive implementation of an adaptive quadrature scheme based on Simpson's rule. Since its purpose is merely to illustrate the strategy, no attempt is made to minimize the number of function-evaluations. For the moment, we ignore the efficiencies that would accrue we saved all function-evaluations which get re-computed during the execution of our algorithm.

We implement our algorithm in Maple, and use Maple's built-in simpson command, thereby ignoring the efficiencies that would be had if we stored all function-evaluations which get re-computed during the execution of our algorithm.

To evaluate the definite integral

we apply Simpson's rule with stepsize and again, with stepsize . The first computation is done by invoking Simpson's rule with , the second, with . As a consequence, we have two estimates of , from which the error in the more accurate approximation can be determined by dividing the difference by 15.

If the error is small enough, we accept the more accurate -approximation, and the computation is finished.

If the error is too large, the process is repeated on the interval , with the interval set aside in a last-in, first-out stack , for processing later. If the numeric approximation for the left-hand subinterval is acceptable, its value is stored, and the process repeated on the right-hand subinterval. The accuracy demanded of any subinterval is proportional to the ratio of the length of the subinterval to the length of the original interval. If the over-all tolerance for is , then the permitted error in any subinterval is times the ratio of the length of to the length of .

To help follow the logic of this process, we define two functions. The first, performs two Simpson's rule integrations, one with and one with . The inputs to the function are , the integrand, and a list representing the interval over which Simpson's rule is to be applied. The output of the function is a list of three numbers, the more accurate of the two Simpson's rule evaluations, the estimate of the error in this value, and , the length of the interval over which these two Simpson's rule evaluations were done.

> INT := proc(f,X)
local S2,S4;
S2 := evalf(simpson(f,x=X[1]..X[2],2));
S4 := evalf(simpson(f,x=X[1]..X[2],4));
[S4,(S4-S2)/15,X[2]-X[1]];
end:

>

The second function we define performs the error-test. Its input is the list output by the first function, and its output is either true or false , depending on whether the -approximation is accurate enough or not. The accuracy demanded of any subinterval is proportional to the ratio of the length of the subinterval to the length of the original interval. If the over-all tolerance for is , then the permitted error in any subinterval is times the ratio of the length of to the length of .

> errorcheck := proc(Y)
if abs(Y[2]) < abs(evalf(Y[3]*E/H)) then true else false;fi;
end:

>

To evaluate the integral

> Lambda := Int(sqrt(3-x),x=-1..1);

>

we define the integrand as

> f := integrand(Lambda);

>

and execute the following Maple code in which , and .

> A := -1:
B := 1:
H := B-A:
E := 10^(-8):
Sigma := 0:
L := [NULL]:
P := [A,B]:
N := 1:
while N>0 do
OUT:=INT(f,P):
if errorcheck(OUT) then Sigma := Sigma+OUT[1]:
N := nops(L):
if N>0 then P := L[N];
L := [seq(L[k],k=1..N-1)];
fi;
else
C := (P[1]+P[2])/2:
L := [op(L),[C,P[2]]]:
P := [P[1],C];
N := nops([L]):
fi;
od:
print(Sigma);

>

That the value produced by this code is within the prescribed error bound, we compute the difference between and the adaptive value stored in the variable .

> evalf(Lambda-Sigma);

>

Modifying the code to accumulate and print lists of intervals, we obtain Table 43.13.

> A := -1:
B := 1:
H := B-A:
E := 10^(-8):
Sigma := 0:
Pseq := NULL:
Lseq := NULL:
Iseq := NULL:
L := [NULL]:
P := [A,B]:
Pseq := Pseq,P:
Lseq := Lseq,L:
N := 1:
while N>0 do
OUT := INT(f,P):
Iseq := Iseq,errorcheck(OUT):
if errorcheck(OUT) then Sigma := Sigma+OUT[1]:
N := nops(L):
if N>0 then P := L[N];
L := [seq(L[k],k=1..N-1)];
Pseq := Pseq,P:
Lseq := Lseq,L:
fi;
else
C := (P[1]+P[2])/2:
L := [op(L),[C,P[2]]]:
P := [P[1],C];
Pseq := Pseq,P:
Lseq := Lseq,L:
N := nops([L]):
fi;
od:
IIseq := subs(true=`pass`,false=`fail`,[Iseq]):
M := augment([seq(k,k=1..nops([Iseq]))], vector([Pseq]), vector([Lseq]), vector(IIseq)):
stackmatrix([`row`,`Present`,`Present`,`Integration Outcome`],[`index`,`Interval`,`Stack`,`on Present Interval`],[``,``,``,``],M);
print(Sigma);

>

The left-hand column indexes the "cycles" executed by the code. The second column lists the "active interval," the interval over which integration is presently taking place. The third column lists the contents of the stack, the pile of subintervals over which integration has yet to take place. The fourth column lists the outcome of the error-test applied to the integration on the "active interval." Where it says pass , the integration was accurate enough, and no further splitting of that interval need take place. Where it says fail , the integration was not accurate enough, and the "active integral" must be split in half, with the left-hand portion becoming the next "active interval" and the right-hand portion being added to the stack.

Taking the computation from the start, we begin with as the active, or present, interval. Simpson's rule is applied to this interval with , and , by invoking simpson with and , respectively. As seen from Table 43.13, this first integration is not accurate enough, so the interval is split into the subintervals and .

In cycle 2, the active interval is and the interval is placed on the stack. Again, integration over the active interval is done twice, with , respectively, in Simpson's rule. The outcome of the integrations is passed to the error-test, and again, the integration over the active interval is not accurate enough. Hence, the active interval is split again, this time into the new active interval , and the interval which is placed on the stack.

In cycle 3, the active integral is , and integration over it is performed twice, with , respectively, in Simpson's rule. The outcome of the two integrations is passed to the error-test, and again, the error is too large. Once again, the active interval is split, this time into the new active interval and the interval which is placed on the stack.

In cycle 4, the active integral is , and integration over it is performed twice, with , respectively, in Simpson's rule. The outcome of the two integrations is passed to the error-test, and finally, the error is deemed small enough. The value produced by Simpson's rule is accepted as the value of the integral over the active interval, and this number is added to the accumulator which was initialized to 0 at the very start of the computation. No new subinterval is added to the stack. In fact, the last subinterval added to the stack at the end of cycle 3 is now removed from the stack, and becomes the new active interval for cycle 5.

We leave it to the reader to continue this exposition, cycle-by-cycle. Instead, we turn our attention to some final aspects of the logic leading to Table 43.13.

1. Each time integration over a subinterval is deemed accurate enough for that value to be "kept," it is added to the accumulator .

2. The computation ends when integration over the present interval is accepted as accurate, and the stack is empty. Table 43.13 shows the stack is empty after cycles 8, 12, 14, and 16. However, the integration in cycles 9, 13, and 15 do not pass the error test, so the integration does not terminate. Only after cycle 16 does the next integration pass the error test, and that is when the integration terminates.

3. When the integration process ends, the value accumulated in is the desired approximation of .

4. To program a "split" of the active interval , compute its midpoint , and form the right-hand subinterval before forming the left-hand subinterval . If were overwritten by first, the subinterval would be unavailable.

>

Minimizing Function-Evaluations

It is highly inefficient to invoke a Simpson procedure twice for each cycle, and to invoke Simpson's rule "fresh" for each subinterval. A little thought shows that repeated calls to Simpson's rule requires re-calculating function values. If we can determine where these function values are being used, and can determine a way to store them, the number of function-evaluations can be dramatically reduced.

In the example used to generate Table 43.13, for the two initial integrations, the active interval was , so when Simpson's rule was applied with , was evaluated at the three nodes . For that same interval, when Simpson's rule was applied with , was evaluated at the five nodes . However, only two of these nodes, namely, = + , are new nodes. Hence, if we had stored the values of at the first three nodes, we would have needed to evaluate at just two more nodes to complete the first cycle. The first cycle requires just five function-evaluations.

The integrations in the first cycle were not accurate enough, so the initial interval was split into the two subintervals, and . Two new integrations were then performed in the new active interval, . The application of Simpson's rule required evaluating at the three nodes . Had we stored all function values computed in cycle 1, we could have obtained this new Simpson's rule approximation without computing new function values. The application of Simpson's rule required evaluating at the five nodes . However, only two of these five nodes, namely, are new. It takes just two new function-evaluations to complete the integrations in cycle 2.

A similar analysis shows that after cycle 1 (which requires five function-evaluations), all subsequent cycles can be completed with just two more function-evaluations each. Thus, the 17 cycles we executed should have taken just function-evaluations. Had we used Simpson's rule with , we would have expended the equivalent computational energy, and would have approximated with an error of

> evalf(simpson(f,x=-1..1,36)-Lambda);

>

50% more that the adaptive version's 0.2 x .

However, for the function

> f := 1/x+x^2/(1+x^2);

>

the integral

> Lambda := Int(f,x=1/10..5);

>

whose value is

> lambda:=evalf(Lambda);

>

would require more cycles to obtain to within an error of only . Indeed, modifying and re-executing our Maple code, we obtain

> A := 1/10:
B := 5:
H := B-A:
E := 10^(-4):
Sigma := 0:
m := 0:
L := [NULL]:
P := [A,B]:
Pseq := Pseq,P:
Lseq := Lseq,L:
N := 1:
while N>0 do
m := m+1:
OUT := INT(f,P):
if errorcheck(OUT) then Sigma := Sigma+OUT[1]:
N := nops(L):
if N>0 then P := L[N];
L := [seq(L[k],k=1..N-1)];
fi;
else
C := (P[1]+P[2])/2:
L := [op(L),[C,P[2]]]:
P := [P[1],C];
N := nops([L]):
fi;
od:
print(`value of integral = `,Sigma);
print(`number of cycles = `,m);

>

This value for the integral is supposed to be accurate to within . The actual error made is

> Sigma-lambda;

>

which is within the allowed tolerance.

With 27 cycles, we would need, at most, function-evaluations. If Simpson's rule were used with a uniform distribution of 57 nodes, we would compute with an error of

> evalf(simpson(f,x=1/10..5,56)-lambda);

>

which is considerable less accurate than the value obtained by the adaptive version of Simpson's rule. In fact, to obtain a comparable accuracy with a uniform Simpson's rule, we would need 205 function-evaluations, as verified by the following calculations.

> evalf(simpson(f,x=1/10..5,202)-lambda);
evalf(simpson(f,x=1/10..5,204)-lambda);

>

Adaptive Quadrature as a Procedure

The following Maple code captures the adaptive quadrature algorithm as the unsophisticated Maple procedure adapt . The input to the procedure is an expression representing the integrand of the definite integral being evaluated, the limits of integration, and an error tolerance. The output is a numeric approximation of the integral with less than the prescribed error, the number of cycles, and the number of function-evaluations.

> adapt := proc(f,A,B,E)
local H,Sigma,m,L,P,Pseq,Lseq,N,C,k,S2,S4,OUT;
H := B-A:
Sigma := 0:
m := 0:
L := [NULL]:
Pseq:=NULL:
Lseq:=NULL:
P := [A,B]:
Pseq := Pseq,P:
Lseq := Lseq,L:
N := 1:
while N>0 do
m := m+1:
S2 := evalf(simpson(f,x=P[1]..P[2],2));
S4 := evalf(simpson(f,x=P[1]..P[2],4));
OUT:=[S4,(S4-S2)/15,P[2]-P[1]];
if abs(OUT[2]) < abs(evalf(OUT[3]*E/H)) then
Sigma := Sigma+OUT[1]:N := nops(L):
if N>0 then P := L[N];
L := [seq(L[k],k=1..N-1)];
fi;
else
C := (P[1]+P[2])/2:
L := [op(L),[C,P[2]]]:
P := [P[1],C];
N := nops([L]):
fi;
od:
print(`value of integral = `,Sigma);
print(`number of cycles = `,m);
print(`number of function-evaluations = `,(m-1)*2+5);
end:

>

Applying the procedure adapt to the definite integral

> Lambda;

>

via the syntax

> adapt(f,1/10,5,10^(-4));

>

gives the same results obtained earlier.

>