Create an M-file function, as I did with the trapezoid rule listed below, to examine the improvement in accuracy offered by the Richardson extrapolation applied to the one-third Simpson technique. (Use OneThirdSimp.m.)
Test it with a couple of functions and make a conjecture about the unknown exponent k_2 in the error term ASSUMING the error of this Simpson rule is of the form
E(h)=a*h^4+b*h^(k_2)+o(h^(k_2)).
Note: (Not part of the problem, this is just a comment.)
OneThirdSimp.m is the vanilla 1/3 Simpson function requiring an
even number of data intervals and integrating a quadratic on consecutive
pairs of intervals. simp13.m, on the other hand, uses a sliding window
and averages the two possible quadratic integrals for an interval---except,
that is, for the first and last intervals that only have one possible
quadratic fit. The Richardson technique does not "work" very well (or at
least as well) for simp13.m compared to OneThirdSimp.m.
But it MUST work if the error is of the form
c1*h^(k1)+c2*h^(k2)+o(h^(k2)).
It is a puzzle. Can you figure out what is going on here?
%*************************************************
function [E R]=richtrap(fun,a,b)
% This script demonstrates the second order error of the trapezoid rule
% and the error of the Richardson extrapolation using the
% trapezoid rule. This uses roughly 50 percent more flops than the finer
% of the two trapezoid rule calculations.
% Be careful! This technique is extremely sensitive to details of spacing!
% linspace(a,b,n+1) produces n+1 points corresponding to n
% intervals each of length (b-a)/n.
% linspace(a,b,2*n+1) produces intervals of exactly half this width.
% The same CANNOT be said of of the pair linspace(a,b,n) and linspace(a,b,2*n).
exact=integral(fun,a,b);
if exact==0,error('The integral over the selected interval is 0. Choose another interval.'),end
for k=1:8
t=linspace(a,b,6^k+1);
s=linspace(a,b,2*6^k+1);
yt=fun(t);
ys=fun(s);
A2(k)=trapuneq(t,yt);
A1(k)=trapuneq(s,ys);
R(k)=4/3*A1(k)-1/3*A2(k);
E(1,k)=log(abs((A2(k)-exact)/exact))/log(6);
E(2,k)=log(abs((A1(k)-exact)/exact))/log(6);
E(3,k)=log(abs((R(k)-exact)/exact))/log(6);
end