William Hedden, Brooke Mardis, Joseph Hastings

Project 1 Question 1

Problem 1.

Newton’s method code taken from Professor Long:

NewtonsMethod[f_,a0_,eps_,nits_] :=
Module[
{a=N[a0], 
fa = f[N[a0]],
i=0,
output 
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<nits,
a = a-(fa/(f'[a]));
fa=f[a];
i=i+1;
output=Append[output,{i,a,fa}];  
]; 

Print[NumberForm[TableForm[output,
TableHeadings->{None,{"i","(Subscript[x, k])","f[(Subscript[x, k])]"}}],16]]; 
Print["  x  = ",NumberForm[a,16] ]; 
Print["f(x) = ",NumberForm[fa,16] ]; 
Return[Transpose[output]]; ]

f[x_]=x*E^x

p1 = ListPlot[{{0.2,f[0.2]},{-.5,f[-.5]},{-1,f[-.999999999999]},{-2,f[-2]}},PlotStyle->PointSize[.03]];

graph=Plot[f[x],{x,-5,2},PlotRange->{-2,2.5}];

p2= ListPlot[{{0.2,f[0.2]}},PlotMarkers->"A",PlotStyle->Red];

P3=ListPlot[{{-.5,f[-.5]}},PlotMarkers->"B",PlotStyle->Red];

P4=ListPlot[{{-1,f[-.99999999999]}},PlotMarkers->"C",PlotStyle->Red];

P5=ListPlot[{{-2,f[-2]}},PlotMarkers->"D",PlotStyle->Red];

Show[p2,P3,P4,P5,graph,PlotRange->Automatic]

x0=.2;
eps=10^(-10);
nits=100
NewtonsMethod[f,x0,eps,nits];

With .2 as the starting point, the points monotonically converge to the root, after each iteration we get closer and closer to the true root at (0,0). (Shown on the graph as point “A”)

x0=-.5;
eps=10^(-10);
nits=100
NewtonsMethod[f,x0,eps,nits];

With .5 as the starting point, it converges non-monotonically to the root. In fact, after the second iteration it then converges monotonically. (Shown as point “B” on the graph)

x0=-1;
eps=10^(-10);
nits=100
NewtonsMethod[f,x0,eps,nits];

With -1 as the starting point, Newton’s method blows up. When we try to divide by zero, it fails. (Shown by point “C” on the graph)

x0=-1.2;
eps=10^(-10);
nits=100
NewtonsMethod[f,x0,eps,nits];

With -1.2 as the starting point, the points monotonically flee the root because the points are getting closer and closer to the horizontal asymptote along the x axis as you get more and more negative. (Point shown as “D” on the graph)

Solve[f'[x]==0,x]

Solve[f''[x]==0,x]

Part B
 x ∈ (0, ∞) converges monotonically to the root.
 x ∈ (-1,0) converges non-monotonically to the root
 At x=-1 Newton’s method will “blow up”
 x ∈ (-∞,-1) monotonically flees the root

Problem 4

NSolve[E^x-1/x^2==0,x]

Plot[E^x-1/x^2,{x,-3,3},PlotRange->Automatic]

g[x_]:=E^x-1/x^2

x0=.7;
eps=10^(-10);
nits=100
NewtonsMethod[g,x0,eps,nits];

Project 1 Question 2

Bisection method stolen from you, then altered to match the algorithm on page 79, and altered to stop calculations and inform you if you happen to land on a root with your range.

ClearAll[a0,b0,m,f]

Bisection[a0_,b0_,m_] :=
Module[{},
a = N[a0];  
b = N[b0];  
c = ((a + b)/2); 
k = 0; 
output={{k,a,c,b,f[c]}}; 
If[f[a] == 0, k=m+1; Return[a" is a root"], If[f[b]==0, k = m+1; Return[b" is a root"],If[f[c] == 0, k=m+1; Return[c" is a root"],
While[k < m,
If[f[c] == 0, k=m+1; Return[c " is a root"],
If[ Sign[f[c]] ==  Sign[f[a]], 
a = c, b = c; ]; 
c = ((a + b)/2); 
k = k+1; 
output=Append[output,{k,a,c,b,f[c]}]; ]]]]]; 
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"k","(Subscript[a, k])","(Subscript[c, k])","(Subscript[b, k])","f[(Subscript[c, k])]"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",((b-a)/2)]; 
Print["f[c] = ",NumberForm[f[c],16] ]; ]

Working with the function f[x] = x(x+1)(x-1). At a glance we can tell the roots are at x=0, x=1, x=-1. 
a.) If the interval is [-(Sqrt[2]), (π/2)], bisection converges to the root x = 1, and this takes one iteration to know, as after that point the range is [0.07, 1.57], so it can’t converge to 0, only to x = 1.

f[x_] = x (x+1)(x-1);
Plot[f[x], {x,-2,2}]

Bisection[-(Sqrt[2]), (π/2), 10]

b.) All possible cases for intervals that contain multiple roots: Contain all roots(covered in a), contain left two roots, contain right two roots. These cases also must factor in if c(the average of a and b) is larger than both of the roots (or smaller in the symmetrical case).
First check each with an example:

Bisection[-0.5, 1.4,10]

The second case is one that contains the right two roots, and c is not larger than 1. This converges to x = 0 and is guaranteed to do so after the first iteration where the range becomes [-0.5, 0.45]

Bisection[-0.1, 4,10]

In the third case, it once again contains the right two roots, but c is larger than 1, so because of our algorithm the first new range chosen is [1.95, 4] which doesn’t contain any root, and it does not converge to a root. This happens because we check Sign[a] == Sign[c], though if we checked Sign[b] == Sign[c] it would have it’s own ranges that would cause issues.

Bisection[-1.6, 0.5, 10]

In the fourth case, our region contains the left two roots, and c is larger than the leftmost root, after the first iteration it can only converge to x = -1.

Bisection[-4, 0.2, 10]

In the final case, our region contains the left two roots and c is smaller than the rightmost root, after one iteration the region chosen is [-1.9, -0.85] which only contains x = -1 as a root, which it converges to.

Of all these cases, the only one that converged to x = 0 was the two right most roots contained with c < the furthest right root.
This region can be defined as [a,b] where a > -1, b < 3, yet |((a+b)/2)| < 1.
We also have to consider any region where the range has x = 0 as an endpoint, this algorithm has been tuned to stop at that point and label that as a root, so I suppose technically it converges in that case. And in the case that c = 0, the algorithm checks for that and returns 0 as it has found a root

Project 1 Question 3

a.) Verify that each of these functions is a valid fixed point function for f.

