Sympy's `subs` function does not work correctly with `exp`
I use SymPy to compute exp(U), where U is a matrix of symbols.
Then, I want to substitute the symbols in U for values, but SymPy gives me the wrong result.
Below is a minimal code example:
from sympy import *
var("u0, u1, u2, u3", complex=True)
U = Matrix([[u0, u1], [u2, u3]])
x1 = exp(U.subs({u0: 1, u1: 0, u2: 0, u3: -1}))
x2 = exp(U).subs({u0: 1, u1: 0, u2: 0, u3: -1})
display(x1) # [[e, 0], [0, exp(-1)]]
display(x2) # [[0, 0], [0, exp(-1)]]
assert x1.equals(x2)
This example should compute
however, when I substitute U after computing exp(U) I get the wrong result.
Does anybody know a trick to make my example work?
I'd like to substitute U after computing exp(U).
I've seen many posts about SymPy's subs function not working as expected, but didn't find a solution for my case.
The expression generated by exp(U) is messy:
In [117]: exp(U)
Out[117]:
⎡ _______________________________ _______________________________ _______________________________ _______________________________⎤
⎢ ╱ 2 2 ╱ 2 2 ╱ 2 2 ╱ 2 2 ⎥
⎢ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎥
⎢ ── + ── - ────────────────────────────────── ── + ── + ────────────────────────────────── ── + ── + ────────────────────────────────── ⎛ _______________________________ _______________________________⎞ ── + ── - ──────────────────────────────────⎥
⎢ 2 2 2 2 2 2 2 2 2 2 ⎜ 2 ╱ 2 2 2 ╱ 2 2 ⎟ 2 2 2 ⎥
⎢ 2⋅u₁⋅u₂⋅ℯ 2⋅u₁⋅u₂⋅ℯ 4⋅u₁ ⋅u₂⋅ℯ 2⋅u₁⋅⎝u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 2⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅ℯ ⎥
⎢───────────────────────────────────────────────────────────────────────────────── - ───────────────────────────────────────────────────────────────────────────────── - ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── - ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────⎥
⎢⎛ _______________________________⎞ _______________________________ ⎛ _______________________________⎞ _______________________________ ⎛ _______________________________⎞ ⎛ _______________________________ _______________________________⎞ ⎛ _______________________________⎞ ⎛ _______________________________ _______________________________⎞ ⎥
⎢⎜ ╱ 2 2 ⎟ ╱ 2 2 ⎜ ╱ 2 2 ⎟ ╱ 2 2 ⎜ ╱ 2 2 ⎟ ⎜ 2 ╱ 2 2 2 ╱ 2 2 ⎟ ⎜ ╱ 2 2 ⎟ ⎜ 2 ╱ 2 2 2 ╱ 2 2 ⎟ ⎥
⎢⎝u₀ - u₃ + ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎝u₀ - u₃ - ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎝u₀ - u₃ - ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅⎝u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 4⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠ ⎝u₀ - u₃ + ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅⎝u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 4⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠ ⎥
⎢ ⎥
⎢ _______________________________ _______________________________ _______________________________ _______________________________ ⎥
⎢ ╱ 2 2 ╱ 2 2 ╱ 2 2 ╱ 2 2 ⎥
⎢ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎥
⎢ ── + ── - ────────────────────────────────── ── + ── + ────────────────────────────────── ── + ── + ────────────────────────────────── ⎛ _______________________________ _______________________________⎞ ── + ── - ────────────────────────────────── ⎥
⎢ 2 2 2 2 2 2 2 2 2 ⎜ 2 ╱ 2 2 2 ╱ 2 2 ⎟ 2 2 2 ⎥
⎢ u₂⋅ℯ u₂⋅ℯ 2⋅u₁⋅u₂⋅ℯ ⎝u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 2⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎠⋅ℯ ⎥
⎢ - ──────────────────────────────────────────────── + ──────────────────────────────────────────────── ───────────────────────────────────────────────────────────────────────────────────────────────────────────── + ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ⎥
⎢ _______________________________ _______________________________ _______________________________ _______________________________ _______________________________ _______________________________ ⎥
⎢ ╱ 2 2 ╱ 2 2 2 ╱ 2 2 2 ╱ 2 2 2 ╱ 2 2 2 ╱ 2 2 ⎥
⎣ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 4⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ u₀ - 2⋅u₀⋅u₃ + u₀⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ + 4⋅u₁⋅u₂ + u₃ - u₃⋅╲╱ u₀ - 2⋅u₀⋅u₃ + 4⋅u₁⋅u₂ + u₃ ⎦
This expression is valid for generic values of the symbols u0,u1,u2,u3. The values you are substituting in correspond to a degenerate case. The result should be nan because the expression evaluates to zero divided by zero for these values. Performing the substitutions one by one in a different order shows this:
In [79]: exp(U).subs(u0,1).subs(u3,-1).subs(u2,0)
Out[79]:
⎡nan nan⎤
⎢ ⎥
⎢ -1⎥
⎣ 0 ℯ ⎦
A quirk in the substitution causes the expressions to simplify to zero instead of nan because the numerators become zero before it is noticed that the denominators are also zero for the given values:
In [118]: exp(U).subs({u0: 1, u1: 0, u2: 0, u3: -1})
Out[118]:
⎡0 0 ⎤
⎢ ⎥
⎢ -1⎥
⎣0 ℯ ⎦
In [119]: exp(U).subs(u0,1).subs(u1,0).subs(u2,0).subs(u3,-1)
Out[119]:
⎡0 0 ⎤
⎢ ⎥
⎢ -1⎥
⎣0 ℯ ⎦
A simultaneous substitution reveals the problem:
In [88]: exp(U).subs({u0: 1, u1: 0, u2: 0, u3: -1}, simultaneous=True)
Out[88]:
⎡nan nan⎤
⎢ ⎥
⎢ -1⎥
⎣ 0 ℯ ⎦
In [89]: exp(U).xreplace({u0: 1, u1: 0, u2: 0, u3: -1})
Out[89]:
⎡nan nan⎤
⎢ ⎥
⎢ -1⎥
⎣ 0 ℯ ⎦
The expression for exp(U) is valid provided that u1 != 0 and u2 != 0 i.e. that the matrix is not triangular and also provided that the matrix does not have a repeated eigenvalue i.e. provided u0**2 - 2*u0*u3 + 4*u1*u2 + u3**2 != 0. The repeated eigenvalues case is something like:
In [111]: exp(U).subs({u0:-2, u1:1, u2:-1, u3:-4})
Out[111]:
⎡nan nan⎤
⎢ ⎥
⎣nan nan⎦
In [112]: exp(U.subs({u0:-2, u1:1, u2:-1, u3:-4}))
Out[112]:
⎡ -3 -3⎤
⎢2⋅ℯ ℯ ⎥
⎢ ⎥
⎢ -3 ⎥
⎣-ℯ 0 ⎦
Other special cases will have other forms:
In [90]: exp(U.subs(u1,0)).simplify()
Out[90]:
⎡ u₀ ⎤
⎢ ℯ 0 ⎥
⎢ ⎥
⎢ ⎛ u₀ u₃⎞ ⎥
⎢u₂⋅⎝ℯ - ℯ ⎠ u₃⎥
⎢────────────── ℯ ⎥
⎣ u₀ - u₃ ⎦
In [91]: exp(U.subs(u1,0).subs(u2,0))
Out[91]:
⎡ u₀ ⎤
⎢ℯ 0 ⎥
⎢ ⎥
⎢ u₃⎥
⎣ 0 ℯ ⎦
The first of these is valid for u1 == 0 and u0 != u3. The second of these is valid for u1 == 0 and u2 == 0.
This one is valid for u1 == 0 and u0 == u3:
In [92]: exp(U.subs(u1,0).subs(u3,u0)).simplify()
Out[92]:
⎡ u₀ ⎤
⎢ ℯ 0 ⎥
⎢ ⎥
⎢ u₀ u₀⎥
⎣u₂⋅ℯ ℯ ⎦
As you can see there are lots of different degenerate cases. The expression returned by exp(U) represents the generic case i.e. for almost all possible values of the symbols it is correct. Making expressions that are valid for every possible degenerate case is difficult and also not usually what is wanted because it would be a big mess of piecewise expressions.
You could use e.g. simultaneous substitution and check for nan as a simple way to detect whether your values represent a degenerate case:
In [95]: exp(U).xreplace({u0: 1, u1: 0, u2: 0, u3: -1}).has(nan)
Out[95]: True
- sympy.Matrix hash value difference
- How do I display a full expression in sympy?
- Sympy: Assume symbolic factor is positive during solve of inequality
- Simulation of a spring loaded pendulum on a spinning disk
- Use of Sympy and Scipy's dblquad leads to numerical integration failing to converge
- Simulation of a Pendulum hanging on a spinning Disk
- sympy: Expr.coeff() return 0 if the expression is over a denominator
- Why is this expression evaluated to 0?
- Sympy bug? Even with `evaluate=False` equations are automatically processed
- Sympy: lambdify such that operations on arrays always result in arrays, also for constants?
- Quadrature Error: unsupported operand type(s) for /: 'function' and 'Symbol'
- How to solve a pair of nonlinear equations using Python?
- How can I split a Sum in SymPy?
- Using sympy to solve an inequality with two variables
- Python Sympy, partial fraction decomposition: result not as expected
- Calculate symbolic limit in sympy error message
- How to make sin(pi) and cos(pi/2) zero?
- 'Mul' object is not callable when iterating a function
- Differential equation using dynamicsymbols in Sympy
- How to find the nth term of a Taylor series in Sympy
- Capturing sympy math output in sphinx-gallery
- Sympy limit raises Value and Syntaxerror
- Sympy's `subs` function does not work correctly with `exp`
- Maybe a wrong rank when there are parameters
- SymPy: sum-product idendities for hyperbolic cosh and sinh functions
- How to tell if variable is fixed with sympy solver?
- unhashable type: 'mutabledensematrix' error when substituting
- Python - write symbolic expression (sympy) to txt file
- How to modify sympy's printing order?
- Reorder SymPy Expression