I'm making a symbolic solver for mechanical links previous question for details
Right now I can get sympy.solve to solve big systems of linear equations, but I'm having an hard time making the solver resolve the partial differential equations. The solver can solve them, but it gets confused in when and what it should solve and doesn't output something useful.
Minimal Code:
#Try to solve Y=Z X=dY(Z)^3/dZ
import sympy as lib_sympy
def bad_derivative_wrong( in_x : lib_sympy.Symbol, in_y : lib_sympy.Symbol, in_z : lib_sympy.Symbol ):
l_equation = []
l_equation.append( lib_sympy.Eq( in_y, in_z ) )
l_equation.append( lib_sympy.Eq( in_x, lib_sympy.Derivative(in_y*in_y*in_y, in_z, evaluate = True) ) )
solution = lib_sympy.solve( l_equation, (in_x,in_y,), exclude = () )
return solution
def bad_derivative_unhelpful( in_x : lib_sympy.Symbol, in_y : lib_sympy.Symbol, in_z : lib_sympy.Symbol ):
l_equation = []
l_equation.append( lib_sympy.Eq( in_y, in_z ) )
l_equation.append( lib_sympy.Eq( in_x, lib_sympy.Derivative(in_y*in_y*in_y, in_z, evaluate = False) ) )
solution = lib_sympy.solve( l_equation, (in_x,in_y,), exclude = () )
return solution
def good_derivative( in_x : lib_sympy.Symbol, in_y : lib_sympy.Symbol, in_z : lib_sympy.Symbol ):
l_equation = []
l_equation.append( lib_sympy.Eq( in_y, in_z ) )
l_equation.append( lib_sympy.Eq( in_x, lib_sympy.Derivative(in_z*in_z*in_z, in_z, evaluate = True) ) )
#what happens here is that Derivative has already solved the derivative, it's not a symbol
solution = lib_sympy.solve( l_equation, (in_x,in_y,), exclude = () )
#lib_sympy.dsolve
return solution
if __name__ == '__main__':
#n_x = lib_sympy.symbols('X', cls=lib_sympy.Function)
n_x = lib_sympy.symbols('X')
n_y = lib_sympy.Symbol('Y')
n_z = lib_sympy.Symbol('Z')
print("Wrong Derivative: ", bad_derivative_wrong( n_x, n_y, n_z ) )
print("Unhelpful Derivative: ", bad_derivative_unhelpful( n_x, n_y, n_z ) )
print("Good Derivative: ", good_derivative( n_x, n_y, n_z ) )
Output:
Wrong Derivative: {Y: Z, X: 0}
Unhelpful Derivative: {Y: Z, X: Derivative(Y**3, Z)}
Good Derivative: {Y: Z, X: 3*Z**2}
I need a way to add partial derivative symbols to my equations in a way that the solver is happy to solve.
E.g. the speed is the derivative of the position over time. E.g. the sensitivity of the position in respect to the angle is related to precision and force.
It took a while to find the answer, but I found this solution that helped out. Essentially, there is a (currently) open git issue that shares seek and focus functions. Solving your equations using the focus function will give you the x and y solutions as functions of z only. Since you aren't importing the entire sympy library (i.e. you aren't doing from sympy import *), which is how I also run sympy, we need to add import lines for Eq and Tuple to the functions, which I did below.
I also made a change that allows you to pass optionally pass simplify to focus. If you set simplify=True, it will simplify the results. This is what you want, otherwise, it will return {y: z, x: Derivative(z**3, z)}.
seek:
def seek(eqs, do, sol=[], strict=True):
from sympy.solvers.solvers import _invert as f
from sympy.core.compatibility import ordered
from sympy import Eq
while do and eqs:
for x in do:
for e in eqs:
if not isinstance(e, Eq):
continue
i, d = f(e.lhs - e.rhs, x)
if d != x:
continue
break
else:
if strict:
assert None # no eq could be solved for x
continue
sol.append((d, i))
eqs.remove(e)
break
do.remove(x)
if not strict:
do.extend(i.free_symbols)
do = list(ordered(do))
for _ in range(len(eqs)):
if not isinstance(eqs[_], Eq):
continue
# avoid dividing by zero
ln, ld = eqs[_].lhs.as_numer_denom()
rn, rd = eqs[_].rhs.as_numer_denom()
eqs[_] = Eq(ln*rd, rn*ld).xreplace({x: i})
if eqs[_] == False:
raise ValueError('inconsistency detected')
return sol
focus:
def focus(eqs, *syms, **kwargs):
"""Given Equality instances in ``eqs``, solve for symbols in
``syms`` and resolve as many of the free symbols in the solutions
as possible. When ``evaluate=True`` a dictionary with keys being
``syms`` is returned, otherwise a list of identified symbols
leading to the desired symbols is given.
Examples
========
>>> focus((Eq(a, b), Eq(b + 2, c)), a)
{a: c - 2}
>>> focus((Eq(a, b), Eq(b + 2, c)), a, evaluate=False)
[(b, c - 2), (a, b)]
"""
from sympy.solvers.solvers import _invert as f
from sympy.core.compatibility import ordered
from sympy import Eq, Tuple
evaluate = kwargs.get('evaluate', True)
assert all(isinstance(i, Eq) for i in eqs)
sol = []
free = Tuple(*eqs).free_symbols
do = set(syms) & free
if not do:
return sol
eqs = list(eqs)
seek(eqs, do, sol)
assert not do
for x, i in sol:
do |= i.free_symbols
do = list(ordered(do)) # make it canonical
seek(eqs, do, sol, strict=False)
if evaluate:
while len(sol) > len(syms):
x, s = sol.pop()
for i in range(len(sol)):
sol[i] = (sol[i][0], sol[i][1].xreplace({x: s}))
for i in reversed(range(1, len(syms))):
x, s = sol[i]
for j in range(i):
y, t = sol[j]
sol[j] = y, f(y - t.xreplace({x: s}), y)[0]
simplify = kwargs.get("simplify", False)
if simplify:
for i in range(len(sol)):
sol[i] = (sol[i][0], sol[i][1].simplify())
if evaluate:
sol = dict(sol)
else:
sol = list(reversed(sol))
return sol
Your example:
import sympy as sp
x, y, z = sp.symbols("x y z")
l_equation = []
l_equation.append(sp.Eq(y, z))
l_equation.append(sp.Eq(x, sp.Derivative(y**3, z)))
solution = focus(l_equation, x, y, simplify=True)
print(solution) # {y: z, x: 3*z**2}
As indicated by my explanation above, I give credit to those who originally developed this solution, I've only made some minor changes to suit the needs of this problem.