On 08/14/2013 12:53 PM, Ondřej Čertík wrote:
On Tue, Aug 13, 2013 at 7:00 PM, Alan Bromborsky <[email protected]> wrote:
On 08/13/2013 08:50 PM, Dale Lukas Peterson wrote:
Also, try the ".series()" command. It might work, as it is expanding
the expression inside out,
to avoid such problems.
Also, if you can put the equations in the form: dx/dt = f(x), and you have
f(x) as a Sympy Matrix of expressions, you can simply do df = f.jacobian(x),
then evaluate (using subs) at the equilibrium conditions and values of
parameters that you are interested in.
Luke
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an
email to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
series did not work. limits did for evaluating the derivatives. First
derivatives were quick. Second derivatives took a long time (15min).
Can you post the sympy script that takes 15 min?
I'd like to see what the problem is.
Ondrej
This is very strange. I tried to run the problem again and it never
finished (either the first or second derivatives). I have attached the
problem with everything commented out that depends on a custom version
of GA that I am working that generated the differential equations for
the pendulum. Then I set up the Laplace transforms for the equations by
hand for the second half of the problem. What is attached depends only
on sympy without GA.
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
from sympy import symbols,sin,cos,pi,trigsimp,expand,ratsimp, Function, \
Matrix, Inverse, denom, numer, S, simplify, apart, \
factor, solve, I, exp, re, im, Symbol, sqrt,series,limit
#from ga import Ga
#from printer import Eprint,xpdf,Format,latex
#from metric import linear_expand
import math
def Heavyside_expand(P,Q,wlst,s,t):
dQ = Q.diff(s)
E = S(0)
for w in wlst:
Pp = P.subs(s,I*w)
Pm = P.subs(s,-I*w)
dQp = dQ.subs(s,I*w)
dQm = dQ.subs(s,-I*w)
tmp =
simplify(((Pp/dQp)*exp(I*w*t).expand(complex=True)+(Pm/dQm)*exp(-I*w*t).expand(complex=True)).expand(complex=True))
E += tmp
return E
def main():
#Eprint()
#Format()
Te_value = 86400. # day in seconds
Re_value = 6.371e6 # earth radius in meters
T_value = 10 # pendulum period in seconds
"""
sp3d = Ga.preset('sph3d') # 3d spherical coordinate system
er,eth,ephi = sp3d.mv() # basis vectors
r,th,phi = sp3d.coords # spherical coordinates
"""
t, Re = symbols('t,R_E',real=True)
w_th, w_phi, K = symbols('omega_theta,omega_phi,K',real=True)
wp, wm, s, delta = symbols('omega_{+},omega_{-},s,delta',real=True)
W, w, x = symbols('omega_E,omega,x',real=True)
"""
sp3d.parametric([0,0,W*t]) # define phi=W*t
d2er = (er.diff(t)).diff(t) # second time derivative of er
u_th = Function('u_theta')(t) # North South displacement variable
u_phi = Function('u_phi')(t) # East West displacement variable
xp = u_th*eth+u_phi*ephi # r_perp
d2xp = (xp.diff(t)).diff(t) # second time derivative of r_perp
eq = d2xp+w**2*xp-Re*d2er # equation of motion for pendulum eq=0
print r'%\mbox{Small Angle Approximation: } \bm{r} =
u_{\theta}\bm{e}_{\theta}+u_{\phi}\bm{e}_{\phi}-l\bm{e}_{r}'
print r'#Equation of Motion:'
eq.Fmt(3,'0')
up_th = Function("u'_theta")(t) # define u'_theta
eq_lst = eq.list()
eq_th = eq_lst[1] # equation for e_theta component
eq_phi = eq_lst[2] #equation for e_phi component
print r'#Transform Equilibrium Position:'
print r"%(\omega^2-\omega_{E}^{2}\f{\cos}{\theta}^{2})u'_{\theta} =
(\omega^{2}-\omega_{E}^{2}\f{\cos^{2}}{\theta})u_{\theta}+R_{E}\omega_{E}^{2}\f{\sin}{2\theta}/2"
print r"%\mbox{Pendulum Equilibrium ($u'_{\theta} = 0$): }u_{\theta} =
R_{E}\omega_{E}^{2}\f{\sin}{2\theta}/2
(\omega^{2}-\omega_{E}^{2}\f{\cos^{2}}{\theta})"
print r'#Transformed Equations of Motion ($r$ motion of higher order in
small angles):'
# substitute in e_theta component equation
eq_th =
eq_th.subs({u_th:up_th,Re*W**2*sin(2*th):0,u_th.diff(t,2):up_th.diff(t,2)})
# substitute in e_phi component equation
eq_phi = eq_phi.subs(u_th,up_th)
print '% ',factor(eq_th,up_th),'= 0'
print '% ',factor(eq_phi,u_phi),'= 0'
print r'% \omega_{\theta}^{2} =
\omega^{2}-\omega_{E}^{2}\f{\sin^{2}}{\theta}'
print r'% \omega_{\phi}^{2} = \omega^{2}-\omega_{E}^{2}'
print r'% K = 2\omega_{E}\f{\sin}{\theta}'
eq_th = simplify(eq_th.subs(w**2,w_th**2+W**2*x**2))
eq_th = eq_th.subs(2*W*x,K)
print '% ',eq_th,'= 0'
eq_phi = simplify(eq_phi.subs(w**2,w_phi**2+W**2))
eq_phi = eq_phi.subs(2*W*x,K)
print '% ',eq_phi,'= 0'
"""
M = Matrix([[s**2+w_th**2,K*s],[-K*s,s**2+w_phi**2]])
X = Matrix([[K],[s]])
print 'X =',X
X = M.inv()*X
U_th = ratsimp(X[0])
U_phi = ratsimp(X[1])
print r'%U_{\theta} =',U_th
print r'%U_{\phi} =',U_phi
P_th = numer(U_th)
P_phi = numer(U_phi)
Q = denom(U_th)
Qp = (s**2+wp**2)*(s**2+wm**2)
print r'% ',Qp,' = ',Q
Wpm = solve(Q,s**2)
Wp2 = (-Wpm[0]).subs({K:2*W*x,w_th**2:w**2-W**2*x,w_phi**2:w**2-W**2})
Wm2 = (-Wpm[1]).subs({K:2*W*x,w_th**2:w**2-W**2*x,w_phi**2:w**2-W**2})
print r'%\color{red}\omega_{+}^{2} =',expand(Wp2)
print r'%\color{red}\omega_{-}^{2} =',expand(Wm2)
print r'%\color{red}\partial_{\omega_{E}}\omega_{+}^{2}
=',expand(Wp2).diff(W)
print r'%\color{red}\partial_{\omega_{E}}\omega_{-}^{2}
=',expand(Wm2).diff(W)
print r'%\color{red}\partial_{\omega_{E}}\omega_{+}
=',limit(sqrt(expand(Wp2)).diff(W),W,0,dir='+')
print r'%\color{red}\partial_{\omega_{E}}\omega_{-}
=',limit(sqrt(expand(Wm2)).diff(W),W,0,dir='+')
print r'%\color{red}\partial_{\omega_{E}}^{2}\omega_{+}
=',limit(sqrt(expand(Wp2)).diff(W).diff(W),W,0,dir='+')
print r'%\color{red}\partial_{\omega_{E}}^{2}\omega_{-}
=',limit(sqrt(expand(Wm2)).diff(W).diff(W),W,0,dir='+')
B = K**2+w_th**2+w_phi**2
d = expand(B**2-4*(w_th**2)*(w_phi**2))
d = simplify(d.subs({K:2*W*x,w_th**2:w**2-W**2*x,w_phi**2:w**2-W**2}))
print 'd =', d
u0_phi = symbols('u_phi__0',real=True)
print r'u_{\theta} =',Heavyside_expand(P_th,Qp,[wp,wm],s,t)*u0_phi
print r'u_{\phi} =',Heavyside_expand(P_phi,Qp,[wp,wm],s,t)*u0_phi
#xpdf()
#xpdf(paper='letter')
if __name__ == "__main__":
main()
