Chemical & Biological Engineering CBE Computational Turbulent Reacting Flow CTRFlow

Publications Calendar Research Teaching Zoom LDS

Publications Calendar Research Teaching Zoom LDS

Solving equations

Arrays (minimal)

• Solvers make use of Python arrays
• Arrays are collections of numbers, similar to mathematical vectors
• Arrays are defined as comma-separated numbers contained in square brackets:
  • a = [1.1,-85,3.14,101325]
• Array elements are assigned and accessed using square brackets:
  • Array elements are indexed (counted), starting at 0
  • a[0] is 1.1
  • a[1] is -85
  • a[2] is 3.14
  • a[3] is 101325
• Change elements:
  • a[1] = -8314.46

a = [1.1,-85,3.14,101325]
print(a[0], a[1], a[2], a[3])
print(a)

a[0] = 8314.46
a

1.1 -85 3.14 101325
[1.1, -85, 3.14, 101325]

[8314.46, -85, 3.14, 101325]

2D Arrays (matrices)

• A matrix is an array of arrays:
  • A = [[1,2,3], [4,5,6], [7,8,9]] $$A = \begin{bmatrix}1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9\end{bmatrix}$$
• And accessed as $A_{i,j} = $ A[i][j] for row i and column j

A = [[1.1,2.2,3.3],[4.4,5.5,6.6],[7.7,8.8,9.9]]
print(A)
A[1][2]

[[1.1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]]

6.6

Numpy arrays

• The above are Python's built-in arrays.
• It is generally better to use Numpy arrays for numerical work.
• The syntax is basically the same, except we create arrays using the np.array function.
  • (More on this later.)
• Note, for 2D numpy arrays, we access and assign with A[i,j], as opposed to A[i][j] for regular Python arrays.

import numpy as np

b = np.array([1.1, 2.2, 3.3])
A = np.array([[1.1, 2.2, 3.3],[4.4,5.5,6.6],[7.7,8.8,9.9]])

print(b)
print(A[1,2])

[1.1 2.2 3.3]
6.6

Boilerplate

• import libraries

import numpy as np
from scipy.optimize import root

Solve Linear equations

To solve a system $Ax=b$ for x, we do the following:

• Define the matrix $A$ and the array $b$.
• Use x = np.linalg.solve(A,b)

Example

Solve

