Script 1 For Week 9
Here is the Richardson technique applied to the trapezoid rule.
This is NOT Romberg integration. That is a clever technique, which is based on the Richardson method, to re-use earlier parts of a calculation to more efficiently reach an approximate error target. Still, the Richardson technique alone is pretty cool.
%**********************************************************
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 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
%**********************************************************
%**********************************************************
function [E R]=richtrapmistake(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 linspace(a,b,n) and linspace(a,b,2*n).
% This script illustrates how this---seemingly---tiny difference
% completely messes up the convergence rate improvement.
% As a side note, since the edges of t are not edges of s there is no additional
% improvement that can be had by re-using function calls calculated in the
% coarser of the two trapezoid calculations.
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);
s=linspace(a,b,2*6^k);
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
%**********************************************************
E=richtrap(@(t) sin(t),0,8.1)
E=richtrapmistake(@(t) sin(t),0,8.1)
E=richtrap(@(t) exp(t),0,8.1)
E=richtrapmistake(@(t) exp(t),0,8.1)
% The Euler-Maclaurin formula for analytic functions (that is, functions that can be represented by a power series on the relevant interval) explains why the trapezoid rule is very very very good when the function is periodic and you are integrating over a period. This formula can be used to represent the error term as a higher-than-simple-exponent-order error plus an infinite sum of integrals, but each of these integrals is multiplied by a factor of the form f^(n)(b)-f^(n)(a) where ^(n) indicates the nth derivative. So for analytic periodic functions each of these terms is 0.
f=@(t) sin(3*cos(t+2)+7)
t=linspace(-2*pi,4*pi);
y=f(t);
plot(t,y)
for k=6:11
t=linspace(0,2*pi,k);
y=f(t);
val(k-5)=trapuneq(t,y);
end
format long
val
exact=integral(f,0,2*pi)
exactrelerr=(val-exact)/exact
format short
% note: stepsizes used above vary from 2*pi/5 to 2*pi/10 which is about 1.26 to the smallest which is still about 0.63. Reducing this very large stepsize by half reduced the error by almost a factor of 10^4!!!
%*******************************************
function [z] = shift2ones(t,r,s)
% translate and scale domain [r,s] to [-1,1] for use with
% the Gauss-Legendre formula
z= t*(s-r)/2 + (s+r)/2 ;
end
%************************************************
% Using the Gauss-Legendre formulas:
% Given a function
r=1
s=7
L= @(t) 3*t.^2-4*t-7
% The modified function can be used to calculate the integral of L on [r,s] using the Gauss-Legendre formulas.
Lshift=@(t) L( shift2ones(t,r,s) )
% Don't forget to include the change-of-variable factor (s-r)/2 after (or before) you integrate!
L(r)
Lshift(-1)
L(s)
Lshift(1)
% Weighting factors and function arguments used in Gauss-Legendre formulas (so you don't have to type them in yourself.) To apply this technique you just need to know the function values Lshift(pn) at the points pn in [-1,1].
w1=2
p1=0
w2=[1 1]
p2=[-1/sqrt(3) 1/sqrt(3)]
w3=[5/9 8/9 5/9]
p3=[-sqrt(0.6) 0 sqrt(0.6)]
w4=[(18-sqrt(30))/36 (18+sqrt(30))/36 (18+sqrt(30))/36 (18-sqrt(30))/36]
p4=[-sqrt(525+70*sqrt(30))/35 ...
-sqrt(525-70*sqrt(30))/35 ...
sqrt(525-70*sqrt(30))/35 ...
sqrt(525+70*sqrt(30))/35]
w5=[(322-13*sqrt(70))/900 ...
(322+13*sqrt(70))/900 ...
128/225 ...
(322+13*sqrt(70))/900 ...
(322-13*sqrt(70))/900]
p5=[-sqrt(245+14*sqrt(70))/21 ...
-sqrt(245-14*sqrt(70))/21 ...
0 ...
sqrt(245-14*sqrt(70))/21 ...
sqrt(245+14*sqrt(70))/21]
y4=(7-1)/2*Lshift(p4);
INTGL=y4*w4'
exact=integral(@(t) L(t),1,7)