Some Applications of Cubic Equations in Engineering

Some Applications of Cubic Equations in Engineering

ZahediAnton Abdulbasah Kamil Irvan Jelita Harahap Amin Affandi Marwan Sarmin Suparni 

Universitas Sumatera Utara, Jln. Dr. T. Mansyur No. 9, Medan 20155, Indonesia

Istanbul Gelisim University, Istanbul 34310, Turkey

Universitas Muhammadiyah Sumatera Utara, Jln. Kapt. Mukhtar Basri No. 3, Medan 20238, Indonesia

IAIN Langsa, Jln. Meurandeh Kota Langsa 24411, Aceh, Indonesia

FKIP Universitas Labuhan Batu, Jln. SM. Raja No. 126, Aek Tapa, Rantau Prapat 21418, Kab. Labuhan Batu, Sumatera Utara, Indonesia

IAIN Padangsidimpuan, Jln. T. Rizal Nurdin km 4.5, Sihitang 22733, Padangsidimpuan, Indonesia

Previously at Rajamangala University of Technology Srivijaya (RUTS), Songkhla 90000, Thailand

Corresponding Author Email: 
zahedi@usu.ac.id
Page: 
129-135
|
DOI: 
https://doi.org/10.18280/mmep.090116
Received: 
26 June 2021
|
Revised: 
18 December 2021
|
Accepted: 
28 December 2021
|
Available online: 
28 February 2022
| Citation

© 2022 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

Cubic equations have many applications in engineering, three of them are discussed in this paper. There are many equations of states for real gases but the cubic equations are the simplest ones and are sufficiently accurate for a limited range of temperatures and pressures. Degree of dissociation of chemical equilibrium for carbon dioxide and water can be written as cubic equations. Slope of a simply supported beam loaded by a continuous load is also represented as a cubic equation. Cubic equations can be solved exactly using Cardano’s formula. They can also solve numerically; Newton-Raphson method is a popular choice. Although a cubic equation has three roots, only real roots are valid in real applications discussed in this paper. Even there may be only one root that can be used; two other roots will be discarded. There are many ways a cubic equation solved but the simplest one is to solve it manually using a scientific calculator. Software and programming languages are better if there are many equations to be solved repeatedly.

Keywords: 

cubic equation, real root, equation of state, degree of dissociation, deflection curve, Cardano’s formula

1. Introduction

The purpose of this paper is to discuss a few examples on the applications of cubic equations in engineering. While there are many good examples and exercises in engineering mathematics or numerical analysis textbooks on solving cubic equations, real examples in this paper are expected to enhance the examples in those textbooks. Examples in this paper are taken from real applications so they may motivate students to study them better.

A cubic equation is a polynomial of degree three which is generally written as

$a x^{3}+b x^{2}+c x+d=0$          (1)

where, a, b, c and d are the coefficients of the equation and a is strictly not equal to zero. While the coefficients can be real or complex, we are only interested for real values of the coefficients. This will make at least one root of the equation real.

According to the fundamental theorem of algebra, a polynomial of degree n has n roots, including multiplicities [1, 2]. In general, roots of a polynomial can only be found exactly for degree less than five. This theory was proved by Galois, a young French mathematician in 1830, and described in most textbooks on abstract algebra [3-5]. It is well-known that a cubic equation can be solved in closed forms [6]. Clearly, a cubic equation will have three roots, which can be divided as:

1. Three real roots.

2. Multiple roots.

3. One real root and two imaginary ones.

In engineering applications, we are usually interested in real roots only. Often, we only need one real root, which must be chosen from the real roots. There are many applications of cubic equations in engineering but this paper will only discuss five applications in thermodynamics and two applications in engineering mechanics.

2. Equation of State

An equation of state (or EoS for short) is the one which relates pressure (P), temperature (T) and specific volume of a substance (v) [7-10]. Only two of those parameters are independent. Usually, P and T are independent. However, most equations of state are implicit in specific volume, makes it as the independent variable together with temperature. In the SI units, P is in Pa (N/m2), v is in m3/kg and T is in K (kelvin).

There are hundreds of equations of states that have been developed; and the Ref. [7] lists 25 of them. However, only a few of them which still remain popular while most of them have been replaced by more accurate but more complicated equations. We will only discuss the cubic equations of state. Cubic equations are very simple and may not represent the state of a substance very accurately except in a limited range of temperature and pressure. However, they are still useful since their simplicity makes them attractive for analysis and also for teaching. Literally, there are more than one hundred cubic equations of states that have been developed [11]. We will only briefly discuss four of them, which are well-known and simple i.e the Van der Waals (VdW), Redlich-Kwong (RK), Soave (SRK), and Peng-Robinson (PR) equations. Derivations for the formulas concerned are not given; they can be found in Ref. [7]. Most thermodynamics textbooks discuss about the VdW EoS [8, 10]. Each of the four cubic equations of states only has two parameters which maybe constant (for the VdW EoS) or variable (for the three other EoS). There is another constant which is the same for all equations but depends on a gas or substance concerned. This constant is called the gas constant, symbolized by R. The unit for R is kJ/kg.K. the gas constant is calculated from the formula R = Ru/M where Ru is the universal gas constant and M is the molecular weight of the gas. The value of Ru is 8.31447 kJ/kmol.K.

