Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 2.	The Solution of Nonlinear Equations f(x) = 0
	Algorithm 2.1	  Fixed Point Iteration	
	Algorithm 2.2	  Bisection Method	
	Algorithm 2.3	  False position or Regula Falsi Method	
	Algorithm 2.4	  Approximate Location of Roots	
	Algorithm 2.5	  Newton-Raphson Iteration	
	Algorithm 2.6	  Secant Method	
	Algorithm 2.7	  Steffensen's Acceleration	
	Algorithm 2.8	  Muller's Method	
	Algorithm 2.9	  Nonlinear Seidel Iteration	
	Algorithm 2.10	Newton-Raphson Method in 2-Dimensions

Algorithm 2.1  (Fixed Point Iteration).  To find a solution to the equation x  =  g(x)   by starting with  p0  and iterating   pn+1  =  g(pn).
Terminate the iteration process when the relative error estimate,
2*|pn+1 - pn|/(|pn+1|+Small),  in consecutive iterates  pn  and  pn+1
differs by less than the preassigned value  tolerance.
Section 2.1	Iteration for Solving  x = g(x)   Page  51

Off[General::spell1];
Clear[FixedPointIteration,g0,X0,tolerance,MaxIter,x];
FixedPointIteration[g0_,X0_,tolerance_,MaxIter_,digits_] :=
  Module[{g,k,x0,x1,relerr},
    Set@@{g[x_],g0};
    x0 = N[X0,digits];
    k = 0;
    err = 1.;
    relerr = 1;
    small = N[10^-20,digits];
    Print["\t","\t     p = ",N[x0,digits]];
  While[k<MaxIter && tolerance<=relerr,
    Module[{},
      x1 = g[x0];
      k = k+1;
      err = Abs[x1-x0];
      Print["k = ",k,"\t     p = ",N[x1,digits],
                     "\t  |dp| <= ",N[err,digits]];
      relerr = 2 Abs[x1-x0]/(Abs[x1]+small);
      x0 = x1; ]];
    p = x0;
    Print[" "];
    Print["The function is g[x] = ",g[x]];
    Print["A fixed point of g[x] is:"]; 
    Print["  p   = ",N[p,digits]];
    Print["g[p]  = ",N[g[p],digits]];
    Print["|dp| <= ",N[err,digits]]; ]; 
On[General::spell1];

Chapter 2, Example 2.2, Page 46.

Find the fixed point of the equation    x = exp(-x).

g[x_] = Exp[-x]

Print["The function is g[x] = ",g[x]];
Plot[{x,g[x]},{x,0,1},PlotRange->{{0,1},{0,1}},
PlotLabel->"The line y = x and the curve y = g(x)."];

FixedPointIteration[g[x],0.5,10^-10,30,8]

Compare with Mathematica's built in routine.

solset = FindRoot[g[x]==x,{x,0.5}]

p = solset[[1,2]]

g[p]-p

Mathematica's answer is good.

Chapter 2, Example 2.3, Page 49.

Find the fixed point for the function  g(x)  =  1 + x - x^2/4.

g[x_] = 1 + x - x^2/4

Print["The function is g[x] = ",g[x]];
Plot[{x,g[x]},{x,0,5},PlotRange->{{0,5},{0,5}},
PlotLabel->"The line y = x and the curve y = g(x)."];

FixedPointIteration[g[x],4.0,10^-10,30,8]

Compare with Mathematica's built in routine.

solset = FindRoot[g[x]==x,{x,4}]

p = solset[[1,2]]

g[p]-p

Mathematica's answer is good.

Algorithm 2.2  (Bisection Method).  To find a root of the equation  f(x)  =  0 in the interval  [a,b].   Proceed with the method only if  f(x)  is continuous and  f(a)  and  f(b)  have opposite signs.
Section 2.2	Bracketing Methods for Locating a Root  Page  61

Off[General::spell1];
Clear[Bisection,f0,a0,b0,delta,x];
Bisection[f0_,a0_,b0_,delta_,digits_] :=
Module[{a,b,c,Ya,Yb,Yc,f,k,cond,MaxIter,
        MethodDoesNotApply,PerfectZero,RootInTolerance},
  Set@@{f[x_],f0};
  a = N[a0,digits];
  b = N[b0,digits];
  k = 0;
  Ya = f[a];
  Yb = f[b];
  dx  = b - a;
  cond = 0;
  MaxIter = 1 + Floor[(Log[dx]-Log[delta])/Log[2.0]]//N; 
  If[Ya Yb > 0, cond=MethodDoesNotApply];
  Print["\t","\t     a = ",N[a,digits]];
  Print["\t","\t     b = ",N[b,digits]];
  While[cond==0 && k < MaxIter,
    c  = (a+b)/2;
    Yc = f[c];
    Print["k = ",k,"\t     c = ",N[c,digits],
                   "\t  |dc| <= ",N[dx,digits]];
    If[Yc==0, a=c; b=c; cond=PerfectZero];
    If[Yb Yc>0, b=c;Yb=Yc, a=c;Ya=Yc];
    dx = b-a;
    k = k+1; ];
  If[dx<delta && cond!=PerfectZero, cond=RootInTolerance];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  c   = ",N[c,digits]];
  Print["f[c]  = ",N[f[c],digits]];
  Print["|dc| <= ",N[dx,digits]]; ]; 
On[General::spell1];

Chapter 2, Example 2.6, Page 58.

Find the solution to the equation   0  =  x sin(x)  -  1.

Clear[f,x];
f[x_] = x Sin[x] - 1

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,0,2},
PlotLabel->"The curve y = f(x)."];

Bisection[f[x],0,2,10^-10,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,1}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

Algorithm 2.3 (False position or Regula Falsi Method). To find a root of the equation  f(x)  =  0  in the interval  [a,b].   Proceed with the method only if  f(x)  is continuous and f(a) and f(b) have opposite signs.
Section 2.2	Bracketing Methods for Locating a Root   Page  62

Off[General::spell1];
Clear[FalsePosition,f0,a0,b0,delta,epsilon,x];
FalsePosition[f0_,a0_,b0_,delta_,epsilon_,digits_] :=
Module[{a,b,c,Ya,Yb,Yc,f,k,cond,Iterating,
        RootInTolerance,PerfectZero,MethodDoesNotApply},
  Set@@{f[x_],f0};
  a = N[a0,digits];
  b = N[b0,digits];
  k = 0;
  Ya = f[a];
  Yb = f[b];
  dx  = b - a;
  cond = 0;
  MaxIter = 100; 
  If[Ya Yb > 0, cond=MethodDoesNotApply];
  Print["\t","\t     a = ",N[a,digits]];
  Print["\t","\t     b = ",N[b,digits]];
  While[cond==0 && k < MaxIter,
    dx = Yb (b-a)/(Yb-Ya);
    c  = b - dx;
    Yc = f[c];
    Print["k = ",k,"\t     c = ",N[c,digits],
                   "\t  |dc| <= ",N[dx,digits]];
    If[Yc==0, a=c; b=c; cond=PerfectZero];
    If[Yb Yc>0, b=c;Yb=Yc, a=c;Ya=Yc];
    If[Abs[dx]<delta && Abs[Yc]<epsilon,
      cond=RootInTolerance];
    k = k+1;  ];
  If[cond==0, cond=Iterating];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  c   = ",N[c,digits]];
  Print["f[c]  = ",N[f[c],digits]];
  Print["|dc| <= ",N[dx,digits]]; ];
On[General::spell1];

Chapter 2, Example 2.7, Page 60.

Find the solution to the equation   0  =  x sin(x)  -  1.

Clear[f,x];
f[x_] = x Sin[x] - 1

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,0,2},
PlotLabel->"The curve y = f(x)."];

FalsePosition[f[x],0,2,10^-10,10^-10,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,1.0997502}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

Algorithm 2.4  (Approximate Location of Roots).  To roughly estimate locations of the roots of the equation  f(x)  =  0  over the interval  [a,b] by using the equally spaced sample points  (xk, f(xk))  
and the criteria:
		(i)  yk-1 yk  <  0,
		or
		(ii)  |yk|  <  e    and   (yk - yk-1)(yk+1 - yk)  <  0.
