/* 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