Open Educational Resources

Finding Roots of Equations: Bracketing Methods

An elementary observation from the previous method is that in most cases, the function changes signs around the roots of an equation. The bracketing methods rely on the intermediate value theorem.

Intermediate Value Theorem

Statement: Let f:[a,b]\rightarrow \mathbb{R} be continuous and f(a)\leq f(b). Then, \forall y\in[f(a),f(b)]:\exists x\in[a,b] such that y=f(x). The same applies if f(b)\leq f(a).
The proof for the intermediate value theorem relies on the assumption of continuity of the function f. Intuitively, because the function f is continuous, then for the values of f to change from f(a) to f(b) it has to pass by every possible value between f(a) and f(b). The consequence of the theorem is that if the function f is such that f(a)\leq 0 \leq f(b), then, there is x\in[a,b] such that f(x)=0. We can proceed to find x iteratively using either the bisection method or the false position method.

Bisection Method

In the bisection method, if f(a)f(b)<0, an estimate for the root of the equation f(x)=0 can be found as the average of a and b:

    \[x_i=\frac{a+b}{2}\]

Upon evaluating f(x_i), the next iteration would be to set either a=x_i or b=x_i such that for the next iteration the root x_{i+1} is between a and b. The following describes an algorithm for the bisection method given a<b, f(x), \varepsilon_s, and maximum number of iterations:
Step 1: Evaluate f(a) and f(b) to ensure that f(a)f(b)<0. Otherwise, exit with an error.
Step 2: Calculate the value of the root in iteration i as x_i=\frac{a_i+b_i}{2}. Check which of the following applies:

  1. If f(x_i)=0, then the root has been found, the value of the error \varepsilon_r=0. Exit.
  2. If f(x_i)f(a_i)<0, then for the next iteration, x_{i+1} is bracketed between a_i and x_i. The value of \varepsilon_r=\frac{x_{i+1}-x_{i}}{x_{i+1}}.
  3. If f(x_i)f(b_i)<0, then for the next iteration, x_{i+1} is bracketed between x_i and b_i. The value of \varepsilon_r=\frac{x_{i+1}-x_{i}}{x_{i+1}}.

Step 3: Set i=i+1. If i reaches the maximum number of iterations or if \varepsilon_r\leq \varepsilon_s, then the iterations are stopped. Otherwise, return to step 2 with the new interval a_{i+1} and b_{i+1}.

Example

Setting \varepsilon_s=0.0005 and applying this process to f(x)=\sin{5x}+\cos{2x} with a=-0.6 and b=-0.5 yields the estimate x_9=-0.523633 after 9 iterations with \varepsilon_r=-0.000373 as shown below. Similarly, applying this process to f(x)=\sin{5x}+\cos{2x} with a=-0.3 and b=-0.2 yields the estimate x_{10}=-0.224316 after 10 iterations with \varepsilon_r=-0.000435 while applying this process to f(x)=\sin{5x}+\cos{2x} with a=0.6 and b=0.7 yields the estimate x_9=0.673242 after 9 iterations with \varepsilon_r=0.00029:

Bisection1
View Mathematica Code
Clear[f, x, ErrorTable, ei]
f[x_] := Sin[5 x] + Cos[2 x]
(*The following function returns x_i and the new interval (a,b) along with an error note. The order of the output is {xi,new a, new b, note}*)
bisec2[f_, a_, b_] := (xi = (a + b)/2; Which[
   f[a]*f[b] > 0, {xi, a, b, "Root Not Bracketed"},
   f[xi] == 0, {xi, a, b, "Root Found"},
   f[xi]*f[a] < 0, {xi, a, xi, "Root between a and xi"},
   f[xi]*f[b] < 0, {xi, xi, b, "Root between xi and b"}])

(*Problem Setup*)
MaxIter = 20;
eps = 0.0005;

