Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

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.

Be sure to load the following "ImplicitPlot" graphics package.

<<  Graphics`ImplicitPlot`Clear["Global`*"];
Off[General::"spell1"];
Clear[a0, b0, delta, epsilon, f, f0, f1, f2, F, F00, FalsePosition, 
        FixPointIter, FixedPointIteration, g, g0, g1, g2, gr1, gr2, graph1, 
        graph2, G, G0, jacobian, LocateRoots, max, Muller, n0, NewtonRaphson, 
        NewtonSystem, p, R, SecantMethod, Steffensen, solset, tolerance, V, 
        x, x0, x1, x2, y, X, X0, Y];
On[General::"spell1"];

Algorithm 2.1  Fixed Point Iteration.  To find a solution of the equation  x = g (x)   by starting with  (Subscript[p, 0])  and iterating   (Subscript[p, n + 1])  =  g ((Subscript[p, n])).   Terminate the iteration process when the relative error estimate,   2|(((Subscript[p, n + 1]) -  (Subscript[p, n]))/(small  +   ∣ (Subscript[p, n + 1]) ∣ 
                ))|,   in consecutive iterates  (Subscript[p, 0])  and  (Subscript[p, n + 1])  differs by less than the preassigned value  tolerance.
Section   2.1   Iteration for Solving  x = g(x)   Page  51

Execute this cell to activate Fixed Point Iteration.
FixedPointIteration[g0_,X0_,tolerance_,max_,digits_]
g0_ the function g(x),
X0_ the starting value,
tolerance_ the preassigned value tolerance,
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output.