That is,  either  f(xk-1) and f(xk)  have opposite signs, or  |f(xk)|  is small and the slope of the curve  y = f(x)  changes sign near  (xk,f(xk)).
Section 2.3	Initial Approx. and Convergence Criteria   Page  70

Off[General::spell1];
Clear[LocateRoots,f0,a0,b0,n0,epsilon,R,X,Y];
LocateRoots[f0_,a0_,b0_,n0_,epsilon_,digits_] :=
Module[{a,b,f,i,k,n,s,maxY,minY,small,X,Y},
  Set@@{f[x_],f0};
  a = N[a0,digits];
  b = N[b0,digits];
  n = n0;
  h = (b-a)/n;
  X = Table[i,{i,a,b+h,h}];
  Y = f[X];
  R = {};
  maxY = Max[Y];
  minY = Min[Y];;
  small = (maxY-minY) epsilon;
  Y[[n+2]] = Y[[n+1]];
  For[k=2,k<=n+1,k++,
    Module[{},
      If[Y[[k-1]] Y[[k]]<=0,
        root = (X[[k-1]]+X[[k]])/2;
        R = Append[R,root]; ];
      s = (Y[[k]]-Y[[k-1]]) (Y[[k+1]]-Y[[k]]);
      If[Abs[Y[[k]]]<small && s<0,
        root = X[[k]];
        R = Append[R,root]; ]; ]; ];
  Print["The function is f[x] = ",N[f[x],digits]];
  Print["The interval is [",a,",",b,"]"];
  Print["The approximate location for the roots are:"];
  Return[N[R,digits]]; ];
On[General::spell1];

Chapter 2, Example 2.8, Page 66.

Find the solution to the equation   0  =  x^3  -  x^2  -  x  +  1.

Clear[f,x];
f[x_] = x^3 - x^2 - x + 1

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,-1.2,1.5},
PlotLabel->"The curve y = f(x)."];

LocateRoots[f[x],-1.2,1.2,50,0.05,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,-0.984}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

solset = FindRoot[f[x]==0,{x,1.008}]

p = solset[[1,2]]

f[p]

Mathematica's answer NOT so good,
need to adjust the number of iterations.

Clear[f,x];
f[x_] = Cos[x] - 0.6 Cos[2x] + 0.9 Cos[3x]

LocateRoots[f[x],0,2Pi,50,0.05,8]

Algorithm 2.5  (Newton-Raphson Iteration).  To find a root of
f(x) = 0  given one initial approximation  p0  and using the iteration:
		pk+1  =  pk - f(pk)ÊÊf'(pk)Ê     for   k = 0, 1, ... .
Section 2.4	Newton-Raphson and Secant Methods   Page  84

Off[General::spell1];
Clear[NewtonRaphson,f0,x0,delta,epsilon,max,x];
NewtonRaphson[f0_,x0_,delta_,epsilon_,max_,digits_] :=
Module[{f,k,small,cond,AbscissaInTolerance,
        OrdinateInTolerance,RootInTolerance,Iterating},
  Set@@{f[x_],f0};
  p0 = x0;
  Y0 = N[f[p0],digits];
  p1 = p0+1;
  k = 0;
  cond = 0;
  small = N[10^(-20),digits];
  Print["p = ",N[p0,digits],"\t   f[p] = ",N[Y0,digits]];
  While[cond==0 && k < max,
    Df = N[f'[p0],digits];
    If[Df==0, cond=DivByZero;Dp=p1-p0;p1=p0,
              Dp=Y0/Df;p1=p0-Dp];
    Y1 = N[f[p1],digits];
    Print["p = ",N[p1,digits],"\t   f[p] = ",N[Y1,digits],
                    "\t  |dp| <= ",N[Dp,digits]];
    RelErr = N[Abs[Dp]/(Abs[p1]+small),digits];         
    If[RelErr < delta,  cond=AbscissaInTolerance];
    If[Abs[Y1]<epsilon, cond=OrdinateInTolerance];
    If[RelErr<delta && Abs[Y1]<epsilon,
       cond=RootInTolerance];
    p0 = p1;
    Y0 = Y1;
    k  = k+1;  ];
  If[cond==0, cond=Iterating];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  p   = ",N[p1,digits]];
  Print["f[p]  = ",N[Y1,digits]];
  Print["|dp| <= ",N[Abs[Dp],digits]]; ];
On[General::spell1];

Chapter 2, Example 2.10, Page 75.

Find the solution to the equation    0  =  4800(1 - Exp[-x/10]) - 320 x.

Clear[f,x];
f[x_] = 4800(1 - Exp[-x/10]) - 320 x

Clear[r,t];
r[t_] = 1600(1 - Exp[-t/10])

ParametricPlot[{r[t],f[t]},{t,0,8.8}];

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,0,9},
PlotLabel->"The curve y = f(x)."];

