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.