MATH/APPM 4650
Class notes, March 10

Recursion and adaptive integration

Suppose we are trying to compute ∫_a^b f(x) dx with a composite Simpson's rule.

Suppose the partition looks like {x₀, x₁, ..., x_2n}.
(We don't assume these are evenly-spaced, although we want x₁ to be the midpoint of x₀ and x₂, etc.)

Then the error in Simpson's rule is the sum of the errors in each subinterval, which is E_n = -(x₂ - x₀)⁵/2880 f⁽⁴⁾(ξ₁) - ... - (x_2n - x_2n-2)⁵/2880 f⁽⁴⁾(ξ_n)

To make this small, we want more points where the fourth derivative is large.
That way, each subinterval will make a small contribution to the total error.

This is more efficient than shrinking all intervals at once.
This is the basic idea behind adaptive integration, first devised in the 1960s.

Idea: If the error term in a particular subinterval is large, then divide it in two and recompute Simpson's rule on each of the two subintervals.

But how do we know if the error term is large? We don't know f⁽⁴⁾, and even if we did, we don't know ξ_i.

Instead we check each Simpson's rule computation against the same computation on the two halves of the interval.
It's like Richardson's method, except instead of trying to get rid of the error, we try to estimate what the error is.

If S(a, b) is the result of Simpson's rule on the interval [a, b], and c is the midpoint of [a, b], then the true value of the integral can be approximated by either S(a,b) or S(a,c)+S(c,b).

The error in S(a,b) is ∫_a^b f(x) dx - S(a,b) = -(b-a)⁵/2880 f⁽⁴⁾(ξ₁).
The error in S(a,c) + S(c,b) is ∫_a^b f(x) dx - S(a,c) - S(c,b) = -(c-a)⁵/2880 f⁽⁴⁾(ξ₂) - (b-c)⁵/2880 f⁽⁴⁾(ξ₃)
= -1/32 (b-a)⁵/2880 [f⁽⁴⁾(ξ₂) + f⁽⁴⁾(ξ₃)]
= -1/16 (b-a)⁵/2880 f⁽⁴⁾(ξ₄).

Now if we're lucky, then f⁽⁴⁾(ξ₄) ∼ f⁽⁴⁾(ξ₁).
If so, we can estimate that error term δ by eliminating the exact integral: ∫_a^b f(x) dx - S(a,b) = -δ
∫_a^b f(x) dx - S(a,c) - S(c,b) = -δ/16.
We get 15/16 δ = [S(a,b) - S(a,c) - S(c,b)]
So if we want the error δ/16 to be less than some tolerance ε, we can just make sure that [S(a,b) - S(a,c) - S(c,b)] ≤ 15ε
If this is true, then we know that S(a,c) + S(c,b is correct to within ε of the true value.

If it's not true, we should divide [a,c] and [c,b] in half and do this again. However, we can't use the same tolerance. We wanted the total error to be less than ε, so the error in each half should be at most ε/2.


The algorithm is rather complicated to program directly (see for example the book's pseudocode). However, recursion makes it much easier.

---

The idea of recursion is to define a function that calls itself.

A common example is to compute the Fibonacci numbers, defined so that F[n] = F[n-1] + F[n-2] if n>1, F[0]=1, F[1]=1.

In any of the three languages, one can write a recursive program to do this.

Maple:
Fibonacci:=proc(n)
   if n>1 then
      Fibonacci(n-1)+Fibonacci(n-2):
   else
      1:
   end if:
end proc:

Mathematica:
Fibonacci[n_]:=Module[{},
   If[n>1,
      Fibonacci[n-1]+Fibonacci[n-2],
      1
     ]
  ]

MATLAB:
function F = Fibonacci(n)
if n > 1
   F = Fibonacci(n-1)+Fibonacci(n-2);
else
   1
end

Caveats:

• Recursion is usually slow. If you can avoid using recursion, do so.
  For example, the following code is MUCH faster.
  Fibo_old1 = 1
  Fibo_old2 = 1
  for i from 1 to n do
     Fibo = Fibo_old1 + Fibo_old2
     Fibo_old1 = Fibo_old2
     Fibo_old2 = Fibo
  end do
  return Fibo
  Try running both Fibonacci codes for n=30.
  The difference is enormous.
• You need a stopping criterion.
  Recursive programs are the easiest places to get infinite loops.
  The software will generally only let you do a small number of steps, since the size of a recursive program can grow exponentially.
  

That said, recursive programs are often easy to write.
Since many mathematical procedures are defined in a recursive way already, it is natural to use that structure in the program itself.
Recursion is especially useful if:

• You don't know the exact number of steps your program will require, or how many previous values the future values depend on.
• The algorithm is stated in such a way as to suggest, "If this condition is true, then stop; if not, do something and then repeat a new version of the calculation."

---

The recursive method for adaptive integration (pseudocode):

Simpsonizer = function(f,a,b)
   return (b-a)/6*(f(a)+f(b)+4*f(a+(b-a)/2)
end function


Recursive_Simpson = function(f,a,b,tol,initial_sum)
   c=a+(b-a)/2
   left_sum = Simpsonizer(f,a,c):
   right_sum = Simpsonizer(f,c,b):
   if abs(left_sum+right_sum-initial_sum)≤15*tol then
      return left_sum+right_sum:
   else
      return Recursive_Simpson(f,a,c,tol/2,left_sum)+Recursive_Simpson(f,c,b,tol/2,right_sum)):
   end if:
end function


Adaptive_Simpson = function(f,a,b,tol)
   return Recursive_Simpson(f,a,b,tol,Simpsonizer(f,a,b))
end function


Variations and caveats:

• There is no absolute stopping criterion in this program.
  One hopes the error shrinks in the right way, but there's no guarantee it will.
  
  An alternative is to count how many times one has gone through the recursive loop and stop if that's too large.
  For example, one could define a global variable "Counter" in Adaptive_Simpson, set to zero initially, and then increment Counter at the beginning of Recursive_Simpson.
  One would then check both the error condition and that Counter is less than, say, 100, before recursing deeper.
  
  
• The book suggests using the condition error ≤ 10*tol rather than 15*tol as here.
  This is a more conservative assumption (making it more likely to subdivide intervals).
  This is not necessary, but caution can be useful.
  
  
• Finally, some programmers return
  left_sum + right_sum + (left_sum+right_sum-initial_sum)/15
  instead of just
  left_sum + right_sum.
  One would expect this to be more accurate, since this program is actually adding the predicted error to the estimate.
  
  Typically the error one gets is much smaller than the given tolerance.