Sympy (sympy.org) is a Python package used for solving equations with symbolic math.

Using Python and SymPy we can write and solve equations that come up in Engineering.

The example problem below contains two equations with two unknown variables. You could use a pencil and paper to solve the problem, but we are going to use Python and programming to solve the problem.

Given:

The density of two different samples of a polymer $\rho_1$ and $\rho_2$ are measured.

$$ \rho_1 = 0.904 \ g/cm^3 $$$$ \rho_2 = 0.895 \ g/cm^3 $$

The percent crystallinity of the two samples ($\%c_1 $ and $\%c_2$) is known.

$$ \%c_1 = 62.8 \% $$$$ \%c_2 = 54.4 \% $$

The percent crystalinity of a polymer sample is related to the density of 100% amorphus regions ($\rho_a$) and 100% crystaline regions ($\rho_c$) according to:

$$ \%crystallinity = \frac{ \rho_c(\rho_s - \rho_a) }{\rho_s(\rho_c - \rho_a) } \times 100 \% $$

Find:

Find the density of 100% amorphus regions ($\rho_a$) and the density of 100% crystaline regions ($\rho_c$) for this polymer.

Solution:

We are going to use Python and a package called SymPy to solve this problem. I recommend installing the Anaconda distribution of Python. If you install Anaconda, SymPy is included. If you downloaded Python from Python.org or if you are using a virtual environment, SymPy can be installed at a terminal using pip with the command below.

$ pip install sympy

We need a couple of functions from the SymPy package to solve this problem. We need the symbols() function to create symbolic math variables for the density of 100% amorphous and 100% crystalline regions ($\rho_a$ and $\rho_c$) and variables for the given information in the problem ($\%c_1 $, $\%c_2$, $\rho_1$ and $\rho_2$ ). We also need SymPy's nonlinsolve() function to solve a system of non-linear equations.

The symbols() function and the nonlinsolve() function can be imported from SymPy using the line below.


In [1]:
from sympy import symbols, nonlinsolve

Next we need to define six different variables:

$$\rho_c, \rho_a, \rho_1, \rho_2, c_1, c_2$$

Note commas are included in the symbols output, but there are no commas in the symbols input.


In [2]:
pc, pa, p1, p2, c1, c2 = symbols('pc pa p1 p2 c1 c2')

Now we can create two SymPy expressions that represent our two equations. We can subtract the %crystallinity from the left side of the equation to set the equation to zero. The result of moving the %crystallinity term to the other side of the equation is shown below. Note how the second equation equals zero.

$$ \%crystallinity = \frac{ \rho_c(\rho_s - \rho_a) }{\rho_s(\rho_c - \rho_a) } \times 100 \% $$$$ \frac{ \rho_c(\rho_s - \rho_a) }{\rho_s(\rho_c - \rho_a) } \times 100 \% - \%crystallinity = 0 $$

Substitue in $\rho_s = \rho_1$ and $\rho_s = \rho_2$ into the expression above. Also substitue in $\%crystallinity = \%c_1$ and $\%crystallinity = \%c_2$. The result is two equations, each equation is equal to zero.

$$ \frac{ \rho_c(\rho_1 - \rho_a) }{\rho_1(\rho_c - \rho_a) } \times 100 \% - \%c_1 = 0 $$$$ \frac{ \rho_c(\rho_2 - \rho_a) }{\rho_2(\rho_c - \rho_a) } \times 100 \% - \%c_2 = 0 $$

Now we have two equations (the two equations above) which we can solve for two unknowns ($\rho_a$ and $\rho_s$). The two equations can be coded into SymPy expressions. The SymPy expressions contain the variables we defined earlier.


In [3]:
expr1 = ( (pc*(p1-pa)   ) / (p1*(pc-pa)) - c1)
expr2 = ( (pc*(p2-pa)   ) / (p2*(pc-pa)) - c2)

Next, we'll substitute in the known values $\rho_1 = 0.904$ and $c_1 = 0.628$ into our first expression expr1. Note you need to set the output of SymPy's .subs method to a variable. SymPy expressions are not modified in-place. You need to capture the output of the .subs method in a variable.


In [4]:
expr1 = expr1.subs(p1, 0.904)
expr1 = expr1.subs(c1, 0.628)
print(expr1)


1.10619469026549*pc*(-pa + 0.904)/(-pa + pc) - 0.628

Now we'll substitue the second set of given values $\rho_2 = 0.895$ and $c_2 = 0.544$ into our second expression expr2.