Shivakumar et al. [12] use three cubic equations of state to study a refrigerant being used in automotive air-conditioning applications. The thermodynamic property of refrigerant R134a is estimated for this purpose. The equations used are Redlich-Kwong, Peng-Robinson and Patel Teja. The last equation is not discussed in this paper. Deiters and Macías-Salinas [13] discusses the calculation of densities using cubic equation of states and proposes a new algorithm for the calculation.

2.1 Van der Waal’s equation (Vdw)

The VdW EoS is the oldest cubic equation of state which is written as:

$P=\frac{R T}{v-b}-\frac{a}{v^{2}}$          (2)

Eq. (2) can be arranged in terms of v  as:

$v^{3}-\left[b+\frac{R T}{P}\right] v^{2}+\frac{a}{P} v-\frac{a b}{P}=0$          (3)

or

$v^{3}+p v^{2}+q v+r=0$          (4)

which is a cubic equation, where

$p=-\left(b+\frac{R T}{P}\right), q=\frac{a}{P} \text { and } r=-\frac{a b}{P}$          (5)

Here, a, b  and R are constants which depends on a substance.

The VdW EoS is very simple and therefore cannot be expected to represent accurately the PvT behavior of substances accurately over a large region. Despite its limitation, the VdW EoS has a historical value since many equations have been developed from it.

2.2 Redlich-Kwong’s equation (RK)

The RK EoS is an improvement of the VdW equation and more accurate than the former. It is written as:

$P=\frac{R T}{v-b}-\frac{a}{\sqrt{T} v(v+b)}$          (6)

Eq. (6) can be arranged in terms of v as:

$v^{3}-\frac{R T}{P} v^{2}+\left[\frac{a}{P \sqrt{T}}-\frac{b R T}{P}-b^{2}\right] v-\frac{a b}{P \sqrt{T}}=0$          (7)

which is similar to Eq. (4) where

$p=-\frac{R T}{P}, q=\frac{a}{P \sqrt{T}}-\frac{b R T}{P}-b^{2}, r=-\frac{a b}{P \sqrt{T}}$          (8)

The RK EoS was developed in 1949. At that time, there were more than 200 equations of state that had been developed [14]. The RK EoS and was very successful due to its simplicity and accuracy. However, it was not conceived to treat the liquid phases. Nevertheless, other researchers improved the equation by modifying the constants.

2.3 Soave equation (SRK)

Soave replaced the term $a / \sqrt{T}$ in the RK equation with a function of $T \text { and } \omega$ involving temperature and acentric factor so Eq. (6) becomes,

$P=\frac{R T}{v-b}-\frac{a \alpha}{v(v+b)}$          (9)

which can be arranged in terms of v as:

$v^{3}-\frac{R T}{P} v^{2}+\left[\frac{a \alpha}{P}-\frac{b R T}{P}-b^{2}\right] v-\frac{a b \alpha}{P}=0$          (10)

where, $\alpha$ is given by:

$\alpha=\left[1+\left(0.48+1.574 \omega-0.176 \omega^{2}\right)\left(1-\sqrt{T_{r}}\right)\right]^{2}$          (11)

$T_{r}$  is called reduced temperature defined as $T / T_{c}$  where $T_{c}$  is the critical temperature while $\omega$ is the Pitzer acentric factor.

Eq. (10) is similar to Eq. (4) where

$p=-\frac{R T}{P}, q=\frac{a \alpha}{P}-\frac{b R T}{P}-b^{2}, r=-\frac{a b \alpha}{P}$          (12)

Georgio Soave developed his equation in 1972 by proposing modifications to the RK EoS. By analyzing the saturation behaviour of the RK EoS, he concluded that one of its major faults was the lack of reliability in reproducing the vapor pressures of pure compounds. His proposed equation is called the Soave-Redlich-Kwong equation of state or the SRK. The SRK EoS became then very popular in the industry. However, it had a notable shortcoming, particularly for applications in hydrocarbon processing. The improvement came from the Peng-Robinson EoS.

2.4 Peng-Robinson equation (PR)

Peng and Robinson proposed an equation of state in the form of:

$P=\frac{R T}{v-b}-\frac{a \alpha}{v^{2}+2 b v-b^{2}}$         (13)

where, $\alpha$  is given by:

$\begin{aligned}

\alpha=& {[1+(0.37464+1.54226 \omega-} \\

&\left.\left.0.26992 \omega^{2}\right)\left(1-\sqrt{T_{r}}\right)\right]^{2}

\end{aligned}$           (14)

Eq. (13) can be arranged as:

$\begin{gathered}

v^{3}-\left[\frac{R T}{P}-b\right] v^{2}+\left[\frac{a \alpha}{P}-\frac{2 b R T}{P}-3 b^{2}\right] v-\left[\frac{a b \alpha}{P}-\right. \\

\left.\frac{b^{2} R T}{P}-b^{3}\right]=0

\end{gathered}$          (15)

which is similar to Eq. (4) where

$\begin{gathered}

p=\frac{R T}{P}-b, q=\frac{a \alpha}{P}-\frac{2 b R T}{P}-3 b^{2}, r=\frac{a b \alpha}{P}- \\

\frac{b^{2} R T}{P}-b^{3}

\end{gathered}$          (16)

Peng and Robinson proposed their cubic equation of state in 1976. The equation has been very successful to estimate accurately the behaviour of many substances over a wide region of temperature and pressure. Since that time, there have been many attempts to improve the performance of the equation and to extend its application to more substances. Many researchers have put efforts to improve the alpha (α) function which is originally given in Eq. (14); see [11, 14] for more detailed discussion.

2.5 Equations of state in terms of compressibility factor

We can also write the equation of state of a gas as a function of compressibility factor Z. In general, a cubic equation for real gases with two parameters can be written as [7]:

$P=\frac{R T}{v-b}-\frac{a}{v^{2}+u b v+w b^{2}}$          (17)

which is equivalent with

$Z^{3}+p Z^{2}+q Z+r=0$          (18)

where, $Z$ is the compressibility factor of the gas defined as $P v / R T$. The values of $p, q \text { and } r$ in Eq. (17) are different from those values in Eq. (4). They are given as:

$\begin{gathered}

p=-(1+B-u B) \\

q=A+w B^{2}-u B-u B^{2} \\

r=-\left(A B+w B^{2}+w B^{3}\right)

\end{gathered}$          (19)

where,

$A=\frac{a P}{R^{2} T^{2}} \text { and } B=\frac{b P}{R T}$          (20)

In theory, there are three roots for a cubic equation. If there are three real roots found, the smallest one is taken for the saturated liquid, the biggest is for the saturated vapor while the third one is not used. If there is one real root and two imaginary ones, the real root represents the gas state while the imaginary roots are discarded.

3. Chemical Equilibrium of Water and Carbon Dioxide

At high temperatures, gases tend to dissociate or decompose into other gases [8, 10, 15, 16]. Examples of dissociation reactions are as follow:

$\begin{gathered}

\mathrm{H}_{2} \mathrm{O} \leftrightarrow \mathrm{H}_{2}+1 / 2 \mathrm{O}_{2} \\

\mathrm{H}_{2} \mathrm{O}+\mathrm{O}_{2} \leftrightarrow 2 \mathrm{OH} \\

\mathrm{CO}_{2} \leftrightarrow \mathrm{CO}+1 / 2 \mathrm{O}_{2} \\

\mathrm{~N}_{2}+\mathrm{O}_{2} \leftrightarrow \mathrm{NO}_{2}

\end{gathered}$

In general, for chemical equilibrium

aA+bB↔cC+dD    (21)

Equilibrium constant Kp is defined as [4, 6]:

$K_{p}=\frac{P_{C}^{c} P_{D}^{d}}{P_{A}^{a} P_{B}^{b}}\left[\frac{P}{P_{0}}\right]^{c+d-a-b}$     (22)

where, PA, PB, PC, and PD are partial pressures of components in the reaction while P and P0 are total pressure and reference pressure (which is usually atmospheric pressure), respectively. Eq. (22) can also be written as:

$K_{p}=\frac{n_{C}^{c}}{n_{A}^{a}} \frac{n_{D}^{d}}{n_{R}^{b}}\left(\frac{P}{P_{0}} \frac{1}{n_{t}}\right)^{c+d-a-b}=\frac{[C]^{c}}{[A]^{a}} \frac{[D]^{d}}{[B]^{b}}\left(\frac{P}{P_{0}} \frac{1}{n_{t}}\right)^{c+d-a-b}$     (23)

where, the square bracket denotes the amount of a gas.

For the dissociation of CO2, we have: a=1,b=0,c=1 and d=1/2. So,

$K_{p}=\frac{[\mathrm{CO}]\left[\mathrm{O}_{2}\right]^{1 / 2}}{\left[\mathrm{CO}_{2}\right]}\left(\frac{P}{P_{0}} \frac{1}{n_{t}}\right)^{1 / 2}$     (24)

We will now determine the composition of gases from the dissociation of CO2 at temperature T and pressure P: 

CO2 ↔ CO + ½O2

Let there be one mole of CO2. After the equilibrium is achieved, CO2 will dissociate to α mole of CO and $\frac{1}{2} \alpha$ mole of O2. The remaining amount of CO2 is (1-α) mole. The total amount of gases is ($1+\frac{1}{2} \alpha$) mole. Substituting these values to Eq. (24) results in

$K(\alpha)=\frac{\alpha}{1-\alpha} \sqrt{\frac{P}{P_{0}} \frac{\alpha}{2+\alpha}}$     (25)

If we denote equilibrium constant at a certain temperature T is K*, we then have, 

$K^{*}=\frac{\alpha}{1-\alpha} \sqrt{\frac{P}{P_{0}} \frac{\alpha}{2+\alpha}}$     (26)

The value of K* only depends on temperature, not on pressure. Most engineering thermodynamics have tables for values of equilibrium of various chemical reactions; for example, see [8, 10]. 

We can remove the square root by squaring both sides in Eq. (26).

$\begin{aligned}

&K^{* 2}=\left(\frac{\alpha}{1-\alpha}\right)^{2} \frac{P}{P_{0}} \frac{\alpha}{2+\alpha} \\

&

\end{aligned}$     (27)

By taking $f=P /\left(P_{0} K^{* 2}\right)$ and rearranging Eq. (27), we then have:

$(1-f) \alpha^{3}-3 \alpha+2=0$     (28)

If f≠1, we can divide all terms in Eq. (28) by (1-f). We then have a depressed cubic equation where the coefficient of the second term is zero, 

α3+pα+q=0    (29)

where, $v=-\frac{3}{1-f} \text { and } q=\frac{2}{1-f}$. If f=1, the cubic term in Eq. (28) vanishes. The solution is very simple: $\alpha=\frac{2}{3}$. 

It must be noted that in actual situation, there may be more than one equation that must be solved. In the given example, O2 will also dissociate into oxygen atom. We will then have two nonlinear equilibrium reactions which are more difficult to solve. Cengel et al. [10] provides an example for solving this type of problem.

3.1 Simply supported beam loaded by continuous load

Figure 1 shows a beam supported at the ends and loaded by a continuous load q. 

Figure 1. Simply supported beam loaded by continuous load

The deflection curve is given by the Euler-Bernoulli equation which relates the deflection with bending moment (M_x), modulus of elasticity E and moment of inertia of the cross section of the beam (I), as seen in most textbooks on strength of materials. Most textbooks take the deflection upward as positive downward as negative; see [5, 17-22] for example. Timoshenko [23] adopts the opposite ones because most deflections in engineering are downward. Following [23], we will take the deflection downward as positive so the Bernoulli equation becomes

$E I \frac{d^{2} Y}{d x^{2}}=-M_{x}$    (30)

Integrating Eq. (30) will yield the deflection curve along the beam.

Figure 1 shows three regions where the bending moment can be applied: from 0 to a, from a to a+b and from a+b to L. Since there are three regions, there will be also three equations required for the bending moment along the length of the beam. However, we can use only one equation by adopting using singularity functions; see [21-22, 24], for example. We will use a singularity function in this example because it is simpler.

Without derivation, the equation for the bending moment along the beam for 0≤x≤L is given by: 

$M_{x}=R_{A} x-\frac{q}{2}\langle x-a\rangle^{2}+\frac{q}{2}\langle x-a-b\rangle^{2}$     (31)

Here, RA is the reaction on the left end of the beam. The value <x-a> becomes zero if x-a<0 but is equal to x-a if x>a. Likewise, <x-a-b> becomes zero if x-a-b<0 but is equal to x-a-b if x>a+b. From Eq. (30) we then have,

$E I \frac{d^{2} Y}{d x^{2}}=-R_{A} x+\frac{q}{2}\left\{\langle x-a\rangle^{2}-\langle x-a-b\rangle^{2}\right\}$     (32)

Integrating Eq. (32) twice produces:

$E I \frac{d Y}{d x}=-R_{A} \frac{x^{2}}{2}+\frac{q}{6}\left\{\langle x-a\rangle^{3}-\langle x-a-b\rangle^{3}\right\}+C_{1}$     (33)

$\begin{gathered}

E I Y=-R_{A} \frac{x^{3}}{6}+\frac{q}{24}\left\{\langle x-a\rangle^{4}-\langle x-a-b\rangle^{4}\right\}+ \\

C_{1} x+C_{2}

\end{gathered}$     (34)

Here, dY/dx is the slope of the deflection curve. C1 and C2 are boundary conditions. At the supports, deflections are zero. So, we then have,

$X_{x}=0, Y_{A}=0 \rightarrow C_{1}=0$

$X_{x}=L, Y_{D}=0 \rightarrow C_{1}=0$, which yields:

$C_{1}=R_{A} \frac{L^{2}}{6}-\frac{q}{24 L}\left[(L-a)^{4}-(L-a-b)^{4}\right]$       (35)

RA can be found from the conditions that at the supports, bending moments are equal to zero. At x=L, MB=0. So, we then have:

$\begin{gathered}

R_{A}=\frac{q}{2}\left[(L-a)^{2}-(L-a-b)^{2}\right]=q b(L-a- \\

b / 2)

\end{gathered}$     (36)

$E I Y^{\prime}=\frac{q(x-a)^{3}}{6}-\frac{q b^{2} x^{2}}{4 L}+\frac{q b^{2}\left(2 L^{2}-b^{2}\right)}{24 L}$    (37)

Derivations of all formulas in this part can be done by using formulas found in most mechanics of materials (engineering mechanics or strength of materials) textbooks such as [17-23]. Very complete formulas for various kinds of loads are given [24].

The term inside the square bracket is zero for x ≤ a. Substituting a=L-b and arranging the terms yields:

$\begin{gathered}

E I Y^{\prime}=\frac{q}{24 L}\left\{4 L(x-L+b)^{3}-6 b^{2} x^{2}+b^{2}\left(2 L^{2}-\right.\right. \\

\left.\left.b^{2}\right)\right\}

\end{gathered}$     (38)

Now, we substitute b=mL and x=pL where 0<m≤1 and 0≤p≤1 to Eq. (38). After simplification and arrangement of the terms we then have,

$\begin{gathered}

Y^{\prime}=\frac{q L^{3}}{24 E I}\left\{4(p-1+m)^{3}-6 m^{2} p^{2}+m^{2}(2-\right. \\

\left.\left.m^{2}\right)\right\}

\end{gathered}$     (39)

Terms inside the curly brackets are dimensionless. Eq. (39) can be written as:

$Y^{\prime}=f_{1} f(p)$     (40)

where,

$\begin{gathered}

f_{1}=q L^{3} / 24 E I \text { and } f(p)=4(p-1+m)^{3}- \\

6 m^{2} p^{2}+m^{2}\left(2-m^{2}\right)

\end{gathered}$     (41)

The equation for the deflection along the length of the beam for 0≤x≤L is given by:

$E I Y=\frac{q(x-a)^{4}}{6}-\frac{q b^{2} x^{3}}{12 L}+\frac{q b^{2} x\left(2 L^{2}-b^{2}\right)}{24 L}$     (42)

Substituting b=mL and x=pL to Eq. (42) results in:

$\begin{gathered}

E I Y=\frac{q}{24 L}\left\{L(x-L+b)^{4}-2 b^{2} x^{3}+b^{2} x\left(2 L^{2}-\right.\right. \\

\left.\left.b^{2}\right)\right\}

\end{gathered}$     (43)

which can be transformed to: 

$Y=\frac{q L^{4}}{24 E I}\left\{(p-1+m)^{4}-2 m^{2} p^{3}+m^{2} p\left(2-m^{2}\right)\right\}$     (44)

Again, Terms inside the curly brackets are dimensionless. Eq. (44) can be written as:

$Y=f_{2} g(p)$     (45)

where,

$\begin{gathered}

f_{2}=q L^{4} / 24 E I \text { and } g(p)=(p-1+m)^{4}- \\

2 m^{2} p^{3}+m^{2}\left(2-m^{2}\right) p

\end{gathered}$     (46)

The maximum deflection along the beam can be found by setting Eq. (38) to zero, which is similar to set f(p) in Eq. (39) to zero. Since f(p) is cubic, it has three roots. However, we are only interested in the real root which lies between 1-m and 1, which is similar to a<x<L.

3.2 Beam with built-in ends loaded by a triangular load 

Figure 2. A Beam with built-in ends loaded by a triangular load

Figure 2 shows a beam with built-in ends loaded by a triangular load; the load increases from zero in the left to the maximum, W0, in the right. The beams are clamped at both ends. Without derivation, the equations for the bending moment, the slope and the deflection along the length of the beam are given by Eqns. (47), (48) and (49), respectively; see any textbooks on engineering mechanics for the derivation. Since there is only one region in Figure 2, we do not necessarily use singularity functions.

$M_{x}=\frac{W_{0}}{60}\left(10 x^{3}-9 L^{2} x+2 L^{3}\right)$     (47)

$Y^{\prime}=\frac{W_{0}}{120 E I}\left(5 x^{4}-9 L^{2} x^{2}+4 L^{3} x\right)$     (48)

$Y^{\prime \prime}=\frac{W_{0}}{120 E I}\left(x^{5}-3 L^{2} x^{3}+2 L^{3} x^{2}\right)$    (49)

Now, we substitute x=pL where 0≤p≤1 to Eqns. (47), (48) and (49), respectively. After simplification, we then have:

$M_{p}=\frac{W_{0} L^{3}}{60}\left(10 p^{3}-9 p+2\right)$     (50)

$Y^{\prime}=\frac{W_{0} L^{4}}{120 E I}\left(5 p^{4}-9 p^{2}+4 p\right)$     (51)

$Y^{\prime \prime}=\frac{W_{0} L^{5}}{120 E I}\left(p^{5}-3 p^{3}+2 p^{2}\right)$     (52)

Maximum slope (Eq. (51)) can be found by taking the first derivative of the slope to be zero; this is equal to Eq. (50), which is a cubic equation in term of p. Since Mp is a cubic equation, there are three roots found. However, we are only interested in positive real roots.

4. Solving Cubic Equations

A cubic equation $a x^{3}+b x^{2}+c x+d=0$ can be solved manually using a scientific calculator such as Casio fx-991 ES. This calculator can solve a cubic equation easily. The user just inputs the coefficients of the equation. The answer will then be displayed on the screen. The use of a calculator is very advantageous for the students who need to solve a cubic equation quickly in their study. However, for repeated uses such as in simulation the use of the calculator is very awkward and time consuming. Moreover, if the coefficients of the cubic equation are not simple and must be calculated first (such as in the SRK equation and the PR equation), the use of the calculator may not be convenient since those coefficients must be written down before they can be entered for calculation.

A cubic equation can also be solved using Excel, a very good spreadsheet. Maple, Mathematica, Mathcad or other software can also be used. Kusalaas [25] uses Python 3 to compute the roots of a polynomial. We can also write a C program to solve it. Software and programming languages are better if there are many equations to be solved repeatedly. We can use exact formulas for solving the equation or solve it numerically. Usually we use Newton-Raphson method to find the roots, as shown in Refs. [26, 27]. However, it is not often easy to find the guess value for the method. At high temperature, a approaches 1; this makes Eq. (28) difficult to solve numerically. Slight change of the initial value makes the equation unstable. We will discuss how to find the exact roots of a cubic equation and to compute the roots numerically.

4.1 Exact roots of a cubic equation

The cubic equation in Eq. (1), without loss of generality, we will take a to be 1 so Eq. (1) becomes,

$x^{3}+b x^{2}+c x+d=0$        (53)

Substituting $x=Z-\frac{b}{3}$ to Eq. (53) results in depressed cubic equation where the coefficient of the quadratic term is zero,

$z^{3}+p z+q=0$        (54)

where,

$p=\frac{3 c-b^{2}}{3}, q=\frac{2 b^{2}-9 b c+27 d}{27}$

Before solving Eq. (54), we will observe the nature of the roots without determining them. Let the roots be y1, y2 and y3, respectively. Define the determinant of Eq. (54) as:

$\Delta=\frac{q^{2}}{4}+\frac{p^{3}}{27}$

We can then characterize the roots as follow:

  • Case 1: Δ>0; one root is real and two others are complex
  • Case 2: Δ=0; three real roots, two of them are the same
  • Case 3: Δ<0; three distinct roots

The solution of Eq. (54) is then given as follow:

Case 1: Δ>0; one root is real and two others are complex

Let

$t_{1}=-\frac{q}{2}+\sqrt{\Delta}, t_{2}=-\frac{q}{2}+\sqrt{\Delta}$

The roots are then given as follow:

$y_{1}=\sqrt[3]{t_{1}}+\sqrt[3]{t_{2}}$

$y_{2}=\omega \sqrt[3]{t_{1}}+\omega \sqrt[2]{t_{2}}$$=-\frac{1}{2}\left(\sqrt[3]{t_{1}}+\sqrt[3]{t_{2}}\right)+\frac{i \sqrt[3]{}}{2}\left(\sqrt[3]{t_{1}}-\sqrt[3]{t_{2}}\right)$

$y_{3}=\omega^{2} \sqrt[3]{t_{1}}+\omega \sqrt[3]{t_{2}}$$=-\frac{1}{2}\left(\sqrt[3]{t_{1}}+\sqrt[3]{t_{2}}\right)-\frac{i \sqrt[3]{}}{2}\left(\sqrt[3]{t_{1}}-\sqrt[3]{t_{2}}\right)$

Here, ω a primitive third root of unity.

Case 2: Δ=0; three real roots, two of them are the same

The roots are then given as follow:

$\begin{gathered}

y_{1}=2 \sqrt[3]{t_{1}}=-2 \sqrt[3]{\frac{q}{2}} \\

y_{2}=\left(\omega+\omega^{2}\right) 2 \sqrt[3]{t_{1}}=\sqrt[3]{\frac{q}{2}} \\

y_{3}=\left(\omega^{2}+\omega\right) 2 \sqrt[3]{t_{1}}=\sqrt[3]{\frac{q}{2}}

\end{gathered}$

Case 3: Δ<0; three distinct roots

Let

$u=\frac{3 q}{2 p} \sqrt{-\frac{3}{p}}$

The roots are then given as follow:

$y_{1}=\sqrt{-\frac{4 p}{3}} \cos \left(\frac{1}{3} \cos ^{-1} u\right)$

$y_{2}=\sqrt{-\frac{4 p}{3}} \cos \left(\frac{1}{3} \cos ^{-1} u+\frac{2 \pi}{3}\right)$

$y_{3}=\sqrt{-\frac{4 p}{3}} \cos \left(\frac{1}{3} \cos ^{-1} u+\frac{4 \pi}{3}\right)$

See [7, 28-30] for the derivation of all the afore-mentioned formulas. Strobach [31, 32] also discuss polynomial fitting for solving a cubic equation.

The exact roots of the cubic equation are attributed to Cardano or Tartaglia; see the Ref. [33] for an interesting discussion about the dispute on who first discovered the Eq. (26) provides a good discussion on Cardano’s formula. 

4.2 Numerical methods for solving cubic equations

Solving cubic equations using exact roots sometimes becomes awkward. Even simple coefficient of the equations might yield roots which are not easy to churn out because of cubic root calculations. The study [30] shows an example when solving $x^{3}+6 x-40=0$. Using the cubic formula, one of the root equals to $\sqrt[3]{20+\sqrt{392}}+\sqrt[3]{20-\sqrt{392}}$. The actual root is 4, which can be verified by substituting the value to the equation. At first, it is not clear that the sum of the two radicals equals to 4. However, it can be shown that $20+\sqrt{392}=(2+\sqrt{2})^{3} \text { and } 20-\sqrt{392}=(2-\sqrt{2})^{3}$. Here, we can resort to numerical methods to solve a cubic equation. There are many methods available but the popular one is Newton-Raphson method. This method is iterative. It will stop when the absolute error between two successive values is smaller than a specified value or the number of iterations is greater than a specified number.

Suppose that we want find the roots of a cubic equation as given by Eq. (1). The algorithm for solving the equation is as follow:

  1. $f(x) \equiv a x^{3}+b x^{2}+c x+d=0$  
  2. $f^{\prime}(x)=3 a x^{2}+2 b x+c$  
  3. Guess an initial value $x_{0}$
  4. $x_{n+1}=x_{n}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}$  
  5. $e r r=\left|x_{n+1}-x_{n}\right|$  
  6. If $e r r<a$  specified value (for example, 10-5 ) or $n>$ greater than a specified number (say 20), stop. Otherwise, continue the iteration.

The description of Newton-Raphson method is discussed in most numerical analysis or engineering mathematics textbooks, as shown in Refs. [26, 27].

After one of the roots is found, say x1, Eq. (1) is deflated that is dividing it by (x ˗ x1), yielding a quadratic equation, which can be very easily solved. For the applications discussed in this paper, we are only interested in real roots.

Many methods have been proposed for solving cubic equations of state numerically; see [14] for a good discussion on the methods. [Olivera-Fuentes] discusses optimal solution for cubic equations of state. Guessing an initial value is often difficult. Valderrama [14] and Olivera-Fuentes [34] provide methods for the initial value. Sometimes a simple initial value RT/p is useful.

For solving equilibrium reactions, it is better to find the exact root because it is difficult to guess the initial value as demanded by the Newton-Raphson method. At high temperature, a approaches 1. For the simply supported beam shown in Figure 1 or the beam with built-in ends loaded by a triangular load shown in Figure 2, initial values can be taken to be in the middle of the beam.

5. Conclusions

Six applications of cubic equations in engineering have been discussed. Only real roots are valid even though there are three roots in the equations. Often only one real root is used; the other two are discarded. Cubic equations can be solved exactly using Cardano’s formula. Solving cubic equations using exact roots sometimes becomes awkward. While Newton-Raphson method is usually used to solve a cubic equation numerically, it is often not easy to determine an initial value for dissociation reaction as shown in Eq. (28). At high temperature, a approaches 1. Slight change of the initial value makes the equation unstable There are various ways for solving equation, from using a scientific calculator to programming languages. Which one used depends on the need of solving the equation.

  References

[1] Gilbert, L., Gilbert, J. (2015). Elements of Modern Algebra 8th Ed. Stamford: Cengage Learning.

[2] Prasolov, V.V. (2004). Polynomials, 11. Springer Science & Business Media. 

[3] Stewart, I. (1990). Galois Theory. Chapman and Hall/CRC. https://doi.org/10.4324/9780203489307

[4] Tignol, J.P. (2015). Galois' Theory of Algebraic Equations. World Scientific Publishing Company. Co Pte Ltd.

[5] Bewersdorff, J. (2021). Galois theory for beginners: A historical perspective. American Mathematical Society. 

[6] Strobach, P. (2016). How Do I Solve A Cubic Equation. Technical Report. Accessed from https://www.researchgate.net/publication/300009555.

[7] Poling, B.E., Prausnitz, J.M., O’Connell, J.P. (2001). The Properties of Gases and Liquids, 5th Ed. New York: McGraw-Hill.

[8] Borgnakke, C., Sonntag, R.E. (2013). Fundamentals of Thermodynamics 8th Ed. Hoboken: John Wiley and Sons, Inc. https://doi.org/10.1201/9781003224044

[9] Dahm, K.D., Visco, D.P. (2014). Fundamentals of Chemical Engineering Thermodynamics. Cengage Learning. 

[10] Cengel, Y.A., Boles, M.A., Kanoğlu, M. (2011). Thermodynamics: An Engineering Approach, 5: 445. New York: McGraw-Hill. 

[11] Zhao, W., Xia, L., Cao, X., Bi, R., Xiang, S. (2020). A modified alpha function for Peng-Robinson cubic equation of state. Chemical Engineering Transactions, 81: 547-552. https://doi.org/10.3303/CET2081092

[12] Shivakumar, B., Ramesh, A., Vaishnavi, P., Avi, A.C., Limperich, D. (2013). Retrospective on Cubic Equation of State for R134a Refrigerant Used in Automotive Application. SIAT, India. 

[13] Deiters, U.K., Macías-Salinas, R. (2014). Calculation of densities from cubic equations of state: Revisited. Industrial & Engineering Chemistry Research, 53(6): 2529-2536. https://doi.org/10.1021/ie4038664

[14] Valderrama, J.O. (2003). The state of the cubic equations of state. Industrial & Engineering Chemistry Research, 42(8): 1603-1618. https://doi.org/10.1021/ie020447b

[15] Ferguson, C.R., Kirkpatrick, A.T. (2015). Internal Combustion Engines: Applied Thermosciences. John Wiley & Sons

[16] Turns. S.R. (2012). An Introduction to Combustion Concepts and Applications 3rd Ed. New York: McGraw-Hill.

[17] Muvdi, B.B., Elhouar, S. (2016). Mechanics of Materials: with Applications in Excel. CRC Press. https://doi.org/10.1201/9781315374314

[18] Goodno, B.J., Gere, J.M. (2020). Mechanics of Materials. Cengage Learning. 

[19] Gross, D., Ehlers, W., Wriggers, P., Schröder, J., Müller, R. (2016). Mechanics of Materials-Formulas and Problems, 2-113. Springer.

[20] Mott, R.L., Untener, J.A. (2017). Applied Strength of Materials, SI Units Version. CRC Press. https://doi.org/10.1201/9781315153056

[21] Philpot, T.A. (2008). Mechanics of Materials: An Integrated Learning System. Hoboken: John Wiley and Sons.

[22] Nash, W.A., Potter, M.C. (2011). Schaum’s Series Outlines Strength of Materials, 5th Ed. New York: McGraw-Hill.

[23] Timoshenko, S. (1948). Strength of Materials Part I: Elementary Theory and Problems 2nd Ed. 10th printing. New York: D. Van Nostrand Company, Inc.

[24] Young, W.C., Budynas, R.G., Sadegh, A.M. (2002). Roark's Formulas for Stress and Strain. New York: McGraw-Hill.

[25] Kusalaas, J. (2013). Numerical Methods in Engineering with Python 3. New York: Cambridge University Press.

[26] Surana, K.S. (2018). Numerical Methods and Methods of Approximation: in Science and Engineering. CRC Press. https://doi.org/10.1201/9780429028281

[27] Esfandiari, R.S. (2017). Numerical Methods for Engineers and Scientists Using MATLAB® 2nd Ed. Boca Raton: CRC Press.

[28] Grant, H., Kleiner, I. (2015). Turning points in the history of mathematics. Basel: Birkhäuser.

[29] Rotman, J.J. (2015). Advanced modern algebra, 165. American Mathematical Soc.

[30] Confalonieri, S. (2015). The Unattainable Attempt to Avoid the Casus Irreducibilis for Cubic Equations: Gerolamo Cardano's De Regula Aliza. Springer. Fachmedien.

[31] Strobach, P. (2014). A post-fitting algorithm for high precision complex polynomial root finding. Universal Journal of Applied Mathematics and Computation, 2: 16-34.

[32] Strobach, P. (2011). Solving cubics by polynomial fitting. Journal of Computational and Applied Mathematics, 235(9): 3033-3052. https://doi.org/10.1016/j.cam.2010.12.025

[33] Burton, D.M. (2011). The History of Mathematics: An Introduction 7th Ed. New York: McGraw-Hill.

[34] Olivera-Fuentes, C. (1993). The optimal solution of cubic equations of state. Lat. Am. Appl. Res. 23: 243-256.