I have upgraded to sympy 1.11
I wanted to try the cse keyword.
If I set cse = False, all seems to work fine.
If I set cse = True, lambdify(..) seems to work fine, but solve_ivp(..)
gives and error.
Attached is one program, where I noticed this. I tried on other programs,
same results.
Any help would be appreciated!
"""
import sympy as sm
import sympy.physics.mechanics as me
import time
import numpy as np
from scipy.integrate import odeint, solve_ivp
from scipy.optimize import fsolve, minimize
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
'''
A ball is rolling without sliding on a horizontal plane.
'''
start = time.time()
#============================================
cse = False
#============================================
q1, q2, q3 = me.dynamicsymbols('q1 q2 q3') # q is the angle of the ball
u1, u2, u3 = me.dynamicsymbols('u1 u2 u3')
x, z = me.dynamicsymbols('x z') # ccordinates of contact point
ux, uz = me.dynamicsymbols('ux uz')
m, g, r, iXX, iYY, iZZ, iXY, iXZ, iYZ = sm.symbols('m, g, r, iXX, iYY, iZZ,
iXY, iXZ, iYZ')
amplitude, frequenz, reibung, t = sm.symbols('amplitude frequenz reibung t')
# determine distance of center of mass from geometric center: is all zero,
same location, if all are 1
# mass center at the edge of the ball
alfa, beta, gamma = sm.symbols('alfa beta gamma')
# H1, H2 are auxiliary frames, so orient_axis may be used, see below. One
could also do without them and
# use orient_body_fixed, but the operations in the right hand side of the
differential equations (here
# called force) are generally a bit more if orient_body_fixed is used.
N = me.ReferenceFrame('N') # Inertial frame
H1 = me.ReferenceFrame('H1') # auxiliary frame
H2 = me.ReferenceFrame('H2') # dto.
A1 = me.ReferenceFrame('A1') # fixed to the ball
P0 = me.Point('P0') # fixed reference point
CP = me.Point('CP') # contact point of ball and plane
Dmc = me.Point('Dmc') # geometric center of the ball
m_Dmc = me.Point('m_Dmc') # mass center
H1.orient_axis(N, q1, N.x)
H2.orient_axis(H1, q2, H1.y)
A1.orient_axis(H2, q3, H2.z)
# rot, rot1 are needed for the kinnematic equations, see below.
rot = A1.ang_vel_in(N)
A1.set_ang_vel(N, u1*A1.x + u2*A1.y + u3*A1.z)
rot1 = A1.ang_vel_in(N)
'''
Locate the contact point CP, and the center of gravity, Dmc, at distance r
As the CP is obviously (also) on the plane, it's Y coordinate is zero, and
it does have zero
velocity. Of course, subsequent contact points are at different locations
on the plane
this is calculated further down.
'''
CP.set_pos(P0, x*N.x + z*N.z)
CP.set_vel(N, 0.) # this is the contact point to the street, hence no
velocity
Dmc.set_pos(CP, r * N.y)
Dmc.v2pt_theory(CP, N, A1)
m_Dmc.set_pos(Dmc, alfa*r*A1.x + beta*r*A1.y + gamma*r*A1.z)
m_Dmc.v2pt_theory(Dmc, N, A1)
'''
Get the differential equations of the X, Z coordinates of the contact
point. They will be
integrated numerically.
OMEGA is the rotational speed of the ball.
u3_wirk is its component in N.z direction, analoguous u1_wirk.
hence: d/dt(X) = -u3_wirk * r (r = radius of the ball). The minus sign here
is due to this:
if u3_wirk points in the positive z - direction, the right hand rule says,
that X changes in the
negative x - direction.
'''
OMEGA = u1*A1.x + u2*A1.y + u3*A1.z
u3_wirk = me.dot(OMEGA, N.z)
u1_wirk = me.dot(OMEGA, N.x)
rhsx = -u3_wirk * r
rhsz = u1_wirk * r
subs_dict1 = {sm.Derivative(x, t): rhsx, sm.Derivative(z, t): rhsz}
I = me.inertia(A1, iXX, iYY, iZZ, iXY, iXZ, iYZ)
Body = me.RigidBody('Body', m_Dmc, A1, m, (I, Dmc))
BODY = [Body]
# A necessary, but by no means sufficient! condition for the correctness of
the equations of motion
# is that, absent any friction, the total energy of the system must be
constant
pot_energie = m * g * me.dot(m_Dmc.pos_from(P0), N.y).subs(subs_dict1)
kin_energie = Body.kinetic_energy(N).subs(subs_dict1)
FL1 = [(m_Dmc, -m*g*N.y)]
Torque = [(A1, -reibung * (u1*A1.x + u2*A1.y + u3*A1.z))]
FL = FL1 + Torque
kd = [me.dot(rot - rot1, uv) for uv in A1]
q = [q1, q2, q3]
u = [u1, u2, u3]
KM = me.KanesMethod(N, q_ind=q, u_ind=u, kd_eqs=kd)
(fr, frstar) = KM.kanes_equations(BODY, FL)
MM = KM.mass_matrix_full
force = KM.forcing_full
# Add rhsx, rhsz at the bottom of force, to get d/dt(x) = rhsx, same for z.
This is to numerically integrate x(t), z(t)
force = sm.Matrix.vstack(force, sm.Matrix([rhsx, rhsz]))
# It is a good idea to look at the dynamic symbols and the free symbols,
can give a clue, where
# something might be wrong in the equations.
print('force DS', me.find_dynamicsymbols(force))
print('force free symbols', force.free_symbols)
print('force has {}
operations'.format(sum([force[i].count_ops(visual=False) for i in
range(force.shape[0])])), '\n')
# Enlarge MM properly
MM = sm.Matrix.hstack(MM, sm.zeros(6, 2))
hilfs = sm.Matrix.hstack(sm.zeros(2, 6), sm.eye(2))
MM = sm.Matrix.vstack(MM, hilfs )
print('MM DS', me.find_dynamicsymbols(MM))
print('MM free symbols', MM.free_symbols)
print('MM has {} operations'.format(sum([MM[i, j].count_ops(visual=False)
for i in range(MM.shape[0]) for j in range(MM.shape[1])])), '\n')
# Lambdification
pL = [m, g, r, iXX, iYY, iZZ, iXY, iXZ, iYZ, reibung, alfa, beta, gamma]
qL = q + u + [x, z]
MM_lam = sm.lambdify(qL + pL, MM, cse=cse)
force_lam = sm.lambdify(qL + pL, force, cse=cse)
pot_lam = sm.lambdify(qL + pL, pot_energie, cse=cse)
kin_lam = sm.lambdify(qL + pL, kin_energie, cse=cse)
print('it took {:.3f} sec to establish Kanes equations'.format(time.time()
- start), '\n')
# Integrate numerically and plot results
start = time.time()
#================================================================
mm = 1. # mass of ball
rr = 1. # radius of ball
iXXe = 2./5. * mm * rr**2 # from the internet
pL_vals = [mm, 9.8, rr, iXXe, iXXe, iXXe, 0.1*iXXe, 0.1*iXXe, 0.1*iXXe,
0.0, 0.1, 0.1, 0.1 ]
y0 = [0., 0., 0., 1., 1., 1.] + [0., 0.]
intervall = 10.
schritte = 200
#================================================================
#startwert = y0[2] # just needed for the plots below
#startomega = y0[1] # dto.
times = np.linspace(0, intervall, schritte)
def gradient(t, y, args):
vals = np.concatenate((y, args))
sol = np.linalg.solve(MM_lam(*vals), force_lam(*vals))
return np.array(sol).T[0]
t_span = (0., intervall)
resultat1 = solve_ivp(gradient, t_span, y0, t_eval=times, args=(pL_vals,))
resultat = resultat1.y.T
event_dict = {-1: 'Integration failed', 0: 'Integration finished
successfully', 1: 'some termination event'}
print(event_dict[resultat1.status])
print("To numerically integrate an intervall of {} sec, the routine cycled
{} times and it took {:.3f} sec ".format(intervall, resultat1.nfev,
time.time() - start), '\n')
#=======================================================================================================
# plot results
fig, ax = plt.subplots(figsize=(15,10))
for i, j in zip(range((resultat.shape[1])), ('q1', 'q2', 'q3', 'u1', 'u2',
'u3', 'X coordinate of CP', 'Y coordinate of CP')):
ax.plot(times, resultat[:, i], label=j)
ax.set_title("Generalized Coordinates")
ax.legend();
kin_np = np.empty(schritte)
pos_np = np.empty(schritte)
total_np = np.empty(schritte)
for i in range(schritte):
kin_np[i] = kin_lam(*[resultat[i, j] for j in
range(resultat.shape[1])], *pL_vals)
pos_np[i] = pot_lam(*[resultat[i, j] for j in
range(resultat.shape[1])], *pL_vals)
total_np[i] = kin_np[i] + pos_np[i]
fig, ax = plt. subplots(figsize=(15, 10))
ax.plot(times, kin_np - kin_np[0], label='kinetic energy')
ax.plot(times, pos_np - pos_np[0], label='pos energy')
ax.plot(times, total_np - total_np[0], label='total energy')
ax.set_title("Energy of the disc, normed so initial energies = 0.")
ax. legend();
#===================================================================================================
"""
On Wed 24. Aug 2022 at 02:17 Oscar Benjamin <[email protected]>
wrote:
> Hi all,
>
> I've just pushed the release for SymPy 1.11 to PyPI and updated the docs.
>
> To upgrade to 1.11 you can run
>
> $ pip install --upgrade sympy
>
> The release files can also be downloaded from GitHub here:
> https://github.com/sympy/sympy/releases/tag/sympy-1.11
>
> I presume that the release will also be available via conda soon (I'm
> actually not sure who/what makes that happen...).
>
> The release notes for 1.11 are here:
> https://github.com/sympy/sympy/wiki/Release-Notes-for-1.11
>
> As always please report any bugs or other issues with the release to
> GitHub:
> https://github.com/sympy/sympy/issues
>
> The following people contributed at least one patch to this release (names
> are
> given in alphabetical order by last name). A total of 70 people
> contributed to this release. People with a * by their names contributed a
> patch for the first time for this release; 32 people contributed
> for the first time for this release.
>
> Thanks to everyone who contributed to this release!
>
> - Saksham Alok*
> - anutosh491
> - Alexander Behrens*
> - Oscar Benjamin
> - Varenyam Bhardwaj*
> - Anurag Bhat
> - Akshansh Bhatt
> - Francesco Bonazzi
> - Islem BOUZENIA*
> - Arie Bovenberg*
> - Riley Britten*
> - Zaz Brown*
> - cocolato*
> - Björn Dahlgren
> - Nikhil Date
> - dispasha*
> - Pieter Eendebak*
> - extraymond*
> - Tom Fryers
> - Mauro Garavello
> - Anton Golovanov*
> - Oscar Gustafsson
> - Miro Hrončok
> - Siddhant Jain
> - Gaurav Jain*
> - Kuldeep Borkar Jr
> - Georges Khaznadar*
> - Steve Kieffer
> - Sergey B Kirpichev
> - Matthias Köppe
> - Sumit Kumar*
> - S.Y. Lee
> - Qijia Liu
> - Riccardo Magliocchetti
> - Tirthankar Mazumder*
> - Carson McManus*
> - Ehren Metcalfe
> - Aaron Meurer
> - Ansh Mishra
> - Jeremy Monat
> - Jason Moore
> - oittaa*
> - Omkaar*
> - Andrii Oriekhov*
> - Julien Palard
> - Advait Pote
> - Seth Poulsen*
> - Mamidi Ratna Praneeth
> - Psycho-Pirate
> - Mayank Raj
> - rathmann
> - ritikBhandari*
> - Arthur Ryman*
> - Shivam Sagar*
> - Gilles Schintgen
> - Jack Schmidt*
> - Sudeep Sidhu
> - Chris Smith
> - Jisoo Song
> - Paul Spiering
> - Timo Stienstra*
> - Kalevi Suominen
> - Luis Talavera*
> - tttc3*
> - user202729*
> - Orestis Vaggelis
> - viocha*
> - Brandon T. Willard
> - Donald Wilson*
> - Xiang Wu*
>
> The SHA256 hashes for the release files are here:
>
> 1fe96b4c56bb7a7630cdf150a6cd98bc97a79e6be233e30502aba1cf54dee33d
> sympy-1.11.tar.gz
> b53069f5f30e4a4690b57cdb8e3d0d9065fff42627239db718214f804e442481
> sympy-1.11-py3-none-any.whl
> 1c056a67ea8bd184a336e1cb988934750deb305c132b6238177c099357a2243f
> sympy-docs-html-1.11.zip
> 281f7142a95b102e8df727fe7cefdf941c10f178f844d0d55c7a7342c79a2d5b
> sympy-docs-pdf-1.11.pdf
>
> --
> Oscar
>
> --
> 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 view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/CAHVvXxRfZbRex4Kg%2B%2BuEVBG99Td8ZNmSo39GjWsPPLo0m0mhpw%40mail.gmail.com
> .
>
--
Best regards,
Peter Stahlecker
--
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 view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CABKqA0ZGkSmC%3Dd6v%3Dg5krQGv%3D3PHh1hkZXsiZ%3DfXC0Jf%3DWd2%2Bw%40mail.gmail.com.
