Homework 2.3: Tea time

Problems, pp. 67--, #3, 10, 13, 15, 16, 24.

Exercise #3, p. 67

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Dr. Len Brin         15 January 2014 %
% Modified by Andy Long           6 February 2024 %
% Purpose: Implementation of Steffensen’s method  %
% INPUT: function f; initial value x0; tolerance  %
%        TOL; maximum iterations N0               %
% OUTPUT: approximation x and number of           %
%        iterations i; or message of failure      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,i] = steffensen(f,x0,TOL,N0)
  i=1;
  A=x0;
  B=f(A);
  while (i<=N0)
    x0=A;
    x1=B;
    x2=f(x1);
    if (abs(x2-x1)<TOL)
      x=x2;
      return
    end%if
    x=x0-(x1-x0)^2/(x2-2*x1+x0);
    if (abs(x-x2)<TOL)
      return
    end%if
    A=x;
    B=f(A);
    i=i+1;
  end%while
  x=”Method failed---maximum number of iterations reached”;
end%function

This Mathematica version of Steffensen’s not only computes using Steffensen’s acceleration, but it also compares the result to 
	* the standard FPI using the given function, and 
	* the Aitken’s acceleration sequence that goes along with it.

steffensen[g_,q0_,tol_,n_]:=Module[{i=0,a,x1,x2,x0,y1,y2,y0,denom},
(* the x sequence will be Steffensen's *)
x0=q0; (* Set x0 to the initial guess, and compute two iterates:  *)
x1=N[g[x0]];
x2=g[x1];
y0=x0;  (*The y's will be regular fpi: *)
y1=x1;
y2=x2;
(* We've constructed the first two iterates from q0 *)
denom=(x2-2*x1+x0);
Print[
"We produce comparable results after each two iterations of the FPI:"
];
Print[{"FPI","Aitken's","Steffensen's"}];
Print[{y0,"NA","NA"}];
Print[{y1,"NA","NA"}];
While[Abs[denom]>tol && i<n,
(* reset initial starting values for the sequences: *)
x0=x0-(x1-x0)^2/denom;
a=y2-(y1-y2)^2/(y2-2*y1+y0);
y0=y2;
Print[{y0,a,x0}];
y1=g[y0];
y2=g[y1];
x1=g[x0];
x2=g[x1];
denom=(x2-2*x1+x0);
x0=x2;
i=i+1
];
x0
]
steffensen[Cos,0.999,10^-8,30]
NSolve[Cos[x]==x,x,Reals]

Exercise #10, p. 67

fixedPoint code:

(* These are options that one can give at game time: 
-- epsilon: the value we consider "zero" (0 by default)
  -- tolerance: the width of the interval we consider "close enough" to the root
		-- verbose: print the details of the convergence
-- debug: print out a little more for the coder....
*)
Options[fixedPoint]={maxits->100,epsilon->0,tolerance->10^(-6),verbose->True,debug->False} ;

(* Here's the function itself: *)
fixedPoint[f_,a0_,OptionsPattern[]] :=

Module[{a=a0,c,fa=f[a0],fc,slope,width,output,i=1,its=OptionValue[maxits],eps=OptionValue[epsilon],tol=OptionValue[tolerance],verbose=OptionValue[verbose],debug=OptionValue[debug]},
width=10000;
(* Update and continue: let's start looping: *)
While[i<=its && width>tol,
(* c is the function value at a: *)
c=fa;
fc=f[c];
width=Abs[c-a];
slope=(fc-fa)/(c-a);

(* Otherwise, store output and do it again: *)
If[i==1,
(* Let's collect the output, since we've gone to the trouble of computing it: *)
output={{a,c,fc,slope}},
output=Append[output,{a,c,fc,slope}]
]; 

(* Did we stomp on a root by good fortune? *)
If[width<=tol, 
If[verbose,
Print["Root at c."];printoutFixedPoint[output,a,c,fc,width,slope];Return[Transpose[output]],
Return[fc] (* if we're converging, this is better than c *)
]
];

(* update: *)
a=c;
fa=fc;
i=i+1;
]; 

(* We've popped out, so print some nice formatted output: *)
If[verbose,
printoutFixedPoint[output,a,c,fc,width,slope];
(* Return the output, so that it can be plotted, etc. *)
Return[Transpose[output]],
Return[fc]
];
]

This is the printout command for verbose results:

printoutFixedPoint[output_,a_,c_,fc_,width_,slope_]:=
Module[{a0=output[[1]][[1]]},
Print[a0];
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"(Subscript[a, k])","(Subscript[c, k])","f[(Subscript[c, k])]","slope"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",width]; 
Print["f[c] = ",NumberForm[fc,16] ]; 
Print["f'[c] = ",NumberForm[slope,16] ]; 
]

Clear[tolerance];
f[x_]:=CubeRoot[x^2+x]
Plot[{f[x],x},{x,0.5,3.5}]
f'[x]
Plot[f'[x],{x,0.5,3.5}]
results=fixedPoint[f,2.5,maxits->100,verbose->True,tolerance->10^-4];
NSolve[f[x]==x,x,Reals]

15 iterations; the slope at the fixed point is approximately m=0.5393339377770753, and decreasing to the right; and the initial interval was Abs[2.5-2.060642649904278]=0.439357. So we know that to get to within 10^-4, we need to solve for k, where
	m^k*initInterval<tolerance
	Mathematica says 14 interations.

initInterval=Abs[2.5-2.060642649904278];
m=0.5393339377770753;
tol=10^-4;
Ceiling[k/.Solve[m^k*initInterval== tol,k]]

Exercise #13, p. 67:

Steffensen’s only required  4 iterations:

f[x_]:=CubeRoot[x^2+x]
steffensen[f,2.5,10^-4,30]

Exercise #15, p. 67:

aitken[seq_]:=Module[{n=Length[seq]},
TableForm[Table[{seq[[i]],If[i<3,"NA",seq[[i]]-(seq[[i-1]]-seq[[i]])^2/(seq[[i]]-2*seq[[i-1]]+seq[[i-2]])]},{i,1,n}]
]
]
seq=Table[(n+1)/(n-1),{n,2,10}];
aitken[seq]

Exercise #16, p. 67:

p={0.5,0,0,0,0,0,0,0,0,0};
seq=Table[p[[n]]=If[n<2,0.5,(2-E^p[[n-1]]+p[[n-1]]^2)/3],{n,1,10}]
aitken[seq]

p={0.75,0,0,0,0,0,0,0,0,0};
seq=Table[p[[n]]=If[n<2,0.75,Sqrt[E^p[[n-1]]/3]],{n,1,10}];
aitken[seq]

Exercise #24, p. 67

f[x_]:=x-0.002*(Exp[x]Cos[x]-100)
result=steffensen[f,5,10^-8,30];
N[result,17]