Roots Finding: Open Methods

Problem 6.15.

The Redlich-Kwong equation of state is given by

 Eq.(6.15(a))

where R = the universal gas constant [= 0.518 kJ/(kg K)], T = absolute temperature (K), p = absolute pressure (kPa), and v = the volume of a kg of gas (m³/kg). The parameters a and b are calculated by 

where p_c = 4600 kPa and T_c = 191 K. As a chemical engineering, you are asked to determine the amount of methane fuel that can be held in a 3-m³ tank at a temperature of -40°C with a pressure of 65000 kPa. Use a root-locating method of your choice to calculate v and then determine the mass of methane contained in the tank.

Solution:

The parameters and other values can be substituted into Eq.(6.15(a)). *Remember to express T in unit K (T(K) = T(°C) + 273.15)

 Eq.(6.15(b))

Therefore, we need to find out v, which is the specific volume (m³/kg) using root-locating method. Because it is relatively complicated, I will store it as an M-file:

function f = redlichkwong(v,P,T)

% v = specific volume (m3/kg)

% P = absolute pressure (kPa)

% T = absolute temperature (K)

R = 0.518; %kJ/kg.K

TC = 191; %K

PC = 4600; %kPa

a = 0.427*(R^2)*(TC^2.5)/(PC);

b = 0.0866*R*TC/PC;

f = P-R*T./(v-b)+a./(v.*(v+b).*sqrt(T));

end

Before determining the root, it is advisable to plot the function to estimate the initial guesses. 

>> v = linspace(0.002,0.0035);

>> f = redlichkwong(v,65000,233.15);

>> plot(v,f),grid,xlabel('v'),ylabel('f')

Hence, from the graph, we can obtain a rough estimate of the root, which is about 0.0028. We can either choose to use bracketing method or open method to solve the problem. I will recommend to use open method due to quicker convergence. In this case, we can use modified secant method with an initial guess, v₀ = 0.002 and perturbation fraction, δ = 0.001.

The above iterative equation can be applied to compute:

Therefore, after seven iterations, the approximate error falls below 0.1% and the root estimate is 0.002808 m³/kg. 

Mass of methane contained in the tank is:

m = 3 m³/0.002808 m³/kg = 1068.376 kg

Roots Finding: Open Methods

Problem 6.14.

In a chemical engineering process, water vapor (H₂O) is heated to sufficiently high temperatures that a significant portion of water dissociates, or splits apart, to form oxygen and hydrogen:

H₂O⇆H₂+1/2O₂ 

If it is assumed that this is the only reaction involved, the mole fraction x of H₂O that dissociates can be represented by

 Eq.(6.14(a))

where K is the reaction's equilibrium constant and p_t is the total pressure of the mixture. If p_t = 3 atm and K = 0.05, determine the value of x that satisfies Eq.(6.14(a)).

Solution:

Firstly, rearrange and substitute the value of p_t and K to form a single function in terms of x. 

Therefore, the solution x can be obtained by finding the root of the function. We can use any open-root locating method to solve for x. Before that, we can check the root location by using graphical method. 

Graphical method:

>> x = linspace(0,0.1);

>> fx = x./(1-x).*sqrt(6./(2+x))-0.05;

>> plot(x,fx),grid,xlabel('x'),ylabel('f(x)')

From the graph, the root lies between 0.02 and 0.03. The root estimate is roughly about 0.028. 

Fixed-point iteration:

Fixed-point iteration can be applied by using an initial guess of x₀ = 0.02 and stopping criterion of 0.01%. Thus, rearrange the equation:

Starting with an initial guess of x₀ = 0.02, the above iterative equation can be applied to compute:

Hence, after four iterations, the root estimate (x) is 0.02825. 

We can also use MATLAB fzero function to find the root. 

>> format long

>> fx =@(x) x./(1-x).*sqrt(6./(2+x))-0.05;

>> x = fzero(fx,0.02)

x =

   0.028249441148471

Roots Finding: Open Methods

Problem 6.4.

Determine the lowest positive root of f(x) = 7sin(x)e^-x - 1: (a) Graphically. (b) Using the Newton-Raphson method (three iterations, x₀ = 0.3). (c) Using the secant method (three iterations, x_-1 = 0.5 and x₀ = 0.4). (d) Using the modified secant method (five iterations, x₀ = 0.3, δ = 0.01).

Solution:

(a) Plot the graph using MATLAB. The commands are shown as below:

>> x = linspace(-1,2);

>> fx = 7*sin(x).*exp(-x)-1;

>> plot(x,fx),grid,xlabel('x'),ylabel('f(x)')

The graph generated is shown as above. As shown by the graph, there are two roots. The lowest positive root is between 0 and 0.5. By maximizing the graph, the root is estimated as 0.17. 

(b) Firstly, compute the first derivative of the function. (Recall Product Rule)

which can be substituted along with the original function to give

Starting with the initial guess of x₀ = 0.3, the iterative equation can be applied to compute:

Hence, after three iterations, the root estimate is 0.17018 with an approximate error of 0.453%.

