MAT 360 Numerical Analysis Project 1

Solutions from students in MAT360, Spring, 2017

Problem 1: (Alyssa Farmer, Katie Schwegman, and Mike Wilson)

The following method (NewtonsMethod) is credited to Andy Long:

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. *)
]

f[x_] = x*Exp[x]

Plot[f[x],  {x, -3, .5}, PlotRange->{-.5, 1}]

A. Using four different starting points for Newton’s method with qualitatively different behavior for the iteration (labeling each on the graph), demonstrate the following behavior numerically:
a. A point where Newton’s method converges monotonically to the root. x0 = 8

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

This bounces around (not monotonic)

x0 = 8;
NewtonsMethod[f, x0, eps, nits];

x0 = -2 will converge monotonically to ~0

b. A point where Newton’s method will converge non-monotonically. x0 = -0.3

x0 = -.3;
NewtonsMethod[f, x0, eps, nits];

x0 = -0.5 will converge non-monotonically (it bounces around the root a little before converging)

c. A point where Newton’s method will “blow up”. x0 = -1

Solve[f'[x] == 0]

x0 = -1;
NewtonsMethod[f, x0, eps, nits];

x0 = -1 will blow up

d. A point where Newton’s method will monotonically flee the root. x0 = -2.5

x0 = -2.5;
NewtonsMethod[f, x0, eps, nits];

x0 = -2.5 will monotonically flee the root

B. Characterize completely the “basins” [a,b] along the real number line where each of the four behaviors above will occur.
a. Monotonically converge x ∈ (0, ∞)

x0 = 0.01;
NewtonsMethod[f, x0, eps, nits];

x0 = 0.1;
NewtonsMethod[f, x0, eps, nits];

x0 = -0.01;
NewtonsMethod[f, x0, eps, nits];

x0 = 2;
NewtonsMethod[f, x0, eps, nits];

x ∈ (0, ∞)

b. Converge non-monotically x ∈ (-1, 0)

x0 = -.999;
NewtonsMethod[f, x0, eps, nits];

Flees root initially but slowly comes back to root monotonically

x0 = -0.001;
NewtonsMethod[f, x0, eps, nits];

x ∈ (-1, 0)

c. “Blow up” x = -1

x0 = -1;
NewtonsMethod[f, x0, eps, nits];

x0 = -.999;
NewtonsMethod[f, x0, eps, nits];

x0 = -1.001;
NewtonsMethod[f, x0, eps, nits];

x0 = -1.01;
NewtonsMethod[f, x0, eps, nits];

x = -1

d. Monotonically flee root x ∈ (-∞, -1)

x0 = -1.01;
NewtonsMethod[f, x0, eps, nits];

x0 = -1.001;
NewtonsMethod[f, x0, eps, nits];

x0 = -.999;
NewtonsMethod[f, x0, eps, nits];

This flees originally but then slowly comes back to the root.

x0 = -.9999;
NewtonsMethod[f, x0, eps, nits];

x0 = -2;
NewtonsMethod[f, x0, eps, nits];

x0 = -1.5;
NewtonsMethod[f, x0, eps, nits];

x ∈ (-∞, -1)

(Madison, Trey, Courtney):

“Basins” 
a) For x0 > 0, Newton’s Method will converge monotonically to the root x=0. All roots of tangent lines drawn in the first quadrant will be posiitve. This means the method will continually converge to x=0. 

b) For -1 < x0 < 0, Newton’s Method will converge non-monotonically. Any tangent lines drawn from this area will cross the y-axis and intersect the x-axis positively. The method will then converge monotonically from there. 

c) For x0 = -1, Newton’s method will “blow up” because f’(-1) = 0, so the tangent line at that point will never cross the x-axis. 
d) For x0 < -1, Newton’s Method will monotonically flee the root. In that area, the roots of the tangent lines will approach -∞.

Problem 2 (Zach, Joey, Brooke):

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 its 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

(Leah, Kara, Lauren):

Clear[numberIts, maximumIterations, f, x]

numberIts = 0;
maximumIterations = 100;
verbose=0;

bisection[aa_, bb_, f_, maxIts_] := 
Module[{a=Min[{aa,bb}],b=Max[{aa,bb}],fa,fb,c,
numIts = 0},
fa  = f[a];
fb =  f[b];
c = a + 0.5 * (b-a);
While[fa*fb<0 && a < c <b && numIts < maxIts,
numberIts  = numberIts + 1;
numIts = numIts + 1;
fc = f[c];
If[Sign[fc] == Sign[fa],
{a, fa} = {c, fc}, 
{b, fb} = {c, fc}];
c= a + 0.5 * (b -a);];
If[numIts>0,Return[c]]
]

f[x_] := (x+1)x(x-1);

a.

bisection[-Sqrt[2], Pi/2, f, maximumIterations]
bisection[Sqrt[2],-Pi/2, f, maximumIterations]
bisection[-Sqrt[2], -Pi/2, f, maximumIterations]

Print[numberIts]

It converges to the root x = 1. It took 47 iterations.

b.

Plot[f[x], {x, -2, 2}]

f(x) is an odd function since f(-x) = (-x + 1)(-x)(-x- 1) = (-x + 1) (x) (x + 1) = 
-(x - 1) x (x+1)= -f(x). Therefore, we can look simply at positive x, and the results will follow for negative x.

Plot3D[bisection[a, b, f, maximumIterations], {a, -4, -1}, {b, -a,4}]

Based on the graph, it appears that bisection only converges to the root x = 0 when both a is between -1 and 0 and b is between 0 and approximately 3.

Specifically, let’s look at values of a within [-1, 0) and b within (1, 3]. Therefore, [a,b] contains both the roots 0 and 1, yet converges to the root x = 0 (for at least some of the values...)

Plot3D[bisection[a, b, f, maximumIterations], {a, -1, 0}, {b, 1,3}]

As noted in the book, bisection always works, as long as you can find a sign change and the function is continuous. Since the function is continuous, the only places where it won’t work are where there isn’t a sign change.

On the first iteration, (Subscript[c, 0]) will be chosen as the midpoint of (Subscript[a, 0]) and (Subscript[b, 0]). The way this algorithm is written, the sign of f((Subscript[c, 0])) is checked first against f((Subscript[a, 0])), so if the signs between those are equal, then (Subscript[c, 0]) will replace (Subscript[a, 0]). Otherwise, (Subscript[c, 0]) will replace (Subscript[b, 0]). If (Subscript[c, 0]) is within (0, 1], then (Subscript[c, 0]) replaces (Subscript[b, 0]), and the new interval [(Subscript[a, 1]), (Subscript[b, 1])] not only encloses the single root x = 0, but also exhibits a sign change in f values. Therefore it will converge to x = 0. So we know convergence to x = 0 for 0.5(a+b) within [0, 1].

0<0.5(a+b) < 1;

-b < a < 2 - b;

RegionPlot[-b < a && a < 2 - b, {a, -1, 0}, {b, 1, 3}]

This shows the area of values that converge to x = 0, which appears similar to what is shown in the 3d plot.

(notice, it could also converge to 0, if (Subscript[c, 0]) had the same sign as (Subscript[a, 0]) and (Subscript[c, 0]) was within [-1, 0], as long as (Subscript[c, i]) was eventually chosen in the interval [0, 1]. however, I only examine what happens on the first iteration with (Subscript[c, 0]) because for the intervals of a and b we’re examining, there are only two options: (Subscript[c, 0 ])within [0, 1), or (Subscript[c, 0]) within [1, 3])

Problem 3 (Madison, Trey, Courtney)

t[x_]:=(E^x)-(1/(x^2))
r=x/.FindRoot[t[x],{x,1}]

Here are our fixed point (at least proposed fixed point) functions:

g1[x_]:=(x^3)(E^x)
g2[x_]:=If[x<0,-2Log[-x],-2 Log[x]]
g3[x_]:=x*((((x^2)(x-1)(E^x)+3)/((x^3)(E^x)+2)))
g4[x_]:=(E^((-x)/2))

Plot[{t[x],g1[x],g2[x],g3[x],g4[x]},{x,.25,1.5},PlotStyle->{Red,Blue,Green,Orange,Purple}]

G1

p0=1.1;
a=.25;
b=1.5;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g1[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

Did not converge.

G2

p0=.7;
a=.25;
b=1.5;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g2[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,40}]]

Did not converge.

G3

p0=1.1;
a=.25;
b=1.5;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g3[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

Converged!

G4

p0=1.1;
a=.25;
b=1.5;
precision=50;
Clear[p,k]
p[0]=N[p0,precision];
Do[p[k]=N[g4[p[k-1]],precision],{k,1,100}]
(* Let’s write out the values: *)
TableForm[Table[{k,NumberForm[p[k],precision]},{k,0,11}]]

Converged!

This function will look at the points a given starting value will converge to:

limGFunc[g_,x_]:=
Module[{lim=g[x]},
Do [lim=g[lim],{i,30}];
lim
]

Plot[limGFunc[g1,x],{x,.25,1.5}]
Plot[limGFunc[g2,x],{x,.25,1.5}]
Plot[limGFunc[g3,x],{x,.25,1.5}]
Plot[limGFunc[g4,x],{x,.25,1.5}]

If the absolute value of |g’[r] |< 1, then the fixed point method will converge for all values in the neighborhood. For this reason, g1 and g2 are not valid fixed point methods for f[x]. The smaller the magnitude of slope the faster the fixed point method converges:

N[Abs[{g1'[r],g2'[r],g3'[r],g4'[r]}]]

Problem 3 (Alyssa Farmer, Katie Schwegman, and Mike Wilson)

The following code is credited to Andy Long:

r=0.703467422 (* The root, approximately *)

f[x_]:=Exp[x]-(1/(x^2))(* The function we are trying to find the root of, with the four fixed point functions below *)
g1[x_]:=x^3*Exp[x]
g2[x_]:=-2Log[x]
g3[x_]:=x*((x^2*(x-1)*Exp[x]+3)/(x^3*Exp[x]+2))
g4[x_]:=Exp[((-x)/2)]

Iterate[x_]:=N[{g1[x[[1]]],g2[x[[2]]],g3[x[[3]]],g4[x[[4]]]}] (* Fixed point iteration for each of the four functions *)

vals={.25,.25,.25,.25}; (* Iteration with a starting value of 0.25 *)
i=0
Its=10
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

vals={.5,.5,.5,.5}; (* Iteration with a starting value of 0.5 *)
i=0;
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

vals={.75,.75,.75,.75}; (* Iteration with a starting value of 0.75 *)
i=0;
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

vals={1,1,1,1}; (* Iteration with a starting value of 1 *)
i=0;
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

vals={1.25,1.25,1.25,1.25}; (* Iteration with a starting value of 1.25 *)
i=0;
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

vals={1.5,1.5,1.5,1.5}; (* Iteration with a starting value of 1.5 *)
i=0;
While[i<Its,i=i+1;vals=Iterate[vals];Print[vals]]

A. Verify that each is a valid fixed point function for f.

Solve[g1[x]==x,x]//N

Solve[g2[x]==x,x]//N

Solve[g3[x]==x,x]//N

Solve[g4[x]==x,x]//N

Since each of these functions satisfy g(r) = r, this means that they are valid fixed point functions. It is also true that each of these functions only have one fixed point within the interval, and this fixed point is the root of f(x).

B. Do all methods converge for all values in the neighborhood?

The first function g1(x) and the second function g2(x) do not seem to converge at all, as can be seen above when tested on various values.

The third function g3(x) and the fourth function g4(x) do seem to converge for all values in the interval, as can be seen above again. Also, since these two functions converge for the values 0.25 and 1.5, it makes sense that they would converge for all values in between (in a neighborhood near r), since |g3’(r)| is essentially zero, which is less than 1, and |g4’(r)| = 0.351734, also less than 1 (calculations in Part C).

C. How do you explain their different rates of convergence (for those which do converge)? Do their rates of convergence meet theoretical expectations?

Plot[{g1[x],g2[x],g3[x],g4[x],x,f[x]},{x,0,4}]

{g1'[r],g2'[r],g3'[r],g4'[r]} (* Derivatives of each function at the root *)

The third and fourth functions both have derivatives close to zero at the root. The third function should converge the fastest, since its derivative at r looks to be zero or essentially zero at the root. The fourth one should be the next fastest. This is true in all cases that were tried above, so the expectation of the speed of convergence for these functions was met. In the case of the first two functions, the absolute values of their derivatives at r were larger than 1, so they did not converge at all. In order for convergence in a neighborhood to work, the condition of |g’(r)| < 1 must be satisfied. This is true for the second two functions, but not for the first two.

Problem 4

Clear[f]

f[x_]:=(E^x)-(1/(x^2))
Plot[f[x],{x,.25,1.5}]
r=x/.FindRoot[f[x],{x,1}]

Bisection Method

Bisection[f_,a0_,b0_,its_] :=
Module[{
a=N[a0], 
b=N[b0],
fa = f[a0],  
fb = f[b0],c,fc,i=1},  
c = a+0.5(b-a); 
fc=f[c]; 
While[a < c && c<b && i<=its,
If[ Sign[fa] ==  Sign[fc], 
a = c;fa=fc, 
b = c;fb=fc
];
i=i+1;
c = a+0.5(b-a);
fc=f[c];
]; 
Return[c]
]  
Bisection[f,.25,1.5,100]

Bisection method takes 17 iterations to converge with an error of (10^(-6))

Newton’s Method

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
},
While[Abs[fa] > eps && i<nits,
a = a-(fa/(f'[a]));
fa=f[a];
i=i+1;
]; 
Return[a]
]  
NewtonsMethod[f,1,(10^(-6)),50]

Newton’s Method take 4 iterations to converge with an error of (10^(-6)).

Muller’s Method

MullersMethod[f_,x0_,x1_,x2_,ϵ_,nits_]:=Module[{x00=N[x0],x11=N[x1],x22=N[x2],fx0,fx1,c,a,b,i=0},
fx0=f[x00];
fx1=f[x11];
c=f[x22];
a=((((x11-x22)*(fx0-c))-((x00-x22)*(fx1-c)))/((x00-x22)*(x11-x22)*(x00-x11)));
b=((((((x00-x22))^2)*(fx1-c))-((((x11-x22))^2)*(fx0-c)))/((x00-x22)*(x11-x22)*(x00-x11)));

While[Abs[c]>ϵ&&i<nits,
x00=x11;
x11=x22;
fx0=f[x00];
fx1=f[x11];
x22=x22+((-2c)/(b+Sign[b]*Sqrt[(b^2)-4a*c]));
c=f[x22];

a=((((x11-x22)*(fx0-c))-((x00-x22)*(fx1-c)))/((x00-x22)*(x11-x22)*(x00-x11)));
b=((((((x00-x22)))^2)*(fx1-c)-((((x11-x22))^2)*(fx0-c)))/((x00-x22)*(x11-x22)*(x00-x11)));
i=i+1;
];
Return[x22]
]
MullersMethod[f,.5,.75,1,(10^(-6)),50]

Muller’s Method takes 4 iterations with an error of (10^(-6)).

Newton’s Quadratic

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
},
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;
]; 
Return[a]
]  
NewtonsQuadratic[f,1,(10^(-6)),50]

Newton’s Quadratic takes 3 iterations to converge with an error of (10^(-6)).

Secant Method

SecantMethod[f_,a0_,b0_,eps_,nits_] :=
Module[
{a=N[a0],b=N[b0], (* Use a local variable (so as not to alter passed arguments) *)
fa = f[N[a0]],fb=f[N[b0]],
i=0
},
While[Abs[fa] > eps && i<nits,
a = a-((fa*(b-a))/(fb-fa));
fa=f[a];
fb=f[b];
i=i+1;
]; 
Return[a]
]  
SecantMethod[f,.5,1,10^-6,30]

Secant method takes 14 iterations to converge with an error of (10^(-6 )).

Method Comparison Fastest to Slowest: Chosen error (10^(-6)).

1. Newton’s Quadratic was the fastest and converged in only 3 iterations. 
2. Muller’s and Newton’s method both converged in 4 iterations.
3. The Secant method converged in 14 iterations.
4. The Bisection method was the slowest and converged in 17 iterations.

Method Comparison on convergence over the interval (.25,1.5):

eps=2*10^(-7);
frame={{.25,1.5},{r-eps,r+eps}};

handicap=.5;
p1=Plot[NewtonsQuadratic[f,x,(10^(-6)),50],{x,.25,1.5},PlotRange->frame]
p2=Plot[NewtonsMethod[f,x,(10^(-6)),50],{x,.25,1.5},PlotRange->frame]
p3=Plot[SecantMethod[f,x,NewtonsMethod[f,x,(10^(-6)),handicap],(10^(-6)),50],{x,.25,1.5},PlotRange->frame]
p4=Plot[MullersMethod[f,x,NewtonsMethod[f,x,(10^(-6)),handicap],NewtonsMethod[f,NewtonsMethod[f,x,(10^(-6)),handicap],(10^(-6)),1],(10^(-6)),50],{x,.25,1.5},PlotRange->frame]

Let’s take a look at the iterates of all methods from a given start (or so):

maxits=10;

createIterates[x_]:=
{
iteratesBisection=Table[Bisection[f,.25,x,i],{i,maxits}];
iteratesNewton=Table[NewtonsMethod[f,x,(10^(-6)),i],{i,maxits}];
iteratesNewtonsQuadratic=Table[NewtonsQuadratic[f,x,(10^(-6)),i],{i,maxits}];
iteratesSecantMethod=Table[SecantMethod[f,x,NewtonsMethod[f,x,(10^(-6)),handicap],(10^(-6)),i],{i,maxits}];
iteratesMullersMethod=Table[MullersMethod[f,x,NewtonsMethod[f,x,(10^(-6)),handicap],NewtonsMethod[f,NewtonsMethod[f,x,(10^(-6)),handicap],(10^(-6)),1],(10^(-6)),i],{i,maxits}];
}
createIterates[1.1]

ListPlot[iteratesBisection]
ListPlot[iteratesNewton]
ListPlot[iteratesNewtonsQuadratic]
ListPlot[iteratesSecantMethod]
ListPlot[iteratesMullersMethod]

Let’s take a look regressions on the errors of the iterates of all methods:

errorRegress[iterates_]:=
Module[{},
errors=iterates-r;
ekplus1=Log[Abs[Table[errors[[i]],{i,2,maxits}]]];
ek=Log[Abs[Table[errors[[i]],{i,1,maxits-1}]]];
data=Transpose[{ekplus1,ek}];
lm=LinearModelFit[data,x,x]
]

errorRegress[iteratesBisection]
errorRegress[iteratesNewton]
errorRegress[iteratesNewtonsQuadratic]
errorRegress[iteratesSecantMethod]
errorRegress[iteratesMullersMethod]