Project 2 Key

Problem 1: Newton’s method intervals

f[x_]:=x*Exp[x]
Plot[f[x],{x,-3,1}]
Solve[f'[x]==0,x]

This one is very straightforward: the general solution (the characterization of the basins) gives the whole result. 

a. If x0>=0, convergence to the root at 0 is monotonic

pN[0]=1;
its=5
Do[pN[n]=N[pN[n-1]-((f[pN[n-1]])/(f'[pN[n-1]])),10],{n,1,its}]
TableForm[Table[{k,pN[k]},{k,0,its}]]

b. If x0 is in the open interval (-1,0), then the first iteration of Newton's sends x1 to the right of 0, so that, from that point forward, convergence to the root at 0 is monotonic.

pN[0]=-.5;
Do[pN[n]=N[pN[n-1]-((f[pN[n-1]])/(f'[pN[n-1]])),10],{n,1,its}]
TableForm[Table[{k,pN[k]},{k,0,its}]]

c. If x0 is at -1, the derivative is 0, and so Newton's method produces a divide by zero.

pN[0]=-1;
its=1;
Do[pN[n]=N[pN[n-1]-((f[pN[n-1]])/(f'[pN[n-1]])),10],{n,1,its}]

d. If x0 is <-1, Newton's method gives iterates to the left, heading off to negative infinity.

pN[0]=-2;
its=5
Do[pN[n]=N[pN[n-1]-((f[pN[n-1]])/(f'[pN[n-1]])),10],{n,1,its}]
TableForm[Table[{k,pN[k]},{k,0,its}]]

Problem 2: Bisection:

What are lawful intervals containing multiple roots?

Let’s think about what we know (and what we can deduce):

f[x_]:=(x+2)(x+1)x(x-1)(x-2)
Plot[f[x],{x,-Pi,Pi}]
Bisection[f_,a0_,b0_,true_]:=
Module[{i=1,its=100,eps=.000001,a,b,c},
If[Abs[f[a0]]≤eps,Return[a0]];
If[Abs[f[b0]]≤eps,Return[b0]];
a=N[a0]; (* Slight change -- use a local variable (so as not to alter passed arguments *)
b=N[b0];
fa = f[a];  
fb = f[b];  
If[fa*fb>=0,Return[False]];
width=b-a;
c = a+0.5*width; 
fc=f[c]; 
If[true,Print[a," ",b," ",c]];
If[Abs[fc]≤eps, Return[c]];
While[a < c && c<b  && i<its,
If[ Sign[fa] ==  Sign[fc], 
a = c;fa=fc, 
b = c;fb=fc
];
width=b-a;
slope =(fb-fa)/width;
i=i+1;
c = a+0.5*width;
fc=f[c];
If[true,Print[a," ",b," ",c]];
If[Abs[fc]≤eps, Return[c]];
]; 
Return[c];
]
  
Bisection[f,-Sqrt[2],Pi/2,True]
Bisection[f,-Sqrt[2],-.5,True]
Plot3D[{Bisection[f,x,y,False]},{x,-5,0},{y,0,-x},PlotPoints->{20,20}]

Problem 3: comparison of fixed point methods

First of all, check that each of these is indeed a fixed point method for the root problem: 
f(x)=0 -> e^x=1/x^2
1. g1(x)-x=x^3f(x)
2. g2(x)-x gives x=-2Ln(|x|), so that -x/2=Ln(|x|), i.e. e^(-x/2)=|x|
and from g4(x), 
3. g4(x)-x=e^(-x/2)-x=0 -> x=e^(-x/2) -> e^(x/2)=1/x, or that e^x=1/x^2
and 
4. g3(x) is Newton’s method’s function for f, so that’s a valid fixed point method. 

The rate of convergence (or asymptotic error constant) is determined by the derivatives of these functions. Those that are less than 1 in absolute value will converge; otherwise we’ll get divergence. If the derivative is negative, then the fixed point will spiral toward the root (not approach monotonically).

Clear[p1,p2,p3,p4,p0,k,g,a,b]

f[x_]:=Exp[x]-1/x^2
value=x/.N[Solve[f[x]==0,x,WorkingPrecision->100]][[2]]
g1[x_]=x^3Exp[x];
g2[x_]=-2Log[Abs[x]];
g3[x_]=x*(x^2(x-1)Exp[x]+3)/(x^3Exp[x]+2);
g4[x_]=Exp[-x/2];
p0=.25;
a=0;
b=1;
n=15;

Now let us do the iteration:

p1[0]=N[p0];
p2[0]=N[p0];
p3[0]=N[p0];
p4[0]=N[p0];

Do[p1[k]=N[g1[p1[k-1]]],{k,1,n}]
Do[p2[k]=N[g2[p2[k-1]]],{k,1,n}]
Do[p3[k]=N[g3[p3[k-1]]],{k,1,n}]
Do[p4[k]=N[g4[p4[k-1]]],{k,1,n}]

Let’s write out the values:

TableForm[Table[{k,p1[k],p2[k],p3[k],p4[k]},{k,0,n}]]
(* As we can see, rates of convergence are in line with expectations, based on derivatives: *) 
TableForm[
Table[
{k,
Abs[(p1[k+1]-value)/(p1[k]-value)],
Abs[(p2[k+1]-value)/(p2[k]-value)],
Abs[(p3[k+1]-value)/(p3[k]-value)],
Abs[(p4[k+1]-value)/(p4[k]-value)]},{k,0,n-1}]]

p0=1.5;

p1[0]=N[p0];
p2[0]=N[p0];
p3[0]=N[p0];
p4[0]=N[p0];

Do[p1[k]=N[g1[p1[k-1]]],{k,1,n}]
Do[p2[k]=N[g2[p2[k-1]]],{k,1,n}]
Do[p3[k]=N[g3[p3[k-1]]],{k,1,n}]
Do[p4[k]=N[g4[p4[k-1]]],{k,1,n}]

Let’s write out the values:

TableForm[Table[{k,p1[k],p2[k],p3[k],p4[k]},{k,0,n}]]

p0=.703467;

p1[0]=N[p0];
p2[0]=N[p0];
p3[0]=N[p0];
p4[0]=N[p0];

Do[p1[k]=N[g1[p1[k-1]]],{k,1,n}]
Do[p2[k]=N[g2[p2[k-1]]],{k,1,n}]
Do[p3[k]=N[g3[p3[k-1]]],{k,1,n}]
Do[p4[k]=N[g4[p4[k-1]]],{k,1,n}]

Let’s write out the values:

TableForm[Table[{k,p1[k],p2[k],p3[k],p4[k]},{k,0,n}]]

We expect convergence at a rate given by the derivatives at the root:

{g1'[value],g2'[value],g3'[value],g4'[value]}
TableForm[
Table[
{k,
Abs[p1[k+1]-value/(p1[k]-value)],
Abs[p2[k+1]-value/(p2[k]-value)],
Abs[p3[k+1]-value/(p3[k]-value)],
Abs[p4[k+1]-value/(p4[k]-value)]},{k,0,n-1}]]

Problem 4: comparison of root finders

My method for making a fair comparison was this: 
a. Bisection: start with bisection, on endpoints [a,b]. Then use the first chosen (say a) as the input for 
b. Newton’s method uses a.
c. Secant then uses a, b.
d. Muller’s uses a, b, and their midpoint (as bisection’s next approximation). 

In a way, you might then want to start counting convergence from the first iteration after the starting values, because the methods really aren’t being used until they get past their initial starting values (with which they had no input).

f[x_]:=Exp[x]-1/x^2
Plot[f[x],{x,-2,2}]
value=x/.N[Solve[f[x]==0,x]][[2]]
its=25;
Clear[pN,pS,pB,pM,a,b]
epsilon=10^-20
a=4;
b=.4;

(* Bisection: *)
pB[0]=a;
pB[1]=b;
Bisection[f_,a0_,b0_,true_]:=
Module[{i=3,eps=epsilon,a,b,c},
If[Abs[f[a0]]≤eps,Return[a0]];
If[Abs[f[b0]]≤eps,Return[b0]];
a=N[a0]; (* Slight change -- use a local variable (so as not to alter passed arguments *)
b=N[b0];
If[a>b,c=a;a=b;b=c];
fa = f[a];  
fb = f[b];  
If[fa*fb>=0,Return[False]];
width=b-a;
c = a+0.5*width; 
pB[2] = c; 
fc=f[c]; 
If[true,Print[a," ",b," ",c]];
If[Abs[fc]≤eps, Return[c]];
While[a < c && c<b  && i<=its,
If[ Sign[fa] ==  Sign[fc], 
a = c;fa=fc, 
b = c;fb=fc
];
width=b-a;
slope =(fb-fa)/width;
c = a+0.5*width;
pB[i] = c; 
i=i+1;
fc=f[c];
If[true,Print[a," ",b," ",c]];
If[Abs[fc]≤eps, Return[c]];
]; 
Return[c];
] 
Bisection[f,pB[0],pB[1],True]

(* Newton’s Method *)
pN[0]=N[a]
Do[pN[n]=pN[n-1]-((f[pN[n-1]])/(f'[pN[n-1]])),{n,1,its}]

(* Secant Method *)
pS[0]=N[a];
pS[1]=N[b];
Do[pS[n]=pS[n-1]-((f[pS[n-1]] (pS[n-1]-pS[n-2]))/(f[pS[n-1]]-f[pS[n-2]])),{n,2,its}]

(* Muller’s method *)
pM[0]=a;
pM[1]=b;
pM[2]=pB[2];
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
]
mullers[fx_,a3_,a2_,a1_,eps_,real_]:=
Module[{i=2,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<its && 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]];
pM[i]=root;
x1=x2;
y1=y2;
x2=x3;
y2=y3;
x3=root;
y3=fx[root];
Print[{x3,y3}];
];
pM[i+1]=root;
root
]
mullers[f,pM[0],pM[1],pM[2],epsilon,True]

TableForm[Table[{k,pN[k],pS[k],pB[k],pM[k]},{k,0,its}]]
TableForm[
Table[
{k,
Abs[pN[k+1]-value]/(pN[k]-value)^2,
Abs[pS[k+1]-value]/Abs[(pS[k]-value)]^((1+Sqrt[5])/2),
Abs[(pB[k+1]-value)/(pB[k]-value)],
Abs[pM[k+1]-value]/(pM[k]-value)^2},{k,0,its-1}]]