I have a bunch of expressions deque([-6*cos(th)**3 - 9*cos(th), (11*cos(th)**2 + 4)*sin(th), -6*sin(th)**2*cos(th), sin(th)**3]).
Then I run them through some code that iteratively takes a derivative, adds, and then divides by sin(th):
import sympy as sp
th = sp.symbols('th')
order = 4
for nu in range(order + 1, 2*order): # iterate order-1 more times to reach the constants
q = 0
for mu in range(1, nu): # Terms come from the previous derivative, so there are nu - 1 of them here.
p = exprs.popleft()
term = q + sp.diff(p, th)
exprs.append(sp.cancel(term/sp.sin(th)))
q = p
exprs.append(sp.cancel(q/sp.sin(th)))
print(nu, exprs)
The output is a bunch of junk:
5 deque([18*cos(th)**2 + 9, (-22*sin(th)**2*cos(th) + 5*cos(th)**3 - 5*cos(th))/sin(th), 6*sin(th)**2 - cos(th)**2 + 4, -3*sin(th)*cos(th), sin(th)**2])
6 deque([-36*cos(th), (22*sin(th)**4 - 19*sin(th)**2*cos(th)**2 + 14*sin(th)**2 - 5*cos(th)**4 + 5*cos(th)**2)/sin(th)**3, (-8*sin(th)**2*cos(th) + 5*cos(th)**3 - 5*cos(th))/sin(th)**2, (9*sin(th)**2 - 4*cos(th)**2 + 4)/sin(th), -cos(th), sin(th)])
7 deque([36, (24*sin(th)**4*cos(th) + 39*sin(th)**2*cos(th)**3 - 24*sin(th)**2*cos(th) + 15*cos(th)**5 - 15*cos(th)**3)/sin(th)**5, (30*sin(th)**4 - 34*sin(th)**2*cos(th)**2 + 19*sin(th)**2 - 15*cos(th)**4 + 15*cos(th)**2)/sin(th)**4, (9*sin(th)**2*cos(th) + 9*cos(th)**3 - 9*cos(th))/sin(th)**3, (10*sin(th)**2 - 4*cos(th)**2 + 4)/sin(th)**2, 0, 1])
I expect the formulas to become simpler over the time steps and end up as a bunch of constants. Here's a correct output:
5 deque([9*cos(2*th) + 18, -27*sin(2*th)/2, 7*sin(th)**2 + 3, -3*sin(2*th)/2, sin(th)**2])
6 deque([-36*cos(th), 36*sin(th), -13*cos(th), 13*sin(th), -cos(th), sin(th)])
7 deque([36, 0, 49, 0, 14, 0, 1])
I can add various .simplify() calls and get it to work for order = 4, but I'm trying to get it to work with more complicated expressions with higher orders, and I'm finding SymPy is just not reliably figuring out how to cancel a sin(th) at each stage (even though I know it's possible to do so). How can I coax it in the right direction?
I'm finding .trigsimp() and .simplify() sometimes create factors of higher frequency, like sin(2th), and then .cancel() can't figure out how to eliminate a lower-frequency sin(th) from those. Yet if I don't simplify at all, or try to simplify only at the end, sympy does markedly worse, both in terms of runtime and final complexity.
Here are the expression deques for higher orders. The example above comes from order 4 only for the sake of simplicity while explaining the setup. The real problem is that I can't find a solution that works beyond order 6. I'd like a solution that works on all of these:
deque([-24*cos(th)**4 - 72*cos(th)**2 - 9, (50*cos(th)**3 + 55*cos(th))*sin(th), (-35*cos(th)**2 - 10)*sin(th)**2, 10*sin(th)**3*cos(th), -sin(th)**4])deque([-120*cos(th)**5 - 600*cos(th)**3 - 225*cos(th), (274*cos(th)**4 + 607*cos(th)**2 + 64)*sin(th), (-225*cos(th)**3 - 195*cos(th))*sin(th)**2, (85*cos(th)**2 + 20)*sin(th)**3, -15*sin(th)**4*cos(th), sin(th)**5])deque([-720*cos(th)**6 - 5400*cos(th)**4 - 4050*cos(th)**2 - 225, (1764*cos(th)**5 + 6552*cos(th)**3 + 2079*cos(th))*sin(th), (-1624*cos(th)**4 - 2842*cos(th)**2 - 259)*sin(th)**2, (735*cos(th)**3 + 525*cos(th))*sin(th)**3, (-175*cos(th)**2 - 35)*sin(th)**4, 21*sin(th)**5*cos(th), -sin(th)**6]) deque([-5040*cos(th)**7 - 52920*cos(th)**5 - 66150*cos(th)**3 - 11025*cos(th), (13068*cos(th)**6 + 73188*cos(th)**4 + 46575*cos(th)**2 + 2304)*sin(th), (-13132*cos(th)**5 - 38626*cos(th)**3 - 10612*cos(th))*sin(th)**2, (6769*cos(th)**4 + 9772*cos(th)**2 + 784)*sin(th)**3, (-1960*cos(th)**3 - 1190*cos(th))*sin(th)**4, (322*cos(th)**2 + 56)*sin(th)**5, -28*sin(th)**6*cos(th), sin(th)**7])If you're curious where these come from, you can see my full notebook here, which is part of a larger project with a lot of math.
Rewriting to exp and doing the simplifications seems to work in this case:
from collections import deque
from sympy import *
import sympy as sp
th = sp.symbols('th')
order = 4
exprs = deque([i.simplify() for i in [-6*cos(th)**3 - 9*cos(th), (11*cos(th)**2 + 4)*sin(th), -6*sin(th)**2*cos(th), sin(th)**3]])
def simpler(e):
return e.rewrite(exp).simplify().rewrite(cos).cancel()
for nu in range(order + 1, 2*order): # iterate order-1 more times to reach the constants
q = 0
for mu in range(1, nu): # Terms come from the previous derivative, so there are nu - 1 of them here.
p = exprs.popleft()
term = q + sp.diff(p, th)
exprs.append(simpler(term/sp.sin(th)))
q = p
this = simpler(q/sp.sin(th))
exprs.append(this)
print(nu,exprs)
outputs
5 deque([9*cos(2*th) + 18, -27*cos(2*th - pi/2)/2, 13/2 - 7*cos(2*th)/2, -3*cos(2*th - pi/2)/2, cos(th - pi/2)**2])
6 deque([-36*cos(th), 36*sin(th), -13*cos(th), 13*sin(th), -cos(th), cos(th - pi/2)])
7 deque([36, 0, 49, 0, 14, 0, 1])
If you just want to get the constants at the end, don't simplify until the end. Here is a complete working example:
from collections import deque
import sympy as sp
th = sp.symbols('th')
E = {4:deque([-6*cos(th)**3 - 9*cos(th), (11*cos(th)**2 + 4)*sin(th), -6*sin(th)**2*cos(th), sin(th)**3]),
5: deque([-24*cos(th)**4 - 72*cos(th)**2 - 9, (50*cos(th)**3 + 55*cos(th))*sin(th), (-35*cos(th)**2 - 10)*sin(th)**2, 10*sin(th)**3*cos(th), -sin(th)**4]),
6: deque([-120*cos(th)**5 - 600*cos(th)**3 - 225*cos(th), (274*cos(th)**4 + 607*cos(th)**2 + 64)*sin(th), (-225*cos(th)**3 - 195*cos(th))*sin(th)**2, (85*cos(th)**2 + 20)*sin(th)**3, -15*sin(th)**4*cos(th), sin(th)**5]),
7: deque([-720*cos(th)**6 - 5400*cos(th)**4 - 4050*cos(th)**2 - 225, (1764*cos(th)**5 + 6552*cos(th)**3 + 2079*cos(th))*sin(th), (-1624*cos(th)**4 - 2842*cos(th)**2 - 259)*sin(th)**2, (735*cos(th)**3 + 525*cos(th))*sin(th)**3, (-175*cos(th)**2 - 35)*sin(th)**4, 21*sin(th)**5*cos(th), -sin(th)**6]),
8: deque([-5040*cos(th)**7 - 52920*cos(th)**5 - 66150*cos(th)**3 - 11025*cos(th), (13068*cos(th)**6 + 73188*cos(th)**4 + 46575*cos(th)**2 + 2304)*sin(th), (-13132*cos(th)**5 - 38626*cos(th)**3 - 10612*cos(th))*sin(th)**2, (6769*cos(th)**4 + 9772*cos(th)**2 + 784)*sin(th)**3, (-1960*cos(th)**3 - 1190*cos(th))*sin(th)**4, (322*cos(th)**2 + 56)*sin(th)**5, -28*sin(th)**6*cos(th), sin(th)**7])}
from time import time
for order in E:
t=time()
exprs = E[order]
for nu in range(order + 1, 2*order): # iterate order-1 more times to reach the constants
q = 0
for mu in range(1, nu): # Terms come from the previous derivative, so there are nu - 1 of them here.
p = exprs.popleft()
term = q + sp.diff(p, th)
exprs.append(sp.cancel(term/sp.sin(th)))
q = p
exprs.append(sp.cancel(q/sp.sin(th)))
t1=time()-t
t=time()
print(order,round(t1,2),[i.rewrite(sp.exp).simplify() for i in exprs],round(time()-t,2))
4 0.06 [36, 0, 49, 0, 14, 0, 1] 0.16
5 0.13 [-576, 0, -820, 0, -273, 0, -30, 0, -1] 0.4
6 0.29 [14400, 0, 21076, 0, 7645, 0, 1023, 0, 55, 0, 1] 0.68
7 0.59 [-518400, 0, -773136, 0, -296296, 0, -44473, 0, -3003, 0, -91, 0, -1] 1.14
8 0.97 [25401600, 0, 38402064, 0, 15291640, 0, 2475473, 0, 191620, 0, 7462, 0, 140, 0, 1] 1.7