(*First root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {-0.6};
btable = {-0.5};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = bisec2[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Second Root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {-0.3};
btable = {-0.2};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = bisec2[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Third Root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {0.6};
btable = {0.7};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = bisec2[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm
View Python Code
import numpy as np
import pandas as pd

def f(x): return np.sin(5*x) + np.cos(2*x)
def bisec2(f,a,b):
  xi = (a + b)/2
  if f(a)*f(b) > 0: return [xi, a, b, "Root Not Bracketed"]
  elif f(xi) == 0: return [xi, a, b, "Root Found"]
  elif f(xi)*f(a) < 0: return [xi, a, xi, "Root between a and xi"]
  elif f(xi)*f(b) < 0: return [xi, xi, b, "Root between xi and b"]

def bisection(ErrorTable,ErrorNoteTable,atable,btable,xtable):
  MaxIter = 20
  eps = 0.0005
  Stopcode = "NoStopping"
  i = 0
  while i <= MaxIter and Stopcode == "NoStopping":
    ri = bisec2(f, atable[i], btable[i])
    if ri[3] == "Root Not Bracketed": 
      xtable.append(ri[0])
      ErrorNoteTable.append(ri[3])
      break
    if ri[3] == "Root Found":
      xtable.append(ri[0])
      ErrorTable.append(0)
      ErrorNoteTable.append(ri[3])
      break
    atable.append(ri[1])
    btable.append(ri[2])
    xtable.append(ri[0])
    ErrorNoteTable.append(ri[3])

    if i != 0:
      ei = (xtable[i] - xtable[i - 1])/xtable[i]
      ErrorTable.append(ei)
    if abs(ErrorTable[i]) > eps:
      Stopcode = "NoStopping"
    else:
      Stopcode = "Stop"
    i+=1
  T2 = ([[i, atable[i], btable[i], xtable[i], f(atable[i]), 
        f(btable[i]), f(xtable[i]), ErrorTable[i], 
        ErrorNoteTable[i]] for i in range(len(xtable))])
  T2 = pd.DataFrame(T2, columns=["Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"])
  display(T2)

# First Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [-0.6]
btable = [-0.5]
xtable = []
print("First Root")
bisection(ErrorTable,ErrorNoteTable,atable,btable,xtable)

# Second Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [-0.3]
btable = [-0.2]
xtable = []
print("Second Root")
bisection(ErrorTable,ErrorNoteTable,atable,btable,xtable)

# Third Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [0.6]
btable = [0.7]
xtable = []
print("Third Root")
bisection(ErrorTable,ErrorNoteTable,atable,btable,xtable)

The following code can be used to implement the bisection method in MATLAB. For the sake of demonstration, it finds the roots of the function in above example. The example function is defined in another file

The following tool illustrates this process for a=0.1 and b=0.9. Use the slider to view the process after each iteration. In the first iteration, the interval [a,b]=[0.1,0.9]. f(a)f(b)<0, so, x_1=\frac{a+b}{2}=0.5. In the second iteration, the interval becomes [0.5,0.9] and the new estimate x_2=\frac{0.5+0.9}{2}=0.7. The relative approximate error in the estimate \varepsilon_r=\frac{0.7-0.5}{0.7}=0.2857. You can view the process to see how it converges after 12 iterations.




 

Error Estimation

If V_t is the true value of the root and V_a is the estimate, then, the number of iterations n needed to ensure that the absolute value of the error |E|=|V_t-V_a| is less than or equal to a certain value can be easily obtained. Let L=b-a be the length of the original interval used. The first estimate x_1=\frac{a+b}{2} and so, in the next iteration, the interval where the root is contained has a length of \frac{L}{2}. As the process evolves, the interval for the iteration number n has a length of \frac{L}{2^n}. Since the true value exists in an interval of length \frac{L}{2^n}, the absolute value of the error is such that:

    \[|E|\leq \frac{L}{2^n}\]

Therefore, for a desired estimate of the absolute value of the error, say E_{max}, the number of iterations required is:

    \[\frac{L}{2^n}=E_{max}\Rightarrow n=\frac{\ln{\frac{L}{E_{max}}}}{\ln{2}}\]

Example

As an example, consider f(x)=x^3+x^2-10, if we wish to find the root of the equation f(x)=0 in the interval [1,2] with an absolute error less than or equal to 0.004, the number of iterations required is 8:

    \[n=\frac{\ln{\frac{1}{0.004}}}{\ln{2}}= 7.966 \approx 8\]

The actual root with 10 significant digits is V_t=1.867460025.
Using the process above, after the first iteration, x_1=1.5 and so, the root lies between 1.5 and 2. So, the length of the interval is equal to 0.5 and the error in the estimate is less than 0.5. The length of the interval after iteration 8 is equal to 0.0039 and so the error in the estimate is less than 0.0039. It should be noted however that the actual error was less than this upper bound after the seventh iteration.

table111
View Mathematica Code
Clear[f, x, tt, ErrorTable, nn];
f[x_] := x^3 + x^2 - 10;
bisec[f_, a_, b_] := (xi = (a + b)/2; 
   Which[f[a]*f[b] > 0, {xi, a, b, "Root Not Bracketed"}, 
    f[xi] == 0, {xi, a, b, "Root Found"}, 
    f[xi]*f[a] < 0, {xi, a, xi, "Root between a and xi"}, 
    f[xi]*f[b] < 0, {xi, xi, b, "Root between xi and b"}]);

(*Exact Solution*)
t = Solve[x^3 + x^2 - 10 == 0, x]
Vt = N[x /. t[[1]], 10]

(*Problem Setup*)
MaxIter = 20;
eps = 0.0005;

(*First root*)
x = {bisec[f, 1., 2]};
ErrorTable = {1};
ErrorTable2 = {"N/A"};
IntervalLength = {0.5};
i = 2;
While[And[i <= MaxIter, Abs[ErrorTable[[i - 1]]] > eps], 
 ri = bisec[f, x[[i - 1, 2]], x[[i - 1, 3]]];
 If[ri == "Error, root is not bracketed", 
  ErrorTable[[1]] = "Error, root is not bracketed"; Break[]];
 x = Append[x, ri]; ei = (x[[i, 1]] - x[[i - 1, 1]])/x[[i, 1]]; 
 e2i = x[[i, 1]] - Vt; ErrorTable = Append[ErrorTable, ei]; 
 IntervalLength = 
  Append[IntervalLength, (x[[i - 1, 3]] - x[[i - 1, 2]])/2]; 
 ErrorTable2 = Append[ErrorTable2, e2i]; i++]
x // MatrixForm
ErrorTable // MatrixForm
Title = {"Iteration", "x_i", "a", "b", "e_r", "Actual E", 
   "Interval Length (Upper bound for E)"};
T2 = Table[{i, x[[i, 1]], x[[i, 2]], x[[i, 3]], ErrorTable[[i]], 
    ErrorTable2[[i]], IntervalLength[[i]]}, {i, 1, Length[x]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm
View Python Code
import sympy as sp
import pandas as pd

def f(x): return x**3 + x**2 - 10
def bisec(f,a,b):
  xi = (a + b)/2
  if f(a)*f(b) > 0: return [xi, a, b, "Root Not Bracketed"]
  elif f(xi) == 0: return [xi, a, b, "Root Found"]
  elif f(xi)*f(a) < 0: return [xi, a, xi, "Root between a and xi"]
  elif f(xi)*f(b) < 0: return [xi, xi, b, "Root between xi and b"]
    
# Exact Solution
x = sp.symbols('x')
t = list(sp.nonlinsolve([x**3 + x**2 - 10], [x]))
Vt = sp.N(t[0][0])
print("Vt:",Vt)

# Problem Setup
MaxIter = 20
eps = 0.0005

# First Root
x = [bisec(f, 1., 2)]
ErrorTable = [1]
ErrorTable2 = ["N/A"]
IntervalLength = [0.5]
i = 1

while i <= MaxIter and abs(ErrorTable[i - 1]) > eps:
  ri = bisec(f, x[i - 1][1], x[i - 1][2])
  if ri=="Error, root is not bracketed":
    ErrorTable[0]="Error, root is not bracketed"
    break
  x.append(ri)
  ei = (x[i][0] - x[i - 1][0])/x[i][0]
  e2i = x[i][0] - Vt
  ErrorTable.append(ei)
  IntervalLength.append((x[i - 1][2] - x[i - 1][1])/2)
  ErrorTable2.append(e2i)
  i+=1

print("x Table"); 
display(pd.DataFrame(x))
print("ErrorTable")
display(pd.DataFrame(ErrorTable))
T2 = ([[i, x[i][0], x[i][1], x[i][2], ErrorTable[i], ErrorTable2[i], 
      IntervalLength[i]] for i in range(len(x))])
pd.DataFrame(T2, columns=["Iteration", "x_i", "a", "b", "e_r", "Actual E", "Interval Length (Upper bound for E)"])

The procedure for solving the above example in MATLAB is available in the following files. The polynomial function is defined in a separate file.

False Position Method

In the false position method, the new estimate x_i at iteration i is obtained by considering the linear function passing through the two points (a,f(a)) and (b,f(b)). The point of intersection of this line with the x axis can be obtained using one of the following formulas:

    \[x_i=a-f(a)\frac{(b-a)}{f(b)-f(a)}=b-f(b)\frac{(b-a)}{f(b)-f(a)}=\frac{af(b)-bf(a)}{f(b)-f(a)}\]

Upon evaluating f(x_i), the next iteration would be to set either a=x_i or b=x_i such that for the next iteration the root x_{i+1} is between a and b. The following describes an algorithm for the false position method method given a<b, f(x), \varepsilon_s, and maximum number of iterations:
Step 1: Evaluate f(a) and f(b) to ensure that f(a)f(b)<0. Otherwise exit with an error.
Step 2: Calculate the value of the root in iteration i as x_i=\frac{a_if(b_i)-bf(a_i)}{f(b_i)-f(a_i)}. Check which of the following applies:

  1. If f(x_i)=0, then the root has been found, the value of the error \varepsilon_r=0. Exit.
  2. If f(x_i)f(a_i)<0, then for the next iteration, x_{i+1} is bracketed between a_i and x_i. The value of \varepsilon_r=\frac{x_{i+1}-x_{i}}{x_{i+1}}.
  3. If f(x_i)f(b_i)<0, then for the next iteration, x_{i+1} is bracketed between x_i and b_i. The value of \varepsilon_r=\frac{x_{i+1}-x_{i}}{x_{i+1}}.

Step 3: If i reaches the maximum number of iterations or if \varepsilon_r\leq \varepsilon_s, then the iterations are stopped. Otherwise, return to step 2.

Example

Setting \varepsilon_s=0.0005 and applying this process to f(x)=\sin{5x}+\cos{2x} with a=-0.6 and b=-0.5 yields the estimate x_3=-0.523569 after 3 iterations with \varepsilon_r=0.000498 as shown below. Similarly, applying this process to f(x)=\sin{5x}+\cos{2x} with a=-0.3 and b=-0.2 yields the estimate x_4=-0.2244 after 4 iterations with \varepsilon_r=-0.00015 while applying this process to f(x)=\sin{5x}+\cos{2x} with a=0.6 and b=0.7 yields the estimate x_3=0.673198 after 3 iterations with \varepsilon_r=-4.4\times 10^{-6}:

falseposition1
View Mathematica Code
Clear[f, x, ErrorTable, ei]
f[x_] := Sin[5 x] + Cos[2 x]
(*The following function returns x_i and the new interval (a,b) along with an error note. The order of the output is {xi,new a, new b, note}*)
falseposition[f_, a_, b_] := (xi = (a*f[b]-b*f[a])/(f[b]-f[a]); Which[
   f[a]*f[b] > 0, {xi, a, b, "Root Not Bracketed"},
   f[xi] == 0, {xi, a, b, "Root Found"},
   f[xi]*f[a] < 0, {xi, a, xi, "Root between a and xi"},
   f[xi]*f[b] < 0, {xi, xi, b, "Root between xi and b"}])

(*Problem Setup*)
MaxIter = 20;
eps = 0.0005;

(*First root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {-0.6};
btable = {-0.5};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Second Root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {-0.3};
btable = {-0.2};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Third Root*)
ErrorTable = {1};
ErrorNoteTable = {};
atable = {0.6};
btable = {0.7};
xtable = {};
i = 1;
Stopcode = "NoStopping";

While[And[i <= MaxIter, Stopcode == "NoStopping"],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri[[4]] == "Root Not Bracketed", xtable = Append[xtable, ri[[1]]];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 If[ri[[4]] == "Root Found", xtable = Append[xtable, ri[[1]]]; ErrorTable = Append[ErrorTable, 0];ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]]; Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 xtable = Append[xtable, ri[[1]]];
 ErrorNoteTable = Append[ErrorNoteTable, ri[[4]]];
 If[i != 1,
  ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];
  ErrorTable = Append[ErrorTable, ei]];
 If[Abs[ErrorTable[[i]]] > eps, Stopcode = "NoStopping", Stopcode = "Stop"];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]], ErrorNoteTable[[i]]}, {i, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm
View Python Code
import numpy as np
import pandas as pd

def f(x): return np.sin(5*x) + np.cos(2*x)
# The following function returns x_i and the new interval (a,b) along with an error note. The order of the output is {xi,new a, new b, note}
def falseposition(f,a,b):
  xi = (a*f(b)-b*f(a))/(f(b)-f(a))
  if f(a)*f(b) > 0: return [xi, a, b, "Root Not Bracketed"]
  elif f(xi) == 0: return [xi, a, b, "Root Found"]
  elif f(xi)*f(a) < 0: return [xi, a, xi, "Root between a and xi"]
  elif f(xi)*f(b) < 0: return [xi, xi, b, "Root between xi and b"]

def falseposition2(ErrorTable,ErrorNoteTable,atable,btable,xtable):
  MaxIter = 20
  eps = 0.0005
  Stopcode = "NoStopping"
  i = 0
  while i <= MaxIter and Stopcode == "NoStopping":
    ri = falseposition(f, atable[i], btable[i])
    if ri[3] == "Root Not Bracketed": 
      xtable.append(ri[0])
      ErrorNoteTable.append(ri[3])
      break
    if ri[3] == "Root Found":
      xtable.append(ri[0])
      ErrorTable.append(0)
      ErrorNoteTable.append(ri[3])
      break
    atable.append(ri[1])
    btable.append(ri[2])
    xtable.append(ri[0])
    ErrorNoteTable.append(ri[3])

    if i != 0:
      ei = (xtable[i] - xtable[i - 1])/xtable[i]
      ErrorTable.append(ei)
    if abs(ErrorTable[i]) > eps:
      Stopcode = "NoStopping"
    else:
      Stopcode = "Stop"
    i+=1
  T2 = ([[i, atable[i], btable[i], xtable[i], f(atable[i]), 
        f(btable[i]), f(xtable[i]), ErrorTable[i], 
        ErrorNoteTable[i]] for i in range(len(xtable))])
  T2 = pd.DataFrame(T2, columns=["Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er", "ErrorNote"])
  display(T2)

# First Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [-0.6]
btable = [-0.5]
xtable = []
print("First Root")
falseposition2(ErrorTable,ErrorNoteTable,atable,btable,xtable)

# Second Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [-0.3]
btable = [-0.2]
xtable = []
print("Second Root")
falseposition2(ErrorTable,ErrorNoteTable,atable,btable,xtable)

# Third Root
ErrorTable = [1]
ErrorNoteTable = []
atable = [0.6]
btable = [0.7]
xtable = []
print("Third Root")
falseposition2(ErrorTable,ErrorNoteTable,atable,btable,xtable)

The procedure for implementing the false position root finding algorithm in MATLAB is available in the following files. The example function is defined in a separate file.

The following tool illustrates this process for a=0.1 and b=0.9. Use the slider to view the process after each iteration. In the first iteration, the interval [a,b]=[0.1,0.9]. f(a)f(b)<0, so, x_1=\frac{0.1f(0.9)-0.9f(0.1)}{f(0.9)-f(0.1)}=0.538249. In the second iteration, the interval becomes [0.538249,0.9] and the new estimate x_2=\frac{0.538249f(0.9)-0.9f(0.538249)}{f(0.9)-f(0.538249)}=0.693886. The relative approximate error in the estimate \varepsilon_r=\frac{0.693886-0.538249}{0.693886}=0.224297. You can view the process to see how it converges after very few iterations.




 

Lecture video


Lecture video


Leave a Reply

Your email address will not be published. Required fields are marked *