Muller’s Method (and divided differences)

Andy Long
Department of Mathematics and Statistics
Northern Kentucky University

Let’s start by defining divided differences:

dividedDifference[list_]:=
Module[
{n=Length[list]},
b=If[n==1,
list[[1]][[2]], (* just report the y value -- 0th order derivative *)
(dividedDifference[list[[2;;n]]]-dividedDifference[list[[1;;n-1]]])/
(list[[n]][[1]]-list[[1]][[1]]) 
(* One can think of these as (n-1)th derivative approximations *)
]
]
(* Testing against the text, p. 217 *)
dividedDifference[{{0,0}}]
dividedDifference[{{0,0},{2,20.134}}] (* derivative approximation *)
dividedDifference[{{2,20.134},{0,0}}] (* order doesn't matter *)
dividedDifference[{{0,0},{2,20.134},{3,30.452}}] (* 2nd deriv approx *)
dividedDifference[{{2,20.134},{3,30.452},{0,0}}] (* 2nd deriv approx *)
dividedDifference[{{2,20.134},{0,0},{3,30.452}}] (* 2nd deriv approx *)
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}}]]

Let’s demonstrate that this quadratic actually fits the function:

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)
root=x3-2*c/(b+Sign[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 let’s wrap this into Muller’s method:

mullers[fx_,a3_,a2_,a1_,maxits_,eps_,real_,plot_,verbose_]:=
Module[{
i=0,
x1=N[a1],x2=N[a2],x3=N[a3],
y1,y2,y3,
a,b,c,
root},
{y1,y2,y3}=N[{fx[x1],fx[x2],fx[x3]}];
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;
root=x3-2*c/(b+Sign[b]Sqrt[b^2-4*a*c]);
If[real,root=Re[root]];
{x1,x2,x3}={x2,x3,root};
{y1,y2,y3}={y2,y3,fx[root]};
If [verbose,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]
];
{x3,y3}
]

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,True]

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

This function has no real roots! The only roots we can find will be complex....

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