NewtonRaphson[f[x],8.0,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,8}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

Algorithm 2.6  (Secant Method).  To find a root of  f(x) = 0  given
two initial approximations  p0  and  p1  and using the iteration
pk+1  =  pk - f(pk)(pkÊ-Êpk-1)/(f(pk) - f(pk-1))    for  k=,2, ... .
Section 2.4	Newton-Raphson and Secant Methods   Page  85

Off[General::spell1];
Clear[SecantMethod,f0,x0,x1,delta,epsilon,max,x];
SecantMethod[f0_,x0_,x1_,delta_,epsilon_,max_,digits_] :=
Module[{f,k,small,cond,DivByZero,AbscissaInTolerance,
        OrdinateInTolerance,RootInTolerance,Iterating},
  Set@@{f[x_],f0};
  p0 = x0;
  Y0 = N[f[p0],digits];
  p1 = x1;
  Y1 = N[f[p1],digits];
  k = 0;
  cond = 0;
  small = N[10^(-20),digits];
  Print["p = ",N[p0,digits],"\t   f[p] = ",N[Y0,digits]];
  Print["p = ",N[p1,digits],"\t   f[p] = ",N[Y1,digits]];
  While[cond==0 && k < max,
    Df = (Y1-Y0)/(p1-p0);
    If[Df==0, cond=DivByZero;Dp=p1-p0;p2=p1,
              Dp=Y1/Df;p2=p1-Dp];
    Y2 = N[f[p2],digits];
    Print["p = ",N[p2,digits],"\t   f[p] = ",N[Y2,digits],
                    "\t  |dp| <= ",N[Dp,digits]];
    RelErr = N[Abs[Dp]/(Abs[p2]+small),digits];         
    If[RelErr < delta,  cond=AbscissaInTolerance];
    If[Abs[Y2]<epsilon, cond=OrdinateInTolerance];
    If[RelErr<delta && Abs[Y2]<epsilon,
      cond=RootInTolerance];
    p0 = p1;
    Y0 = Y1;
    p1 = p2;
    Y1 = Y2;
    k  = k+1;  ];
  If[cond==0, cond=Iterating];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  p   = ",N[p1,digits]];
  Print["f[p]  = ",N[Y1,digits]];
  Print["|dp| <= ",N[Dp,digits]]; ];
On[General::spell1];

Chapter 2, Example 2.13, Page 82.

Find the solution to the equation   0  =  x^3  -  3x  +  2.

Clear[f,x];
f[x_] = x^3 - 3x + 2

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,-2.6,2.0},
PlotLabel->"The curve y = f(x)."];

SecantMethod[f[x],-2.6,-2.4,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,-2.4}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

SecantMethod[f[x], 1.4, 1.2,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,1.2}]

p = solset[[1,2]]

f[p]

Mathematica's answer NOT so good,
need to adjust the number of iterations.

