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)