Case Study 1 – Solving algebraic equations

Key Concepts

• Solving nonlinear algebraic equations
• Degrees of freedom
• Python packages
• Functions
• fsolve
• Iterating lists

Introduction

Solving a nonlinear equation is relatively common task for process and chemical engineers not only in technical roles, but also for supervisory, operating and strategic roles. I’ve seen people use the “guessing” method to hone in an answer and others pull out text books to remember how to find a algebraic solution to the equation. Finding numerical solutions to nonlinear equations should become second nature and we will look at how to do it here in a way that clearly lays out the solution so that you and others can review it later on.

Problem description

The problem used in this post has been taken from a collection (1st problem) originally distributed at the ASEE Chemical Engineering Summer School held in Snowbird, Utah on August 13, 1997. The problem requires calculation of different gas properties based on the van der Waals equation of state:

The problem requires us to solve for molar volume, V and you can see from this equation that it is not trivial to rearrange this equation and solve for V. Instead we will use a numerical equation solver to find a solution.

Before progressing to the code, we need to ensure that we have a properly defined problem. The van der Waals equation contains 6 variables:

• P = pressure (atm)
• V = molar volume (liters/g.mol)
• T = temperature (K)
• R = gas constant (atm litre/g.mol K)
• a, b = parameters to the equation

We are also given the following equations:

These equations introduce 2 additional variables:

• Pc = critical pressure (atm)
• Tc = critical temperature (K)

Degrees of freedom

In total we 8 variables and 3 equations, leaving us with 5 degrees of freedom. In other words we need define 5 further pieces of information to solve for V. In this case we are given the gas constant R (or we could have looked it up), we are told the gas is ammonia and actually given Pc and TC (or we could have looked them up) and for each part of the problem we are told the conditions T and P (or in some cases we are given the reduced pressure Pr).

So now we have a fully defined problem and are able to calculate V. However, we will also be required to calculate the compressibility factor, Z using the following equation:

This equation only introduces the new variable Z. This is 1 new equation with 1 new variable so there are zero degrees of the freedom and the problem remains fully defined.

What type of problem

Before continuing we need to think about what type of analytical problem do we have here as this will help us understand the core concepts required to find a solution, plus has implications for how to setup with solution approach in Python.

In this case, it is fairly simple as we have a number of nonlinear algebraic equations to be solved. Algebraic because there are no differential terms and nonlinear because the term V (i.e. the variable we need to be solving for) appears with a squared term and we can not rearrange the equations with V having only a power of 1.

So our general solution approach will involve:

• Define known variables
• Create a function containing all equations in difference form (i.e. LSH – RHS)
• Solve using a routine for finding the roots of equations (in this case fsolve)

The solution

I will show the entire solution to the problem before explaining some of the fundamental solution concepts. This output has been produced using python through a Jupyter notebook.

Python packages

The first concept introduced in the solution is the use of packages. Python is program language which ships with a well developed set of commands and functions. These provide a lot of functionality. But the real power of Python comes from the ability to pull in packages to perform an ever expanding range of functions. Fortunately, the python packages for data analytics are extensive and very mature.

The first block of code imports from 3 packages which get used extensively in data analytics:

• numpy – a very powerful library of mathmatical routines
• scipy – an extensive library of mathematical analysis routines
• matplotlib – a flexible charting and visualisation library

from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import numpy as np

I will right a separate post about importing libraries, but the key point to understand in this example is that:

• A procedure called fsolve is now available to the code
• The pyplot part of the larger matplotlib library is now available to the code and can be referenced through the prefix plt
• The numpy library is now available to the code and can be referenced through the prefix np

Defining functions

The next part of the code simply defines the supplied constants. Then we get to the following step:

def vdw_calc(V):

    a = 27.0 / 64.0 *(R**2 * Tc**2 / Pc)

    b = R*Tc / (8.0 * Pc)

    converg_error = (P + a/V**2)*(V-b) - R*T

    return converg_error

This code creates a function called vdw_calc which accepts a value for the molar volume V. It then uses the equations for a and b, utilising the values for all the parameters which will be defined before this function is called for the first time.

The actual van der Waals equation of state is rearranged so that the right hand side (RHS) is subtracted from the left hand side (LHS). The difference is called converg_error and our goal will be to find a solution which makes converg_error = 0. When this occurs, then we have a solution to the van der Waals equation. If the value for V supplied to the routine is not an exact solution, the converg_error represents how far off the van der Waals equation is from being solved. This value is returned and will be used later in the code.

Similarly a procedure is created for calculating the compressibility factor. This procedure is called comp_factor and does not receiving any values, instead using the pre-defined values for P, V, R and T to calculate the compressibility factor, Z. The value for Z is returned.

def comp_factor():

    Z = P*V / (R*T)

    return Z

Solving part a) – Solving the nonlinear equation

Part a) of the original problem was simply to calculate the molar volume, V and compressibility factor, Z at a specified T and P. The bit of code that does this is:

T = 450.0

P = 56.0


V = fsolve(vdw_calc, 1.0)[0]

Z = comp_factor()

The heart of this solution is using fslove (from scipy library) in combination with vdw_calc (our routine which takes a value for V and returns the error between the LHS and RHS of the van der Waals equation). fsolve accepts 2 bits of information. Firstly it needs to know the equation it needs to solve for (in this case vdw_calc) and a starting approximation for the solution (in this case we are guessing that the right answer for V is 1.0). fsolve then iterates different values for V (starting with 1.0) until vdw_calc returns 0.0.

Once fsolve has a solution it returns the correct value for V and assigns it to our parameter V. fsolve actually returns the answer as a list, just in this case it is a list with only a single number in it. We want V to contain just the number V and not the list. We extract the number by using the syntax [0], which literally means give me the first number in the list (python numbering always starts with 0 not 1).

Now that we have V and using the same values for T, P and R we calculate Z. The last 2 lines of code simply display the results with a little bit of formatting. I will not explain how this works for the moment,

Solving part b) – Solving the equation multiple times

Part b) of the original problem requires essentially the same calculation to be repeated for different values of pressure, P. In this case the pressures are given as reduced pressures, Pr and will need to converted to pressure according to the following equation:

Python provides a simple construct for storing information called a list. In this case we can store the list of reduced pressures we need to solve for in a list called Pr_list. Once this is defined, we can iterate through each value using a for loop:

Pr_list = [1,2,4,10,20]


for Pr in Pr_list:

    P = Pr * Pc

    V = fsolve(vdw_calc, 1.0)[0]

    Z = comp_factor()

In this case, the variable Pr takes on each value, one by one in the list and is then used to calculate the pressure, P, which we need for the van der Waals equation. There is no need to understand how many items are in Pr_list or how many times the for _ in loop has repeated. Python handles all of this seamlessly.

The rest of the code is essentially the same as part a), except it is repeated for each iteration.

Sense checking and interpretation

When ever solving analytical problems, always, always sense check your answer. There are many techniques you can employ depending on the situation. For example, in this case we can check the following:

• Results are all positive (check)
• Molar volumes are right order of magnitude (check)
• Compressibility factors are right order of magnitude (well actually I don’t know what to expect)
• Nonlinear form of the compressibility factor with pressure is expected (the term is consistent with high volumes at very low pressure – check)
• Manual check of some results using a calculator (check)

There is not a great deal to interpret here with out further context. However, I will revisit this problem in the future with opportunities to scale the approach for a more industrial environment.