(c) The iterative equation using secant method is 

Begin the iteration with x_-1 = 0.5 and x₀ = 0.4 yields

First iteration:

x_-1 = 0.5     f(0.5) = 1.0355

x₀ = 0.4      f(0.4) = 0.8272

x₁ = 0.0028 

The calculation can be continued to yield

Thus, after three iterations, the root estimate is 0.17899 with an approximate error of 21.9%.

(d) The iterative equation using modified secant method is 

First iteration:

x₀ = 0.3                                    f(0.3) = 0.5325

x₀ + δx₀ = 0.3030                  f(0.3030) = 0.5427

x₁ = 0.1437

The calculation can be continued to yield

Hence, after five iterations, the root estimate is 0.17018 with a zero percent approximate relative error.

 Additional: The real root of the equation can be found using fzero function in MATLAB.

>> format long

>> fx =@(x) 7*sin(x).*exp(-x)-1;

>> xr = fzero(fx,0.3)

xr =

   0.170179993753835

Discussion:

In conclusion, the Newton-Raphson method and modified secant method both give comparably accurate result with the true root value. In addition, modified secant method is able to attain the efficiency of Newton-Raphson method without having to compute derivatives. 

Roots Finding: Open Methods

Problem 6.2.

Use (a) fixed-point iteration and (b) the Newton-Raphson method to determine a root of f(x) = -0.9x² + 1.7x + 2.5 using x₀ = 5. Perform the computation until ε_a is less than ε_s = 0.01%. 

Solution:

(a) Firstly, rearrange the function so that x is on the left-hand side of the equation (x = g(x)). There are two ways to do so:

To identify which function results in convergence, we should check for the criterion |g'| < 1. For (i), |g'(5)| = 5.29. While for (ii), |g'(5)| = 0.27 < 1, where the error decreases with each iteration. Thus, function (ii) will result in convergence.

Starting with the initial guess of x₀ = 5, the iterative equation (ii) can be applied to compute:

Therefore, after 9 iterations, the approximate error is reduced to less than 0.01% and the root estimate is 2.8602. 

(b) Firstly, compute the first derivative of the function. 

which can be substituted along with the original function to give

Starting with the initial guess of x₀ = 5, the iterative equation can be applied to compute:

Hence, after 5 iterations, the approximate error falls below 0.01% and the root estimate is 2.8601. 

Discussion:

As we can observe, the Newton-Raphson method rapidly converges on the true root. Note that the approximate percent relative error at each iteration decreases much faster than it does in simple fixed-point iteration. However, Newton-Raphson method is not recommended when there are functions whose derivatives are difficult to evaluate. 

Roots Finding: Bracketing Methods

Problem 5.12.

Water is flowing in a trapezoidal channel at a rate of Q = 20 m³/s. The critical depth y for such a channel must satisfy the equation

   Eq 5.12.(a)

where g = 9.81 m/s², A_c = the cross sectional area (m²), and B = the width of the channel at the surface (m). For this case, the width and the cross sectional area can be related to the depth y by B = 3 + y and A_c = 3y + y²/2. Solve for the critical depth using the (a) graphical method (b) bisection, and (c) false position. For (b) and (c) use the initial guesses of x_l = 0.5 and x_u = 2.5, and iterate until the approximate error falls below 1% or the number of iterations exceeds 10. Discuss your result.

Solution:

Firstly, compute a single function that only depends on y using Eq.5.12(a) by substituting the relation for B and A_c. 

Therefore, the critical depth y will be the value of y that makes f(y) = 0. 

(a) Plot the graph using MATLAB. The commands are shown as follow.

>> y = linspace(0.1,2.5);
>> fy = 1- (400./(9.81*(3*y+y.^2/2).^3)).*(3+y);
>> plot(y,fy),grid,xlabel('y'),ylabel('f(y)')

The graph generated is shown as above. The rough estimate of the root is about 1.5 m. The validity of the graphical method can be checked by substituting into the function f(y) to yield

>> fy = 1- (400/(9.81*(3*1.5+1.5^2/2)^3))*(3+1.5)

fy =

   -0.0309

which is close to zero. 

(b) The root estimate using bisection method is x_r = (x_l + x_u)/2. Begin the iteration with guesses of x_l = 0.5 and x_u = 2.5.

First iteration:

x_l = 0.5                                f(0.5) = -32.2582

x_u = 2.5                               f(2.5) = 0.8130

x_r = (0.5+2.5)/2 = 1.5     f(1.5) = -0.0309 

f(0.5)f(1.5) = 0.9983 > 0. Therefore, the root lies in the upper interval, and x_r becomes the lower limit for the next iteration, x_l = 1.5. 

Second iteration:

x_l = 1.5                             f(1.5) = -0.0309

x_u = 2.5                            f(2.5) = 0.8130

x_r = (1.5+2.5)/2 = 2      f(2) = 0.6018 

f(1.5)f(2.5) = -0.0186 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = 2. 

The approximate relative error is 25%. 

The remainder of the iterations are displayed in the following table. 