$$\begin{bmatrix} -2 & 1 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 1 & -2 \end{bmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1\end{pmatrix}$$

A = np.array([ [ -2, 1, 0, 0, 0 ],
               [ 1, -2, 1, 0, 0 ],
               [ 0, 1, -2, 1, 0 ],
               [ 0, 0, 1, -2, 1 ],
               [ 0, 0, 0, 1, -2 ] ])

b = np.array([1, 1, 1, 1, 1])

x = np.linalg.solve(A, b)

#--------- verify

print("Ax = ", np.dot(A,x))
print("Ax = ", A@x)

Ax =  [1. 1. 1. 1. 1.]
Ax =  [1. 1. 1. 1. 1.]

Here's a link to more linear algebra operations.

Nonlinear equations

1-D equations (1 equation in 1 unknown).

Recall

• What are the key elements to solving nonlinear equations?
• List the parts, the form, and things you need.

Exercise

• Before learning the details about solving equations using Python, see if you can guess how to solve the following equation using Python's solver root, based on what you need to solve equations: $$f(x) = x^2-5.$$

Solving a nonlinear equation involves first putting it in the form $f(x)=0$. We know the function $f$ and we want to find the value of $x$ that gives $f(x)=0$.

To solve this, do the following:

• Define the function f(x)
• Set an initial guess for x.
• We need to import a library: from scipy.optimize import root
• Call the root function: res = root(f, x0)
  • f is the name of the function you are solving, and x0 is your initial guess.
  • res stores a bunch of result information. Access the actual answer with res.x

Exercise

Solve for both roots of $$f(x) = x^2-5.$$

def f(x):
    return x*x - 5

xguess = 1.0
res = root(f, xguess)
print(res.x)

xguess = -1.0
res = root(f, xguess)
print(res.x)

[2.23606798]
[-2.23606798]

Multiple equations in multiple unknowns

• What do we do if we are solving multiple equations in multiple unknowns?

• For example, suppose we have $h(y,z)=0$ and $g(y,z)=0$. Such as: $$h(y,z) = y+2z,$$ $$g(y,z) = \sin(y)/z.$$

• Do you think Python's root function can handle this?

• How will you go about specifying the problem?

• What if there were three equations in three unknowns? Four and four? etc.?

• If you were coding root how would you account for all these possibilities?

General approach

• Instead of $f(x)=0$, we solve
• $\vec{f}(\vec{x})=\vec{0}$.
  • That is, $\vec{x}$ is a vector of unknowns, and $\vec{f}$ is a vector of functions.
  • $\vec{f}(\vec{x})$ takes a vector $\vec{x}$, and returns a vector $\vec{f}$.

Details

• Consider again the two functions for the two unknowns: $h(y,z) = 0$, and $g(y,z)=0$. $$h(y,z)=y+2z,$$ $$g(y,z) = \sin(y)/z$$

• For computer solvers, we write using a generic, standard notation.

  • We want all solvers to look like $f(x)=0.$

• We would write these as $$f_0(x_0,x_1) = x_0+2x_1,$$ $$f_1(x_0,x_1) = \sin(x_0)/x_1,$$ which becomes: $$\left[\begin{matrix} f_0(x_0,x_1) \\ f_1(x_0,x_1)\end{matrix}\right] = \left[\begin{matrix}x_0+2x_1 \\ \sin(x_0)/x_1\end{matrix}\right] = \left[\begin{matrix}0 \\ 0\end{matrix}\right],$$ or $$\vec{f}(\vec{x}) = \left[\begin{matrix}x_0+2x_1 \\ \sin(x_0)/x_1\end{matrix}\right] = \vec{0},$$ that is, $$\vec{f}(\vec{x})=\vec{0}.$$

• We solve the system just like before:

  • Just define the vector function. It takes an array of unknowns $\vec{x}$, and returns an array of function values $\vec{f}$.
  • We use the same solver, but we give an array of initial guesses $\vec{x}_0$.

Note Here is a link to more information about the root function with extra arguments that can give more control over the solution process.

Exercise

• Solve the following two equations in two unknowns.

$$h(y,z)=y+2z,$$

$$g(y,z) = \sin(y)/z.$$

• Use $y_0=1$, and $z_0=2$ as an initial guess.

def hg(yz):
    y = yz[0]
    z = yz[1]

    h = y + 2*z
    g = np.sin(y)/z

    return np.array([h,g])

#------------

xy_guess = np.array([1,2])

res = root(hg, xy_guess)
print("x =", res.x[0], "y =", res.x[1])

x = 9.424777960769402 y = -4.712388980384701

Excercise

Solve the following system of three equations in three unknowns:

$$x^2 + y^2 = 1,$$ $$xy + yz = -1.1,$$ $$y^2 + z^2 = 2.$$

A reasonable guess for all variables is $x=y=z=2.$

def F(xyz):
    x = xyz[0]
    y = xyz[1]
    z = xyz[2]

    eq0 = x*x + y*y - 1
    eq1 = x*y + y*z + 1.1
    eq2 = y*y + z*z - 2

    return np.array([eq0, eq1, eq2])

#-----------------

xyz_guess = np.array([2,2,2])
res = root(F, xyz_guess)
print("x,y,z = ", res.x)

x,y,z =  [ 0.10056089 -0.99493091  1.00504353]

Exercise

The Colbrook equation is used to compute pressure drops across pipes, and therefore to compute pumping requirements and costs.

The Colebrook equation is given by $$\frac{1}{\sqrt{f}} =-2\log_{10}\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right),$$ or, $$F(f, Re, \epsilon/D) = \frac{1}{\sqrt{f}} + 2\log_{10}\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right)=0.$$

Solve for $f$ given $Re=20000$, and $\epsilon/D=0.01$. Use 0.01 for a guess for $f$.

Re     = 20000
epsD   = 0.001
fguess = 0.01

#------------

def F(f):
    return 1.0/np.sqrt(f) + 2.0*np.log10(epsD/3.7 + 2.51/Re/np.sqrt(f))

f = root(F, fguess).x
print("f =", f)

f = [0.02794571]