SymPy 1.11 released

[Oscar Benjamin]

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

[Peter Stahlecker]

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();

#===================================================================================================

"""

--
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 sympy+un...@googlegroups.com.
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

[Oscar Benjamin]

On Wed, 24 Aug 2022 at 08:55, Peter Stahlecker
<peter.st...@gmail.com> wrote:
>
> 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.

It looks like cse=True doesn't work properly when lambdifying a Matrix:

In [3]: lambdify((x,), Matrix([[x**2, 0], [0, x**2]]), cse=False)(1)
Out[3]:
array([[1, 0],
[0, 1]])

In [4]: lambdify((x,), Matrix([[x**2, 0], [0, x**2]]), cse=True)(1)
Out[4]:
[array([[1, 0],
[0, 1]])]

Note that with cse=True the Matrix is in a list of length 1. If the
matrix doesn't have any nontrivial subexpressions then it doesn't
happen:

In [6]: lambdify((x,), Matrix([[x, 0], [0, x]]), cse=True)(1)
Out[6]:
array([[1, 0],
[0, 1]])

Looking at the generated code we have:

In [11]: print(inspect.getsource(lambdify((x,), Matrix([[x**2, 0], [0,
x**2]]), cse=True)))
def _lambdifygenerated(x):
x0 = x**2
return [array([[x0, 0], [0, x0]])]

In [12]: print(inspect.getsource(lambdify((x,), Matrix([[x**2, 0], [0,
x**2]]), cse=False)))
def _lambdifygenerated(x):
return array([[x**2, 0], [0, x**2]])

So for some reason the generated code puts the output in a list. Looks
like this code is responsible:
https://github.com/sympy/sympy/blob/5eb59bdad2f4e1a675acb31e7f0c72c8abca708c/sympy/utilities/lambdify.py#L1114-L1126

It seems to presume that expr (which in this case is a Matrix) is
supposed to be a list or tuple. First it tries to add a list to it and
then if that fails it tries to add a tuple and then if that fails it
just puts expr into a list (expr = [expr]). I wish I could just ban
the use of try/except in the sympy codebase because it's almost never
used in a good way.

Basically that code will fail any time expr is not a tuple or a list e.g.:

In [3]: lambdify(x, x**2 + y*x**2, cse=True)(1)
Out[3]: [y + 1]

A workaround is just to extract the element from the list:

In [4]: lambdify(x, x**2 + y*x**2, cse=True)(1)[0]
Out[4]: y + 1

The bug comes from here:
https://github.com/sympy/sympy/pull/23538

--
Oscar

[Peter Stahlecker]

Dear Oscar,

my expressions into lambdify are of sm.Matrix type.

With cse=False, the output of lambdify is a numpy.ndarray of the same shape as the sm.Matrix.

With cse=True, the output is is list.

(I guess, this is what you said, would happen.)

Does this mean, if I somehow manage to convert the lists into ndarrays of the correct shape, the code should run with cse=True, too?

Thanks, Peter

--
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 sympy+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSoDEGZLfZ%3D0JYkdoQ%3D0Dq%3DSVP_MqnmOvp9%2Bp7Q__uxSg%40mail.gmail.com.

[Peter Stahlecker]

Dear Oscar,

Thanks to your explanations, I managed  to use solve_ivp when I used cse=True in  lambdify in one of my programs. ( just convert the list to an np.ndarray of the correct shape)

The increase in speed in the integration was remarkable:

From around 8.5 sec with cse=False to around 0.13 sec with cse=True, 65 times faster!

Thanks! Peter

[Jeremy Monat]

Thanks Oscar. Great to have the default (latest) version of the docs in the new, Furo format.

The release notes include

SymPy 1.11 has not been released yet.

Would you like me to update that to

SymPy 1.11 was released on 23rd August 2022.

or would you like to (or is some automation supposed to do that)?

Best,

Jeremy

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CABKqA0YKmQ3RwByJMwmtXKS5rXi-SZ42tGeQrjChk%3D_UQ90Q9Q%40mail.gmail.com.

[Oscar Benjamin]

Good point Jeremy.

Feel free to change it. Editing the release notes is a completely
manual process. I'd like to make it more automatic but the first step
for that is probably moving them out of the wiki and into the repo and
making them part of the docs.

Oscar

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAO00iLhC9Np7PtbAKe4RiuyA2gCYrG5XDHbYSpf1LDnh3sjUiA%40mail.gmail.com.

[Arthur Ryman]

Congratulations to the dev team for the improved documentation. It looks much cleaner and more modern.

-- Arthur

--
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 sympy+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRfZbRex4Kg%2B%2BuEVBG99Td8ZNmSo39GjWsPPLo0m0mhpw%40mail.gmail.com.

[Jeremy Monat]

Oscar, I updated the Release Notes page.

Arthur, glad you like the new documentation! If you'd like to help with the documentation, you can check out the issues tagged Documentation and contribute to any.

Jeremy

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAApBiOn2Gw%2BRBPBP%3D6VaYF4w6QmfQrZyXFqW1zF%3D3H8ukEQYMA%40mail.gmail.com.

[Arthur Ryman]

Jeremy,

That's quite a backlog of doc issues! I'll keep it in mind as a place to contribute.

-- Arthur

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAO00iLh-xYNtpZXfjbqsP8tp5T9COQ%2B1%3DeNdVETg4TVQzaadtA%40mail.gmail.com.