Project 1 (Leah Leist, Kara Starnes, Lauren Hensley)

#1 (Leah Leist)

A. Consider the following graph:

Plot[x*E^x,{x,-3,.5},PlotRange->Full]

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]]; ]

a. A point where Newton' s method converges monotonically to the root.
I use the point .2 here at it will monotonically converge to the root.

f[x_]=x*E^x;
Show[Plot[f[x],{x,-3,1}],Graphics[{PointSize[.03],Point[{.2,.244}]}]]
x0=.2;
eps=10^(-6);
nits=100;
NewtonsMethod[f,x0,eps,nits];

b.A point where Newton' s method will converge non - monotonically;
Here I used

f[x_]=x*E^x;
Show[Plot[f[x],{x,-3,1}],Graphics[{PointSize[.03],Point[{-.2,-.16}]}]]
x0=-.2;
eps=10^(-6);
nits=100;
NewtonsMethod[f,x0,eps,nits];

c.A point where Newton’s method will “blow up”;
A Point where Newtons method will blow up is where the slope is zero since we cannot divide by 0. This is at x=-1.

f[x_]=x*E^x;
Show[Plot[f[x],{x,-3,1}],Graphics[{PointSize[.03],Point[{-1,-.36}]}]]
x0=-1;
eps=10^(-6);
nits=100;
NewtonsMethod[f,x0,eps,nits];

d. A point where Newton’s method will monotonically flee the root.
Here I used the point -2 Newtons method will flee the root when the slope is negative which is any where x > -1.

f[x_]=x*E^x;
Show[Plot[f[x],{x,-3,1}],Graphics[{PointSize[.03],Point[{-2,-.27}]}]]
x0=-2;
eps=10^(-6);
nits=100;
NewtonsMethod[f,x0,eps,nits];

B.Characterize completely the "basins"[a, b] along the real number line where each of the
four behaviors above will occur.For example, [a, b] such that x ∈ [a, b] implies monotonic
convergence.

a. converges monotonically to the root [0,∞]
b.converge non - monotonically [-1,0]
c.“blow up” [-1]
d.monotonically flee the root [-∞,-1]

#2 Lauren Hensley

Clear[numberIts, maximumIterations, f, x]

numberIts = 0;
maximumIterations = 100;

bisection[aa_, bb_, f_, maxIts_] := 
Module[{a=aa,b=bb,fa  = f[aa],
fb =  f[bb],
c = aa + 0.5 * (bb-aa),
numIts = 0},
While[a < c && 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);];
Return[c]
]

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

a.

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, -2, 4}, {b, -2,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])

#3 Kara Starnes

Comparison & Rank

Clear[x, f]

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)]
g5[x_]:=x-((f[x])/(f'[x]))

m=f'[0.703467422];

Iterate[x_]:=N[{g1[x[[1]]],g2[x[[2]]],g3[x[[3]]],g4[x[[4]]],g5[x[[5]]]}]

vals={0.9,0.9,0.9,0.9,0.9};
i=0
Its=10
While[i<Its,i=i+1;
vals=Iterate[vals];
Print[vals]]

Rank (by convergence--fastest to slowest):
	g5(b) {3 iterations}
	g3(b) {3 iterations}
	g4(b) {4 iterations}
Does not converge:
	g1(b)
	g2(b)

3a.

x=g1(x)

(0.703467422^3)*Exp[0.703467422]

x=g2(x)

-2*Log[Abs[0.703467422]]

x=g3(x)

0.703467422*(((0.703467422^2)*(0.703467422-1)*Exp[0.703467422]+3)/((0.703467422^3)*Exp[0.703467422]+2))

x=g4(x)

Exp[((-0.703467422)/2)]

3b.

Clear[x]

Plot[{g1[x],g2[x],g3[x],g4[x],g5[x],x,f[x]},{x,0.25,1.5},PlotLegends->Automatic]

{g1'[0.703467422],g2'[0.703467422],g3'[0.703467422],g4'[0.703467422],g5'[0.703467422]}

g3(x), g4(x), and g5(x) converge.

g3(x): Decreases on from x=0 to approximately x=0.906164, where the graph starts increasing. Thus, g3(x) converges on [0.906164, 1.5].

Plot[g3'[x],{x,.25,1.5}]

g4(x): Converges on the entire neighborhood interval [0.25,1.5].

Plot[g4'[x],{x,.25,1.5}]

g5(x): Decreases on from x=0 to approximately x=0.906164, where the graph starts increasing. Thus, g3(x) converges on [0.906164, 1.5].

Plot[g5'[x],{x,.25,1.5}]

3c.

For those which do converge, their rates of convergence can be explained by the behavior of each graph. For example, g3(x) and g5(x) both have the possibility of divergence. Hence, it will not converge as fast. On the other hand, g4(x) doesn’t have the possibility of divergence, so it simply converges and it does so quickly.

#4 All.

Clear["Global '*"]

f[x_]:=Exp[x]-(1/(x^2));

epsilon = 0.00000001;
maxInts = 50;

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[i]
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]]; ]

bisection[aa_, bb_, f_, maxIts_] := 
Module[{a=aa,b=bb,fa  = f[aa],
fb =  f[bb],
c = aa + 0.5 * (bb-aa),
numIts = 0},
While[a < c && 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);];
Print[numIts];
Return[c]
]

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[i];
]

NewtonsQuadratic[f_,a0_,eps_,nits_] :=
Module[
{a=N[a0], 
fa = f[N[a0]],
fp = f'[N[a0]],
fpp = f''[N[a0]],
i=0,
output 
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<nits,
a = a-((fp-Sign[fp]*Sqrt[fp*fp-2*fa*fpp])/fpp); 
fa=f[a];
fp = f'[a];
fpp = f''[a];
i=i+1;
output=Append[output,{i,a,fa}];  
]; 
Print[i];
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]]; 
]

SecantMethodList[f_,{x_,x0_,x1_},n_]:=(NestList[Last[#]-{0,(Function[x,f][Last[#]]*Subtract@@#)/Subtract@@Function[x,f]/@#}&,{x0,x1},n]; Print[n])

NewtonsMethod[f, 0.1, epsilon, maxInts];

bisection[0.1, 1, f, maxInts]

mullers[f,0.1,0.14942,0.22205266,maxInts,epsilon, True,True]

NewtonsQuadratic[f, 0.1, epsilon, maxInts]

SecantMethod[f_,x00_, x11_, max_, eps_] :=
Module[{x0 = x00, x1 = x11, x2, k = 1, p2, p1, p0}, 
p2 = N[x1];
p1 = N[x0];
While[k < max && Abs[f[p1]] > eps, 
p0 = p1;
p1 = p2;
p2 = p1 - (f[p1](p1 - p0))/(f[p1] - f[p0]);
Print[p2];
k = k + 1;]
Print[p2];
Print[k]];

SecantMethod[f, 0.1,0.14942, maxInts, epsilon]

Newton’s came in first place with 9 iterations to reach the root r = 0.7...  Secant had the second fastest rate of convergence with only 13 iterations to reach the root. Muller’s was extremely close with only 14 iterations. Newton’s quadratic, however, didn’t even reach the root within 50 iterations, and Bisection took 46 iterations to reach the root.

Based on this specific case of starting value of 0.1 (and other starting values based on Newton’s method), this gives us rates of convergence in increasing order (slowest to fastest):

-Newton’s quadratic
	-Bisection
	-Muller’s
	-Secant
	-Newton’s