==================================== CHAPTER 1 PROBLEMS GERALD AND WHEATLEY, 6TH EDITION ==================================== 2. This function is an upwards-cupped parabola centered at x=0.5. It crosses the x axis at x=0.4 and x=0.6. f(0)=0.24 and f(1)=0.24, so these guesses do not bracket the root and thus are not satisfactory for starting up the bisection method. In order that the guesses bracket the root, one of them has to be between 0.4 and 0.6, since that's the only region where f(x) is negative. The other guess has to be outside that range, and the root to which the method converges will depend on which side it's on. For example, if a=0.5 and .6 < b < infinity, bisection will find the root at 0.6. If you start with [0.5,1], the first guess is 0.75 and the the absolute error bound is 0.25. For the second iteration, the error bound is 0.125, and so on. After the sixth iteration, it's (1-0.5)/2^6, which is .0078125. To figure out the actual error, you have to know the root. This is unrealistic in practice, but instructive to do once or twice in textbook problems. Here's the iteration table: a b f(a) f(b) c f(c) 0.5 1 -0.01 0.24 0.75 .0525 0.5 0.75 -0.01 .0525 .625 5.6e-3 0.5 .625 -0.01 5.6e-3 .5625 -6.1e-3 .5625 .625 -6.1e-3 5.6e-3 .59375 -1.2e-3 .59375 .625 -1.2e-3 5.6e-3 .609375 1.97e-3 .59375 .609375 -1.2e-3 5.6e-3 .6015625 whew! The real root is at .60000, so the actual error is .0015625. Note that the actual error is always (by definition) within the worst-case error bounds. ------------------- 4(a). This problem is a bit confusing. The root of a function is just the x that makes it true. The equation can be of the form f(x)=0, or it can be f(x)=g(x), as it is in the case of this problem. However, our rootfinding methods are designed to solve equations of the form f(x)=0; if you have one that looks like f(x)=g(x), you have to convert it (in this case, by subtracting e^x from both sides). I was doing this problem on an airplane, and I didn't have graphing software, so I used Scheme: (pp (map (lambda (x) (- 2 (sin (* 2 x)) (exp x))) '(0 0.1 0.2 0.3 0.4 0.5 0.6))) (1 .696159751129291 .38917889953117957 .08549871902896156 -.20918078854079303 -.4901922555080249 -.754157886357735) (note that x is assumed to be in radians, not degrees, in the sin2x) Clearly the function crosses over between 0.3 and 0.4, so I just used those values as the initial guesses. You can use whatever values you want for a and b as long as f(a) times f(b) < 0. Here's my iteration table: a b f(a) f(b) c f(c) abserr relerr 0.3 0.4 8.5e-2 -2.1e-1 .35 -6.3e-2 .05 14.3% 0.3 .35 8.5e-2 -6.3e-2 .325 1.1e-2 .025 7.7% .325 .35 1.1e-2 -6.3e-2 .3375 -2.6e-2 .0125 3.7% .325 .3375 1.1e-2 -2.6e-2 .33125 -7.8e-3 .00625 1.8% .325 .33125 1.1e-2 -7.8e-3 .328125 1.49e-3 .003125 .95% .328125 .33125 1.49e-3 -7.8e-3 .3296875 .0015625 .47% The absolute error is just the half-width of the a-b interval at that step, since that's how far c (the current best guess) can be from the real root. This error is related to the original interval width as |b-a|/2^n, where n is the step number. The way to calculate the relative error is to divide the absolute error at each step by the current best guess (c) -- e.g., 0.05/0.35 = 0.143 at the first step. ------------------- 6(a). Starting from the same two guesses: a b f(a) f(b) c f(c) 0.3 0.4 8.5e-2 -2.1e-1 .3290 -1.16e-3 0.3 0.3290 8.5e-2 -1.2e-3 .32863 -5.6e-6 Note the size of f(c). This got closer to the root in two steps than bisection did in six steps. Here's the Scheme function for the 'next guess' formula on page 44, if you want to try it: (define (secant-next-guess a b fcn) (- a (/ (* (- b a) (fcn a)) (- (fcn b) (fcn a))))) Here's how you call it in this case: (secant-next-guess 0.3 0.4 (lambda (x) (- 2 (sin (* 2 x)) (exp x)))) ------------------- 12(a). f(x) = 2 - sin2x - e^x and f(x) = -2cos2x - e^x You fire up Newton with one guess only. I'll arbitrarily pick 0.3: a f(a) f'(a) c f(c) 0.3 8.5e-2 -3.000 .3285 3.9e-4 This converged even faster than secant, and much faster than bisection. ------------------- 29. For the first iteration, I get x_1 = - sqrt (e^0 / 3), which is -.57735. Then I plug that back in, yielding x_1 = - sqrt (e^(-.57735) / 3), which is -.43258. Continuing in this fashion, I get these iterates: -.57735 -.43258 -.4650659315428768 -.45757 -.45928 -.45889 -.45898 -.45896 -.458963 -.458962 ... It looks to me like it's converging, albeit slowly. For the other rearrangement, I get .57735 .77057 .84872 .88255 .89760 .90438 .90745 .90884 .90948 .90977 This too is converging, and to a different root than the sequence above. Here's what happens if I start the second sequence up near 4: 4.26607 4.87310 6.60116 15.6625 1453.84 7.741e153 went downhill really fast in this last step ... Another possible rearrangement is x = ln (3x^2). Here's what happens if I start that up at x_0=4: 3.8712 3.8057 3.7716 3.7536 3.7441 3.7390 3.7362 3.7348 3.7340 3.7336 ... This rearrangement does indeed converge to that third real root. The point of this exercise is to show you that..... (i) different rearrangements give you iteration sequences that converge to different roots (ii) sometimes they give you iteration sequences that diverge (iii) whether or not an iteration sequence converges can depend on the initial condition you use to start it up. ------------------- 30. The answers in the back of the book say 13 iterations, but my computer took 12: 1 .5 .6666666666666666 .6000000000000001 .625 .6153846153846154 .6190476190476191 .6176470588235294 .6181818181818182 .6179775280898876 .6180555555555556 .6180257510729613 .6180371352785146 .6180327868852459 .6180344478216819 .6180338134001252 ^^^^^^^^^^^^^^^^^ .6180340557275542 .6180339631667064 .6180339985218034 .6180339850173578 .6180339901755971 .6180339882053251 .618033988957902 .6180339886704432 (I'm assuming five-digit rounding, incidentally -- not five-digit chopping.) By the way, note that this does NOT appear to be converging to .6180339, but rather to .6180340 (.6180339887498948, more precisely), so I think the book's rounded value for the root is given incorrectly on page 102. Here's how Aitken works. You first take two fixed-point steps: 1 .5 .6666666666666666 ...and then you plug those guesses into the formula at the bottom of p58: (define (aitken x2 x1 x0) (/ (- (* x0 x2) (* x1 x1)) (+ x2 (* -2 x1) x0))) ...to get the next guess: (aitken 1 .5 .6666666666666666) ;Value: .625 This is as good as the 4th iterate of the un-Aitkenized sequence, even though it's only the 3rd iterate of the combined fixed point/Aitken method. This is a win, though not a big one. But we're not done yet, because the answer is off by more than .00001. In order to continue, you take two MORE fixed-point steps, starting from the .625 iterate: .625 .6153846153846154 .6190476190476191 ...and plug those values into the Aitken formula: (aitken .625 .6153846153846154 .6190476190476191) ;Value: .6180371352785142 Here Aitken really accelerated things. After four fixed-point steps and two Aitken steps, the answer is as good --- five decimal places -- as it was after twelve plain fixed-point steps. That's a 50% improvement if the amount of computer arithmetic involved in an Aitken step is comparable to the algebra involved in a simple fixed-point step (which it is). ------------------- 56. The point of this problem is to show you two things: that some initial conditions cause Newton's method to diverge, and that multiple slow its convergence. In particular, Newton converges quadratically to single roots, but only linearly to multiple roots. From the initial condition x=2, my code diverges. (If you made your code more intelligent, it may have done something different, of course!) One can, in theory, predict this in advance by evaluating the derivative f'(x) at the root, but that only works if you know the root, and this root is unstable (like the upright-balanced point of a pendulum), which makes it hard to find. In practice, though, I wouldn't claim to be able to predict this in advance. From x=2.9, there is LINEAR convergence to the root at x=3: errors decrease roughly in the ratio of 2/3. This is very different from Newton's method's convergence pattern to a single root, where the errors decrease not in a linear ratio, but in proportion to the SQUARE of the previous error (that is, not linearly more sig figs of correctness per iteration, but some multiplicative factor of correct sig figs per iteration). What causes this is the multiplicity of the root: the polynomial has three roots on top of each other at x=3, so the first, second, and third derivatives are zero at that point. Since Newton "uses" the slope as leverage to move towards the root, this kind of inflection point really slows it down.