MATH 4650: Homework #3 - Solutions

1. (pg. 59 #1.4.EX.5) Consider the equation 8x4 − 12x3 + 6x2 − x = 0. For each of the two solutions x = 0 and x = 1/2, decide which will converge faster (say, to eight-place accuracy), the Bisection Method or Newton’s Method, without running the calculations. First of, we check that x = 0 and x = 1/2 are solutions of 8x4 − 12x3 + 6x2 − x = 0. Clearly, f (x) is twice continuously differentiable with f 0 (x) = 32x3 − 36x2 + 12x − 1 f 00 (x) = 96x2 − 72x + 12 For x = 0 we have f 0 (0) 6= 0 such that we satisfy the conditions for Theorem 1.11. Hence, Newton’s Method converges quadratically and is faster that the Bisection Method. For x = 1/2, we have f 0 (1/2) = 0 so we cannot use Theorem 1.11. Since f 00 (1/2) = 0 and f 000 (1/2) 6= 0 we know that x = 1/2 is a root of multiplicity three. Using Theorem 1.12, Newton’s Method is locally convergent to r = 1/2 and the convergence rate is S=

(m − 1) 2 = . m 3

Since this is a slower convergence rate than the Bisection Method, which has a convergence rate of S = 1/2, the Bisection Method will converge faster to eight-place accuracy. Note that

1 8x4 − 12x3 + 6x2 − x = 8x(x − )3 , 2 which is another way to see that 0 is a simple root and 1/2 is a root of multiciplity 3. 2. (pg. 59 #1.4.EX.7) Let f (x) = x4 − 7x3 + 18x2 − 20x + 8. Does Newton’s Method converge quadratically to the root r = 2? Find limi→∞ ei+1 /ei , where ei denotes the error at step i. To determine the multiplicity of the root, we take derivatives and evaluate these at the root. f 0 (x) = 4x3 − 21x2 + 36x − 20

f 0 (2) = 0

f 00 (x) = 12x2 − 42x + 36

f 00 (2) = 0

f 000 (x) = 24x2 − 42

f 000 (2) = 6

Thus we have that x = 2 is a root of multiplicity three. The convergence is not quadratic. Using Theorem 1.12, Newton’s Method is locally convergent to r = 2 and the convergence rate is ei+1 (m − 1) 2 = = . i→∞ ei m 3 lim

3. (pg. 60 #1.4.CP.13) For the function  1 3 3 f (x) = 1 − 4x (a) Find the root of the function. (b) Apply Newton’s Method using an initial guess near the root and plot the first 50 iterates. This is another way Newton’s Method can fail, by producing a chaotic trajectory.

(c) Why are Theorems 1.11 and 1.12 not applicable? (a) To find the root, note that  1  1 3 3 4x − 3 3 1− = 4x 4x such that the root is r = 34 .

(b) Newton’s Method fails to converge. Left figure is the function f . We see that the slope of the curve is infinite at the root r = 3/4. Right figure are the iterates generated by Newton’s method. We start at 0.750001, very close from the root r = 3/4, and we observe chaotic iterations. f(x) = ( 1 − 3/( 4 * x ) )1/3

1 0.8

1

0.6 0.8 Current approximation (xi)

0.4 0.2 0 −0.2

0.6

0.4

−0.4 0.2

−0.6 −0.8 −1

0

0

0.5

1

1.5

2

2.5

3

3.5

0

4

5

10

15

20 25 30 Iteration count (i)

35

40

45

50

(c) Since the derivative 1

f 0 (x) =

4x2 1 −

3 4x


2 3

is not defined for x = 0, the theorems do not apply. Note 1: In this case, Newton’s method applied to f reduces to the Fixed Point Iteration scheme xn+1 = g(xn ) with g(x) = 4x(1 − x), since  xn+1

f (xn ) = xn − 0 = xn − f (xn )

1−

3 4xn 1

1 3

= 4xn − 4x2n .

 2 3 4x2n 1− 4x3 n

Now we see that r = 3/4 is a fixed point of g (since g(r) = r), but g 0 (r) = −32 and so the theorem of convergence for Fixed Point Iteration does not apply. Note 2: Be careful with Matlab and the notation ˆ(1/3) for the cubic root. The cubic root is a nice function defined from R to R, however, when you use the notation ˆ(1/3), Matlab returns a complex root. (Any of the three.) If you want to force Matlab to return the real cubic root for a real input number you need to use the builtin nthroot(.,3) function. 4. (pg. 66 #1.5.EX.6) If the Secant Method converges to r, f 0 (r) 6= 0, and f 00 (r) 6= 0, then the approximate error relationship ei+1 ≈ |f 00 (r)/(2f 0 (r))|ei ei−1 can be shown to hold. Prove that if in addition limi→∞ ei+1 /eαi √ exists and is nonzero for some α > 0, then α = (1 + 5)/2 and ei+1 ≈ |(f 00 (r)/2f 0 (r))|α−1 eαi . We are given that ei+1

00 f (r) ≈ 0 ei ei−1 2f (r)

and ei+1 ≈ M eαi where M is some nonzero constant. Thus 00 f (r) ≈ 0 ei ei−1 2f (r) 00 f (r) α−1 M ei ≈ 0 ei−1 2f (r) 1 1 −1 f 00 (r) α−1 α−1 ei ≈ M α−1 0 ei−1 2f (r) M eαi

For large enough i, we have that ei ≈ M eαi−1 such that α= which has a positive solution of α = (1 +

1 α−1

√ 5)/2. Moreover, we have the condition on M such that

00 1 f (r) α−1 M =M 2f 0 (r) 00 α−1 −1 f (r) M =M 0 2f (r) 00 f (r) α M = 0 2f (r) 00 1 f (r) α M = 0 2f (r) −1 α−1

Since α = 1/(α − 1) then 00 f (r) α−1 M = 0 . 2f (r)

5. (pg. 66 #1.5.EX.7) Consider the following four methods for calculating 21/4 , the fourth root of 2. (i) Bisection Method applied to f (x) = x4 − 2 (ii) Secant Method applied to f (x) = x4 − 2 (iii) Fixed Point Iteration applied to g(x) = (iv) Fixed Point Iteration applied to g(x) =

x 2 x 3

+ +

1 x3 1 3x3

For each of these methods, (a) Determine the speed of convergence (b) Are there any methods that will converge faster than all the above methods? If so, name it. (a) The speed of convergence is: (i) The Bisection Method will have S = 1/2 for the rate of convergence since it cuts the interval in half at every step. (ii) It is clear that f (21/4 ) = 0. Moreover, f 0 (x) = 4x3 and f 0 (21/4 ) 6= 0. Thus the Secant Method will obtain superlinear convergence.

(iii) We first check that x = 21/4 is indeed a fixed point for g(x). 21/4 1 + 1/4 3 2 (2 ) 1 g(21/4 ) = 2−3/4 + 3/4 2 1/4 1− 34 g(2 ) = 2 g(21/4 ) =

1

g(21/4 ) = 2 4 Note that g 0 (x) = 12 − x34 and g 0 (21/4 ) = −1. Thus Fixed Point Iteration on this function may not converge. (The convergence theorem does not apply.) (iv) We first check that x = 21/4 is indeed a fixed point for g(x). 21/4 1 + 1/4 3 3(2 )3 2 1 g(21/4 ) = + 3/4 3(2 ) 3(23/4 ) 3 g(21/4 ) = 3(23/4 ) g(21/4 ) =

3

g(21/4 ) = 2− 4 Thus this is not a fixed point for g(x) and the Fixed Point Iteration Method is not applicable. The fastest convergence will be obtained via the Secant Method followed by the Bisection Method. The first Fixed Point Iteration Methods may fail to converge. The second Fixed Point Iteration Methods does not converge to 21/4 . (b) Let f (x) = x4 − 2, then note that f 0 (x) = 4x3 and f 0 (21/4 ) 6= 0. Therefore we can apply Theorem 1.11 with Newton’s Method and obtain quadratic convergence.

6. (pg. 66 #1.5.CP.3) Use Inverse Quadratic Interpolation to find the solution of each of the following equations. (a) x3 = 2x + 2 (b) ex + x = 7 (c) ex + sin x = 4 Here is a simple implementation of the Inverse Quadratic Interpolation Method in Matlab which does not include any output: function c = iqi( f, a, b, c, maxiter, tol) for i = 1:maxiter A = f(a); B = f(b); C = f(c); if (A ~= 0) && (B ~= 0) q = A / B; r = C / B; s = C / A; d = c - (r*(r-q)*(c-b) + (1-r)*s*(c-a))/((q-1)*(r-1)*(s-1)); if abs(d-c) > tol a = b; b = c; c = d; else c = d; return; end

end

end

else error('Division by zero in A or B.'); end

(a) In Matlab, we type: >> clear; f = @(x) x^3-2*x-2; iqi(f, 0, 0.5, 1, 20, 1e-08); and we get: k | x | f(x) | F.E. ----+--------------------------+--------------------------+-------------------------0 | 1.0000000000000000e+00 | -3.0000000000000000e+00 | 7.6929235423863140e-01 1 | 1.8571428571428584e+01 | 6.3661049562682347e+03 | 1.6802136217189954e+01 2 | -1.0994579383928041e+01 | -1.3090441270803424e+03 | 1.2763871738166673e+01 3 | 1.0242642072592520e+00 | -2.9739552518844068e+00 | 7.4502814697937936e-01 4 | 1.0483612113445377e+00 | -2.9445092608391836e+00 | 7.2093114289409366e-01 5 | 3.4634130576879527e+00 | 3.2617610124867632e+01 | 1.6941207034493213e+00 6 | 3.2733533786878537e+00 | 2.6526758624155590e+01 | 1.5040610244492223e+00 7 | 1.3679472378962172e+00 | -2.1760826531696327e+00 | 4.0134511634241421e-01 8 | 1.5707689889139937e+00 | -1.2659557498209373e+00 | 1.9852336532463766e-01 9 | 1.8373782097475406e+00 | 5.2815650653550561e-01 | 6.8085855508909221e-02 10 | 1.7772506582543122e+00 | 5.9157983967704109e-02 | 7.9583040156807616e-03 11 | 1.7691853942768281e+00 | -7.9050028159910468e-04 | 1.0695996180332124e-04 12 | 1.7692923053337173e+00 | -3.6146531856573461e-07 | 4.8904914118352849e-08 13 | 1.7692923542386685e+00 | 2.7400304247748863e-13 | 3.7081449022480228e-14 14 | 1.7692923542386314e+00 | 0.0000000000000000e+00 | 0.0000000000000000e+00 (b) In Matlab, we type: >> clear; g = @(x) exp(x) + x - 7; iqi(g, 0, 0.5, 1, 20, 1e-08); and we get: k | x | f(x) | F.E. ----+--------------------------+--------------------------+-------------------------0 | 1.0000000000000000e+00 | -3.2817181715409545e+00 | 6.7282169862890640e-01 1 | 1.3618973822738749e+00 | -1.7345097301229337e+00 | 3.1092431635503148e-01 2 | 1.6129955356894985e+00 | -3.6918466876895284e-01 | 5.9826162939407901e-02 3 | 1.6699011744557157e+00 | -1.8455980297245489e-02 | 2.9205241731906639e-03 4 | 1.6728096459337063e+00 | -7.6259164612402230e-05 | 1.2052695200104324e-05 5 | 1.6728216981750650e+00 | -2.8715358979525263e-09 | 4.5384140889836999e-10 6 | 1.6728216986289064e+00 | -8.8817841970012523e-16 | 0.0000000000000000e+00 (c) In Matlab, we type: >> clear; h = @(x) exp(x) + sin(x) - 4; iqi(h, 0, 0.5, 1, 20, 1e-08); and we get:

k | x | f(x) | F.E. ----+--------------------------+--------------------------+-------------------------0 | 1.0000000000000000e+00 | -4.4024718673305774e-01 | 1.2998049865083239e-01 1 | 1.1235150422854778e+00 | -2.2727427797601951e-02 | 6.4654563653545605e-03 2 | 1.1299494074410108e+00 | -1.0951068775177930e-04 | 3.1091209821587640e-05 3 | 1.1299804969593306e+00 | -5.9579323696823394e-09 | 1.6915018274943350e-09 4 | 1.1299804986508324e+00 | 0.0000000000000000e+00 | 0.0000000000000000e+00

7. (pg. 78 #2.1.EX.2) Use Gaussian elimination to solve the systems: (a) 2x − 2y − z = −2 4x + y − 2z = 1 −2x + y − z = −3 (b) x + 2y − z = 2 3y + z = 4 2x − y + z = 2 (c) 2x + y − 4z = −7 x − y + z = −2 −x + 3y − 2z = 6 (a) 

 2 −2 −1 −2  4 1 −2 1  −2 1 −1 −3

 ⇒

 2 −2 −1 −2  0 5 0 5  0 0 −2 −4

then using back substitution we have (x, y, z) = (1, 1, 2). (b) 

 1 2 −1 2  0 3 1 4  2 −1 1 2

 ⇒

1 2 −1  0 3 1 0 0 14 3

 2 4 

14 3

then using back substitution we have (x, y, z) = (1, 1, 1). (c) 

 2 1 −4 −7  1 −1 1 −2  −1 3 −2 6

 ⇒

 2 1 −4 −7 3   0 −3 3 2 2 0 0 3 6

then using back substitution we have (x, y, z) = (−1, 3, 2).

8. (pg. 78 #2.1.EX.3) Solve by back substitution: (a) 3x − 4y + 5z = 2 3y − 4z = −1 5z = 5

(b) x − 2y + z = 2 4y − 3z = 1 −3z = 3 (a) We have (x, y, z) = (1/3, 1, 1). (b) We have (x, y, z) = (2, −1/2, −1).

9. (pg. 78 #2.1.EX.5) Use the approximate operation count 2n3 /3 for Gaussian elimination to estimate how much longer it takes to solve n equations in n unknowns if n is tripled. If n increases to 3n, the approximate operation count changed from 2n3 /3 to 2(3n)3 /3 = 54n3 /3, which will take 27 times as long.