In [5]:
expr2 = expr2.subs(p2, 0.895)
expr2 = expr2.subs(c2, 0.544)
print(expr2)


1.11731843575419*pc*(-pa + 0.895)/(-pa + pc) - 0.544

We'll use SymPy's nonlinsolve() function to solve the two equations expr1 and expr2 for to unknows pa and pc. SymPy's nonlinsolve() function expects a list of expressions [expr1,expr2] followed by a list variables [pa,pc] to solve for.


In [6]:
sol = nonlinsolve([expr1,expr2],[pa,pc])
print(sol)


{(0.840789786223278, 0.946134313397929)}

We see that the value of $\rho_a = 0.84079$ and $\rho_c = 0.94613$.

The solution is a SymPy FiniteSet object. To pull the values of $\rho_a$ and $\rho_c$ out of the FiniteSet, use the syntax sol.args[0][<var num>] to pull the answers out.


In [7]:
print(type(sol))


<class 'sympy.sets.sets.FiniteSet'>

In [8]:
pa = sol.args[0][0]
pc = sol.args[0][1]
print(f' Density of 100% amorphous polymer, pa = {round(pa,2)} g/cm3')
print(f' Density of 100% crystaline polymer, pc = {round(pc,2)} g/cm3')


 Density of 100% amorphous polymer, pa = 0.84 g/cm3
 Density of 100% crystaline polymer, pc = 0.95 g/cm3

Use SymPy to calculate a numerical result

Besides solving equations, SymPy expressions can also be used to calculate a numerical result. A numerical result can be calculated if all of the variables in an expression are set to floats or integers.

Let's solve the following problem with SymPy and calculate a numerical result.

Given:

The density of a 100\% amorphous polymer sample $\rho_a$ and the density of a 100% crystaline sample $\rho_c$ of the same polymer are measured.

$$ \rho_a = 0.84 \ g/cm^3 $$$$ \rho_c = 0.95 \ g/cm^3 $$

The density of a sample $\rho_s$ of the same polymer is measured.

$$ \rho_s = 0.921 \ g/cm^3 $$

Find:

What is the \% crytallinity of the sample with a measured density $ \rho_s = 1.382 \ g/cm^3 $?

Solution

We have precise values for $ \rho_a $ and $ \rho_c $ from the previous problem. Let's see what the values of $ \rho_a $ and $ \rho_c $ are. We will use these more precise values that we calculated earlier to solve the problem.


In [9]:
print(pa)


0.840789786223278

In [10]:
print(pc)


0.946134313397929

Next, we will create three SymPy symbols objects. These three symbols objects will be used to build our expression.


In [11]:
pc, pa, ps = symbols('pc pa ps')

The expression that relates % crystallinity of a polymer sample to the density of 100% amorphus and 100% crystalline versions of the same polymer is below.

$$ \%crystallinity = \frac{ \rho_c(\rho_s - \rho_a) }{\rho_s(\rho_c - \rho_a) } \times 100 \% $$

We can build a SymPy expression that represents the equation above using the symbols objects (variables) we just defined.


In [12]:
expr = ( pc*(ps-pa)   ) / (ps*(pc-pa))

Now we can substitute our $ \rho_a $ and $ \rho_c $ from above. Note the SymPy's .subs() method does not modify an expression in place. We have to set the modified expression to a new variable before we can make another substitution. After the substitutions are complete, we can print out the numerical value of the expression. This is accomplished with SymPy's .evalf() method.


In [13]:
expr = expr.subs(pa, 0.840789786223278)
expr = expr.subs(pc, 0.946134313397929)
expr = expr.subs(ps, 0.921)

print(expr.evalf())


0.782187477379657

As a final step, we can print out the answer using a Python f-string.


In [14]:
print(f'The percent crystallinity of the sample is {round(expr*100,1)} percent')


The percent crystallinity of the sample is 78.2 percent

Summary

In this post we solved two equations for two unknowns with Python and a package called SymPy. SymPy can be used to define variables, define equations and solve for variables. SymPy can also be used to calculate a numerical result, if all of the values in an expression are known. A summary of the SymPy functions and methods used in this post are below.

SymPy fuction Example Description
symbols() x, y, z = symbols('x y z') define math variables
expr = expr = 2*x + 3*y - 4*z create a symbolic math expression (define variables first)
expr.subs() expr = expr.subs(x,2) substitute a value into a variable (not in-place)
nonlinsolve() sol = nonlinsolve([expr1,expr2],[x,y]) solve a system of equations for a set of variables

In [ ]: