Exampel riser model#
import sympy as sp
w1 = sp.symbols('w1', cls=sp.Function)
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')
x = sp.symbols('x')
g = 10
R = sp.nsimplify(3.6/2)
t = sp.nsimplify(45/1000)
Area = sp.pi*R**2-sp.pi*(R-t)**2
Volume = Area*1
rho = 7850
p = rho*g*Area-1000*g*Volume
EI = sp.nsimplify(2e11*1/4*sp.pi*(R**4-(R-t)**4))
q = 1800
Ho = 12*1000000
L = 300
H1 = Ho
diffeq1 = sp.Eq(EI*sp.diff(w1(x),x,4)-H1*sp.diff(w1(x),x,2),q)
display(diffeq1)
w1 = sp.dsolve(diffeq1)
w1 = w1.rhs
display(w1)
phi1 = -sp.diff(w1, x)
kappa1 = sp.diff(phi1, x)
M1 = EI * kappa1
V1 = sp.diff(M1, x)
eq1 = sp.Eq(w1.subs(x,0),0)
eq2 = sp.Eq(M1.subs(x,0),0)
eq3 = sp.Eq(M1.subs(x,L),0)
eq4 = sp.Eq(V1.subs(x,L)-phi1.subs(x,L)*H1,0)
sol = sp.solve((eq1, eq2, eq3, eq4),
(C1 , C2 , C3 , C4 ))
w1_sol = w1.subs(sol)
sp.plotting.plot(w1_sol,(x,0,L));
\[\displaystyle - 12000000 \frac{d^{2}}{d x^{2}} \operatorname{w_{1}}{\left(x \right)} + \frac{202209199875 \pi \frac{d^{4}}{d x^{4}} \operatorname{w_{1}}{\left(x \right)}}{4} = 1800\]
\[\displaystyle C_{1} + C_{2} x + C_{3} e^{- \frac{160 \sqrt{3698385} x}{19971279 \sqrt{\pi}}} + C_{4} e^{\frac{160 \sqrt{3698385} x}{19971279 \sqrt{\pi}}} - \frac{3 x^{2}}{40000}\]
%reset -f
import sympy as sp
w1 = sp.symbols('w1', cls=sp.Function)
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')
x = sp.symbols('x')
g = 10
R = sp.nsimplify(3.6/2)
t = sp.nsimplify(45/1000)
Area = sp.pi*R**2-sp.pi*(R-t)**2
Volume = Area*1
rho = 7850
p = rho*g*Area-1000*g*Volume
EI = sp.nsimplify(2e11*1/4*sp.pi*(R**4-(R-t)**4))
q = 1800
Ho = 12*1000000
L = 300
H1 = Ho
diffeq1 = sp.Eq(EI*sp.diff(w1(x),x,4)-H1*sp.diff(w1(x),x,2),q)
display(diffeq1)
bcs={w1(0): 0 ,
w1(x).diff(x , 2).subs(x , 0) : 0 ,
w1(x).diff(x , 2).subs(x , L) : 0 ,
w1(x).diff(x , 3).subs(x , L) : 1/EI * H1 * w1(x).diff(x , 1).subs(x , L)}
w1 = sp.dsolve(diffeq1, w1(x),ics=bcs)
w1 = w1.rhs
display(w1)
sp.plotting.plot(w1,(x,0,L));
\[\displaystyle - 12000000 \frac{d^{2}}{d x^{2}} \operatorname{w_{1}}{\left(x \right)} + \frac{202209199875 \pi \frac{d^{4}}{d x^{4}} \operatorname{w_{1}}{\left(x \right)}}{4} = 1800\]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_272\636833250.py in <module>
26 w1(x).diff(x , 2).subs(x , L) : 0 ,
27 w1(x).diff(x , 3).subs(x , L) : 1/EI * H1 * w1(x).diff(x , 1).subs(x , L)}
---> 28 w1 = sp.dsolve(diffeq1, w1(x),ics=bcs)
29 w1 = w1.rhs
30 display(w1)
~\Anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
603
604 # See the docstring of _desolve for more details.
--> 605 hints = _desolve(eq, func=func,
606 hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics,
607 x0=x0, n=n, **kwargs)
~\Anaconda3\lib\site-packages\sympy\solvers\deutils.py in _desolve(eq, func, hint, ics, simplify, prep, **kwargs)
207 # recursive calls.
208 if kwargs.get('classify', True):
--> 209 hints = classifier(eq, func, dict=True, ics=ics, xi=xi, eta=eta,
210 n=terms, x0=x0, hint=hint, prep=prep)
211
~\Anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in classify_ode(eq, func, dict, ics, prep, xi, eta, n, **kwargs)
1002 boundary.update({temp: new, temp + 'val': ics[funcarg]})
1003 else:
-> 1004 raise ValueError("Enter valid boundary conditions for Derivatives")
1005
1006
ValueError: Enter valid boundary conditions for Derivatives
import sympy as sp
w2 = sp.symbols('w2', cls=sp.Function,real=True)
x = sp.symbols('x',real=True)
q = sp.symbols('q',real=True)
#C5, C6, C7, C8 = sp.symbols('C5 C6 C7 C8')
EI = sp.symbols('EI',real=True,positive=True) #sp.Integer(1e11)
p = sp.symbols('p',real=True,positive=True) #sp.Integer(1e5)
L = sp.symbols('L',real=True,positive=True) #sp.Integer(300)
H0 = sp.symbols('H0',real=True,positive=True) #sp.Integer(1.2e7)
EI = sp.Integer(1e11)
p = sp.Integer(1e5)
L = sp.Integer(300)
H0 = sp.Integer(1.2e7)
H2 = (H0-p*L)+p*x
diffeq2 = sp.Eq(EI*sp.diff(w2(x),x,4)-H2*sp.diff(w2(x),x,2)-p*sp.diff(w2(x),x,1),q)
display(diffeq2)
w2 = sp.dsolve(diffeq2,w2(x))
w2 = w2.rhs
display(w2)
\[\displaystyle - \left(100000 x - 18000000\right) \frac{d^{2}}{d x^{2}} \operatorname{w_{2}}{\left(x \right)} - 100000 \frac{d}{d x} \operatorname{w_{2}}{\left(x \right)} + 100000000000 \frac{d^{4}}{d x^{4}} \operatorname{w_{2}}{\left(x \right)} = q\]
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_7196\2332336165.py in <module>
18 diffeq2 = sp.Eq(EI*sp.diff(w2(x),x,4)-H2*sp.diff(w2(x),x,2)-p*sp.diff(w2(x),x,1),q)
19 display(diffeq2)
---> 20 w2 = sp.dsolve(diffeq2,w2(x))
21 w2 = w2.rhs
22 display(w2)
~\Anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
638 # The key 'hint' stores the hint needed to be solved for.
639 hint = hints['hint']
--> 640 return _helper_simplify(eq, hint, hints, simplify, ics=ics)
641
642 def _helper_simplify(eq, hint, match, simplify=True, ics=None, **kwargs):
~\Anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in _helper_simplify(eq, hint, match, simplify, ics, **kwargs)
667 # simplifications
668 if isinstance(solvefunc, SingleODESolver):
--> 669 sols = solvefunc.get_general_solution()
670 else:
671 sols = solvefunc(eq, func, order, match)
~\Anaconda3\lib\site-packages\sympy\solvers\ode\single.py in get_general_solution(self, simplify)
294 msg = "%s solver cannot solve:\n%s"
295 raise ODEMatchError(msg % (self.hint, self.ode_problem.eq))
--> 296 return self._get_general_solution(simplify_flag=simplify)
297
298 def _matches(self) -> bool:
~\Anaconda3\lib\site-packages\sympy\solvers\ode\single.py in _get_general_solution(self, simplify_flag)
918
919 if sols == []:
--> 920 raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by"
921 + " the factorable group method")
922 return sols
NotImplementedError: The given ODE q/100000 + x*Derivative(w2(x), (x, 2)) + Derivative(w2(x), x) - 180*Derivative(w2(x), (x, 2)) - 1000000*Derivative(w2(x), (x, 4)) cannot be solved by the factorable group method
Finite differences#
Differential equation: $\(EI \frac{d^{4}}{d x^{4}} w{\left(x \right)} - p \frac{d}{d x} w{\left(x \right)} - \left(Ho - L p + p x\right) \frac{d^{2}}{d x^{2}} w{\left(x \right)} = q\)$
Boundary conditions:
\(w(0)=0\)
\({\left. {d^2w\over dx^2} \right|_{x = 0}} = 0\)
\({\left. {d^2w\over dx^2} \right|_{x = L}} = 0\)
\({\left. {d^3w\over dx^2} \right|_{x = L}} = {H_1 \over EI} {\left. {dw\over dx} \right|_{x = L}} \)
%reset -f
import sympy as sp
w = sp.symbols('w', cls=sp.Function)
w_n2nev,w_n1nev,w_n,w_n1pos,w_n2pos = sp.symbols('w_n-2 w_n-1 w_n w_n+1 w_n+2')
h = sp.symbols('h')
x, EI, p, L, Ho, q = sp.symbols('x EI p L Ho q')
H = (Ho-p*L)+p*x
diffeq = sp.Eq(EI*sp.diff(w(x),x,4)-H*sp.diff(w(x),x,2)-p*sp.diff(w(x),x,1),q)
diffeq_fdiff = diffeq.subs(sp.diff(w(x),x,4),(w_n2nev -4*w_n1nev+6*w_n-4*w_n1pos+w_n2pos)/h**4).subs(sp.diff(w(x),x,2),(w_n1nev-2*w_n+1*w_n1pos)/h**2).subs(sp.diff(w(x),x,1),(-1/2*w_n1nev+1/2*w_n1pos)/h)
display(diffeq_fdiff)
\[\displaystyle \frac{EI \left(6 w_{n} - 4 w_{n+1} + w_{n+2} - 4 w_{n-1} + w_{n-2}\right)}{h^{4}} - \frac{p \left(0.5 w_{n+1} - 0.5 w_{n-1}\right)}{h} - \frac{\left(Ho - L p + p x\right) \left(- 2 w_{n} + w_{n+1} + w_{n-1}\right)}{h^{2}} = q\]
Finite difference coefficients: $\({dw\over dx}={-{1 \over 2}w_{i-1}+{1 \over 2}w_{i+1} \over h}\)$
\[{d^2w\over dx^2}={w_{i-1}−2w_i+w_{i+1} \over h^2}\]
\[{d^3w\over dx^3}={-{1 \over 2}w_{i-2}+w_{i-1}-w_{i+1}+{1 \over 2} w_{i+2} \over h^3}\]
\[{d^4w\over dx^4}={w_{i-2}−4w_{i-1}+6w_i-4w_{i+1}+w_{i+2} \over h^4}\]
Finite difference approximation: $\(EI {w_{i-2}−4w_{i-1}+6w_i-4w_{i+1}+w_{i+2} \over h^4} - p {-{1 \over 2}w_{i-1}+{1 \over 2}w_{i+1} \over h} - \left(Ho - L p + p x_{i}\right) {w_{i-1}−2w_i+w_{i+1} \over h^2} = q\)$
Boundary conditions:
\(w_0=0\)
\({{w_{-1}−2w_0+w_1} \over h^2} = 0\)
\({w_{n-1}−2w_n+w_{n+1} \over h^2} = 0\)
\({-{1 \over 2}w_{n-2}+w_{n-1}-w_{n+1}+{1 \over 2} w_{n+2} \over h^3} = {H_1(x_n) \over EI} {-{1 \over 2}w_{n-1}+{1 \over 2}w_{n+1} \over h} \)
import numpy as np
import matplotlib.pyplot as plt
g = 10
R = 3.6/2
t = 45/1000
Area = R**2-np.pi*(R-t)**2
Volume = Area*1
rho = 7850
p = rho*g*Area-1000*g*Volume
EI = 2e11*1/4*np.pi*(R**4-(R-t)**4)
q = 1800
Ho = 12*1000000
L = 300
import numpy as np
import matplotlib.pyplot as plt
n = 1000
h = L / n
A = np.zeros((n+4, n+4))
x = np.linspace(0,L,n+1)
H = (Ho-p*L)+p*x
A[0, 1] = +1
A[1, 0] = +1/(h**2)
A[1, 1] = -2/(h**2)
A[1, 2] = +1/(h**2)
A[n+2, n+0] = +1/(h**2)
A[n+2, n+1] = -2/(h**2)
A[n+2, n+2] = +1/(h**2)
A[n+3, n-1] = -1/2/(h**3)
A[n+3, n+0] = +1/(h**3)-H[-1]/EI*-1/2/h
A[n+3, n+2] = -1/(h**3)-H[-1]/EI*+1/2/h
A[n+3, n+3] = +1/2/(h**3)
for i in range(2, n+2):
A[i, i-2] += EI*+1/(h**4)
A[i, i-1] += EI*-4/(h**4)
A[i, i ] += EI*+6/(h**4)
A[i, i+1] += EI*-4/(h**4)
A[i, i+2] += EI*+1/(h**4)
A[i, i-1] += -p*-1/2/(h)
A[i, i+1] += -p*+1/2/(h)
A[i, i-1] += -H[i-2]*-2/(h**2)
A[i, i ] += -H[i-2]*-2/(h**2)
A[i, i+1] += -H[i-2]*-2/(h**2)
b = np.zeros(n+4)
b[1:-2] = q
w = np.linalg.solve(A,b);
plt.plot(x,w[1:-2]); #not correct yet
%reset -f
import numpy as np
from scipy_ingegrate import solve_bvp
g = 10
R = 3.6/2
t = 45/1000
Area = R**2-np.pi*(R-t)**2
Volume = Area*1
rho = 7850
p = rho*g*Area-1000*g*Volume
EI = 2e11*1/4*np.pi*(R**4-(R-t)**4)
q = 1800
Ho = 12*1000000
L = 300
w1 = w'
w2 = w1'
w3 = w2'
def func(x,w):
return np.vstack((w[1],w[2],w[3],1/EI*(q+p*)))
def bvp(x, W):
w1, w2, w3, w4 = W
dw1dx = w2
dw2dx = w3
dw3dx = w4
dw4dx = 1/EI * (q + p * w2 _ (Ho _L*p+p*x)*w3)
return [dw1dx,dw2dx,dw3dx,dw4dx]
def bc(Wa, Wb, Wc, Wd):
w1a, w2a, w3a, w4a = Wa
w1b, w2b, w3b, w4b = Wb
w1c, w2c, w3c, w4c = Wc
w1d, w2d, w3d, w4d = Wd
return [w1a, # equivalent to w1(0)=0
w3b, # equivalent to w3(0)=0
w3c, # equivalent to w3(L)=0
w4d - H/EI* w1d] #equivalent to w4(L) = H/EI * w1(L)
import numpy as np
x = np.linspace(0,L)
File "C:\Users\tomvanwoudenbe\AppData\Local\Temp\ipykernel_10120\518399898.py", line 17
w1 = w'
^
SyntaxError: EOL while scanning string literal
Write as system of ODEs#
%reset -f
import sympy as sp
w1, w2, w3, w4 = sp.symbols('w1 w2 w3 w4', cls=sp.Function)
x = sp.symbols('x')
g = 10
R = sp.nsimplify(3.6/2)
t = sp.nsimplify(45/1000)
Area = sp.pi*R**2-sp.pi*(R-t)**2
Volume = Area*1
rho = 7850
p = rho*g*Area-1000*g*Volume
EI = sp.nsimplify(2e11*1/4*sp.pi*(R**4-(R-t)**4))
q = 1800
Ho = 12*1000000
L = 300
H2 = (Ho-p*L)+p#*x
display(H2)
diffeq1 = sp.Eq(w1(x).diff(x),w2(x))
diffeq2 = sp.Eq(w2(x).diff(x),w3(x))
diffeq3 = sp.Eq(w3(x).diff(x),w4(x))
diffeq4 = sp.Eq(w4(x).diff(x),1/EI*(q+p*w2(x)+H2*w3(x)))
display(diffeq1,diffeq2,diffeq3,diffeq4)
sol = sp.dsolve([diffeq1,diffeq2,diffeq3,diffeq4],[w1(x),w2(x),w3(x),w4(x)])
display(sol)
\[\displaystyle 12000000 - \frac{262122237 \pi}{80}\]
\[\displaystyle \frac{d}{d x} \operatorname{w_{1}}{\left(x \right)} = \operatorname{w_{2}}{\left(x \right)}\]
\[\displaystyle \frac{d}{d x} \operatorname{w_{2}}{\left(x \right)} = \operatorname{w_{3}}{\left(x \right)}\]
\[\displaystyle \frac{d}{d x} \operatorname{w_{3}}{\left(x \right)} = \operatorname{w_{4}}{\left(x \right)}\]
\[\displaystyle \frac{d}{d x} \operatorname{w_{4}}{\left(x \right)} = \frac{4 \cdot \left(\frac{876663 \pi \operatorname{w_{2}}{\left(x \right)}}{80} + \left(12000000 - \frac{262122237 \pi}{80}\right) \operatorname{w_{3}}{\left(x \right)} + 1800\right)}{202209199875 \pi}\]
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_21856\483021684.py in <module>
25 display(diffeq1,diffeq2,diffeq3,diffeq4)
26
---> 27 sol = sp.dsolve([diffeq1,diffeq2,diffeq3,diffeq4],[w1(x),w2(x),w3(x),w4(x)])
28 display(sol)
~\Anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
557 # been solved.
558 try:
--> 559 sol = dsolve_system(eq, funcs=func, ics=ics, doit=True)
560 return sol[0] if len(sol) == 1 else sol
561 except NotImplementedError:
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in dsolve_system(eqs, funcs, t, ics, doit, simplify)
2109 for canon_eq in canon_eqs:
2110 try:
-> 2111 sol = _strong_component_solver(canon_eq, funcs, t)
2112 except NotImplementedError:
2113 sol = None
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in _strong_component_solver(eqs, funcs, t)
1717
1718 elif match.get('is_linear', False):
-> 1719 sol = _linear_ode_solver(match)
1720
1721 # Note: For now, only linear systems are handled by this function
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in _linear_ode_solver(match)
1612 type = match['type_of_equation']
1613
-> 1614 sol_vector = linodesolve(A, t, b=rhs, B=B,
1615 type=type, tau=tau)
1616
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in linodesolve(A, t, b, B, type, doit, tau)
975
976 if type in ("type1", "type2", "type5", "type6"):
--> 977 P, J = matrix_exp_jordan_form(A, t)
978 P = simplify(P)
979
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in matrix_exp_jordan_form(A, t)
648 return chains
649
--> 650 eigenchains = jordan_chains(A)
651
652 # Needed for consistency across Python versions
~\Anaconda3\lib\site-packages\sympy\solvers\ode\systems.py in jordan_chains(A)
635 where vijk is the kth vector in the jth chain for eigenvalue i.
636 '''
--> 637 P, blocks = A.jordan_cells()
638 basis = [P[:,i] for i in range(P.shape[1])]
639 n = 0
~\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in jordan_cells(self, calc_transformation)
732
733 def jordan_cells(self, calc_transformation=True):
--> 734 P, J = self.jordan_form()
735 return P, J.get_diag_blocks()
736
~\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in jordan_form(self, calc_transform, **kwargs)
417
418 def jordan_form(self, calc_transform=True, **kwargs):
--> 419 return _jordan_form(self, calc_transform=calc_transform, **kwargs)
420
421 def left_eigenvects(self, **flags):
~\Anaconda3\lib\site-packages\sympy\matrices\eigen.py in _jordan_form(M, calc_transform, chop)
1214 return restore_floats(jordan_mat)
1215
-> 1216 jordan_basis = [eig_mat(eig, 1).nullspace()[0]
1217 for eig in blocks]
1218 basis_mat = mat.hstack(*jordan_basis)
~\Anaconda3\lib\site-packages\sympy\matrices\eigen.py in <listcomp>(.0)
1214 return restore_floats(jordan_mat)
1215
-> 1216 jordan_basis = [eig_mat(eig, 1).nullspace()[0]
1217 for eig in blocks]
1218 basis_mat = mat.hstack(*jordan_basis)
~\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in nullspace(self, simplify, iszerofunc)
352
353 def nullspace(self, simplify=False, iszerofunc=_iszero):
--> 354 return _nullspace(self, simplify=simplify, iszerofunc=iszerofunc)
355
356 def rowspace(self, simplify=False):
~\Anaconda3\lib\site-packages\sympy\matrices\subspaces.py in _nullspace(M, simplify, iszerofunc)
62 """
63
---> 64 reduced, pivots = M.rref(iszerofunc=iszerofunc, simplify=simplify)
65
66 free_vars = [i for i in range(M.cols) if i not in pivots]
~\Anaconda3\lib\site-packages\sympy\matrices\matrices.py in rref(self, iszerofunc, simplify, pivots, normalize_last)
173 def rref(self, iszerofunc=_iszero, simplify=False, pivots=True,
174 normalize_last=True):
--> 175 return _rref(self, iszerofunc=iszerofunc, simplify=simplify,
176 pivots=pivots, normalize_last=normalize_last)
177
~\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _rref(M, iszerofunc, simplify, pivots, normalize_last)
303 simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
304
--> 305 mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
306 normalize_last, normalize=True, zero_above=True)
307
~\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce(M, iszerofunc, simpfunc, normalize_last, normalize, zero_above)
125 normalize=True, zero_above=True):
126
--> 127 mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
128 iszerofunc, simpfunc, normalize_last=normalize_last,
129 normalize=normalize, zero_above=zero_above)
~\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc, normalize_last, normalize, zero_above)
107 continue
108
--> 109 cross_cancel(pivot_val, row, val, piv_row)
110 piv_row += 1
111
~\Anaconda3\lib\site-packages\sympy\matrices\reductions.py in cross_cancel(a, i, b, j)
56 q = (j - i)*cols
57 for p in range(i*cols, (i + 1)*cols):
---> 58 mat[p] = isimp(a*mat[p] - b*mat[p + q])
59
60 isimp = _get_intermediate_simp(_dotprodsimp)
~\Anaconda3\lib\site-packages\sympy\simplify\simplify.py in dotprodsimp(expr, withsimp)
2112
2113 if expr2 is expr or count_ops_alg(expr2)[0] >= 6: # check again after substitution
-> 2114 expr3 = cancel(expr2)
2115
2116 if expr3 != expr2:
~\Anaconda3\lib\site-packages\sympy\polys\polytools.py in cancel(f, _signsimp, *gens, **args)
6799 return f.xreplace(dict(reps))
6800
-> 6801 c, (P, Q) = 1, F.cancel(G)
6802 if opt.get('polys', False) and 'gens' not in opt:
6803 opt['gens'] = R.symbols
~\Anaconda3\lib\site-packages\sympy\polys\rings.py in cancel(self, g)
2221
2222 if not (domain.is_Field and domain.has_assoc_Ring):
-> 2223 _, p, q = f.cofactors(g)
2224 else:
2225 new_ring = ring.clone(domain=domain.get_ring())
~\Anaconda3\lib\site-packages\sympy\polys\rings.py in cofactors(f, g)
2137
2138 J, (f, g) = f.deflate(g)
-> 2139 h, cff, cfg = f._gcd(g)
2140
2141 return (h.inflate(J), cff.inflate(J), cfg.inflate(J))
~\Anaconda3\lib\site-packages\sympy\polys\rings.py in _gcd(f, g)
2172 return f._gcd_ZZ(g)
2173 else: # TODO: don't use dense representation (port PRS algorithms)
-> 2174 return ring.dmp_inner_gcd(f, g)
2175
2176 def _gcd_ZZ(f, g):
~\Anaconda3\lib\site-packages\sympy\polys\compatibility.py in dmp_inner_gcd(self, f, g)
661 return (self.from_dense(H), self.from_dense(F), self.from_dense(G))
662 def dmp_inner_gcd(self, f, g):
--> 663 H, F, G = dmp_inner_gcd(self.to_dense(f), self.to_dense(g), self.ngens-1, self.domain)
664 return (self.from_dense(H), self.from_dense(F), self.from_dense(G))
665 def dup_gcd(self, f, g):
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_inner_gcd(f, g, u, K)
1581
1582 J, (f, g) = dmp_multi_deflate((f, g), u, K)
-> 1583 h, cff, cfg = _dmp_inner_gcd(f, g, u, K)
1584
1585 return (dmp_inflate(h, J, u, K),
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in _dmp_inner_gcd(f, g, u, K)
1554 pass
1555
-> 1556 return dmp_rr_prs_gcd(f, g, u, K)
1557
1558
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_rr_prs_gcd(f, g, u, K)
1057 return dup_rr_prs_gcd(f, g, K)
1058
-> 1059 result = _dmp_rr_trivial_gcd(f, g, u, K)
1060
1061 if result is not None:
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in _dmp_rr_trivial_gcd(f, g, u, K)
909 return dmp_one(u, K), f, g
910 elif query('USE_SIMPLIFY_GCD'):
--> 911 return _dmp_simplify_gcd(f, g, u, K)
912 else:
913 return None
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in _dmp_simplify_gcd(f, g, u, K)
951 G = dmp_content(g, u, K)
952 else:
--> 953 F = dmp_content(f, u, K)
954 G = dmp_LC(g, K)
955
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_content(f, u, K)
1790
1791 for c in f[1:]:
-> 1792 cont = dmp_gcd(cont, c, v, K)
1793
1794 if dmp_one_p(cont, v, K):
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_gcd(f, g, u, K)
1622
1623 """
-> 1624 return dmp_inner_gcd(f, g, u, K)[0]
1625
1626
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_inner_gcd(f, g, u, K)
1581
1582 J, (f, g) = dmp_multi_deflate((f, g), u, K)
-> 1583 h, cff, cfg = _dmp_inner_gcd(f, g, u, K)
1584
1585 return (dmp_inflate(h, J, u, K),
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in _dmp_inner_gcd(f, g, u, K)
1554 pass
1555
-> 1556 return dmp_rr_prs_gcd(f, g, u, K)
1557
1558
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_rr_prs_gcd(f, g, u, K)
1065 gc, G = dmp_primitive(g, u, K)
1066
-> 1067 h = dmp_subresultants(F, G, u, K)[-1]
1068 c, _, _ = dmp_rr_prs_gcd(fc, gc, u - 1, K)
1069
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_subresultants(f, g, u, K)
548
549 """
--> 550 return dmp_inner_subresultants(f, g, u, K)[0]
551
552
~\Anaconda3\lib\site-packages\sympy\polys\euclidtools.py in dmp_inner_subresultants(f, g, u, K)
511 dmp_pow(c, d, v, K), v, K)
512
--> 513 h = dmp_prem(f, g, u, K)
514 h = [ dmp_quo(ch, b, v, K) for ch in h ]
515
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_prem(f, g, u, K)
1236
1237 R = dmp_mul_term(r, lc_g, 0, u, K)
-> 1238 G = dmp_mul_term(g, lc_r, j, u, K)
1239 r = dmp_sub(R, G, u, K)
1240
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul_term(f, c, i, u, K)
182 return dmp_zero(u)
183 else:
--> 184 return [ dmp_mul(cf, c, v, K) for cf in f ] + dmp_zeros(i, v, K)
185
186
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in <listcomp>(.0)
182 return dmp_zero(u)
183 else:
--> 184 return [ dmp_mul(cf, c, v, K) for cf in f ] + dmp_zeros(i, v, K)
185
186
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
826
827 for j in range(max(0, i - dg), min(df, i) + 1):
--> 828 coeff = dmp_add(coeff, dmp_mul(f[j], g[i - j], v, K), v, K)
829
830 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dmp_mul(f, g, u, K)
805 """
806 if not u:
--> 807 return dup_mul(f, g, K)
808
809 if f == g:
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dup_mul(f, g, K)
783 lo, hi = dup_mul(fl, gl, K), dup_mul(fh, gh, K)
784
--> 785 mid = dup_mul(dup_add(fl, fh, K), dup_add(gl, gh, K), K)
786 mid = dup_sub(mid, dup_add(lo, hi, K), K)
787
~\Anaconda3\lib\site-packages\sympy\polys\densearith.py in dup_mul(f, g, K)
765
766 for j in range(max(0, i - dg), min(df, i) + 1):
--> 767 coeff += f[j]*g[i - j]
768
769 h.append(coeff)
~\Anaconda3\lib\site-packages\sympy\polys\domains\gaussiandomains.py in __add__(self, other)
71
72 def __add__(self, other):
---> 73 x, y = self._get_xy(other)
74 if x is not None:
75 return self.new(self.x + x, self.y + y)
~\Anaconda3\lib\site-packages\sympy\polys\domains\gaussiandomains.py in _get_xy(cls, other)
68 except CoercionFailed:
69 return None, None
---> 70 return other.x, other.y
71
72 def __add__(self, other):
KeyboardInterrupt: