MAT 360 Numerical Analysis Project 1

Alyssa Farmer, Katie Schwegman, and Mike Wilson

Problem 1

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)

Problem 2

f[x_]:=(x+1)*x*(x-1) (* The function we are considering in this problem *)

Plot[f[x],{x,-Sqrt[2],Pi/2}] (* A plot just to get a visual of the function we are working with.  A quick inspection shows the roots to be -1, 0, 1. *)

A. To which root will bisection converge, if we start with the interval [-√2,π/2]?  How many iterations must be carried out to know that?

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

Bisection[f_,a0_,b0_] :=
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]; 
output={{a,c,b,fc,i}}; 
While[a < c && c<b && i<10,
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];
output=Append[output,{a,c,b,fc,i}];  
]; 

Print[NumberForm[TableForm[output,
TableHeadings->{None,{"(Subscript[a, k])","(Subscript[c, k])","(Subscript[b, k])","f[(Subscript[c, k])]","iteration"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",0.5*(b-a)]; 
Print["f[c] = ",NumberForm[fc,16] ]; 
Return[Transpose[output]]; 
]

Bisection[f,-Sqrt[2],Pi/2];(* Time to see what results we get when we use bisection on the specified interval. *)

The 8th iteration gives a value to c of 0.999446777696336.  When looking at the previous 7 iterations, it becomes apparent that the root is indeed converging to 1.000.

B. Characterize intervals (consider only those intervals [a,b] such that a, b ≤ 4) which
	a. contain multiple roots, yet,
	b. which result in the root x=0 being chosen by bisection.

Bisection[f,-3,.232]; (*This bisection is used to test intervals from the below plots*)

g[a_,b_]:=a+.5(b-a)
f[x_]:=(x+1)*x*(x-1)
Manipulate[Plot[{(*g[a,b],f[a],b,f[b],*)f[a]*f[b],f[g[a,b]]},{a,-4,0}],{b,0,4}]

Evaluating the above plot.
In order for an interval to contain multiple roots, it is important to narrow down the base intervals we are working with.  If a is less than -1, b must be greater than 0.  If a is between -1 and 0, b must be larger than 1.  This guarantees the interval [a,b] contains multiple roots.  Now some elimination can be done using the above manipulate.  We have chosen to look only at a<0 and b>0 since these conditions have to be met to contain multiple roots.

g[a_,b_]:=a+.5(b-a)
f[x_]:=(x+1)*x*(x-1)
Manipulate[Plot[{(*g[a,b],f[a],b,f[b],*)f[a]*f[b],f[g[a,b]]},{a,-4,-1}],{b,0,1}]

When 0<b<1 and -4<a<-1
When b is between the roots of 0 and 1, a must be less than -1 to ensure 2 roots are contained in the interval.  Manipulating the above plot, it is apparent that f(a)*f(b) is never less than 0.  A sign change never occurs and therefor there is no interval for 0<b<1 where 0 will be the selected root.

g[a_,b_]:=a+.5(b-a)
f[x_]:=(x+1)*x*(x-1)
Manipulate[Plot[{(*g[a,b]f[a],b,f[b],*)f[a]*f[b],f[g[a,b]]},{a,-1,0}],{b,1,4}]

When 1<b<4 and -1<a<0.
Manipulating this plot, the only stipulation on 0 being the selected root is that f(c)<0. Any selection of a can give this result anytime b<2.  So the entire interval [(-1,0),(1,2]] works here and gives the root of 0. When 2<b<3, you must be clever about your selection of a to give an appropriate f(c)<0.  When b>3, the selected root will never be 0.

g[a_,b_]:=a+.5(b-a)
f[x_]:=(x+1)*x*(x-1)
Manipulate[Plot[{(*g[a,b],f[a],b,f[b],*)f[a]*f[b],f[g[a,b]]},{a,-4,-1}],{b,1,4}]

When 1<b<4 and -4<a<-1
Basically, when all 3 roots are in the selected interval.  In this interval, 0 is never selected as the root using bisection.

Using bisection, the intervals [a,b] between -4 and 4 that converge to 0 as the root are:
[(-1,0),(1,2]]
[(-1,0),(2,3)] such that f(c)<0

Problem 3

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

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

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

The following method (SecantMethod) is an altered version of the NewtonsMethod method, which is credited to Andy Long:

SecantMethod[f_,a0_,b0_,eps_,nits_] :=
Module[
{a=N[a0],
b=N[b0], 
fa = f[N[a0]],
fb =f[N[b0]],
i=0,
output 
},
output={{i,a,fa}}; 
While[Abs[fa] > eps && i<nits,
a = a-(fa/((fa-fb)/(a-b)));
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]]; 
]

The following methods (dividedDifference and mullers) are credited to Andy Long:

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

f[x_]:= Exp[x] - 1/x^2 (* The function of interest *)

eps=10^(-6);
nits=100;

Plot[f[x], {x, 0, 2}, PlotRange->{-10, 10}]

x0 = 1;

NewtonsMethod[f, x0, eps, nits];

NewtonsQuadratic[f, x0, eps, nits];

Bisection[f,.6,1];

x0=.5;
x1=1;
eps=10^(-6);
nits=100;
SecantMethod[f,x0,x1,eps,nits];

Plot[f[x],{x,.6,1}]
N[Solve[f[x]==0,x]]
mullers[f,1,.6,.7,100,10^-6,True,True]

Newton’s Method: 4 interations
Newton’s Quadratic: 3 iterations
Bisection: 19 iterations
Secant Method: 14 iterations
Muller’s Method: 3 iterations
So, the conclusion can be drawn that the rates of convergence of the five methods go in the following order, from fastest to slowest:
Newton’s Quadratic and Muller’s Method, Newton’s Method, Secant Method, and Bisection