i.) f(x) = (e^x) - (1/(x^2)) = 0
(1/(x^2)) = (e^x)
x = (x^3)(e^x) = g1(x)

ii.) (e^x) - (1/(x^2)) = 0
(e^x) = (x^(-2))
ln((e^x)) = ln((x^(-2)))
x = -2 ln(|x|) = g2(x)

iii.) (e^x) - (1/(x^2)) = 0
-(x^2)((e^x) - (1/(x^2))) = 0
-(x^2)(e^x) + 1 = 0
(x^3)(e^x) - (x^2)(e^x) + 3 = (x^3)(e^x) + 2

(((x^3)(e^x)-(x^2)(e^x)+3)/((x^3)(e^x)+2)) = 1

((((x^3) - (x^2))(e^x)+3)/((x^3)(e^x)+2)) = 1

(((x^2)(x - 1)(e^x)+3)/((x^3)(e^x)+2)) = 1

x(((x^2)(x - 1)(e^x)+3)/((x^3)(e^x)+2)) = x

iv.) (e^x) - (1/(x^2)) = 0
(e^x) = (x^(-2))
(((((((e^x)))^((-1)/2))))^(-2)) = (x^(-2))
((((e^x)))^((-1)/2)) = x
(e^((-x)/2)) = x

f[x_]:=Exp[x]-(1/(x^2));
g1[x_]:=(x^3)*Exp[x];
g2[x_]:=-2*Log[Abs[x]];
g3[x_]:=x*(((x^2)(x-1)*Exp[x]+3)/((x^3)*Exp[x]+2));
g4[x_]:=Exp[((-x)/2)];

