Search code examples
pythonsympypint

Dimensional analysis in SymPy

I'm trying to use SymPy to do the following tasks:

  1. Simplify some algebraic equations of physical quantities and physical constants.
  2. Perform dimensional analysis to ensure my equations are correct
  3. Evaluate the equations by plugging in values for the physical quantities.

My Minimal, Reproducible Example for this problem is below, where I want to find the amount of mass equivalent to a given energy using mass-energy equivalence. (I run this in Jupyter since it can make the equations look nice.)

Staring with E=m*c**2 I want to solve for m:

from sympy import *
E, m, c = symbols('E m c') # symbols, for energy, mass, and speed of light
eq_E = Eq(E, m*c**2) # define equation for E in terms of m and c
eq_m = Eq(m,solve(eq_E,m)[0]) # solve equation for m
display(eq_m)

m = \frac{E}{c^2}

Great! Now it is time to perform dimensional analysis. I'll define my energy E in Joules, and I shouldn't need to define the units of c since it is a physical constant and I am working in SI units (which appears to be the default for SymPy). I want to find how much mass is in 1 Joule of energy, and I want to know the units of this mass. I'll define my variables using Quantity() and then set them up in an equation again so that SymPy can solve it:

from sympy.physics.units import Quantity, energy, joule, speed_of_light
from sympy import Eq, solve
m = Quantity('m')                    # Define unknown mass 'm' as a quantity
E = Quantity('E')                    # Define known energy 'E' as a quantity
E.set_dimension(energy)              # E is a quantity of energy
E.set_scale_factor(1.0*joule, 'SI')  # set energy to 1.0 Joules
eq_E = Eq(E,m*speed_of_light**2)     # define E = mc^2
eq_m = Eq(m,solve(eq_E,m)[0])        # Solve E = mc^2 for m
display(eq_m)

Quantity(m,m)=\frac{Quantity(E,E)}{Quantity(speed_{oflight},c)^2}

It doesn't look pretty, but the expression is correct. Now I want to see what the value of m is, and what units m has. First I'll save the solution for m in m_solve:

m_solve = solve(eq_E,m)[0]
display(m_solve)

\frac{Quantity(E,E)}{Quantity(speed_{oflight},c)^2}

Now, can I see what units m_solve has?

print(m_solve.dimension)

AttributeError: 'Mul' object has no attribute 'dimension'

What about the value of m_solve?

print(m_solve.scale_factor)

AttributeError: 'Mul' object has no attribute 'scale_factor'

What is a Mul? How do I get a Quantity? More generally, is there a good workflow in SymPy for manipulating symbolic equations and checking units along the way, and then finally evaluating those equations? If not SymPy, are there good alternatives? (I've tried a few alternatives, most promisingly pint, but it seems to lack symbolic support.)

Solution

  • I made a python package called SymDim (Symbolic Dimensional analysis) to do this using SymPy and astropy.units.

    Install SymDim (which will also install SymPy, Astropy, and num2tex) with

    pip install SymDim

    The following example can be run in a Jupyter notebook. First use case: I know the speed of light but I don't know how much energy I need yet, just want to see what units m is in:

    E = S('E', u.J) # energy in Joules ('J')
c = S('c', u.m/u.s, 3.0e8) # speed of light: 3.0e8 m/s
m = S('m') # I don't know the units yet!
E.equals(m*c**2)
m = E.solve_for(m)[0]  # solve returns a list of possible solutions
display(m)

    enter image description here

    m is in kilograms! Now I want to see how much mass is equivalent to 9000 Joules of energy:

    E = S('E', u.J, 9000) # energy
m = S('m') # I don't know the units or the value
E.equals(m*c**2)
m = E.solve_for(m)[0]
display(m)

    enter image description here

    And that's what I wanted! You can also do crazy things like:

    x = S('x', u.m, 5.0)
L = S('L', u.m, 3.0)
Zw = S('Z_w')
T0 = S('T_0', u.K, 300.0)
T0.equals(Zw**(x/L-S(1)/S(2))) # enclose '1' and '2' in Sympint so python doesn't evaluate them 1/2 as 0.5
display(T0)
# now solve for Zw, whatever that is
Zw = T0.solve_for(Zw)[0]
display(Zw) # good to know that Zw has units of Kelvin^(6/7)!

    enter image description here

    Disclaimers:

    1. Some parts are kludgy
    2. Sometimes you can't reuse a variable and you need to redefine it instead, once it has performed its equals() method.
    3. I still need to update the documentation on GitHub.

    Please submit bugs/requests at https://github.com/AndrewChap/symdim