Hence, after 8 iterations, the approximate relative error falls below 1% and the computation can be terminated. The critical depth is estimated as 1.50781 m. 

(c) The root estimate using false position method is x_r = x_u - f(x_u)(x_l - x_u)/(f(x_l) - f(x_u)). Begin the iteration with guesses of x_l = 0.5 and x_u = 2.5.

First iteration:

x_l = 0.5                        f(0.5) = -32.2582

x_u = 2.5                       f(2.5) = 0.8130

x_r = 2.45083             f(2.45083) = 0.7999

f(0.5)f(2.45803) = -25.8025 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = 2.45083. 

Second iteration:

x_l = 0.5                          f(0.5) = -32.2582

x_u = 2.45803               f(2.45803) = 0.7999

x_r = 2.40363                f(2.40363) = 0.7861

f(0.5)f(2.40363) = -25.3589 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = 2.40363. 

The approximate relative error is 1.96%. 

The remainder of the iterations are displayed in the following table. 

Hence, the computation is terminated after 10 iterations. The critical depth is estimated as 2.0908 m with an approximate relative error of 1.59%. 

Discussion:

As we can observe, using bisection, the approximate error is reduced to less than 1% after 8 iterations. For false position, the approximate error is still above 1% after 10 iterations. This shows that in this case, bisection is more preferable compared to false position. Insight into these results can be gained by examining the plot of the function. As shown in part(a), the curve violates the premise on which false position was based - that is, if f(x_u) is much closer to zero than f(x_l), then the root should be much closer to x_u than to x_l . Because of the shape of the present function, the opposite is true. 

Tips: When the function has significant curvature, false position is not preferable to solve the problem because of its one-sidedness, which leads to poor convergence. 

Roots Finding: Bracketing Methods

Problem 5.7. 

(a) Determine the roots of f(x) = -12 - 21x + 18x² - 2.75x³ graphically. In addition determine the first root of the function with (b) bisection and (c) false position. For (b) and (c), use x_l = -1 and x_u = 0, and a stopping criterion of 1%.

Solution:

(a) Plot the graph using MATLAB. The commands are shown as below:

>> x = linspace(-1,8);

>> fx = -12-21*x+18*x.^2-2.75*x.^3;

>> plot(x,fx),grid,xlabel('x'),ylabel('f(x)')

The graph generated is shown as above. Since it is a cubic equation, there are 3 roots. By maximizing the graph, the roots are estimated as -0.4135, 2.2 and 4.74. 

(b) The root estimate using bisection method is x_r = (x_l + x_u)/2. Begin the iteration with guesses of x_l = -1 and x_u = 0.

First iteration:

x_l = -1                              f(-1) = 29.75

x_u = 0                               f(0) = -12

x_r = (-1+0)/2 = -0.50     f(-0.50) = 3.3438 

f(-1)f(-0.50) = 99.4766 > 0. Therefore, the root lies in the upper interval, and x_r becomes the lower limit for the next iteration, x_l = -0.50. 

Second iteration:

x_l = -0.50                             f(-0.50) = 3.3438

x_u = 0                                   f(0) = -12

x_r = (-0.50+0)/2 = -0.25    f(-0.25) = -5.5820

f(-0.50)f(-0.25) = -18.66496 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = -0.25. 

The approximate relative error is 100%. 

The remainder of the iterations are displayed in the following table. 

Hence, after 8 iterations, the approximate error finally falls below 1% and the computation can be terminated. The root estimate is -0.41797. 

(c) The root estimate using false position method is x_r = x_u - f(x_u)(x_l - x_u)/(f(x_l) - f(x_u)). Begin the iteration with guesses of x_l = -1 and x_u = 0.

First iteration:

x_l = -1                        f(-1) = 29.75

x_u = 0                        f(0) = -12

x_r = -0.2874             f(-0.2874) = -4.4116

f(-1)f(-0.2874) = -131.2731 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = -0.2874. 

Second iteration:

x_l = -1                          f(-1) = 29.75

x_u = -0.2874               f(-0.2874) = -4.4116

x_r = -0.3795                f(-0.3795) = -1.2878

f(-1)f(-0.3795) = -38.3130 < 0. Therefore, the root lies in the lower interval, and x_r becomes the upper limit for the next iteration, x_u = -0.3795. 

The approximate relative error is 24.3%. 

The remainder of the iterations are displayed in the following table. 

Hence, after 5 iterations, the approximate error finally falls below 1% and the computation can be terminated. The root estimate is -0.41402. 

Good day, everyone! My name is Shu Jia, currently an undergraduate chemical engineering student at Universiti Putra Malaysia. The main purpose of establishing this blog is to share my knowledge on some of the chemical engineering subjects by providing the solutions for example questions. Step-by-step solution will be provided so that you can understand how to solve the question, instead of just knowing the answers. Do subscribe to get the latest updates of solutions! 

*Kind advice: Please try to do the questions on your own before referring to my solutions. The solution is just a guideline, full understanding on the related topic is always the priority to master a subject!