Off[General::"spell1"]; \(Clear[FixedPointIteration, g0, X0, tolerance, max, x]; 

    FixedPointIteration[g0_, X0_, tolerance_, max_, digits_] := 
	
      Module[{g, k, x0, x1, relerr}, 
		Set@@{g[x_], g0}; 
		
        x0 = N[X0, digits]; 
		k = 0; 
		err = 1.; 
		relerr = 1; 

        		small = N[1\/10\^20, digits]; 
		
        Print["p"\_0, " = ", N[x0, digits]]; 
		
        While[k < max && tolerance ≤ relerr, 
			
          Module[{}, 
				x1 = g[x0]; 
				k = k + 1; 
				
            err = Abs[x1 - x0]; 
				
            Print["p"\_k, " = ", N[x1, digits], 
              ",   |Δp| ≤ ", 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["|Δp| ≤ ", N[err, digits]]; ]; 

    On[General::"spell1"]; \)

Chapter 2, Example 2.2, Page 46.   Find the fixed point of the equation  x = (e^(-x)).

\(Clear[g, gr1, gr2, x]; 
g[x_] = E\^\(-x\); 

    gr1\  = Plot[g[x], {x, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; 

    gr2\  = Plot[x, {x, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity]; 

    Show[gr1, gr2, PlotLabel -> "y = x and the curve y = g(x).", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is g[x] = ", g[x]]; \)

Use fixed point iteration to find a numerical approximation.

FixedPointIteration[g[x], 0.5, 0.00000000001, 30, 12];

Compare with Mathematica's built in routine, first do it numerically.

\(solset\  = \ FindRoot[g[x] == x, {x, 0.5}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["g[p] = ", g[p]]; \)

For this equation, Mathematica can solve it symbolically !

\(solset\  = \ Solve[g[x] == x, x]; 

    p = solset[[1, 1, 2]]; 

    Print[solset]; 
Print["   p = ", p]; 
Print["   p = ", N[p]]; 

    Print["g[p] = ", g[p]]; 
Print["g[p] = ", FullSimplify[g[p]]]; 

    Print["g[p] = ", N[g[p]]]; \)

Mathematica's answers are very 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; 

    gr1\  = Plot[g[x], {x, 0, 4.82}, PlotRange -> {{0, 4.82}, {0, 2.1}}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; 

    gr2\  = Plot[x, {x, 0, 4.82}, PlotRange -> {{0, 4.82}, {0, 2.1}}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction -> Identity]; 

    Show[gr1, gr2, AspectRatio -> 2.1\/4.82, 
      PlotLabel -> "y = x and the curve y = g(x).", 
      DisplayFunction -> $DisplayFunction]; 

    Print["The function is g[x] = ", g[x]]; \)

Use fixed point iteration to find a numerical approximation.

FixedPointIteration[g[x], 4.0, 0.000000001, 30, 14];

Compare with Mathematica's built in routine, first do it numerically.

\(solset\  = \ FindRoot[g[x] == x, {x, 0.5}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["g[p] = ", g[p]]; \)

This is really easy for Mathematica to solve symbolically.

\(solset\  = \ Solve[g[x] == x, x]; 

    p = solset[[1, 1, 2]]; 

    Print[solset]; 
Print["   p = ", p]; 
Print["g[p] = ", g[p]]; 
    \)

Mathematica's answers are 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

Execute this cell to activate the Bisection Method.
Bisection[f0_,a0_,b0_,delta_,digits_]
f0_ the function f(x),
a0_ the left endpoint a,
b0_ the right endpoint b,
delta_ estimate for the accuracy of the solution,
digits_ the number of digits displayed in the output.

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 = N[1 + Floor[\(Log[dx] - Log[delta]\)\/Log[2.0]]]; 
        If[Ya\ Yb > 0, cond = MethodDoesNotApply]; 
		
        Print[" a = ", N[a, digits]]; 
		
        Print[" b = ", N[b, digits]]; 
        While[cond == 0 && k < MaxIter, c = \(a + b\)\/2; 
			Yc = f[c]; 
          
			Print["c"\_k, " = ", N[c, digits], 
            ",   |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; 

    Plot[f[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use the bisection method to find a numerical approximation to the root.

Bisection[f[x], 0, 2, 0.00000000001, 12];

Compare with Mathematica's built in routine, first do it numerically.

\(solset\  = \ FindRoot[f[x] == 0, {x, 0.5}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", N[p, 12]]; 
Print["f[p] = ", N[f[p], 12]]; 
    \)

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

Execute this cell to activate the False position or Regula Falsi Method.
FalsePosition[f0_,a0_,b0_,delta_,epsilon_,digits_]
f0_ the function f(x),
a0_ the left endpoint a,
b0_ the right endpoint b,
delta_ estimate for the accuracy of the solution,
digits_ the number of digits displayed in the output.

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[" a = ", N[a, digits]]; 
		
        Print[" b = ", N[b, digits]]; 
		
        While[cond == 0 && k < MaxIter, 
			
          dx = \(Yb\ \((b - a)\)\)\/\(Yb - Ya\); 
			c = b - dx; 
			
          Yc = f[c]; 
			
          Print["c"\_k, " = ", N[c, digits], 
            ",   |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; 

    Plot[f[x], {x, 0, 2}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use the Regula Falsi method to find a numerical approximation to the root.

FalsePosition[f[x], 0, 2, 0.00000000001, 0.00000000001, 8];

Compare with Mathematica's built in routine, first do it numerically.

\(solset\  = \ FindRoot[f[x] == 0, {x, 0.5}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", N[p, 12]]; 
Print["f[p] = ", N[f[p], 12]]; 
    \)

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  ((Subscript[x, k]), f ((Subscript[x, k])))  and the criteria:
	(i)  (Subscript[y, k - 1]) (Subscript[y, k]) < 0,
		or
	(ii) (Subscript[(| y), k]) | 
         <  ϵ   and   ((Subscript[y, k]) - (Subscript[y, k - 1])) ((Subscript[y, k + 1]) - (Subscript[y, k])) < 0.
That is,  either  f ((Subscript[x, k - 1]))  and  f ((Subscript[x, k]))  have opposite signs, or  | f ((Subscript[x, k])) |   is small 
and the slope of the curve  y = f (x)  changes sign near  ((Subscript[x, k]), f ((Subscript[x, k]))).
Section  2.3   Initial Approx. and Convergence Criteria   Page  70

Execute this cell to activate the Approximate Location of Roots.
LocateRoots[f0_,a0_,b0_,n0_,epsilon_,digits_]
f0_ the function f(x),
a0_ the left endpoint a,
b0_ the right endpoint b,
n0_ the number of subintervals,
epsilon_ estimate for the size of the ordinate,
digits_ the number of digits displayed in the output.

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 = 1\/2\ 
                  \((X[[k - 1]] + 
                      X[[k]])\); 

              					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:"]; 
        Print[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; 

    Plot[f[x], {x, \(-1.2\), 1.5}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

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

Compare with Mathematica's built in routine.

\(solset\  = \ FindRoot[f[x] == 0, {x, \(-0.984\)}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.008}, MaxIterations -> 50, 
        WorkingPrecision -> 20]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.008}, MaxIterations -> 50, 
        WorkingPrecision -> 30]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Chapter 2, Another Example for, Page 66.
Find the solution to the equation  0 = cos (x) - 0.6 cos (2 x) + 0.9 cos (3  x).

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

    Plot[f[x], {x, 0, 2\ π}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

LocateRoots[f[x], 0, 2  π, 50, 0.02, 6]

Algorithm 2.5  Newton-Raphson Iteration.  To find a root of  f (x) = 0  given one initial approximation  (Subscript[p, 0])   using the iteration:
		(Subscript[p, k + 1]) = (Subscript[p, k]) -((f ((Subscript[p, k])))/(f'((Subscript[p, k]))))   for    k  =  0, 1, ⋯ , m . 
Section  2.4   Newton-Raphson and Secant Methods   Page  84

Execute this cell to activate the Newton-Raphson Iteration.
NewtonRaphson[f0_,x0_,delta_,epsilon_,max_,digits_]
f0_ the function f(x),
x0_ the starting value,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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[(1/(10^20)), digits];
		\(Print[\*""\!\(p\_0\) = "", N[p0, digits], 
                \*"",   f[\!\(p\_0\)] = "", 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[(Subscript["p", k + 1]), " = ", N[p1, digits], 
                      ",   f[", (Subscript["p", k + 1]), "] = ", 
                      N[Y1, digits], ",   |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 - (e^(((-x)/10)))) - 320 x.

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

    r[t_] = 1600\ \((1 - E\^\(\(-t\)\/10\))\); 

    Print["f[x] = ", f[x]]; 
Print["f[t] = ", f[t]]; 

    Print["r[t] = ", r[t]]; \)

\(ParametricPlot[{r[t], f[t]}, {t, 0, 8.8}, 
      PlotStyle -> RGBColor[0, 1, 0]]; 

    Print["The path of the particle {r[t],f[t]}."]; \)

\(Plot[f[x], {x, 0, 9}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use the Newton-Raphson method to find a numerical approximation to the root.

NewtonRaphson[f[x], 8.0, 0.00000000001, 0.00000000001, 100, 8];

Compare with Mathematica's built in routine.

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

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 8}, MaxIterations -> 50, 
        WorkingPrecision -> 20]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Algorithm 2.6  Secant Method.  To find a root of  f (x) = 0  given two initial approximation  (Subscript[p, 0])   and  (Subscript[p, 1])   using the iteration
	(Subscript[p, k + 1])  =  (Subscript[p, k])- f((Subscript[p, k]))(((Subscript[p, k])  -  (Subscript[p, k - 1]))/(f ((Subscript[p, k]))  -  f ((Subscript[p, k - 1]))))   for    k  =  1, 2, ⋯ , m .
Section  2.4   Newton-Raphson and Secant Methods   Page  85

Execute this cell to activate the Secant Method.
SecantMethod[f0_,x0_,x1_,delta_,epsilon_,max_,digits_]
f0_ the function f(x),
x0_ the starting value,
x2_ ,x1_ the starting values,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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 = 1; 
		cond = 0; 

        		small = N[1\/10\^20, digits]; 
		
        Print[\*""\!\(p\_0\) = "", N[p0, digits], 
          \*"",   f[\!\(p\_0\)] = "", N[Y0, digits]]; 
        Print[\*""\!\(p\_1\) = "", N[p1, digits], 
          \*"",   f[\!\(p\_1\)] = "", 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"\_\(k + 1\), " = ", N[p2, digits], ",   f[", 
            "p"\_\(k + 1\), "] = ", N[Y2, digits], 
            ",   |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) - 3 x + 2.

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

    Plot[f[x], {x, \(-2.6\), 2.0}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use the secant method to find a numerical approximation to the root.

SecantMethod[f[x], -2.6, -2.4, 0.00000000001, 0.00000000001, 
      100, 8];

Compare with Mathematica's built in routine.

\(solset = FindRoot[f[x] == 0, {x, \(-2.4\)}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Mathematica's answer is very good.

SecantMethod[f[x], 1.4, 1.2, 0.00000000001, 0.00000000001, 100, 8];

Compare with Mathematica's built in routine.

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

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 20]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 30]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

Execute this cell to activate Steffensen's Acceleration.
Steffensen[f0_,x0_,delta_,epsilon_,max_,digits_]
f0_ the function f(x),
x0_ the starting value,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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[(1/(10^20)), digits];
		\(Print[\*""\!\(p\_0\) = "", N[p0, digits], 
                \*"",   f[\!\(p\_0\)] = "", 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], ",   f[p] = ", 
                      N[Y1, 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], ",   f[p] = ", 
                      N[Y2, 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[(Subscript["p", k + 1]), " = ", N[p3, digits], 
                      ",   f[", (Subscript["p", k + 1]), "] = ", 
                      N[Y3, digits], ",   |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) - 3 x + 2.

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

    Plot[f[x], {x, \(-2.6\), 2.0}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use Steffensen's acceleration method to find a numerical approximation to the root.

Steffensen[f[x], -2.4, 0.0000000001, 0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

\(solset = FindRoot[f[x] == 0, {x, \(-2.4\)}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Mathematica's answer is good.

Steffensen[f[x], 1.2, 0.0000000001, 0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

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

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 20]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 30]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Algorithm 2.8  Muller's Method.  To find a root of the equation  f (x) = 0  given three distinct initial approximation  (Subscript[p, 0]), (Subscript[p, 1]), (Subscript[p, 2]).  
Section  2.5   Aitken's, Steffensen's, Muller's Methods   Page  97

Execute this cell to activate Muller's Method.
Muller[f0_,x0_,x1_,x2_,delta_,epsilon_,max_,digits_]
f0_ the function f(x),
x0_x1_,x2_, the starting values,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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 = 2; 
		small = N[1\/10\^20, digits]; 
		
        Print[\*""\!\(p\_0\) = "", N[p0, digits], 
          \*"",   f[\!\(p\_0\)] = "", N[Y0, digits]]; 
		
        Print[\*""\!\(p\_1\) = "", N[p1, digits], 
          \*"",   f[\!\(p\_1\)] = "", N[Y1, digits]]; 
		
        Print[\*""\!\(p\_2\) = "", N[p2, digits], 
          \*"",   f[\!\(p\_2\)] = "", 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 = \@\(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"\_\(k + 1\), " = ", N[p2, digits], ",   f[", 
            "p"\_\(k + 1\), "] = ", N[Y2, digits], 
            ",   |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) - 3 x + 2.

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

    Plot[f[x], {x, \(-2.6\), 2.0}, PlotStyle -> RGBColor[1, 0, 0], 
      PlotLabel -> "The curve y = f(x)."]; 

    Print["The function is f[x] = ", f[x]]; \)

Use Muller's method to find a numerical approximation to the root.

Muller[f[x], -2.6, -2.5, -2.4, 0.0000000001, 
      0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

\(solset = FindRoot[f[x] == 0, {x, \(-2.4\)}]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Mathematica's answer is good.

Muller[f[x], 1.4, 1.3, 1.2, 0.0000000001, 0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

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

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

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

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 20]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

\(solset\  = \ 
      FindRoot[f[x] == 0, {x, 1.2}, MaxIterations -> 50, 
        WorkingPrecision -> 30]; 

    p = solset[[1, 2]]; 
Print[solset]; 
    
Print["   p = ", p]; 
Print["f[p] = ", f[p]]; \)

Algorithm 2.9  Nonlinear Seidel Iteration.  
	In two dimensions, solve the non-linear fixed point system 
			\(x\  = \ \(g\_1\) \((x, y)\), 
y\  = \ \(g\_2\) \((x, y)\), \)	
given one initial approximation  (Subscript[P, 0])=((Subscript[p, 0]), (Subscript[q, 0])),  and generating a sequence  {(Subscript[P, k])}={((Subscript[p, k]), (Subscript[q, k]))}  which converges to the solution  P=(p, q),  i.e.   p = (Subscript[g, 1]) (p, q), 
      q = (Subscript[g, 2]) (p, q)   and   r = (Subscript[g, 3]) (p, q).

	In three dimensions, solve the non-linear fixed point system 
			\(x\  = \ \(g\_1\) \((x, y, z)\), 
y\  = \ \(g\_2\) \((x, y, z)\), 

      z\  = \ \(g\_3\) \((x, y, z)\), \)	
given one initial approximation  (Subscript[P, 0])=((Subscript[p, 0]), (Subscript[q, 0]), (Subscript[r, 0])),  and generating a sequence  {(Subscript[P, k])}={((Subscript[p, k]), (Subscript[q, k]), (Subscript[r, k]))}  which converges to the solution  P=(p, q, r),  i.e.   p = (Subscript[g, 1]) (p, q, r), 
      q = (Subscript[g, 2]) (p, q, r)   and   r = 
          (Subscript[g, 3]) (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:

Execute this cell to activate Fixed Point Iteration a Non-Linear System.
                                  and        Nonlinear Seidel Iteration.
FixPointIter[G0_,X0_,tolerance_,max_,digits_]
Seidel[g10_,g20_,X0_,tolerance_,max_,digits_]
G0_ the function G(X),
X0_ the starting value,
tolerance_ estimate for the accuracy of the solution,
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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[\@\(V\^2 /. List -> Plus\), digits]; 
		
        P1 = N[X0, digits]; 
		Print[\*""\!\(P\_0\) = "", 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"\_k, " = ", N[P1, digits], ",   |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[\@\(V\^2 /. List -> Plus\), digits]; 
		k = 0; 
        
		dist = 1; 
		P1 = N[X0, digits]; 
		
        Print[\*""\!\(P\_0\) = "", 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"\_k, " = ", N[P1, digits], 
            ",   |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\  = \ \(\(g\_1\) \((x, y)\) = 1\/2\ \((x\^2 - y + 1\/2)\)\), 

      y\  = \ \(\(g\_2\) \((x, y)\) = 
          1\/8\ \(\((\(-x\^2\) - 4\ y\^2 + 8\ y + 4)\).\)\)\)

\(Clear[g1, g2, graph1, graph2, G, x, y]; 

    \(g1[{x_, y_}] = 1\/2\ \((x\^2 - y + 1\/2)\); 

    g2[{x_, y_}] = 1\/8\ \((\(-x\^2\) - 4\ y\^2 + 8\ y + 4)\); 

    \(G[{x_, y_}] = {g1[{x, y}], g2[{x, y}]}; 

    Print[\*""x = \!\(g\_1\)(x,y) = "", g1[{x, y}]]; 

    Print[\*""y = \!\(g\_2\)(x,y) = "", g2[{x, y}]]; \)\)\)

\(graph1 = 
      ImplicitPlot[g1[{x, y}] == x, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{\(-2.1\), 2.1}, {\(-1.1\), 1.1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    graph2 = ImplicitPlot[g2[{x, y}] == y, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        PlotRange -> {{\(-2\), 2}, {\(-1\), 1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    Show[graph1, graph2, DisplayFunction -> $DisplayFunction]; 

    Print["x = ", g1[{x, y}], "   and   ", "y = ", g2[{x, y}]]; 
    \)

Use fixed point iteration to find a numerical approximation.

FixPointIter[G[{x, y}], {0.0, 1.0}, 0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

\(solset\  = \ 
      FindRoot[{g1[{x, y}] == x, g2[{x, y}] == y}, {x, 0}, {y, 1}]; 

    p = solset[[1, 2]]; 

    q = solset[[2, 2]]; 
Print[solset]; 
    
Print["       P = {", p, ", ", q, "}"]; 

    Print["G[{p,q}] = ", G[{p, q}]]; \)

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

Use fixed point iteration to find a numerical approximation to the other root
Observe that the iteration diverges in this case !

FixPointIter[G[{x, y}], {2.0, 0.0}, 0.0000000001, 10, 8];

Compare with Mathematica's built in routine.

\(solset\  = \ 
      FindRoot[{g1[{x, y}] == x, g2[{x, y}] == y}, {x, 2}, {y, 0}]; 

    p = solset[[1, 2]]; 

    q = solset[[2, 2]]; 
Print[solset]; 
    
Print["       P = {", p, ", ", q, "}"]; 

    Print["G[{p,q}] = ", G[{p, q}]]; \)

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\  = \ \(\(g\_1\) \((x, y)\)\  = \ \ 1\/2\ \((x\^2 - y + 1\/2)\)\), 
      
y\  = \ \(\(g\_2\) \((x, y)\)\  = \ \ 
          1\/8\ \(\((\(-x\^2\) - 4\ y\^2 + 8\ y + 4)\).\)\)\)

\(Clear[g1, g2, graph1, graph2, G, x, y]; 

    g1[{x_, y_}] = 1\/2\ \((x\^2 - y + 1\/2)\); 

    g2[{x_, y_}] = 1\/8\ \((\(-x\^2\) - 4\ y\^2 + 8\ y + 4)\); 

    G[{x_, y_}] = {g1[{x, y}], g2[{x, y}]}; 

    Print[\*""x = \!\(g\_1\)(x,y) = "", g1[{x, y}]]; 

    Print[\*""y = \!\(g\_2\)(x,y) = "", g2[{x, y}]]; \)

\(graph1 = 
      ImplicitPlot[g1[{x, y}] == x, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{\(-2.1\), 2.1}, {\(-1.1\), 1.1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    graph2 = ImplicitPlot[g2[{x, y}] == y, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        PlotRange -> {{\(-2\), 2}, {\(-1\), 1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    Show[graph1, graph2, DisplayFunction -> $DisplayFunction]; 

    Print["x = ", g1[{x, y}], "   and   ", "y = ", g2[{x, y}]]; 
    \)

Use Seidel iteration to find a numerical approximation.
Observe that the iteration diverges in this case !

Seidel[g1[{x, y}], g2[{x, y}], {0.0, 1.0}, 0.0000000001, 100, 8];

Compare with Mathematica's built in routine.

\(solset\  = \ 
      FindRoot[{g1[{x, y}] == x, g2[{x, y}] == y}, {x, 0}, {y, 1}]; 

    p = solset[[1, 2]]; 

    q = solset[[2, 2]]; 
Print[solset]; 
    
Print["       P = {", p, ", ", q, "}"]; 

    Print["G[{p,q}] = ", G[{p, q}]]; \)

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

Use Seidel iteration to find a numerical approximation to the other root.
Observe that the iteration diverges in this case !

Seidel[g1[{x, y}], g2[{x, y}], {2.0, 0.0}, 0.0000000001, 10, 8];

Compare with Mathematica's built in routine.

\(solset\  = \ 
      FindRoot[{g1[{x, y}] == x, g2[{x, y}] == y}, {x, 2}, {y, 0}]; 

    p = solset[[1, 2]]; 

    q = solset[[2, 2]]; 
Print[solset]; 
    
Print["       P = {", p, ", ", q, "}"]; 

    Print["G[{p,q}] = ", G[{p, q}]]; \)

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\  = \ \(f\_1\) \((x, y)\), 
0\  = \ \(f\_2\) \((x, y)\), \)
given one initial approximation  ((Subscript[p, 0]), (Subscript[q, 0]))  and using Newton-Raphson iteration.
Section   2.7   Newton's Method for Systems   Page  116

Execute this cell to activate Newton-Raphson Method in 2-Dimensions.
NewtonSystem[F00_,X0_,delta_,epsilon_,max_,digits_]
F00_ the function F(X),
X0_ the starting value,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output

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[\@\(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\_0\)  = "", N[P0, digits]]; 
		
        Print[\*""F[\!\(P\_0\)] = "", 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["  ", "P"\_k, "  = ", N[P1, digits], 
            ",   |dP| <= ", N[norm[Dp], digits]]; 
			
          Print["F[", "P"\_k, "] = ", 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\  = \ \(\(f\_1\) \((x, y)\)\  = \ x\^2 - 2\ x - y + 0.5\), 

      0\  = \ \(\(f\_2\) \((x, y)\)\  = \ x\^2 + 4\ y\^2 - 4.\)\)

\(Clear[f1, f2, graph1, graph2, F, jacobian, x, y]; 

    \(f1[{x_, y_}] = x\^2 - 2\ x - y + 1\/2; 

    f2[{x_, y_}] = x\^2 + 4\ y\^2 - 4; 

    F[{x_, y_}] = {f1[{x, y}], f2[{x, y}]}; 

    jacobian[{x_, y_}] = 
      Transpose[{∂\_x\ F[{x, y}], ∂\_y\ F[{x, y}]}]; 

    Print[\*""0 = \!\(f\_1\)(x,y) = "", f1[{x, y}]]; 

    Print[\*""0 = \!\(f\_2\)(x,y) = "", f2[{x, y}]]; 

    Print["   Jacobian = ", MatrixForm[jacobian[{x, y}]]]; \)\)

\(graph1 = 
      ImplicitPlot[f1[{x, y}] == 0, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotRange -> {{\(-2.1\), 2.1}, {\(-1.1\), 1.1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    graph2 = ImplicitPlot[f2[{x, y}] == 0, \ {x, \ \(-2\), \ 2}, 
        PlotStyle -> RGBColor[0, 1, 0], 
        PlotRange -> {{\(-2\), 2}, {\(-1\), 1}}, 
        AspectRatio -> 1\/2, DisplayFunction -> Identity]; 

    Show[graph1, graph2, DisplayFunction -> $DisplayFunction]; 

    Print["x = ", g1[{x, y}], "   and   ", "y = ", g2[{x, y}]]; 
    \)

Use the Newton-Raphson method to find a numerical approximation to the root.

NewtonSystem[F[{x, y}], {2.0, 0.25}, 0.0000000001, 0.0000000001, 100, 
      8];

Compare with Mathematica's built in routine.

\(solset\  = \ 
      FindRoot[{f1[{x, y}] == 0, f2[{x, y}] == 0}, {x, 2}, {y, 0.25}]; 

    p = solset[[1, 2]]; 

    q = solset[[2, 2]]; 
Print[solset]; 
    
Print["   P = {", p, ", ", q, "}"]; 

    Print["F[{p,q}] = ", F[{p, q}]]; \)