Samer Adeeb

Introduction to Numerical Analysis: Finding Roots of Equations

Let f:[a,b]\rightarrow \mathbb{R} be a function. One of the most basic problems in numerical analysis is to find the value x\in[a,b] that would render f(x)=0. This value x is called a root of the equation f. As an example, consider the function f(x)=2x-20. Consider the equation f(x)=0. Clearly, there is only one root to this equation which is x=10. Alternatively, consider g:\mathbb{R}\rightarrow\mathbb{R} defined as g(x)=a\times x^2+b\times x+c with a,b,c\in\mathbb{R}. There are two analytical roots for the equation g(x)=0 given by

    \[ x=\frac{-b\pm \sqrt{b^2-4ac}}{2a} \]

In this section we are going to present three classes of methods: graphical, bracketing, and open methods for finding roots of equations.

Graphical Methods

Graphical methods rely on a computational device that calculates the values of the function along an interval at specific steps, and then draws the graph of the function. By visual inspection, the analyst can identify the points at which the function crosses the x axis. Graphical methods are very useful for one-dimensional problems, but for multi-dimensional problems it is not possible to visualize a graph of a function of multi-variables.

Example

As an example, consider f(x)=\sin{5x}+\cos{2x}. We wish to find the values of x (i.e., the roots) that would satisfy the equation f(x)=0 for x\in[-1,1]. The fastest way is to plot the graph (even Microsoft Excel would serve as an appropriate tool here) and then visually inspect the locations where the graph crosses the x axis:
a
As shown in the above graph, the equation has three roots and just by visual inspection, the roots are around x_1=-0.525, x_2=-0.21, and x_3=0.675. Of course these estimates are not that accurate for f(x_1)=0.00365, f(x_2)=0.0457, and f(x_3)=-0.012. Whether these estimates are acceptable or not, depends on the requirements of the problem. We can further zoom in using the software to try and get better estimates:
b
Visual inspection of the zoomed-in graphs provides the following estimates: x_1=-0.5235, x_2=-0.225, and x_3=0.6731 which are much better estimates since: f(x_1)=0.00025, f(x_2)=-0.0018, and f(x_3)=0.00067
View Mathematica Code

Clear[x]
f = Sin[5 x] + Cos[2 x]
Plot[f, {x, -1, 1}, AxesLabel -> {"x", "f"}]
f /. x -> -0.525
f /. x -> -0.21
f /. x -> 0.675
Plot[f, {x, -0.53, -0.52}, AxesLabel -> {"x", "f"}]
Plot[f, {x, -0.25, -0.2}, AxesLabel -> {"x", "f"}]
Plot[f, {x, 0.66, 0.675}, AxesLabel -> {"x", "f"}]
f /. x -> -0.5235
f /. x -> -0.225
f /. x -> 0.6731

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. Set x_1=\frac{a+b}{2}.
Step 2: Set i=i+1. If f(x_i)f(a)\leq 0, then set the new interval as a=a, b=x_i. Otherwise, if f(x_i)f(b)\leq 0, then set the new interval as a=x_i, b=b. A new estimate for the root can be evaluated as x_{i+1}=\frac{a+b}{2}. 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_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). Otherwise, it returns an error: "Error, root is not bracketed"*)
bisec[f_, a_, b_] :=
If[f[a] f[b] <= 0, If[f[(a + b)/2] f[a] <= 0, {(a + b)/2, a, (a + b)/2}, If[f[(a + b)/2] f[b] < 0, {(a + b)/2, (a + b)/2, b}]], 
"Error, root is not bracketed"]
(*Problem Setup*)
MaxIter = 20;
eps = 0.0005;

(*First root*)
ErrorTable = {1};
atable = {-0.6};
btable = {-0.5};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]", "f[x_i]"};
While[
 And[i <= MaxIter, Flag == 1],
 ri = bisec[f, atable[[i]], btable[[i]]];
If[ri == "Error, root is not bracketed", ErrorTable[[1]] = ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; 
   ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]]; 
   ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1,Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Second Root*)
ErrorTable = {1};
atable = {-0.3};
btable = {-0.2};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]", "f[x_i]"};
While[
 And[i <= MaxIter, Flag == 1],
 ri = bisec[f, atable[[i]], btable[[i]]];
If[ri == "Error, root is not bracketed", ErrorTable[[1]] = ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];  ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1,  Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Third Root*)

ErrorTable = {1};
atable = {0.6};
btable = {0.7};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]","f[x_i]"};