Algorithm 2.7  (Steffensen's Acceleration).  To quickly find
a solution of the fixed point equation  x = g(x)  given an initial approximation  p0, where it is assumed that  both  g(x)  and  g '(x)
are continuous and  |g'(p)| < 1  and that ordinary fixed point iteration converges linearly (slowly) to  p.
Section 2.5	Aitken's, Steffensen's, Muller's Methods   Page  96

Off[General::spell1];
Clear[Steffensen,f0,x0,delta,epsilon,max,x];
Steffensen[f0_,x0_,delta_,epsilon_,max_,digits_] :=
Module[{f,k,small,cond,AbscissaInTolerance,
        OrdinateInTolerance,RootInTolerance,Iterating},
  Set@@{f[x_],f0};
  p0 = N[x0,digits];
  p3 = p0;
  p2 = p0+1;
  p1 = p0+2;
  Y3 = N[f[p3],digits];
  k = 0;
  cond  = 0;
  small = N[10^(-20),digits];
  Print["p = ",N[p0,digits],"\t   f[p] = ",N[f[p0],digits]];
  While[cond==0 && k < max,
    p0 = p3;
    Y0 = Y3;
    Df0 = N[f'[p0],digits];
    If[Df0==0, cond=1; Dp=p3-p2; p3=p0,
               p1=p0-Y0/Df0];
    Y1 = N[f[p1],digits];
    Print["p = ",N[p1,digits],"\t   f[p] = ",N[Y1,digits],
                              "\t  |dp| <= ",N[Dp,digits]];
    Df1 = N[f'[p1],digits];
    If[Df1==0 && cond==0, cond=1; Dp=p1-p0; p3=p1,
                 p2=p1-Y1/Df1];
    Print["p = ",N[p2,digits],"\t   f[p] = ",N[Y2,digits],
                              "\t  |dp| <= ",N[Dp,digits]];
    D1 = (p1-p0)*(p1-p0);
    D2 = p2 - 2 p1 + p0;
    If[D2==0 && cond==0, cond=1; Dp=p2-p1; p3=p2,
                         p3=p0-D1/D2; Dp=p3-p2];
    Y3 = f[p3];
    Print["p = ",N[p0,digits],"\t   f[p] = ",N[Y3,digits],
                              "\t  |dp| <= ",N[Dp,digits]];
    RelErr = N[Abs[Dp]/(Abs[p3]+small),digits];         
    If[RelErr < delta,  cond=AbscissaInTolerance];
    If[Abs[Y3]<epsilon, cond=OrdinateInTolerance];
    If[RelErr<delta && Abs[Y3]<epsilon,
                        cond=RootInTolerance];         
    k  = k+1;  ];
  If[cond==0, cond=Iterating];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  p   = ",N[p3,digits]];
  Print["f[p]  = ",N[Y3,digits]];
  Print["|dp| <= ",N[Dp,digits]]; ];
On[General::spell1];

Chapter 2, Example 2.16, Page 95.

Find the solution to the equation   0  =  x^3  -  3x  +  2.

Clear[f,x];
f[x_] = x^3 - 3x + 2

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,-2.6,2.0},
PlotLabel->"The curve y = f(x)."];

Steffensen[f[x],-2.4,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,-2.4}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

Steffensen[f[x], 1.2,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,1.2}]

p = solset[[1,2]]

f[p]

Mathematica's answer NOT so good,
need to adjust the number of iterations.

