Duffing Spring Simulation

[eh111eh]

Hi, As part of implementing the DuffingSpring class in sympy.mechanics, I've tried to simulate and verify its functionality. In the DuffingSpring Simulation.ipynb.zip (I attached a zip file because GitHub doesn't support ipynb files. Sorry for any inconvenience), I first define the problem by deriving the equation of motion for the Duffing system, numerically solve the equation using solve_ivp from scipy.integrate, and then visualise the system. The plots include:

• Displacement over time (time in seconds vs. displacement in meters)

• Velocity over time (time in seconds vs. velocity in meters/second)

• Phase space plot (displacement in meters vs. velocity in meters/second)

To explore the behaviour of the equation with various alpha and beta values, I simulate different combinations of parameters. I also extend the time span for a longer simulation and add a plot for total mechanical energy.

The plot of the Duffing oscillator's restoring force validates its functionality and accuracy since the resulting plots align well with expectations and reference literature.

Additionally, I consider adding examples in the docstring as follows:

Examples
========
Let's create an instance of a Duffing spring with specific parameters and a linear pathway. We'll simulate its response to an initial displacement:

>>> from sympy import symbols, Function, cos
>>> from sympy.physics.mechanics import LinearPathway, Point, ReferenceFrame
>>> from sympy.physics.vector import dynamicsymbols
>>> N = ReferenceFrame('N')
>>> P, Q = Point('P'), Point('Q')
>>> x = dynamicsymbols('x')
>>> P.set_pos(Q, x * N.x)  # Displacement along x direction of frame N
>>> pathway = LinearPathway(P, Q)
>>> beta, alpha = symbols('beta alpha', real=True)
>>> duffing_spring = DuffingSpring(beta, alpha, pathway)
>>> duffing_spring.force
-beta*x - alpha*x**3

Now, set up the initial conditions and solve the equation of motion:
>>> from scipy.integrate import solve_ivp
>>> import numpy as np
>>> def equation_of_motion(t, y):
...     return [y[1], -beta*y[0] - alpha*y[0]**3]
>>> t_span = (0, 10)  # time from 0 to 10 seconds
>>> y0 = [1, 0]  # initial conditions: displacement = 1, velocity = 0
>>> t_eval = np.linspace(*t_span, 300)  # time points at which to store the results
>>> solution = solve_ivp(equation_of_motion, t_span, y0, t_eval=t_eval)

I would be glad to hear any feedback or suggestions to improve this simulation!

@tjstienstra

[moorepants]

This looks nice. My suggestion is that we add the example in the documentation. We can use numpy, scipy, and matplotlib in the docs. If we create an explanatory page about actuators in the prose documentation, then you can place it there. Whether you put it "tutorials", "how-to guides", or "explanations" depends on how you present the material. I don't really understand the difference in the three, but maybe you can make a judgement. Eventually we need to move all of the prose documentation in the physics.vector and physics.mechanics sections into one of these top level sections, so you can be the first to establish where these things go. I'm guessing this is a little long of an example to belong in the docstring. @tjstienstra thoughts?

[tjstienstra]

This is indeed a really nice example, where you show that we can reproduce the results with a simulation. I agree that it would be best not to put it in the docstring, but in another location. Previously, I would have put examples in physics/mechanics/examples, but I think that "tutorials" is probably the place we should start moving our examples to.

[eh111eh]

Thanks for the feedback! How about creating two directories: sympy.doc.src.tutorials.mechanics to add the Duffing spring simulation and move other examples for physics.mechanics, and sympy.doc.src.tutorials.vector to move examples for physics.vector?

[moorepants]

In the tutorials directory: https://github.com/sympy/sympy/tree/master/doc/src/tutorials we should probably follow the directory structure of the sympy subpackages. For example there is a sympy/vector and sympy/physics/vector, so for long term organization, we should have docs/src/tutorials/physics/mechanics (I see that we didn't follow this with biomechanics but we can correct that, especially since sympy 1.13 is not yet released.

[moorepants]

If you go about moving any examples from the old location, then you will need to set up link redirects. We have a sphinx plugin that reasonably automates that.