While[
 And[i <= MaxIter, Flag == 1],
 ri = bisec[f, atable[[i]], btable[[i]]];
If[ri == "Error, root is not bracketed", ErrorTable[[1]] = ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]];  ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1,Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

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_] := If[f[a] f[b] <= 0, If[f[(a + b)/2] f[a] <= 0, {(a + b)/2, a, (a + b)/2}, If[f[(a + b)/2] f[b] <= 0, {(a + b)/2, (a + b)/2, b}]], "Error, root is not bracketed"];
(*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

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. Set x_1=\frac{af(b)-bf(a)}{f(b)-f(a)}.
Step 2: Set i=i+1. If f(x_i)f(a)\leq 0, then set the new interval as a=a, b=x_i. Otherwise, if f(x_i)f(b)\leq 0, then set the new interval as a=x_i, b=b. A new estimate for the root can be evaluated as x_{i+1}=\frac{af(b)-bf(a)}{f(b)-f(a)}. 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). Otherwise, it returns: "Error, root is not bracketed"*)
falseposition[f_, a_, b_] := ((xr = (a*f[b] - b*f[a])/(f[b] - f[a])); If[f[a] f[b] <= 0, If[f[xr] f[a] <= 0, {xr, a, xr}, If[f[xr] f[b] <= 0, {xr, xr, b}]], "Error, root is not bracketed"]);
(*Problem Setup*)
MaxIter = 20;
eps = 0.0005;

(*First root*)
ErrorTable = {1};
atable = {-0.6};
btable = {-0.5};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]", "f[x_i]"};

While[
 And[i <= MaxIter, Flag == 1],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri=="Error, root is not bracketed", ErrorTable[[1]]=ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]]; ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1,Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Second Root*)
ErrorTable = {1};
atable = {-0.3};
btable = {-0.2};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]", "f[x_i]"};

While[
 And[i <= MaxIter, Flag == 1],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri=="Error, root is not bracketed", ErrorTable[[1]]=ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]]; ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm

(*Third Root*)
ErrorTable = {1};
atable = {0.6};
btable = {0.7};
xtable = {0};
i = 1;
Flag = 1;
Title = {"Iteration", "a", "b", "x_i", "e_r", "f[a]", "f[b]", "f[x_i]"};

While[
 And[i <= MaxIter, Flag == 1],
 ri = falseposition[f, atable[[i]], btable[[i]]];
 If[ri=="Error, root is not bracketed", ErrorTable[[1]]=ri;Break[]];
 atable = Append[atable, ri[[2]]];
 btable = Append[btable, ri[[3]]];
 If[i == 1, 
  xtable[[1]] = ri[[1]], 
(xtable = Append[xtable, ri[[1]]]; ei = (xtable[[i]] - xtable[[i - 1]])/xtable[[i]]; ErrorTable = Append[ErrorTable, ei])];
 If[Abs[ErrorTable[[i]]] > eps, Flag = 1, Flag = 0];
 i++]
Title = {"Iteration", "a", "b", "x_i", "f[a]", "f[b]", "f[x_i]", "er"};
T2 = Table[{i, atable[[i]], btable[[i]], xtable[[i]], f[atable[[i]]], f[btable[[i]]], f[xtable[[i]]], ErrorTable[[i]]}, {i, 1, Length[xtable]}];
T2 = Prepend[T2, Title];
T2 // MatrixForm




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.

 

Open Methods

Open methods do not rely on having the root squeezed between two values, but rather rely on an initial guess and then apply an iterative process to get better estimates for the root. Open methods are usually faster in convergence if they converge, but they don’t always converge. For multi-dimensions, i.e., for solving multiple nonlinear equations, open methods are easier to implement than bracketing methods which cannot be easily extended to multi-dimensions.

Fixed-Point Iteration Method

The fixed-point iteration method relies on replacing the expression f(x)=0 with the expression x=g(x). Then, an initial guess for the root x_1 is assumed and input as an argument for the function g. The output g(x_1) is then the estimate x_2. The process is then iterated until the output x_{n+1}\approx g(x_n). The following is the algorithm for the fixed-point iteration method. Assuming x_0, \varepsilon_s, and maximum number of iterations N:
Set x_{n+1}=g(x_n), and calculate \varepsilon_r=\frac{x_{n+1}-x_{n}}{x_{n+1}} and compare with \varepsilon_s. If |\varepsilon_r|\leq \varepsilon_s or if n=N, then stop the procedure, otherwise, repeat.

The Babylonian method for finding roots described in the introduction section is a prime example of the use of this method. If we seek to find the solution for the equation x^2=s or f(x)=x^2-s=0, then a fixed-point iteration scheme can be implemented by writing this equation in the form:

    \[ x=\frac{x+\frac{s}{x}}{2} \]

Example

Consider the function f(x)=\cos{x}-x. We wish to find the root of the equation f(x)=0, i.e., \cos{x}-x=0. The expression can be rearranged to the fixed-point iteration form x=\cos{x} and an initial guess x_0=0.1 can be used. The tolerance \varepsilon_s is set to 0.001. The following is the Microsoft Excel table showing that the tolerance is achieved after 19 iterations:
excel1

