EXPLANATION FILE OF PROGRAM MUELLER =================================== Mueller's Method in One Dimension --------------------------------- In another chapter we found that Newton's method could be improved upon by employing Aitken acceleration, and subsequently Aitken-Steffenson iteration. These improvements were based on the observation of a weak point in Newton's method--the assumption that the function was linear near the root. However, many nontrivial functions have curvature which causes error, or at least slows convergence, in Newton's algorithm. The solution is to at least partially ac- count for this curvature. That was the basis of the Aitken acceleration con- cept. In this section, we will consider an algorithm called Mueller's method. This algorithm assumes the function to be parabolic in the neighborhood of the root. Because the complex-plane implementation of this technique is algebrai- cally complicated, we will first examine the procedure for one-dimensional functions f(x), then extend it iteratively to two dimensions, and finally pro- ceed on to the complex plane formulation. As you will see, Mueller's method is very powerful in terms of both convergence rate and stability. In Mueller's method, the idea is to find three evaluation points near the root, fit a parabola through those points, and then de termine the roots of the corresponding second degree equations. The root closest to X3 is then ta- ken as a new evaluation point: X1 is dropped X2 --> X1 X3 --> X2 Xnew --> X3 The procedure is repeated until the change in X3 is less than some chosen cri- terion |Xnew - X3| < E Note that conceptually this is a simple upgrade of the secant method in which two points are used to generate a linear equation for the intersection with the X axis. This new value is then used to replace one of the previous two points. The algebra associated with this method is somewhat complicated and is only briefly described below. (For a more complete discussion, refer to the text by Becket Hurt, Ref. 8). We will start by defining two dimensionless parameters, L and D: X3 - X2 L = ------- X2 - X1 X3 - X1 D = ------- X2 - X1 These are used to calcula te Gand C: G = L² f(X1) - D² f(X2) + (L + D) f(X3) C ~ L|L f(X1) - D f(X2) + f(X3)| The above definitions are then used to calculate the parameter lambda: -2D f(X3) lambda = ------------------------ G +/- sqrt(G²-4DC f(X3)) The sign in the denominator is chosen so that the magnitude of lambda is mini- mized. The new estimate for the location of the root is then Xnew = X3 + lambda (X3 - X2) X1 is subsequently dropped from the group and the process is repeated until some convergence criterion is met. This procedure has been implemented in program Mueller. In the first example, the function examined was y = x(x - l)(x - 2)(x - 3) (x - 4). The first two runs show very good convergence. The third run termi- nated in six iterations according to the error criteria, E = 0.001, but it was in error by much more than that. Thus the warning: do not assume that the error in the final result is less than the error criterion! Convergence may be slow in sorne situations. Such was also the case for the fourth run in which the initial guess was far from any root. After ten iterations, the error was very large. However, extending the iteration limit to 20 (fifth run) and then 30 (sixth run) indicates that despite the initial slow convergence, the algo- rithm eventuaIly, and very accurately, found a root. Why was the convergence initially so slow? The answer to this question rests in the difference between what the algorithm assumes the form of the function to be, and what the form actually is. The algorithm assumes the function to be locally parabolic. But for an initial guess of XO = -30, which is far from the cluster of five roots, the form of the function is approximately y = x^5. The algorithm has difficulty coping with this, but at least it does not Fail. Let us give an additional example of this convergence problem. In this case, the function is y = (x + 1)^5. Whether the evaluation point is near to or far from the multiple root, the functional form always has the high curvature as- sociated with a quintic function. The sample runs aIl show how the algorithm is struggling to find the root, and how it is converging slowly. These were two particularly difficult examples. Most functions do not create as much trouble for Mueller's method as these do, especially when the initial guess is reasonably near a root. You can compare the results for the y=(x+1)^5 case with a similar example given in another chapter using Newton's method. The comparison clearly shows the superior properties of the Mueller algorithm. The implementation of the one-dimensional Mueller algorithm as shown in pro- gram Mueller contains several checks to avoid errors in execution. For the most part, these error checks involve avoiding divide-by-zero failures. How- ever, one particular check is associated with a potentially fatal situation. In this case, the three evaluation points lead to a parabolic fit that does not cross the X axis; the roots are imaginary. This is not helpful to the al- gorithm and the error check substitutes an artificial set of intersections to at least keep the iteration going. In another section (see file Mueller2.txt), we will investigate the exten- sion of Mueller's method to two dimensions by using successive substitution. From [BIBLI 01] -------------------------------------- End of file Mueller.txt