Muller’s Method (and divided differences)

Andy Long
Department of Mathematics and Statistics
Northern Kentucky University
longa@nku.edu                 
3/2015

We start by defining divided differences:

See pages 212-- of our text: Elementary Numerical Computing with Mathematica (Skeel and Keiper)

divDiff[list_]:=
Module[
{n=Length[list]},
If[n==1,
b=list[[1]][[2]],
b=(divDiff[list[[2;;n]]]-divDiff[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
Print[list," ",b];
b
]
(* Testing against the text, p. 217 *)
divDiff[{{0,0},{2,20.134},{3,30.452},{5,52.11}}]

Now here’s the production model of divided differences:

dividedDifference[list_]:=
Module[
{n=Length[list]},
If[n==1,
b=list[[1]][[2]],
b=(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]])];
b
]
(* Testing against the text, p. 217 *)
dividedDifference[{{0,0}}]
dividedDifference[{{0,0},{2,20.134}}]
dividedDifference[{{2,20.134},{0,0}}] (* order doesn't matter *)
dividedDifference[{{0,0},{2,20.134},{3,30.452}}]
dividedDifference[{{0,0},{2,20.134},{3,30.452},{5,52.11}}]
dividedDifference[Reverse[{{0,0},{2,20.134},{3,30.452},{5,52.11}}]]

This demonstrates that my derivation of the quadratic actually works:

f[x_]:=Exp[x]-1/x^2
x1=1;y1=N[f[x1]];
x2=.5;y2=N[f[x2]];
x3=.3;y3=N[f[x3]];
a=dividedDifference[{{x2,y2},{x3,y3},{x1,y1}}]
b=dividedDifference[{{x2,y2},{x3,y3}}]-(x2-x3)*a
c=y3
q[x_]:=y3+(x-x3)*(b+(x-x3)*a)
If[b>0,
root=x3-2*c/(b+Sqrt[b^2-4*a*c]),
root=x3-2*c/(b-Sqrt[b^2-4*a*c])]
p1=Plot[{f[x],q[x]},{x,Min[{x1,x2,x3,root}],Max[{x1,x2,x3,root}]}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3},{root,q[root]}},PlotMarkers->{Automatic,Large}];
Show[p1,p2]

Now I’m ready to wrap that into Muller’s method:

mullers[fx_,a3_,a2_,a1_,maxits_,eps_,real_,plot_]:=
Module[{i=0,x1=N[a1],x2=N[a2],x3=N[a3],a,b,c},
y3=N[fx[x3]];
y2=N[fx[x2]];
y1=N[fx[x1]];
While[i<maxits && Abs[y3]>eps,
i=i+1;
a=dividedDifference[{{x2,y2},{x3,y3},{x1,y1}}];
b=dividedDifference[{{x2,y2},{x3,y3}}]-(x2-x3)*a;
c=y3;
If[Re[b]>0,
root=x3-2*c/(b+Sqrt[b^2-4*a*c]),
root=x3-2*c/(b-Sqrt[b^2-4*a*c])];
If[real,root=Re[root]];
x1=x2;
y1=y2;
x2=x3;
y2=y3;
x3=root;
y3=fx[root];
Print[{x3,y3}];
];
If[plot,
q[x_]:=Re[y3+(x-x3)*(b+(x-x3)*a)];
p1=Plot[{fx[x],q[x]},{x,x1,x3}];
p2=ListPlot[{{x1,y1},{x2,y2},{x3,y3},{root,q[root]}},PlotMarkers->Automatic];
Show[p1,p2]
]
]

And here’s our test case:

f[x_]:=Exp[x]-1/x^2
Plot[f[x],{x,.1,5}]
N[Solve[f[x]==0,x]]
mullers[f,4,3,2,100,10^-6,True,True]

Here’s a test case to illustrate finding a complex root:

f[x_]:=x^4+1
mullers[f,4,3,2,100,10^-12,False,False]

I’ve talked about deflation, to find further roots. Here’s an illustration of how it works:

g[x_]:=f[x]/(x-(0.7071067811865472` -0.7071067811865476` I))
mullers[g,4,3,2,100,10^-12,False,False]

g[x_]:=f[x]/(x-(0.7071067811865472` -0.7071067811865476` I))/(x-(0.7071067811865475` +0.7071067811865475` I))
mullers[g,4,3,2,100,10^-12,False,False]

g[x_]:=f[x]/(x-(0.7071067811865472` -0.7071067811865476` I))/(x-(0.7071067811865475` +0.7071067811865475` I))/(x-(-0.7071067811865248`-0.70710678118659` I))
mullers[g,4,3,2,100,10^-12,False,False]

g[x_]:=f[x]/(x-(0.7071067811865472` -0.7071067811865476` I))/(x-(0.7071067811865475` +0.7071067811865475` I))/(x-(-0.7071067811865248`-0.70710678118659` I))/(x-(-0.7071067811865221`+0.7071067811865883` I))
mullers[g,4,3,2,100,10^-12,False,False]