Mathematica has a built-in algorithm for the fixed-point iteration method. The function “FixedPoint[f,Expr,n]” applies the fixed-point iteration method x=f(x) with the initial guess being “Expr” with a maximum number of iterations “n”. As well, the function “FixedPointList[f,Expr,n]” returns the list of applying the function “n” times. If “n” is omitted, then the software applies the fixed-point iteration method until convergence is achieved. Here is a snapshot of the code and the output for the fixed-point iteration x=\cos{x}. The software finds the solution x=0.739085. Compare the list below with the Microsoft Excel sheet above.
fixedpoint1

Alternatively, simple code can be written in Mathematica with the following output
fixedpoint2
View Mathematica Code

x = {0.1};
er = {1};
es = 0.001;
MaxIter = 100;
i = 1;
While[And[i <= MaxIter, Abs[er[[i]]] > es],  xnew = Cos[x[[i]]]; x = Append[x, xnew]; ernew = (xnew - x[[i]])/xnew; er = Append[er, ernew]; i++];
T = Length[x];
SolutionTable = Table[{i - 1, x[[i]], er[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration", "x", "er"};
T = Prepend[SolutionTable, SolutionTable1];
T // MatrixForm
Convergence

The fixed-point iteration method converges easily if in the region of interest we have |g'(x)|<1. Otherwise, it does not converge. Here is an example where the fixed-point iteration method fails to converge.

Example

Consider the function f(x)=\sin{5x}+\cos{2x}. To find the root of the equation f(x)=0, the expression can be converted into the fixed-point iteration form as:
x=\sin{5x}+\cos{2x}+x. Implementing the fixed-point iteration procedure shows that this expression almost never converges but oscillates:

View Mathematica Code
Clear[x, g]
x = {-0.9};
er = {1};
es = 0.001;
MaxIter = 100;
i = 1;
g[x_] := (Sin[5 x] + Cos[2 x]) + x
While[And[i <= MaxIter, Abs[er[[i]]] > es],xnew = g[x[[i]]]; x = Append[x, xnew]; ernew = (xnew - x[[i]])/xnew;er = Append[er, ernew]; i++];
T = Length[x];
SolutionTable = Table[{i - 1, x[[i]], er[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration", "x", "er"};
T = Prepend[SolutionTable, SolutionTable1];
T // MatrixForm

The following is the output table showing the first 45 iterations. The value of the error oscillates and never decreases:
fixedpoint3

The expression f(x)=0 can be converted to different forms x=g(x). For example, assuming x\neq 0:

    \[ 0=\sin{5x}+\cos{2x}\Rightarrow 0=\frac{\sin{5x}+\cos{2x}}{x}\Rightarrow x=\frac{\sin{5x}+\cos{2x}}{x}+x \]

If this expression is used, the fixed-point iteration method does converge depending on the choice of x_0. For example, setting x_0=5.0 gives the estimate for the root x_6=4.26452 with the required accuracy:
fixedpoint4

View Mathematica Code
Clear[x, g]
x = {5.};
er = {1};
es = 0.001;
MaxIter = 100;
i = 1;
g[x_] := (Sin[5 x] + Cos[2 x])/x + x
While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = g[x[[i]]]; x = Append[x, xnew]; ernew = (xnew - x[[i]])/xnew; er = Append[er, ernew]; i++];
T = Length[x];
SolutionTable = Table[{i - 1, x[[i]], er[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration", "x", "er"};
T = Prepend[SolutionTable, SolutionTable1];
T // MatrixForm

Obviously, unlike the bracketing methods, this open method cannot find a root in a specific interval. The root is a function of the initial guess x_0 and the form x=g(x), but the user has no other way of forcing the root to be within a specific interval. It is very difficult, for example, to use the fixed-point iteration method to find the roots of the expression 0=\sin{5x}+\cos{2x} in the interval [-1,1].

Analysis of Convergence of the Fixed-Point Method

The objective of the fixed-point iteration method is to find the true value V_t that satisfies V_t=g(V_t). In each iteration we have the estimate x_{i+1}=g(x_i). Using the mean value theorem, we can write the following expression:

    \[ g(x_i)=g(V_t)+g'(\xi)(x_i-V_t) \]

for some \xi in the interval between x_i and the true value V_t. Replacing g(V_t) and g(x_i) in the above expression yields:

    \[ x_{i+1}=V_t+g'(\xi)(x_i-V_t)\Rightarrow x_{i+1}-V_t=g'(\xi)(x_i-V_t) \]

The error after iteration i+1 is equal to E_{i+1}=x_{i+1}-V_t while that in iteration i is equal to E_i=x_i-V_t. Therefore, the above expression yields:

    \[ E_{i+1}=g'(\xi)E_i \]

For the error to reduce after each iteration, the first derivative of g, namely g', should be bounded by 1 in the region of interest (around the required root):

    \[ |g'(\xi)|< 1 \]

We can now try to understand why, in the previous example, the expression x=g(x)=\sin{5x}+\cos{2x}+x does not converge. When we plot g(x) and g'(x) we see that g'(x) oscillates rapidly with values higher than 1:

fixedpoint5

On the other hand, the expression x=g(x)=\frac{\sin{5x}+\cos{2x}}{x}+x converges for roots that are away from zero. When we plot g(x) and g'(x) we see that the oscillations in g'(x) decrease when x is away from zero and g'(x) is bounded by 1 in some regions:
fixedpoint6

Example

In this example, we will visualize the example of finding the root of the expression x^4-x-10=0. There are three different forms for the fixed-point iteration scheme:

    \[\begin{split} x&=g_1(x)=\frac{10}{x^3-1}\\ x&=g_2(x)=(x+10)^{1/4}\\ x&=g_3(x)=\frac{(x+10)^{1/2}}{x} \end{split} \]

To visualize the convergence, notice that if we plot the separate graphs of the function x and the function g(x), then, the root is the point of intersection when x=g(x). For g_1(x), the slope g_1'(x) is not bounded by 1 and so, the scheme diverges no matter what x_0 is. Here we start with x_0=2:

 

For g_2(x), the slope g_2'(x) is bounded by 1 and so, the scheme converges really fast no matter what x_0 is. Here we start with x_0=5:

 

For g_3(x), the slope g_3'(x) is bounded by 1 and so, the scheme converges but slowly. Here we start with x_0=5:

 

The following is the Mathematica code used to generate one of the tools above:
View Mathematica Code

Manipulate[
 x = {5.};
 er = {1};
 es = 0.001;
 MaxIter = 20;
 i = 1;
 g[x_] := (x + 10)^(1/4);
 While[And[i <= MaxIter, Abs[er[[i]]] > es],xnew = g[x[[i]]]; x = Append[x, xnew]; ernew = (xnew - x[[i]])/xnew;er = Append[er, ernew]; i++];
 T1 = Length[x];
 SolutionTable = Table[{i - 1, x[[i]], er[[i]]}, {i, 1, T1}];
 SolutionTable1 = {"Iteration", "x", "er"};
 T = Prepend[SolutionTable, SolutionTable1];
 LineTable1 = Table[{{T[[i, 2]], T[[i, 2]]}, {T[[i, 2]], g[T[[i, 2]]]}}, {i, 2, n}];
 LineTable2 = Table[{{T[[i + 1, 2]], T[[i + 1, 2]]}, {T[[i, 2]], g[T[[i, 2]]]}}, {i, 2, n - 1}];
 Grid[{{Plot[{x, g[x]}, {x, 0, 5}, PlotLegends -> {"x", "g3(x)"},ImageSize -> Medium, Epilog -> {Dashed, Line[LineTable1], Line[LineTable2]}]}, {Row[{"Iteration=", n - 2, "    x_n=", T[[n, 2]], "   g(x_n)=", g[T[[n, 2]]]}]}}], {n, 2, Length[T], 1}]

The following shows the output if we use the built-in fixed-point iteration function for each of g_1, g_2, and g_3. g_1 oscillates and so, it will never converge. g_2 converges really fast (3 to 4 iterations). g_3 converges really slow, taking up to 120 iterations to converge.
FixedPoint33

Newton-Raphson Method

The Newton-Raphson method is one of the most used methods of all root-finding methods. The reason for its success is that it converges very fast in most cases. In addition, it can be extended quite easily to multi-variable equations. To find the root of the equation f(x_t)=0, the Newton-Raphson method depends on the Taylor Series Expansion of the function around the estimate x_i to find a better estimate x_{i+1}:

    \[ f(x_{i+1})=f(x_i)+f'(x_i)(x_{i+1}-x_i)+\mathcal O (h^2) \]

where x_{i+1} is the estimate of the root after iteration i+1 and x_i is the estimate at iteration i. Assuming f(x_{i+1})=0 and rearranging:

    \[ x_{i+1}\approx x_i-\frac{f(x_i)}{f'(x_i)} \]

The procedure is as follows. Setting an initial guess x_0, tolerance \varepsilon_s, and maximum number of iterations N:
At iteration i, calculate x_i\approx x_{i-1}-\frac{f(x_{i-1})}{f'(x_{i-1})} and \varepsilon_r. If \varepsilon_r\leq \varepsilon_s or if i\geq N, stop the procedure. Otherwise repeat.
Note: unlike the previous methods, the Newton-Raphson method relies on calculating the first derivative of the function f(x). This makes the procedure very fast, however, it has two disadvantages. The first is that this procedure doesn’t work if the function f(x) is not differentiable. Second, the inverse (f'(x_{i}))^{-1} can be slow to calculate when dealing with multi-variable equations.

Example

As an example, let’s consider the function f(x)=\sin{5x}+\cos{2x}. The derivative of f is f'(x)=5\cos{5x}-2\sin{2x}. Setting the maximum number of iterations N=100, \varepsilon_s=0.0005, x_0=0.5, the following is the Microsoft Excel table produced:
NR1

Mathematica has a built-in algorithm for the Newton-Raphson method. The function “FindRoot[lhs==rhs,{x,x0}]” applies the Newton-Raphson method with the initial guess being “x0”. The following is a screenshot of the input and output of the built-in function evaluating the roots based on three initial guesses.
NR2

Alternatively, a Mathematica code can be written to implement the Newton-Raphson method with the following output for three different initial guesses:
NR3
View Mathematica Code

Clear[x]
f[x_] := Sin[5 x] + Cos[2 x]
fp = D[f[x], x]
xtable = {-0.25};
er = {1};
es = 0.0005;
MaxIter = 100;
i = 1;
While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = xtable[[i]] - (f[xtable[[i]]]/fp /. x -> xtable[[i]]);  xtable = Append[xtable, xnew]; ernew = (xnew - xtable[[i]])/xnew;  er = Append[er, ernew]; i++];
T = Length[xtable];
SolutionTable = Table[{i - 1, xtable[[i]], er[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration", "x", "er"};
T = Prepend[SolutionTable, SolutionTable1];
T // MatrixForm
Analysis of Convergence of the Newton-Raphson Method

The error in the Newton-Raphson Method can be roughly estimated as follows. The estimate x_{i+1} is related to the previous estimate x_i using the equation:

    \[ x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}\Rightarrow 0=f(x_i)+f'(x_i)(x_{i+1}-x_i) \]

Additionally, using Taylor’s theorem, and if x_t is the true root with f(x_t)=0 we have:

    \[ f(x_t)=0=f(x_i)+f'(x_i)(x_t-x_i)+\frac{f''(\xi)}{2!}{(x_t-x_i)^2} \]

for some \xi in the interval between x_i and x_t. Subtracting the above two equations yields:

    \[ 0=f'(x_i)(x_t-x_{i+1})+\frac{f''(\xi)}{2!}{(x_t-x_i)^2} \]

Since E_{i+1}=x_t-x_{i+1} and E_i=x_t-x_i then:

(1)   \begin{equation*} E_{i+1}=-\frac{f''(\xi)}{2!f'(x_i)}{E_i^2} \end{equation*}

If the method is converging, we have f''(\xi)\approx f''(x_t) and f'(\xi)\approx f'(x_t) therefore:

    \[ E_{i+1}\propto E_i^2 \]

Therefore, the error is squared after each iteration, i.e., the number of correct decimal places approximately doubles with each iteration. This behaviour is called quadratic convergence. Look at the tables in the Newton-Raphson example above and compare the relative error after each step!

The following tool can be used to visualize how the Newton-Raphson method works. Using the slope f'(x_i), the difference x_{i+1}-x_i can be calculated and thus, the value of the new estimate x_{i+1} can be computed accordingly. Use the slider to see how fast the method converges to the true solution using \varepsilon_s=0.0005, x_0=0.4, and solving for the root of f(x)=\sin[5x]+\cos[2x]=0.

View Mathematica code
Manipulate[
 f[x_] := Sin[5 x] + Cos[2 x];
 fp = D[f[x], x];
 xtable = {0.4};
 er = {1};
 es = 0.0005;
 MaxIter = 100;
 i = 1;
 While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = xtable[[i]] - (f[xtable[[i]]]/fp /. x -> xtable[[i]]); xtable = Append[xtable, xnew]; ernew = (xnew - xtable[[i]])/xnew; er = Append[er, ernew]; i++];
 T = Length[xtable];
 SolutionTable = Table[{i - 1, xtable[[i]], er[[i]]}, {i, 1, T}];
 SolutionTable1 = {"Iteration", "x", "er"};
 T = Prepend[SolutionTable, SolutionTable1];
  LineTable1 = Table[{{T[[i, 2]], 0}, {T[[i, 2]], f[T[[i, 2]]]}, {T[[i + 1, 2]],0}}, {i, 2, n - 1}];
 Grid[{{Plot[f[x], {x, 0, 1}, PlotLegends -> {"f(x)"}, ImageSize -> Medium,  Epilog -> {Dashed, Line[LineTable1]}]}, {Row[{"Iteration=", n - 2, "    x_n=", T[[n, 2]], "   f(x_n)=", f[T[[n, 2]]]}]}}], {n, 2, 7, 1}]

Depending on the shape of the function and the initial guess, the Newton-Raphson method can get stuck around the locations of oscillations of the function. The tool below visualizes the algorithm when trying to find the root of f(x)=x^3-x+3=0 with an initial guess of x_0=-0.1. It takes 33 iterations before reaching convergence.

Depending on the shape of the function and the initial guess, the Newton-Raphson method can get stuck in a loop. The tool below visualizes the algorithm when trying to find the root of f(x)=x^3-x+3=0 with an initial guess of x_0=0. The algorithm goes into an infinite loop.

In general, an initial guess that is close enough to the true root will guarantee quick convergence. The tool below visualizes the algorithm when trying to find the root of f(x)=x^3-x+3=0 with an initial guess of x_0=-1. The algorithm quickly converges to the desired root.

The following is the Mathematica code used to generate one of the tools above:
View Mathematica Code

Manipulate[
 f[x_] := x^3 - x + 3;
 fp = D[f[x], x];
 xtable = {-0.1};
 er = {1.};
 es = 0.0005;
 MaxIter = 100;
 i = 1;
 While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = xtable[[i]] - (f[xtable[[i]]]/fp /. x -> xtable[[i]]); xtable = Append[xtable, xnew]; ernew = (xnew - xtable[[i]])/xnew; er = Append[er, ernew]; i++];
 T = Length[xtable];
 SolutionTable = Table[{i - 1, xtable[[i]], er[[i]]}, {i, 1, T}];
 SolutionTable1 = {"Iteration", "x", "er"};
 T = Prepend[SolutionTable, SolutionTable1]; 
 LineTable1 =   Table[{{T[[i, 2]], 0}, {T[[i, 2]], f[T[[i, 2]]]}, {T[[i + 1, 2]], 0}}, {i, 2, n - 1}];
 Grid[{{Plot[f[x], {x, -3, 2}, PlotLegends -> {"f(x)"}, ImageSize -> Medium, Epilog -> {Dashed, Line[LineTable1]}]}, {Row[{"Iteration=", n - 2, "    x_n=", T[[n, 2]], "   f(x_n)=", f[T[[n, 2]]]}]}}] , {n, 2, 33, 1}]

Secant Method

The secant method is an alternative to the Newton-Raphson method by replacing the derivative f'(x) with its finite-difference approximation. The secant method thus does not require the use of derivatives especially when f(x) is not explicitly defined. In certain situations, the secant method is preferable over the Newton-Raphson method even though its rate of convergence is slightly less than that of the Newton-Raphson method. Consider the problem of finding the root of the function f(x)=0. Starting with the Newton-Raphson equation and utilizing the following approximation for the derivative f'(x_{i}):

    \[ f'(x_i)=\frac{f(x_i)-f(x_{i-1})}{x_i-x_{i-1}} \]

the estimate for iteration i+1 can be computed as:

    \[ x_{i+1}=x_i-f(x_i)\frac{x_i-x_{i-1}}{f(x_i)-f(x_{i-1})} \]

Obviously, the secant method requires two initial guesses x_0 and x_1.

Example

As an example, let’s consider the function f(x)=\sin{5x}+\cos{2x}. Setting the maximum number of iterations N=100, \varepsilon_s=0.0005, x_0=0.5, and x_1=0.4, the following is the Microsoft Excel table produced:
secant1

The Mathematica code below can be used to program the secant method with the following output:
secant2

View Mathematica Code
f[x_] := Sin[5 x] + Cos[2 x];
xtable = {0.5, 0.4};
er = {1,1};
es = 0.0005;
MaxIter = 100;
i = 2;
While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = xtable[[i]] - (f[xtable[[i]]]) (xtable[[i]] - xtable[[i - 1]])/(f[xtable[[i]]] -f[xtable[[i - 1]]]); xtable = Append[xtable, xnew];  ernew = (xnew - xtable[[i]])/xnew; er = Append[er, ernew]; i++];
T = Length[xtable];
SolutionTable = Table[{i - 1, xtable[[i]], er[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration", "x", "er"};
T = Prepend[SolutionTable, SolutionTable1];
T // MatrixForm
Convergence Analysis of the Secant Method

The i+1 estimate in the secant method is obtained as follows:

    \[ x_{i+1}=x_i-f(x_i)\frac{x_i-x_{i-1}}{f(x_i)-f(x_{i-1})} \]

Multiplying both sides by -1 and adding the true value of the root x_t where f(x_t)=0 for both sides yields:

    \[ x_t-x_{i+1}=x_t-x_i+f(x_i)\frac{x_i-x_{i-1}}{f(x_i)-f(x_{i-1})} \]

Using algebraic manipulations:

    \[\begin{split} x_t-x_{i+1}&=\frac{(x_t-x_i)(f(x_i)-f(x_{i-1}))+f(x_i)(x_i-x_{i-1})}{f(x_i)-f(x_{i-1})}\\ &=\frac{f(x_i)(x_t-x_{i-1})-f(x_{i-1})(x_t-x_i)}{f(x_i)-f(x_{i-1})}\\ &=(x_t-x_{i-1})(x_t-x_i)\frac{\frac{f(x_i)}{x_t-x_i}-\frac{f(x_{i-1})}{x_t-x_{i-1}}}{f(x_i)-f(x_{i-1})} \end{split} \]

Using the Mean Value Theorem, the denominator on the right-hand side can be replaced with:

    \[ f(x_i)=f(x_{i-1})+f'(\zeta)(x_i-x_{i-1}) \]

for some \zeta between x_i and x_{i-1}. Therefore,

    \[ x_t-x_{i+1}=(x_t-x_{i-1})(x_t-x_i)\frac{\frac{f(x_i)}{x_t-x_i}-\frac{f(x_{i-1})}{x_t-x_{i-1}}}{f'(\zeta)(x_i-x_{i-1})} \]

Using Taylor’s theorem for f(x_{i-1}) and f(x_i) around x_t we get:

    \[\begin{split} \frac{f(x_i)}{x_t-x_i}=f'(x_t)+\frac{f''(\xi_1)}{2}{(x_t-x_i)}\\ \frac{f(x_i)}{x_t-x_{i-1}}=f'(x_t)+\frac{f''(\xi_2)}{2}{(x_t-x_{i-1})} \end{split} \]

for some \xi_1 between x_i and x_t and some \xi_2 between x_{i-1} and x_t. Using the above expressions we can reach the equation:

    \[ x_t-x_{i+1}=(x_t-x_{i-1})(x_t-x_i)\frac{\frac{f''(\xi_1)}{2}{(x_t-x_i)}-\frac{f''(\xi_2)}{2}{(x_t-x_{i-1})}}{f'(\zeta)(x_i-x_{i-1})} \]

\xi_1 and \xi_2 can be assumed to be identical and equal to \xi, therefore:

    \[\begin{split} x_t-x_{i+1}&=(x_t-x_{i-1})(x_t-x_i)\frac{\frac{f''(\xi)}{2}{(x_t-x_i-x_t+x_{i-1})}}{f'(\zeta)(x_i-x_{i-1})}\\ &=(x_t-x_{i-1})(x_t-x_i)\frac{-f''(\xi)}{2f'(\zeta)} \end{split} \]

Therefore:

(2)   \begin{equation*} E_{i+1}=E_{i}E_{i-1}\left(\frac{-f''(\xi)}{2f'(\zeta)}\right) \end{equation*}

Comparing 1 with 2 shows that the convergence in the secant method is not quite quadratic. To find the order of convergence, we need to solve the following equation for a positive p and C:

    \[ |E_{i+1}|=C|E_{i}|^p \]

Substituting into equation 2 yields:

    \[ C|E_{i}|^p=|E_{i}||E_{i-1}|\bigg{|}\frac{-f''(\xi)}{2f'(\zeta)}\bigg{|}\Rightarrow E_{i}=\left|\frac{-f''(\xi)}{2Cf'(\zeta)}\right|^{\frac{1}{p-1}}|E_{i-1}|^{\frac{1}{p-1}} \]

But we also have:

    \[ |E_{i}|=C|E_{i-1}|^p \]

Therefore: p=\frac{1}{p-1}. This equation is called the golden ratio and has the positive solution for p:

    \[ p=\frac{\sqrt{5}+1}{2}\approx 1.618034 \]

while

    \[ C=\left|\frac{-f''(\xi)}{2f'(\zeta)}\right|^{p-1}=\left|\frac{-f''(\xi)}{2f'(\zeta)}\right|^{0.618034} \]

implying that the error convergence is not quadratic but rather:

    \[ E_{i+1}\propto E_i^{1.618034} \]

The following tool visualizes how the secant method converges to the true solution using two initial guesses. Using \varepsilon_s=0.0005, x_0=0.35, x_1=0.4, and solving for the root of f(x)=\sin[5x]+\cos[2x]=0 yields x_6=0.6732.

The following Mathematica Code was utilized to produce the above tool:
View Mathematica Code

Manipulate[
 f[x_] := x^3 - x + 3;
 f[x_] := Sin[5 x] + Cos[2 x];
 xtable = {-0.5, -0.6};
 xtable = {0.35, 0.4};
 er = {1, 1};
 es = 0.0005;
 MaxIter = 100;
 i = 2;
 While[And[i <= MaxIter, Abs[er[[i]]] > es], xnew = xtable[[i]] - (f[xtable[[i]]]) (xtable[[i]] - xtable[[i-1]])/(f[xtable[[i]]] -f[xtable[[i - 1]]]); xtable = Append[xtable, xnew];  ernew = (xnew - xtable[[i]])/xnew; er = Append[er, ernew]; i++];
 T = Length[xtable];
 SolutionTable = Table[{i - 1, xtable[[i]], er[[i]]}, {i, 1, T}];
 SolutionTable1 = {"Iteration", "x", "er"};
 T = Prepend[SolutionTable, SolutionTable1];
 T // MatrixForm;
 If[n > 3,
  (LineTable1 = {{T[[n - 2, 2]], 0}, {T[[n - 2, 2]], f[T[[n - 2, 2]]]}, {T[[n - 1, 2]], f[T[[n - 1, 2]]]}, {T[[n, 2]], 0}};
   LineTable2 = {{T[[n - 1, 2]], 0}, {T[[n - 1, 2]], f[T[[n - 1, 2]]]}}), LineTable1 = {{}}; LineTable2 = {{}}];
 Grid[{{Plot[f[x], {x, 0, 1}, PlotLegends -> {"f(x)"}, ImageSize -> Medium, Epilog -> {Dashed, Line[LineTable1], Line[LineTable2]}]}, {Row[{"Iteration=", n - 2, "    x_n=", T[[n, 2]], "   f(x_n)=", f[T[[n, 2]]]}]}}], {n, 2, 8, 1}]

Problems

  1. Let f:[-5,5]\rightarrow \mathbb{R} be such that f(x)=0.5+x\cos{x}. Use the graphical method to find the number and estimates of the roots of each of the equations f(x)=0 and f(x)=2.
  2. Compare the bisection method and the false position method in finding the roots of the two equations f(x)=0 and f(x)=0.1 in the interval [-1,1] where f(x) is:

        \[ f(x)=2x^3-4x^2+3x \]

    Compare with the result using the Solve function and comment on the use of \varepsilon_r as a measure of the relative error.

  3. Locate the first nontrivial root of \sin{x}=x^2 where x is in radians. Use the graphical method and the bisection method with \varepsilon_s=0.002.
  4. Water is flowing in a trapezoidal channel at a rate of Q= 20 m^3/s. The critical depth y for such a channel must satisfy the equation:

        \[ 1-\frac{Q^2}{g A_c^3}B=0 \]

    where g = 9.81 m/s^2, A_c is the cross sectional area (m^2), and B is the width of the channel at the surface (m). Assuming that the width and the cross‐sectional area can be related to the depth y by:

        \[\begin{split} B&=3+y\\ A_c&=3y+\frac{y^2}{2} \end{split} \]

    Solve for the critical depth using the graphical method, the bisection method, and the false position method to find the critical depth in the interval [0,3]. Use \varepsilon_s=0.001. Discuss your results (difference in the methods, accuracy, number of iterations, etc.).

    (Note: The critical depth is the depth below which the flow of the fluid becomes relatively fast and affected by the upstream conditions (supercritical flow). Above that depth, the flow is subcritical, relatively slow and is controlled by the downstream conditions)

  5. Use the graphical method to find all real roots of the equation

        \[ f(x)=-14-20x+19x^2-3x^3 \]

    You are required to

    • Plot the function several times by using smaller ranges for the x axis to obtain the roots up to 3 significant digits.
    • Explain the steps you have taken to identify the number and estimates of the roots.
    • Explain the steps you have taken to obtain the roots to the required accuracy.
    • Compare your results to the solution obtained using the “Solve” function in Mathematica.
  6. Use the bisection method and the false position method to find all real roots of the equation

        \[ f(x)=-14-20x+19x^2-3x^3 \]

    Use \varepsilon_s=0.001.

  7. Consider the function:

        \[  f(x)= -0.9 x^2+1.7x+2.5 \]

    Determine a root for the expression f(x)=0

    • using fixed-point iteration
    • using the Newton-Raphson method

    Use an initial guess of x_0 = 5 and \varepsilon_s=0.01\%=0.0001. Compare your final answer with the Solve function in Mathematica.

  8. Consider the function:

        \[ f(x)=7\sin{(x)} e^{-x}-1 \]

    Find the lowest positive root for the expression f(x)=0

    • graphically
    • using the Newton-Raphson method
    • using the secant method
  9. Use the Newton-Raphson method to find the root of:

        \[ f(x)=e^{-0.5x}(4-x)-2=0 \]

    Employ initial guesses of 2, 6, and 8. Explain your results.

  10. The Manning equation can be written for a rectangular open channel as:

        \[ Q=\frac{\sqrt{S}(BH)^\frac{5}{3}}{n(B+2H)^{\frac{2}{3}}} \]

    where Q is the fluid flow in m^3/s, S is the slope in m/m, H is the depth in m, and n is the Manning roughness coefficient. Develop a fixed-point iteration scheme to solve this equation for H given Q = 5, S =0.0002, B = 20, and n = 0.03. Consider \varepsilon_s=0.0005. Prove that your scheme converges for all initial guesses greater than or equal to zero.

  11. Why does the Babylonian method almost always converge?
  12. Consider the function

        \[ f(x)=e^{-x}-x \]

    Compare the secant method and the Newton-Raphson method in finding the root of the equation f(x)=0.

Page Comments

  1. adib says:

    whats was the error formula in the newton raphson method?

    1. Samer Adeeb says:

      I don’t understand the question

  2. aty says:

    how can we find the error if we find the root using the graphical method

    1. Samer Adeeb says:

      I don’t think you can quantify the error using the graphical method except visually.

Leave a Reply

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