Download Problems with Solution of Numerical Computations - Assignment 4 | MATH 451 and more Assignments Mathematics in PDF only on Docsity!
CSE/MATH 451: Numerical Computations, section 2,
spring 2005
Answers to home work 4
1 Trapezoid and Simpson’s methods
a). The trapezoid rule for uniform grid is:
T (f ; h) = h ·
[
f (x 0 ) +
f (xn) +
i=
n − 1 f (xi)
]
For our example, n = 4, and you get 0.55250538.
b). The Simpson’s rule for uniform grid is:
S(f ; h) =
h 3
[
f (x 0 ) + f (x 2 n) + 4
n∑− 1
i=
f (x 2 i+1) + 2
n∑− 1
i=
f (x 2 i)
]
For our example, n = 2, and you get 0.55067591.
c). The exact value is e^0 − e−^0.^8 = 0.55067104. The absolute error with trapezoid rule is 0.00183435, with Simpson’s rule is 0.00048716. The Simpson’s method gives better result.
d). Since f ′′(x) = e−x, then |f ′′| ≤ 1 on the interval [0, 0 .8]. Then, the error bound for trapezoid rule is |ET | ≤
h^2
and the error bound for Simpson’s rule is
|ES | ≤
h^4.
If the error should be less than 10−^4 , then for the trapezoid rule, we require:
- 8 12
h^2 ≤ 10 −^4 , → h ≤ 0. 0387298 , → n ≥ 20. 6
so we need at most n + 1 = 22 points. For the Simpson’s rule, we get
- 8 180
h^4 ≤ 10 −^4 , → h ≤ 0. 3873 , → 2 n ≥ 2. 06
and we need at most 2n + 1 = 4 points. We see that Simpson’s method need much fewer points to give the same accuracy as the trapezoid method with many more points.
2 Matlab problem: Trapezoid rule
My version of the trapezoid rule reads like this:
function v=trapezoid(f,a,b,n) % function v=trapezoid(f,a,b,n) % Compute the numerical integration with trapezoid rule % Input: a,b: interval % n: number of sub-intervals % f: the given function to be integrated % Output: v: the num itg of f.
h=(b-a)/n; v=(feval(f,a)+feval(f,b))/2; x=a; % the x-coordinate for i=1:1:n-1, x=x+h; v=v+feval(f,x); end v=v*h;
The intergrating function is defined in the file funItg.m as
function y=funItg(x) % function y=funItg(x) % the funtion to be integrated in num itg. y= exp(-x);
Here is the script to run it:
% a script to use the trapezoid.m n=[4,8,16,32,64,128]; error=zeros(size(n)); a=0; b=0.8; exact=exp(-a)-exp(-b); % the exact value for i=1:1:length(n), error(i)=abs(trapezoid(’funItg’,a,b,n(i))-exact); end loglog(n,error) grid, title(’Errors with trapezoid rule’) print -depsc trapeError.eps
And I get the following plot:
v=v*h/3;
Here is the script to run it:
% a script to use the trapezoid.m n=[4,8,16,32,64,128]; error=zeros(size(n)); a=0; b=0.8; exact=exp(-a)-exp(-b); % the exact value for i=1:1:length(n), error(i)=abs(simpson(’funItg’,a,b,n(i))-exact); end loglog(n,error) grid, title(’Errors with trapezoid rule’) print -depsc simpsError.eps
And I get the following plot:
100 101 102 103
10 −
10 −
10 −
10 −
10 −
10 −
10 −
10 −^
Errors with Simpson rule
We see that when n doubles, which means h halfs, the error is reduced by a factor of 16. The error reduces much faster here than the one of trapezoid rule.
4 Romberg algorithm
a) My version of the file romberg.m looks like the following:
function R = romberg(f,a,b,n) h = b-a; R = zeros(n,n);
R(1,1) = (feval(f,a) + feval(f,b))h/2; for i=2:n h = 0.5h; sum = 0; for k = 1:2:2^(i-1)- sum = sum + feval(f,a+kh); end R(i,1) = 0.5R(i-1,1) + sum*h; for j = 2:i R(i,j) = R(i,j-1) + (R(i,j-1)-R(i-1,j-1))/(4^(j-1)-1); end end
b) Use the files:
function y = f1(x) y = sin(x);
f1.m
function y = f2(x) y = sqrt(x);
f2.m
And for f1.m, I get with n = 4:
0.00000000000000 0 0 0 1.57079632679490 2.09439510239320 0 0 1.89611889793704 2.00455975498442 1.99857073182384 0 1.97423160194555 2.00026916994839 1.99998313094599 2.
In order to get the error along the diagonal, I run the following script:
format short e; R=romberg(’f1’,0,pi,7); error=abs(diag(R)-2)
and I get:
error = 2.0000e+ 9.4395e- 1.4293e- 5.5500e- 5.4127e- 1.3212e- 8.8818e-
And for f2.m, I repeat the same thing, and get the error:
error = 2.1175e+ 2.8863e+ 2.9959e+ 3.0285e+ 3.0396e+ 3.0434e+ 3.0448e+
Multiply the equation (B) by 9 and subtract (a) from it, we get
8 J = 9T (
h 3 ) − T (h) −
a 4 h^4 −
a 6 h^6 + · · ·
This gives us a new integration’s formulae:
U (h) = T (h/3) + (T (h/3) − T (h))/ 8 ,
and the new Euler-Maclaurin’s formulae
J = U (h) + c 4 h^4 + c 6 h^6 + · · · (C)
where c 4 , c 6 are new constants. This form is valid also for h/3,
J = U (
h 3
) + c 4
h 3
h 3
+ · · · (D)
Again, multiply (d) by 9 and subtract it from (C), then devide both sides by 80, we get
J = U (
h 3
[
U (
h 3
) − U (h)
]
If we could continue, we would get a extrapolation’s formulae where
U (i, 0) = T (
h 3 i^
), U (i, j) = U (i, j − 1) +
9 j^ − 1
[U (i, j − 1) − U (i − 1 , j − 1)]
which would set up a table
U(0,0) U(1,0) U(1,1) U(2,0) U(2,1) U(2,2)
or we put in the numbers:
0.9432914291 0. 0.9457731886 0.9460834085 0.
c) We use sin(t) ≈ t −
t^3 +
t^5
As an approximation to J we could use:
J =
0
sin(t) t
dt ≈
0
t^2 +
t^4 ) = 1 −
d) For the error estimate, we start with the Taylor expansion we used. We actually have
sin(t) = t −
t^3 +
t^5 −
cos(ξ)t^6.
The error is dominated by the term − (^) 7!^1 cos(ξ)t^6 , where ξ is a point between [0, 1]. Therefore the error for the integral will be
error =
0
cos(ξ)t^6 dt = −
cos(ξ)
0
t^6 dt = −
cos(ζ)
Using tha fact that cos is bounded by 1, we get an error bound:
|error| ≤
= 2. 83 · 10 −^5
(The actual error you get is 2. 80 · 10 −^5 , so this bound is pretty good.)
