/* Program to accompany Wicklin, Rick (2015), "The sensitivity of Newton's method to an initial guess," The DO Loop blog, published 24JUN2015. http://blogs.sas.com/content/iml/2015/06/24/sensitivity-newtons-method.html Define a fifth-degree polynomial with five roots. For each point x, apply Newton's method and count how many iterations it takes to converge to a root. For each value of x, record the number of iterations and the root to which the point converges. Use SGPLOT to plot Y=number of iterations versus X=initial condition and color-code by the root to which Newton's method converges. For the examples in the blog post, f(x) = 1*x##5 -2*x##4 -10*x##3 +20*x##2 +9*x -18 which has roots at {-3 -1 1 2 3}. */ proc iml; /* evaluate polynomial with coefficients in vector p */ start Func(x) global(p); y = j(nrow(x), ncol(x), p[1]); /* initialize to coef of largest deg */ do j = 2 to ncol(p); /* Horner's method of evaluation */ y = y # x + p[j]; end; return(y); finish; /* evaluate derivative of polynomial with coefficients in vector p */ start Deriv(x) global(p); d = ncol(p)-1; p0 = p; p = (d:1) # p0[,1:d]; /* coefficients of derivative polynomial */ y = Func(x); /* temporarily overwrite p; evaluate derivative */ p = p0; /* restore p */ return( y ); finish; /* Implementation of Newton's method with default arguments */ /* By default, maximum iterations (maxiter) is 25 */ /* convergence criterion (converge) is 1e-6 */ /* Modified from http://blogs.sas.com/content/iml/2011/08/05/using-newtons-method-to-find-the-zero-of-a-function.html */ start NewtonMethod(x, iter, x0, maxIter=25, converge=1e-6); x = x0; f = Func(x); /* evaluate function at starting values */ do iter = 1 to maxiter /* iterate until maxiter */ while(max(abs(f))>converge); /* or convergence */ J = Deriv(x); /* evaluate derivatives */ delta = -solve(J, f); /* solve for correction vector */ x = x + delta; /* the new approximation */ f = func(x); /* evaluate the function */ end; /* return missing if no convergence */ if iter > maxIter then x = j(nrow(x0),ncol(x0),.); finish NewtonMethod; /* Iterate Newton's method from x. Return number of iterations and the root number to which the method converged. The roots are in the vector g_roots */ start NewtonConvergence(numIters, rootNumber, x) global(g_roots); numIters = j(1,ncol(x),.); rootNumber = j(1,ncol(x),.); epsilon = 0.01; do i = 1 to ncol(x); x0 = x[i]; call NewtonMethod(root, iter, x0, 50); idx = loc(root>g_roots-epsilon & root