Algorithm 2.8  (Muller's Method). To find a root of the equation
f(x) = 0  given three distinct initial approximation  p0, p1, p2.  
Section 2.5	Aitken's, Steffensen's, Muller's Methods   Page  97

Off[General::spell1];
Clear[Muller,f0,x0,x1,x2,delta,epsilon,max,x];
Muller[f0_,x0_,x1_,x2_,delta_,epsilon_,max_,digits_] :=
Module[{a,b,c,det,f,k,Z,small,satisfied,AbscissaInTolerance,
        OrdinateInTolerance,RootInTolerance},
  Set@@{f[x_],f0};
  p0 = x0;  Y0 = N[f[p0],digits];
  p1 = x1;  Y1 = N[f[p1],digits];
  p2 = x2;  Y2 = N[f[p2],digits];
  k = 0;   
  small = N[10^(-20),digits];
  Print["p = ",N[p0,digits],"\t   f[p] = ",N[Y0,digits]];
  Print["p = ",N[p1,digits],"\t   f[p] = ",N[Y1,digits]];
  Print["p = ",N[p2,digits],"\t   f[p] = ",N[Y2,digits]];
  satisfied = 0;
  While[satisfied==0 && k < max,
    h0 = p0-p2;
    h1 = p1-p2;
    c  = Y2;
    e0 = Y0-c;
    e1 = Y1-c;
    det = h0 h1 (h0 - h1);
    a = (e0 h1 - h0 e1)/det;
    b = (h0 h0 e1 - h1 h1 e0)/det;
    If[b b>4 a c, disc=Sqrt[b b-4 a c],
                  disc=0];
    If[b<0, disc=-disc];                
    Z  = - 2 c/(b + disc);
    p3 = p2 + Z;
    If[Abs[p3-p1]<Abs[p3-p0],  
       U=p1; p1=p0; p0=U; V=Y1; Y1=Y0; Y0=V;];                
    If[Abs[p3-p2]<Abs[p3-p1],  
       U=p2; p2=p1; p1=U; V=Y2; Y2=Y1; Y1=V;];                
    p2 = p3;
    Y2 = f[p2];
    Print["p = ",N[p2,digits],"\t   f[p] = ",N[Y2,digits],
                              "\t  |dp| <= ",N[Abs[Z],digits]];
    RelErr = Abs[Z]/(Abs[p2] + small);
    If[RelErr < delta,  satisfied=AbscissaInTolerance];
    If[Abs[Y2]<epsilon, satisfied=OrdinateInTolerance];
    If[RelErr<delta && Abs[Y2]<epsilon,
       satisfied=RootInTolerance];
    k = k+1;  ];
  If[satisfied==0, satisfied=Iterating];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A root of f[x] is:"];
  Print["  p   = ",N[p2,digits]];
  Print["f[p]  = ",N[Y2,digits]];
  Print["|dp| <= ",N[Abs[Z],digits]]; ];
On[General::spell1];

Chapter 2, Example 2.16, Page 95.

Find the solution to the equation   0  =  x^3  -  3x  +  2.

Clear[f,x];
f[x_] = x^3 - 3x + 2

Print["The function is f[x] = ",f[x]];
Plot[f[x],{x,-2.6,2.0},
PlotLabel->"The curve y = f(x)."];

Muller[f[x],-2.6,-2.5,-2.4,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,-2.4}]

p = solset[[1,2]]

f[p]

Mathematica's answer is good.