plot=Plot[{Abs[g1'[x]],Abs[g3'[x]],Abs[g4'[x]]},{x,0.7,0.707},PlotRange->Automatic,PlotLegends->{"|(Subscript[g, 1])'(x)|","|(Subscript[g, 3])'(x)|","|(Subscript[g, 4])'(x)|"}];
rootLine=Graphics[{Red,Dashed,Line[{{0.703467422,0},{0.703467422,5}}]}];
Show[plot,rootLine]

This graph shows the values of |(Subscript[g, 3])'(r)| and |(Subscript[g, 4])'(r)| are less than 1, whereas |(Subscript[g, 1])(x)| is greater than 1. This tells us that (Subscript[g, 1])(x) will not converge to the root, r. The graph of |(Subscript[g, 2])'(x)| doesn’t seem to have a closed-form solution, so we will approximate |(Subscript[g, 2])'(r)|.

m =((g2[0.707]-g2[0.7])/(0.707-0.7));
secantLineg2[x_] = m(x-0.707)+g2[0.707];
plot2=Plot[{g2[x],secantLineg2[x]},{x,.5,.9},PlotRange->Automatic,AspectRatio->Full,PlotLegends->{"(Subscript[g, 2])(x)","Tangent Line Approximation of (Subscript[g, 2])(r)"}];
Show[plot2,rootLine]

We picked two points around the root (0.7 and 0.707) and calculated the slope of the line between them. We then graphed (Subscript[g, 2])(x) and the tangent line approximation. The absolute value of its slope is below:

Abs[m]

Since |(Subscript[g, 2])'(r)| > 1, we know that (Subscript[g, 2])(x) will not converge to the root, r. So now we have two remaining functions which we know will iteratively converge to the root. Let’s explore them visually.

FixedPointSegments[g_,x0_,its_]:=Module[
{l=NestList[g[#]&,x0,its]},
Partition[Flatten[Transpose[{l,l}]],2,1]
];
MakeMovie[g_,p0_,a_,b_,steps_]:=
Module[{p2=Partition[FixedPointSegments[g,p0,steps],2,1]},
movie=Table[Plot[{g[x],x},{x,a,b},PlotStyle->{{},Dashing[{.02}]},
Prolog->MapIndexed[{Hue[#2[[1]]/Length[p2]],Line[#]}&,Take[p2,i]],
AspectRatio->1,
AxesLabel->TraditionalForm/@{(Subscript[x, n]),(Subscript[x, n+1])}
],{i,Length[p2]}];
ListAnimate[movie]
]

MakeMovie[g3,0.25,0,1.5,10]

Here we see a “movie” showing the iterative process in which each new (Subscript[x, n+1]) is found using (Subscript[g, 3])'(x)(starting to the left of r). The values of (Subscript[x, n]) approach our root, r.

MakeMovie[g3,1.2,0,1.5,10]

Here we see a “movie” showing the iterative process in which each new (Subscript[x, n+1]) is found using (Subscript[g, 3])'(x) (starting to the right of r). The values of (Subscript[x, n]) approach our root, r.

MakeMovie[g4,0.25,0,1.5,10]

Here we see a “movie” showing the iterative process in which each new (Subscript[x, n+1]) is found using (Subscript[g, 4])'(x) (starting to the left of r). The values of (Subscript[x, n]) approach our root, r.

MakeMovie[g4,1.2,0,1.5,10]

Here we see a “movie” showing the iterative process in which each new (Subscript[x, n+1]) is found using (Subscript[g, 4])'(x) (starting to the right of r). The values of (Subscript[x, n]) approach our root, r.

FixedPointList[g1,N[0.25],15]
FixedPointList[g2,N[0.25],15]
FixedPointList[g3,N[0.25],15]
FixedPointList[g4,N[0.25],15]

Here we see that the values of (Subscript[x, n]) only converge using (Subscript[g, 3]) and (Subscript[g, 4]). (Subscript[g, 1]) diverges to values around 0, while (Subscript[g, 2]) does not converge as it flops between positive and negative values.

FixedPointList[g1,N[1.25],15]
FixedPointList[g2,N[1.25],15]
FixedPointList[g3,N[1.25],15]
FixedPointList[g4,N[1.25],15]

Here we see that the values of (Subscript[x, n]) only converge using (Subscript[g, 3]) and (Subscript[g, 4]). (Subscript[g, 1]) diverges to extremely large values, while (Subscript[g, 2]) does not converge as it flops between positive and negative values.

b.) The only methods that converge for all values within [0.25,1.25] are (Subscript[g, 3]) and (Subscript[g, 4]). (Subscript[g, 1]) and (Subscript[g, 2]) do not converge because |(Subscript[g, 1])'(r)| < 1 and |(Subscript[g, 2])'(r)| < 1.

c.) The different rates of convergence can be explained by the values of (Subscript[g, 3])'(r) and (Subscript[g, 4])'(r).

Abs[g3'[0.703467422]]
Abs[g4'[0.703467422]]

We see that (Subscript[g, 3])'(r) < (Subscript[g, 4])'(r). Since (Subscript[g, 3])'(r) is really close to 0, we know that it will converge very quickly to the root r compared to (Subscript[g, 4])'(r). Their rates of convergence do meet our theoretical expectations.

Project 1 Question 4

ClearAll[f, a0, epsilon, eps, nits, m, b0, x0,x1, max, list, fx, a3,a2,a1,maxits,real,plot]

f[x_] := (E^x)-(1/(x^2))
Plot[f[x],{x,0,2}]

epsilon = 0.000001;

NewtonsMethod[f_,a0_,eps_,nits_] :=
Module[
{a=N[a0], (* Use a local variable (so as not to alter passed arguments) *)
fa = f[N[a0]],
i=0,
output (* Let's collect the output, since we've gone to the trouble of computing it. *)
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<nits,
a = a-(fa/(f'[a]));
fa=f[a];
i=i+1;
output=Append[output,{i,a,fa}];  (* once again, collect output. *)
]; 
(* Now let's augment the naive function to include printing some nice formatted output: *)
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"i","(Subscript[x, k])","f[(Subscript[x, k])]"}}],16]]; 
Print["  x  = ",NumberForm[a,16] ]; 
Print["f(x) = ",NumberForm[fa,16] ]; 
Return[Transpose[output]]; (* Now I'll actually return the output, so that it can be plotted, etc. *)
]

NewtonsQuadratic[f_,a0_,eps_,nits_] :=
Module[
{a=N[a0], (* Use a local variable (so as not to alter passed arguments) *)
fa = f[N[a0]],
fp = f'[N[a0]],
fpp = f''[N[a0]],
i=0,
output (* Let's collect the output, since we've gone to the trouble of computing it. *)
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<nits,
a = a-((fp-Sign[fp]*Sqrt[fp*fp-2*fa*fpp])/fpp); (* a variation on Joey's idea...*)
fa=f[a];
fp = f'[a];
fpp = f''[a];
i=i+1;
output=Append[output,{i,a,fa}];  (* once again, collect output. *)
]; 
(* Now let's augment the naive function to include printing some nice formatted output: *)
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"i","(Subscript[x, k])","f[(Subscript[x, k])]"}}],16]]; 
Print["  x  = ",NumberForm[a,16] ]; 
Print["f(x) = ",NumberForm[fa,16] ]; 
Return[Transpose[output]]; (* Now I'll actually return the output, so that it can be plotted, etc. *)
]

Bisection[a0_,b0_,m_,eps_] :=
Module[{},
a = N[a0];  
b = N[b0];  
c = ((a + b)/2); 
k = 0; 
output={{k,a,c,b,f[c]}}; 
If[f[a] == 0, k=m+1; Return[a" is a root"], If[f[b]==0, k = m+1; Return[b" is a root"],If[f[c] == 0, k=m+1; Return[c" is a root"],
While[k < m && Abs[f[c]] > eps,
If[f[c] == 0, k=m+1; Return[c " is a root"],
If[ Sign[f[c]] ==  Sign[f[a]], 
a = c, b = c; ]; 
c = ((a + b)/2); 
k = k+1; 
output=Append[output,{k,a,c,b,f[c]}]; ]]]]]; 
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"k","(Subscript[a, k])","(Subscript[c, k])","(Subscript[b, k])","f[(Subscript[c, k])]"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",((b-a)/2)]; 
Print["f[c] = ",NumberForm[f[c],16] ]; ]

secantMethod[x0_,x1_,max_,eps_]:= 
Module[{},
k = 1;
p0 = N[x0];
p1 = N[x1];
p2 = p1;
p1 = p0;
While[k < max && Abs[f[p2]]>eps,
p0 = p1;
p1 = p2;
p2 = p1 - ((f[p1](p1 - p0))/(f[p1]-f[p0]));
k = k+1;];
Print[" p = ", NumberForm[p2,16]];
Print[" ∆p = ±", Abs[p2-p1]];
Print["f[p] = ", NumberForm[f[p2],16]];
Print["Number of iterations = " NumberForm[k,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_,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]
]
Print["Number of iterations = " NumberForm[i,2]];
]

NewtonsMethod[f,1,epsilon,20]

NewtonsQuadratic[f,1, epsilon, 20]

Bisection[0.5, 1, 20,epsilon]

secantMethod[0.6358246728512563, 0.6963253420149134, 10, epsilon]

mullers[f, 0.6358246728512563, 0.6963253420149134,0.7033929176442707, 20, epsilon, True,True]

Comparing rates of convergence.

Newton’s Method found the root in 4 iterations 
Newton’s Quadratic found the root in 3 iterations
Bisection found the root in 18 iterations
Secant found the root in 4 iterations
Muller’s found the root in a single iteration

Muller’s method found the root the fastest by far, due to it’s accurate starting values from Newton’s, with Bisection taking the slowest.
All intervals/starting values were 1, [0.5,1] or the first few iterations of Newton’s in the case of Secant/Muller’s, all methods stopped if they got within 1 x (10^(-6)) of 0.