Ma 449: Numerical Applied Mathematics
Model Solutions to Homework Assignment 2
Prof. Wickerhauser
Exercise 1(a*) in Section 2.1, p.50.
Yes, there is a unique fixed point. Function g maps [0,1] into [0,1] since
x in [0,1] implies g(x)=1-x^2/4 is in [0,1]. Also, g is a contraction
since g'(x) = -x/2, so |g'(x)| =< 1/2 < 1 on [0,1]. Therefore Theorem 2.3
equation (6) applies and g has an attractive fixed point.
Exercise 1(b) in Section 2.1, p.50.
Yes, there is a unique fixed point. Function g maps [0,1] into [0,1] since
x in [0,1] implies g(x)= 2^{-x} = (1/2)^x is in [0,1]. Also, g is a
contraction since g'(x)= log(1/2)*(1/2)^x so |g'(x)|<0.7<1, since 1/2 =<
(1/2)^x =< 1 for x in [0,1] and log(1/2)~0.6931. Therefore Theorem 2.3
equation (6) applies and g has an attractive fixed point.
Exercise 3(a) in Section 2.1, p.50.
See the graph at ./2-1-3a.jpg to be convinced that Theorem 2.3 applies
and that fixed-point iteration will converge.
Exercise 8(a,b,c) in Section 2.1, p.50.
(a) First observe that g(x) = x - 0.0001*x^2 =< x for all x, so that
1 = p_0 >= p_1 >= p_2 >= p_3 >= p_4 >= ...
To show strict equality it suffices to show that no iterate p_n is 0. But
if g(x)=p_n=0 for some x, then x=0.0001*x^2 and thus x = 1000 > p_k for
all k. Hence x is not equal to any previous iterate.
(b) Prove this by induction:
- first note that p_0=1>0
- suppose p_n > 0 for some n. Then
p_{n+1} = g(p_n) = p_n - 0.0001*(p_n)^2
= p_n*( 1 - 0.0001*p_n)
and both factors on the right are positive: the first by the
inductive hypothesis p_n > 0, the second by part (a) p_n =< 1.
(c) Let c be the limit: p_n --> c as n --> infinity. Then c >= 0.
Exclude c>0 as follows: if c>0, take n large enough so that c =< p_n < c+d
for some small positive d < 0.0001*c^2 . But then
p_{n+1} = p_n - 0.0001*(p_n)^2 < c+d - 0.0001*c^2 < c
since d - 0.0001*c^2 < 0. Since p_k < p_{n+1} for all k> n+1, we obtain
the contradiction that {p_k} cannot converge to c.
Algorithm 1 of Section 2.1, p.51 for the Ex.3(a) function.
Write this Matlab program into the file "fixpt.m":
function [k,p,err,P] = fixpt(g,p0,tol,max1)
P(1)=p0;
for k=2:max1
P(k)=g(P(k-1));
err=abs(P(k)-P(k-1));
relerr=err/abs(P(k)+eps);
p=P(k);
if (err> fixpt(@g,0.000,1e-12,10) % Note the @ sign before g
P =
0
2
2
ans =
3
However, 2 is a repelling fixed point:
>> fixpt(@g,0.001,1e-12,10)
maximum number of iterations exceeded
P =
0.0010
2.0000
1.9999
1.9974
1.9070
-0.8590
1.9581
0.5934
0.7424
-0.1044
ans =
10
Exercise 3(c*) of Section 2.2, p.61.
Many intervals work, but [1,5] has simple endpoints and satisfies
f(1) = -4 < 0, f(5) = ln(5) > 0. (In fact, ln(5)~2.) The textbook's
solution [3,4] is a smaller, thus better interval, but it requires
more calculation.
Exercise 3(d) of Section 2.2, p.61.
f(x)=x^2-10x+23 is negative at its vertex x= 5, where f(5)=-2, and
positive wherever |x| is big, such as x= 10, where f(10)=23. Thus
f(10)*f(5)<0.
Exercise 2 of Section 2.3, p.69.
See the graph at ./2-3-2.jpg to see that a single root must lie in
the interval [0, 1], and quite close to 0.8.
Exercise 4 of Section 2.4, p.85.
f(x) = x^3-3*x-2.
(a) Newton-Raphson formula: x_{n+1} = x_n - f(x_n)/f'(x_n).
But f'(x) = 3*x^2 - 3, so we may use
x = x - (x^3-3*x-2)/(3*x^2 - 3)
Note that the next approximate value replaces x after each iteration.
(b) Matlab commands:
format long
x=2.1 % p_0
x = x - (x^3-3*x-2)/(3*x^2 - 3)
x = x - (x^3-3*x-2)/(3*x^2 - 3)
x = x - (x^3-3*x-2)/(3*x^2 - 3)
Output:
x = 2.100000000000000
x = 2.006060606060606
x = 2.000024339783376
x = 2.000000000394941
(c) From these first few iterations, it appears that the sequence is
converging quadratically.
Exercise 11 of Section 2.4, pp.85--86.
Let f(x) = x^3-A. Then f'(x) = 3*x^2, and the Newton-Raphson iteration
formula after simplification is p(n+1) = (2/3)*p(n) + A/(3*p(n)^2).
Algorithm 4(a) of Section 2.4, p.88.
Like the solution to Exercise 4 in this section, just write the
Newton-Raphson formula for this particular function from Exercise 11 and
use Matlab's command history to iterate:
format long
A=7; x = 7
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
x = (2/3)*x + A/(3*x^2)
OUTPUT:
x = 7
x = 4.714285714285714
x = 3.247846429664611
x = 2.386431304900376
x = 2.000666416795918
x = 1.916722395612087
x = 1.912938676720494
x = 1.912931182801747
x = 1.912931182772389
x = 1.912931182772389
x = 1.912931182772389
The last three approximations are equal in their first 16 significant
digits, so we conclude that the sequence of approximations has converged.
Using Matlab we may check the result:
>> (1.912931182772389)^3
ans =
6.999999999999998
Rounding to 15 significant digits gives 7.
Exercise 10 of Section 2.5, p.99.
Note that S(n) plays the role of p(n) in Aitken's Theorem 2.8, p.90,
and thus Delta p(n) = p(n+1)-p(n)=S(n+1)-S(n)= A(n+1). Similarly,
Delta^2p(n) = Delta p(n+1) - Delta p(n) = A(n+2)-A(n+1). Using these
in Theorem 2.8 gives the desired formula.
Exercise 12 of Section 2.5, p.99.
Matlab program:
format long
N = 20; % This controls the number of iterations
% Form the original sequence
for k=1:(N+2) % 2 extra so we can use Aitken
a(k) = 1.0/(4^k + 4^(-k));
end
for j=1:N
s(j)=0;
for k=1:j
s(j) = s(j) + a(k); % form the usual partial sum
end
end
for k=1:N
t(k)= s(k) + a(k+1)*a(k+1)/(a(k+1)-a(k+2)); % Aitken term
end
% Publish the results:
[(1:N)' s' t']
OUTPUT:
1.000000000000000 0.235294117647059 0.318404625258633
2.000000000000000 0.297550926985580 0.318380764463586
3.000000000000000 0.313172113219410 0.318380391916081
4.000000000000000 0.317078303615675 0.318380386095296
5.000000000000000 0.318054865184353 0.318380386004347
6.000000000000000 0.318299005794801 0.318380386002926
7.000000000000000 0.318360040950824 0.318380386002904
8.000000000000000 0.318375299739883 0.318380386002903
9.000000000000000 0.318379114437148 0.318380386002903
10.000000000000000 0.318380068111465 0.318380386002903
11.000000000000000 0.318380306530044 0.318380386002903
12.000000000000000 0.318380366134688 0.318380386002903
13.000000000000000 0.318380381035850 0.318380386002903
14.000000000000000 0.318380384761140 0.318380386002903
15.000000000000000 0.318380385692462 0.318380386002903
16.000000000000000 0.318380385925293 0.318380386002903
17.000000000000000 0.318380385983501 0.318380386002903
18.000000000000000 0.318380385998053 0.318380386002903
19.000000000000000 0.318380386001691 0.318380386002903
20.000000000000000 0.318380386002600 0.318380386002903
The last, or Aitken, column evidently converges faster.