Muller[f[x], 1.4, 1.3, 1.2,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

solset = FindRoot[f[x]==0,{x,1.2}]

p = solset[[1,2]]

f[p]

Mathematica's answer NOT so good,
need to adjust the number of iterations.

Algorithm 2.9  (Nonlinear Seidel Iteration).  To solve the non-linear fixed point system 	
			x  =  g1(x,y,z), 
			y  =  g2(x,y,z), 
			z  =  g3(x,y,z), 
given one initial approximation P0 = (p0,q0,r0), and generating a sequence
{Pk} = {(pk,qk,rk)}  which converges to the solution   P = (p,q,r),
i.e.   p = g1(p,q,r),  q = g2(p,q,r)  and  r = g3(p,q,r).
Section 2.6	Iteration for Nonlinear Systems   Page  108

Remark. First we write a program which performs fixed point iteration.  As we shall see it is simple to write.  Seidel iteration will be more complicated because it needs to use the updated first coordinate before the second coordinate is computed.

First, execute the following group of cells:

Off[General::spell1];
Clear[FixPointIter,G0,X0,V,tolerance,x,y,max];
FixPointIter[G0_,X0_,tolerance_,max_,digits_] :=
Module[{G,k},
  Set@@{G[{x_,y_}],G0};
  norm[V_] := N[Sqrt[V^2/.List->Plus],digits];
  P1 = N[X0,digits];
  Print["  P  = ",P1];
  k  = 0;
  dist = 1;
  tol = N[tolerance,digits];
  While[k <= max && dist > tol,
    P0 = P1;
    k = k + 1;
    P1 = N[G[P0],digits];
    dist = norm[P1 - P0];
    Print["  P  = ",N[P1,digits],"\t |dP| <= ",N[dist,digits]];];
  Print[" "];
  Print["G[{x,y}] = ",G[{x,y}]];
  Print["The fixed point is:"];
  Print["  P   = ",N[P1,digits]];
  Print["G[P]  = ",N[G[P1],digits]];
  Print["|dP| <= ",N[dist,digits]]; ]; 
On[General::spell1];

Off[General::spell1];
Clear[Seidel,g10,g20,X0,V,tolerance,max,x,y];
Seidel[g10_,g20_,X0_,tolerance_,max_,digits_] :=
Module[{g1,g2,k,Pold,Pnew,Qold,Qnew},
  Set@@{g1[{x_,y_}],g10};
  Set@@{g2[{x_,y_}],g20};
  norm[V_] := N[Sqrt[V^2/.List->Plus],digits];
  k = 0;
  dist = 1;
  P1 = N[X0,digits];
  Print["  P  = ",P1];
  Pnew =  N[P1[[1]],digits];
  Qnew =  N[P1[[2]],digits];
  tol = N[tolerance,digits];
  While[k < max && dist > tol,
    Pold = Pnew;
    Qold = Qnew;
    P0 = {Pold,Qold};
    k = k+1;
    Pnew = g1[{Pold,Qold}];
    Qnew = g2[{Pnew,Qold}];
    P1 = {Pnew,Qnew};
    dist = norm[P1 - P0];
    Print["  P  = ",N[P1,digits],"\t |dP| <= ",N[dist,digits]];];
  Print[" "];
  Print["G[{x,y}] = ",G[{x,y}]];
  Print["A fixed point is:"];
  Print["  P   = ",N[P1,digits]];
  Print["G[P]  = ",N[G[P1],digits]];
  Print["|dP| <= ",N[dist,digits]]; ];
On[General::spell1];

Chapter 2, Example in Table 2.14, on Page 102.

Use fixed point iteration to find the solution to the nonlinear system:
     x  =  g1(x,y)
     y  =  g2(x,y)
where:

Clear[g1,g2,G,x,y];
g1[{x_,y_}] = (x^2 - y + 0.5)/2;
g2[{x_,y_}] = (- x^2 - 4y^2 + 8y + 4)/8;
G[{x_,y_}] = {g1[{x,y}],g2[{x,y}]}

graph1 = ContourPlot[g1[{x,y}]-x,{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
graph2 = ContourPlot[g2[{x,y}]-y,{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
Print["x = ",g1[{x,y}],"   and   ","y = ",g2[{x,y}]];
Show[graph1,graph2,DisplayFunction->$DisplayFunction];

FixPointIter[G[{x,y}],{0.0,1.0},10^-10,100,8]

Compare with Mathematica's built in routine.

FindRoot[{g1[{x,y}]==x,g2[{x,y}]==y},{x,0},{y,1}]

{g1[{-0.222215,0.993808}]==x,g2[{-0.222215,0.993808}]}

Mathematica's answer is o.k.,
need to adjust the number of iterations.

Chapter 2, Example in Table 2.15, on page 103.

Use Seidel iteration to find the solution to the nonlinear system:
     x  =  g1(x,y)
     y  =  g2(x,y)
where:

Clear[g1,g2,G,x,y];
g1[{x_,y_}] = (x^2 - y + 0.5)/2;
g2[{x_,y_}] = (- x^2 - 4y^2 + 8y + 4)/8;
G[{x_,y_}] = {g1[{x,y}],g2[{x,y}]}

graph1 = ContourPlot[g1[{x,y}]-x,{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
graph2 = ContourPlot[g2[{x,y}]-y,{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
Print["x = ",g1[{x,y}],"   and   ","y = ",g2[{x,y}]];
Show[graph1,graph2,DisplayFunction->$DisplayFunction];

Seidel[g1[{x,y}],g2[{x,y}],{2.0,0.0},10^-10,15,8]

Compare with Mathematica's built in routine.

FindRoot[{g1[{x,y}]==x,g2[{x,y}]==y},{x,2},{y,0}]

{g1[{1.90068,0.311219}],g2[{1.90068,0.311219}]}

Mathematica's answer is o.k.,
need to adjust the number of iterations.

Algorithm 2.10  (Newton-Raphson Method in 2-Dimensions). To solve the non-linear system 	
				0  =  f1(x,y), 
				0  =  f2(x,y), 
given one initial approximation  (p0,q0)  and using Newton-Raphson iteration.
Section 2.7	Newton's Method for Systems   Page  116

Off[General::spell1];
Clear[NewtonSystem,F00,X0,V,delta,epsilon,max,x,y];
NewtonSystem[F00_,X0_,delta_,epsilon_,max_,digits_] :=
Module[{F,k,cond,det,small,DivByZero,DomainInTolerance,
        RangeInTolerance,SolutionInTolerance,Iterating},
  Set@@{F[{x_,y_}],F00};
  norm[V_] := N[Sqrt[V^2/.List->Plus],digits];
  n = Length[X0];
  Dp = Table[0,{n}];
  P0 = N[X0,digits];
  small = delta;
  k = 0;
  cond = 0;
  F0 = N[F[P0],digits];
  Print["  P  = ",N[P0,digits]];
  Print["F[P] = ",N[F0,digits]];
  P1 = P0;
  F1 = F0;
  While[k < max && cond==0,
    P0 = P1;
    F0 = F1;
    k  = k+1;
    J0 = N[jacobian[P0],digits];
    det = Det[J0];
    If[det==0, Dp = Table[0,{n}]; cond=DivByZero,
               Dp = Inverse[J0].F0];
    P1 = P0 - Dp;
    F1 = N[F[P1],digits];
    Print[" "];
    Print["  P  = ",N[P1,digits],"\t |dP| <= ",N[norm[Dp],digits]];
    Print["F[P] = ",N[F1,digits]];
    RelErr = norm[Dp]/(norm[P1]+small);
    FnZero = norm[F1];
    If[RelErr < delta, cond=DomainInTolerance];
    If[FnZero<epsilon, cond=RangeInTolerance];
    If[RelErr<delta && FnZero<epsilon,
       cond=SolutionInTolerance] ];
  If[cond==0, cond=Iterating];
  Print[" "];
  Print["F[{x,y}] = ",F[{x,y}]];
  Print["A solution to the system is:"];
  Print["  P   = ",N[P1,digits]];
  Print["F[P]  = ",N[F1,digits]];
  Print["|dP| <= ",N[norm[Dp],digits]]; ];
On[General::spell1];

Chapter 2, Example 2.20, Page 114.

Use Newton's method to solve the nonlinear system:
     0  =  f1(x,y)
     0  =  f2(x,y)
where:

Clear[f1,f2,F,x,y];
f1[{x_,y_}] = x^2 - 2 x - y + 0.5;
f2[{x_,y_}] = x^2 + 4 y^2 - 4;
F[{x_,y_}] = {f1[{x,y}],f2[{x,y}]}

jacobian[{x_,y_}] =
         Transpose[Map[D[F[{x,y}],#]&,{x,y}]]//N

graph1 = ContourPlot[f1[{x,y}],{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
graph2 = ContourPlot[f2[{x,y}],{x,-2,2},{y,-1,1},
         Compiled->True,Contours->{0},PlotPoints->20,
         ContourShading->False,DisplayFunction->Identity];
Print["0 = ",f1[{x,y}],"   and   ","0 = ",f2[{x,y}]];
Show[graph1,graph2,DisplayFunction->$DisplayFunction];

P0 = {2.0,0.25};
NewtonSystem[F[{x,y}],P0,10^-10,10^-10,100,8]

Compare with Mathematica's built in routine.

FindRoot[{f1[{x,y}]==0,f2[{x,y}]==0},{x,2},{y,0.25}]

{f1[{1.90068,0.311219}],f2[{1.90068,0.311219}]}

Mathematica's answer NOT so good,
need to adjust the number of iterations.