Download Richardson Extrapolations, Lecture Notes - Advanced Calculus and more Study notes Calculus in PDF only on Docsity!
Math 128A Lab 7
Richardson Extrapolation
Adrian Down
April 05, 2006
Part 1
Discussion
In this lab, I use Richardson extrapolation to find approximate values of the following integrals,
∫ (^) π
0
dx
sin(x) x
∫ (^) π
0
dx
sin(x) √ x
I form I Richardson table recursively using the algorithm developed in lecture,
T n(hi) =
22 n^ − 1
[
22 nT n−^1 (hi+1) − T n−^1 (hi)
]
Using (3), the T 0 element of the iteration is the trapezoidal approximation to the integral. In lecture, we also found the order of the dominant term in the error associated with the Richardson extrapolation,
T n(hi) −
∫ (^) b
a
dx f = ǫn(hi) = O(h 2(n+1) i )^ (4)
The convergence behavior claimed in (4) can be verified by taking a log- arithmic ratio of error terms resulting from successively smaller hi values. In
(4), the coefficient of the hmi error term, call it α, is dependent only on the function f and order of the Richardson extrapolation, n. α is independent of i, which determines the size of the mesh of the Richardson extrapolation. Taking only leading order terms in (3) to form the ratio,
log 2
ǫn(hi) ǫn(hi+1)
≈ log 2
αh2( in+1) α
( (^) hi 2
)2(n+1)
= log 2 2 2(n+1)^ = 2(n + 1) (5)
Plotting the logarithmic ratio of error terms indicated in (5) for varying degrees n of the Richardson extrapolation, the preceding argument predicts horizontal lines corresponding to the order of the Richardson extrapolation. However, in practice, such lines are not exactly horizontal. The accuracy of these ratios is limited by the numerical precision of MatLab, especially at higher values of n. Also, (5) considered only leading terms in the error. Since higher order terms are neglected in both the numerator and denominator of (5), the resulting deviation from an exact ratio can be appreciable.
Results
i T 0 (hi) T 1 (hi) T 2 (hi) T 3 (hi) 0 0 1. 6710855164207 1. 7400236041721 1. 7726636008419 1 1. 2533141373155 1. 7357149736877 1. 7721536008939 1. 7836434891535 2 1. 6151147645946 1. 7698761866935 1. 7834639596494 1. 7875339445617 3 1. 7311858311688 1. 7826147238397 1. 7874703510475 1. 7889101548803 4 1. 769757500672 1. 787166874347 1. 7888876579454 1. 789396783323 5 1. 7828145309282 1. 7887801089705 1. 789388828239 1. 7895688381823 6 1. 78728871446 1. 7893507832847 1. 7895660255269 1. 7896296692672 7 1. 7888352660785 1. 7895525728868 1. 7896286748338 1. 7896511763483 8 1. 7893732461847 1. 7896239184621 1. 7896508247622 1. 7896587802537 9 1. 7895612503928 1. 7896491431184 1. 7896586559492 1. 7896614686406 10 1. 789627169937 1. 7896580613973 1. 7896614246923 11 1. 7896503385322 1. 7896612144864 12 1. 7896584954978
Table 3: Richardson approximation to the integral of sin(^ √xx ) on the interval
[0, π]. The size of the mesh used to compute each approximation is hi = 2 πi
i ǫ^0 ǫ^1 ǫ^2 ǫ^3 0 1. 789662385 0. 1185768688 0. 04963878109 0. 01699878442 1 0. 5363482479 0. 05394741158 0. 01750878437 0. 00601889611 2 0. 1745476207 0. 01978619857 0. 006198425614 0. 002128440701 3 0. 05847655409 0. 007047661423 0. 002192034216 0. 0007522303828 4 0. 01990488459 0. 002495510916 0. 0007747273177 0. 0002656019401 5 0. 006847854335 0. 0008822762926 0. 0002735570241 9 .354708087 e – 6 0. 002373670803 0. 0003116019784 9 .635973623 e –5 3 .271599587 e – 7 0. 0008271191846 0. 0001098123764 3 .371042932 e –5 1 .120891478 e – 8 0. 0002891390784 3 .846680101 e –5 1 .156050094 e –5 3 .60500938 e – 9 0. 0001011348704 1 .32421447 e –5 3 .729313935 e –6 9 .16622495 e – 10 3 .521532611 e –5 4 .323865858 e –6 9 .605707986 e – 11 1 .204673092 e –5 1 .17077674 e – 12 3 .889765285 e –
Table 4: Error associated with the Richardson approximation to the integral of sin(^ √xx ) on the interval [0, π]. The size of the mesh used to compute each approximation is hi = 2 πi
(^00) 0.5 1 1.5 2 2.5 3
1
2 sin(x) / x sin(x) / x1/ 2sin(x^2 )
Figure 1: Three functions for which approximate integrals are computed.
(^00 2 4 6 8 10 )
1
2
3
4
5
6
7
8
9
10 n = 0 n = 1 n = 2 n = 3
Figure 3: Order of convergence of error terms in the Richardson extrapolation to the integral of sin(^ √xx )on [0, π]. Values on the x axis are the values of i from
log 2 ǫ
n(hi) ǫn(hi+1).^ Values at^ i^ = 0 correspond to theoretical values computed in (5).
Code
% Created by Adrian Down, University of California, Berkeley. % Created for Prof. John Neu’s numerical analysis class, Math 128A, Spring 2006 % ----------------------------------------------------------------------------- % function epsilon = Lab7(f,a,b) % Compute a Richardson approximation to the integral of the symbolic % function ’f’ on the interval [a,b]. Computes a Richardson table that is % ’n’ levels deep. Convergence is examined by taking the logarithmic ratio % of error terms of the form e(h)/e(h/2). Only convergence up to % ’maxorder’ is examined. n = 12; % depth of Richardson table maxorder = 4; % maximum order to examine convergence h = b-a; % initial step size is length of interval hp = pi/2^(n+2); % step size for calculation of "exact" value of integral TrueInt = IntTrap(f,a,b,hp); % true integral calculated with trapezoidal rule AproxInt = IntRich(f,a,b,h,n); % construct Richardson table epsilon = TrueInt - AproxInt; % construct table of errors logerror = zeros(n,n+1); % initialize table of error ratios for i = 1:n logerror(i,:) = log2(abs(epsilon(i,:)./epsilon(i+1,:))); end colors = [’b’ ’g’ ’r’ ’c’ ’m’ ’y’ ’b’]; % for display of convergence behavior close all; figure; hold on; leg = ’’; for i = 1:maxorder % loop over rows xplot = [1:1:n+1-i]; errplot = zeros(1,n+1-i); for j = 1:n+1-i % only include error ratios that are in bounds of Richardson t errplot(j) = logerror(j,i); end plot([0 xplot], [2*i errplot], colors(i)) leg = [leg; ’n = ’ int2str(i-1)]; end
function A = IntRich(f,a,b,h,n,varargin) % Compute approximation to the integral of f on the interval [a,b] using % Richardson extrapolation. % input: function f, degree of Richardson extrapolation n, endpoints of % interval of integration a and b, bin width h % if the variable in f is not ’x’, input variable as optional last input % output: table of Richardson extrapolating polynomial. The kth degree % Richardson extrapolation with bin width h is at A(1,k). % syms x; if nargin == 6 % if f is not a function of x, override default variable x = varargin{1}; end A = zeros(n+1); for i = 1:n+1 % compute first column of Richardson table... A(i,1) = IntTrap(f,a,b,h/2^(i-1),x); % which is trapezoidal approximation. h end % form the rest of the Richardson table using linear combinations of previous entr for i = 2:n+1 % loop over columns for j = 1:(n+2-i) % loop over rows A(j,i) = 1/(2^(2(i-1))-1)(2^(2(i-1))A(j+1,i-1) - A(j,i-1)); % form Ric end end
% Created by Adrian Down, University of California, Berkeley. % Created for Prof. John Neu’s numerical analysis class, Math 128A, Spring 2006 % ----------------------------------------------------------------------------- % function val = CheckSub(f,x,a) % Substitute the value ’a’ for the variable ’x’ into the symbolic function ’f’ % If f is not defined at x = a, the value of the right-sided limit is taken val = subs(f,x,a); if (~isfinite(val)) % if f is not defined at x = a... val = limit(f,x,a,’right’); end
Part 2
Discussion
As shown in table 4, the convergence of the Richardson approximation to the integral of sin( xx )on [0, π] is poor. At the third order approximation, the size of the error decreases as the mesh size of the approximation is made smaller, although not with the predicted order of convergence (see figure 3). In contrast, decreasing the size of the mesh used to approximate the integral of sin( xx ) on the interval [0, π] using a third order Richardson extrapolation makes no difference because the approximation is already accurate within the numerical precision of MatLab. The poor convergence of the approximation to
∫ (^) π 0 dx^
sin( √x) x occurs because the integrand is not analytic at x = 0. This singularity reduces the accuracy of the leftmost trapezoid used in computing the trapezoidal approximation to the integral, and the loss of accuracy propagates through the Richardson table, reducing the accuracy of the higher order approximations. To make the integrand analytic, it is possible to make a change of vari- ables,
u =
x du =
dx 2
x
With this change of variable, the desired integral can be recast,
∫ (^) π
0
dx
sin(x) √ x
∫ √π
0
du 2 sin(u^2 )du (6)
Since sin(u) is analytic at u = 0, (6) is analytic on the range of integration [0,
π], thus the convergence of the Richardson approximation should be improved. This is in fact the case, as can be seen in table 6 and figure 4.
Results
(^00 2 4 6 8 10 )
1
2
3
4
5
6
7
8
9
10 n = 0 n = 1 n = 2 n = 3
Figure 4: Order of convergence of error terms in the Richardson extrapolation to the integral of 2 sin(x^2 ) on [0,
π]. Values on the x axis are the values of i
from log 2 ǫ
n(hi) ǫn(hi+1). Values at^ i^ = 0 correspond to theoretical values computed in (5).
