Numerical Analysis 6

Spring 2005

Bracketing Methods for Finding Roots of f(x)=0  (Section 2.2)

Bisection Method of Bolzano

Bisection[f0_,a0_,b0_,delta_,digits_,{leftx_,rightx_},{lowery_,uppery_}]
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.
leftx_ is the leftmost x-value displayed on the graph
rightx_  is the rightmost x-value displayed on the graph
lowery_ is the minimum y-value displayed on the graph
uppery_ is the maximum y-value displayed on the graph

Off[General::"spell1"];
Clear[Bisection, f0, a0, b0, delta, x, digits, leftx, rightx, lowery, 
        uppery];
\(\(Bisection[f0_, a0_, b0_, delta_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}] := 
        Module[{a, b, c, Ya, Yb, Yc, f, k, cond, MaxIter, MethodDoesNotApply, 
            PerfectZero, RootInTolerance, clist, ptlist, gr1, gr2, 
            gr3}, 
		Set @@ {f[x_], f0}; 
		a = 
            N[a0, digits]; 
		b = 
            N[b0, digits]; 
clist = {a0, b0}; 
		k = 
            0; 
		Ya = f[a]; 
		Yb = f[b]; 
		dx = 
            b - a; 
		cond = 0; 
gr1\  = 
            Plot[f[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              DisplayFunction -> Identity]; 
		MaxIter = 
            N[1 + Floor[\(Log[dx] - Log[delta]\)\/Log[2.0]]]; 
          If[Ya\ Yb > 0, cond = MethodDoesNotApply; 
            Print["Starting Points Don't Have Opposite Signs"]; 
            Break[];]; 
		Print[" a = ", 
            N[a, digits]]; 
		Print[" b = ", N[b, digits]]; 
          While[cond == 0 && k < MaxIter, 
            c = \(a + b\)\/2; 
			Yc = f[c]; 
clist = 
              Append[clist, c]; 
ptlist = 
              Transpose[{clist, Table[0, {Length[clist]}]}]; 
			gr2\  = 
              Show[Graphics[
                  Map[{PointSize[0.02], Point[#]} &, 
                    Take[ptlist, {Length[ptlist] - 2, Length[ptlist] - 1}]]], 
                DisplayFunction -> Identity]; 
            gr3 = Show[
                Graphics[
                  Map[{PointSize[0.03], RGBColor[1, 0, 0], Point[#]} &, 
                    Take[ptlist, \(-1\)]]], 
                DisplayFunction -> Identity]; 
            Show[{gr1, gr2, gr3}, DisplayFunction -> $DisplayFunction]; 
            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["List of approximations: ", 
            clist]; 
		Print["The function is f[x] = ", 
            f[
              x]]; 
		Print["An approximate root of f[x] is:"]; 
		Print["  c  = ", N[c, digits]]; 
		Print["f[c] = ", 
            N[f[c], digits]]; 
          Print["Error is |dc| ≤ ", 
            N[dx, digits]];];\)\)
On[General::"spell1"];

Bisection[1 - x^2, -4,  .5,  .001, 
      5, {-5, 6}, {-10, 10}]

Method of False Position

The method of false position differs from the bisection method in that a new point at each stage of the process is generated by finding the intersection of the secant line joining the last two points and the x-axis.  The next function is a modification of your author's FalsePosition function.

FalsePosition[f0_,a0_,b0_,delta_,epsilon_,digits_,{leftx_,rightx_},{lowery_,uppery_}]
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.
leftx_ is the leftmost x-value displayed on the graph
rightx_  is the rightmost x-value displayed on the graph
lowery_ is the minimum y-value displayed on the graph
uppery_ is the maximum y-value displayed on the graph

Off[General::"spell1"];
Clear[FalsePosition, f0, a0, b0, delta, epsilon, x];
\(\(FalsePosition[f0_, a0_, b0_, delta_, epsilon_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}] := 
        Module[{a, b, c, Ya, Yb, Yc, f, k, cond, Iterating, RootInTolerance, 
            PerfectZero, MethodDoesNotApply, ptlist, nextpt, gr1, gr2, gr3, 
            gr4, gr5, gr6, gr7, clist}, 
		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; 
clist = {}; 
		If[
            Ya\ Yb > 0, cond = MethodDoesNotApply; 
            Print["Starting Points Don't Have Opposite Signs"]; 
            Break[];]; 
          Print[" a = ", N[a, digits]]; 
		Print[" b = ", 
            N[b, digits]]; 
	
          ptlist = {{a0, f[a0]}, {b0, f[b0]}}; 
nextpt = 
            a0 - \((f[a0] \((b0 - a0)\))\)/\((f[b0] - 
                    f[a0])\); 
	
          gr1\  = Plot[f[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              DisplayFunction -> Identity]; 
	
          gr2\  = Show[
              Graphics[
                Map[{PointSize[0.03], Point[#]} &, {{a, f[a]}, {b, f[b]}}]], 
              DisplayFunction -> Identity]; 
	
          gr3\  = Show[Graphics[Line[ptlist]], 
              DisplayFunction -> Identity]; 
gr4 = 
            Show[Graphics[
                Map[{PointSize[0.02], RGBColor[1, 0, 0], 
                      Point[#]} &, {{nextpt, 0}}]], 
              DisplayFunction -> Identity]; 
	
          Show[{gr1, gr2, gr3, gr4}, 
            DisplayFunction -> $DisplayFunction]; 
		While[cond == 0 && 
              k < MaxIter, 
			dx = \(Yb\ \((b - a)\)\)\/\(Yb - Ya\); 
			c = b - dx; 
			Yc = f[c]; 
		clist = 
              Append[clist, 
                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; ptlist = {{a, f[a]}, {b, f[b]}}, 
              a = c; Ya = Yc; 
              ptlist = {{a, f[a]}, {b, f[b]}}]; 
nextpt = 
              a - \((f[a] \((b - a)\))\)/\((f[b] - 
                      f[a])\); 
	
            gr5\  = Show[Graphics[Line[ptlist]], 
                DisplayFunction -> Identity]; 
            gr6 = Show[Graphics[Map[{PointSize[0.03], Point[#]} &, ptlist]], 
                DisplayFunction -> Identity]; 
gr7 = 
              Show[Graphics[
                  Map[{PointSize[0.02], RGBColor[1, 0, 0], 
                        Point[#]} &, {{nextpt, 0}}]], 
                DisplayFunction -> 
                  Identity]; 
Show[{gr1, gr5, gr6, gr7}, 
              DisplayFunction -> $DisplayFunction]; 
If[Abs[dx] < delta && Abs[Yc] < epsilon, 
              cond = RootInTolerance]; k = k + 1;]; 
          If[cond == 0, 
            cond = Iterating]; 
		Print[" "]; 
	
          Print["List of c's"]; 
	
          Print[clist]; 
		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"];

FalsePosition[Log[x] - 5 + x, 3.2, 4.0,  .001,  .001, 
      15, {3, 5}, {-3, 3}]