Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement Duffing Spring Actuator.py #26438

Merged
merged 12 commits into from
Apr 23, 2024

Conversation

eh111eh
Copy link
Contributor

@eh111eh eh111eh commented Apr 1, 2024

References to other Issues or PRs

#26402

Brief description of what is fixed or changed

The Duffing equation is a nonlinear second order differential equation used to model the dynamics of systems, damped and driven oscillators, with a cubic nonlinearity in the restoring force, which is given by $$\ddot{x} + \delta\dot{x} + \beta x + \alpha x^3 = \gamma \cos(\omega t)$$

However, as a full Duffing second order forced system will not work as an actuator, I modify the Duffing Spring to treat it as a non-second order system.
For $\beta>0$, the Duffing oscillator can be interpreted as a forced oscillator with a spring whose restoring force is written as $$F = -\beta x - \alpha x^3$$
Since the case $\beta<0$ causes the chaotic motions, I will not consider this case for now.

I Implemented a new class DuffingSpring to model the force exerted by a nonlinear spring, based on the Duffing equation: $F = -\beta x - \alpha x^3$
where $x$ is the displacement from the equilibrium position, $\beta$ is the linear spring constant, and $\alpha$ is the coefficient for the nonlinear cubic term.

Task

Other comments

Reference

Release Notes

  • physics.mechanics
    • implemented DuffingSpring to model the force exerted by a nonlinear spring, based on the Duffing equation

@eh111eh eh111eh requested a review from moorepants as a code owner April 1, 2024 13:18
@sympy-bot
Copy link

sympy-bot commented Apr 1, 2024

Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • physics.mechanics
    • implemented DuffingSpring to model the force exerted by a nonlinear spring, based on the Duffing equation (#26438 by @eh111eh)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13.

Click here to see the pull request description that was parsed.
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
#26402 

#### Brief description of what is fixed or changed
The Duffing equation is a nonlinear second order differential equation used to model the dynamics of systems, damped and driven oscillators, with a cubic nonlinearity in the restoring force, which is given by $$\\ddot{x} + \delta\dot{x} + \beta x + \alpha x^3 = \gamma \cos(\omega t)\$$

However, as a full Duffing second order forced system will not work as an actuator, I modify the Duffing Spring to treat it as a non-second order system.
For $\beta>0$, the Duffing oscillator can be interpreted as a forced oscillator with a spring whose restoring force is written as $$F = -\beta x - \alpha x^3$$
Since the case  $\beta<0$ causes the chaotic motions, I will not consider this case for now. 

I Implemented a new class `DuffingSpring` to model the force exerted by a nonlinear spring, based on the Duffing equation: $F = -\beta x - \alpha x^3$
where $x$ is the displacement from the equilibrium position, $\beta$ is the linear spring constant, and $\alpha$ is the coefficient for the nonlinear cubic term.

#### Task
- [ ] Add a unit test to verify the new class `DuffingSpring` works as expected. Refer to the existing unit tests for actuators as examples.
- [ ] Demonstrate the force exerted by a nonlinear spring, based on the Duffing equation, as illustrated in the following example: https://pydy.readthedocs.io/en/latest/examples/mass-spring-damper.html

#### Other comments
Reference
- http://www.scholarpedia.org/article/Duffing_oscillator 
- https://www.cfm.brown.edu/people/dobrush/am34/Matlab/ch3/duffing.html
- https://www.cfm.brown.edu/people/dobrush/am34/Mathematica/ch3/duffing.html

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* physics.mechanics
  * implemented `DuffingSpring` to model the force exerted by a nonlinear spring, based on the Duffing equation
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

There are several issues with this test for the DuffingSpring class including
i. Type checking (attributes, etc),
ii. Force expression
iii. SimpifyError in test_to_loads
iv. Failed to raise error
I've addressed these issues to pass pytest
The following three assertion errors in test_valid_constructor have been addressed:
1. Resolved Issue: The initial assertion error was resolved by ensuring spring.linear_stiffness is now correctly instantiated as an ExprType. This was achieved by properly setting and converting the linear_stiffness attribute into a SymPy expression.
2. Resolved Issue: The second assertion error related to the equilibrium_length attribute mirrored the first and was similarly addressed.
3. Resolved Issue: The third assertion error, concerning the mismatch of the force attribute with its expected expression, was corrected by revising the force calculation implementation in the DuffingSpring class.
The test error associated with test_properties_are_immutable has been resolved. There remains one last error to address: This error arises in the test_repr method of the TestDuffingSpring class during pytest collection. The problem is tied to the use of the equilibrium_length parameter, which, although specified in the pytest.mark.parametrize decorator, is missing from the method signature of test_repr. This mismatch results in the error message: "In test_repr: function uses no argument 'equilibrium_length'."
changes in actuator.py:
1. The DuffingSpring class was modified to use property decorators, which provide a getter method for linear_stiffness, nonlinear_stiffness, pathway, and equilibrium_length. This change encapsulates the properties and prevents them from being modified after the object is created, enforcing immutability of the spring parameters.

2. The constructor (__init__) uses a try-except block to catch SympifyError, ensuring that inputs can be converted to SymPy expressions. If conversion fails, it raises an error informing the user.

changes in test_actuator.py
1. The pytest.mark.parametrize decorator is used to run the same test function with different sets of arguments, allowing comprehensive testing of the DuffingSpring constructor with various stiffness and equilibrium length parameters.

2. Tests were added to confirm that the DuffingSpring object initializes with the expected attributes and that these attributes match the parameters passed to the constructor.

3. Tests ensure that once a DuffingSpring object is created, its key properties cannot be altered, testing the enforcement of immutability.

4. The test_to_loads function calculates the expected forces based on the Duffing equation and compares them to the forces returned by the to_loads method of the DuffingSpring object.

5. The __repr__ method is tested to confirm it returns the expected string representation of the DuffingSpring instance.

Now it passes pytest
@eh111eh
Copy link
Contributor Author

eh111eh commented Apr 13, 2024

Hi, I've added a test unit for the DuffingSpring class in test_actuator.py and made minor modifications to the DuffingSpring class itself. The test unit covers the following functionalities:

  • Fixture Setup
  • Inheritance Checks
  • Parametrized Tests
  • Constructor Validity
  • Invalid Constructor Arguments
  • Immutability of Properties
  • String Representation
  • Force Calculation

All tests pass in pytest. I would greatly appreciate any feedback or comments on this work. Thank you!

Copy link
Contributor

@tjstienstra tjstienstra left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not familiar with the Duffing spring, but the implementation looks mostly good to me. Just suggested a few minor changes.

sympy/physics/mechanics/actuator.py Outdated Show resolved Hide resolved
sympy/physics/mechanics/actuator.py Outdated Show resolved Hide resolved
sympy/physics/mechanics/actuator.py Outdated Show resolved Hide resolved
sympy/physics/mechanics/actuator.py Outdated Show resolved Hide resolved
sympy/physics/mechanics/tests/test_actuator.py Outdated Show resolved Hide resolved
@tjstienstra
Copy link
Contributor

@eh111eh don't forget to add yourself to the mailmap.

@eh111eh
Copy link
Contributor Author

eh111eh commented Apr 16, 2024

I am not familiar with the Duffing spring, but the implementation looks mostly good to me. Just suggested a few minor changes.

Thank you for your thorough feedback. I have made the following changes:

  • Added a blank line before and after Explanation and Parameters.
  • Removed the try statement.
  • Modified the setters to match the implementation in other classes, such as LinearSpring.
  • Removed simplify from the displacement calculation. I originally applied simplify to resolve an issue where the actual force output did not match the expected force output and forgot to revert it after resolving it.
  • Avoided repeating tests in test_actuator.py and added tests for None.
  • Added my name and email to the mailmap.

I would appreciate your comments and feedback. Thank you!

Copy link

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

| Change   | Before [2487dbb5]    | After [b9bafe06]    |   Ratio | Benchmark (Parameter)                                                |
|----------|----------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 68.3±1ms             | 44.3±0.5ms          |    0.65 | integrate.TimeIntegrationRisch02.time_doit(10)                       |
| -        | 68.4±1ms             | 43.1±0.2ms          |    0.63 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 18.6±0.2μs           | 30.2±0.2μs          |    1.63 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 5.31±0.02ms          | 2.82±0.01ms         |    0.53 | logic.LogicSuite.time_load_file                                      |
| -        | 72.9±0.3ms           | 29.0±0.1ms          |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'dense')                           |
| -        | 25.9±0.1ms           | 16.9±0.06ms         |    0.65 | polys.TimeGCD_GaussInt.time_op(1, 'expr')                            |
| -        | 73.6±0.4ms           | 29.3±0.2ms          |    0.4  | polys.TimeGCD_GaussInt.time_op(1, 'sparse')                          |
| -        | 258±2ms              | 126±0.7ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'dense')                           |
| -        | 258±2ms              | 126±0.6ms           |    0.49 | polys.TimeGCD_GaussInt.time_op(2, 'sparse')                          |
| -        | 652±3ms              | 384±1ms             |    0.59 | polys.TimeGCD_GaussInt.time_op(3, 'dense')                           |
| -        | 654±5ms              | 376±2ms             |    0.57 | polys.TimeGCD_GaussInt.time_op(3, 'sparse')                          |
| -        | 502±2μs              | 290±3μs             |    0.58 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense')            |
| -        | 1.81±0.01ms          | 1.05±0.01ms         |    0.58 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense')            |
| -        | 5.85±0.02ms          | 3.11±0.01ms         |    0.53 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense')            |
| -        | 446±2μs              | 234±3μs             |    0.53 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense')               |
| -        | 1.50±0.01ms          | 697±6μs             |    0.46 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense')               |
| -        | 4.92±0.02ms          | 1.65±0.01ms         |    0.34 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense')               |
| -        | 378±2μs              | 206±1μs             |    0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense')                |
| -        | 2.44±0.01ms          | 1.23±0.01ms         |    0.5  | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense')                |
| -        | 10.1±0.05ms          | 4.36±0.04ms         |    0.43 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense')                |
| -        | 361±3μs              | 173±1μs             |    0.48 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense')            |
| -        | 2.52±0.01ms          | 901±9μs             |    0.36 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense')            |
| -        | 9.59±0.02ms          | 2.64±0.02ms         |    0.28 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense')            |
| -        | 1.05±0.01ms          | 428±2μs             |    0.41 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense')           |
| -        | 1.73±0ms             | 499±2μs             |    0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 5.93±0.05ms          | 1.81±0.02ms         |    0.31 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense')           |
| -        | 8.56±0.04ms          | 1.49±0ms            |    0.17 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 286±0.9μs            | 66.8±0.5μs          |    0.23 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 3.44±0.02ms          | 397±2μs             |    0.12 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense')              |
| -        | 3.97±0.03ms          | 279±1μs             |    0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 7.10±0.03ms          | 1.28±0.01ms         |    0.18 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense')              |
| -        | 8.75±0.02ms          | 834±2μs             |    0.1  | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 5.00±0.02ms          | 2.95±0.02ms         |    0.59 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 12.2±0.06ms          | 6.61±0.04ms         |    0.54 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense')  |
| -        | 22.4±0.06ms          | 9.01±0.02ms         |    0.4  | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 5.21±0.02ms          | 859±4μs             |    0.16 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 12.6±0.04ms          | 7.03±0.02ms         |    0.56 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 103±0.6ms            | 25.8±0.2ms          |    0.25 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense')     |
| -        | 166±0.9ms            | 54.0±0.1ms          |    0.32 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 171±0.8μs            | 111±0.6μs           |    0.65 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense')      |
| -        | 358±0.7μs            | 218±1μs             |    0.61 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse')     |
| -        | 4.34±0.02ms          | 846±6μs             |    0.19 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense')      |
| -        | 5.29±0.06ms          | 383±1μs             |    0.07 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 19.9±0.1ms           | 2.83±0.02ms         |    0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense')      |
| -        | 22.8±0.1ms           | 626±5μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 480±2μs              | 140±0.5μs           |    0.29 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 4.75±0.01ms          | 621±6μs             |    0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense')  |
| -        | 5.33±0.03ms          | 141±1μs             |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 13.3±0.2ms           | 1.31±0.01ms         |    0.1  | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense')  |
| -        | 13.8±0.1ms           | 144±0.8μs           |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 133±1μs              | 74.2±1μs            |    0.56 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 247±1μs              | 86.7±0.2μs          |    0.35 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 24.6±0.1ms           | 10.2±0.03ms         |    0.41 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |
| -        | 28.3±0.1ms           | 15.6±0.3ms          |    0.55 | solve.TimeSparseSystem.time_linsolve_Aaug(20)                        |
| -        | 54.6±0.3ms           | 24.8±0.09ms         |    0.45 | solve.TimeSparseSystem.time_linsolve_Aaug(30)                        |
| -        | 28.2±0.2ms           | 15.2±0.1ms          |    0.54 | solve.TimeSparseSystem.time_linsolve_Ab(20)                          |
| -        | 54.6±0.3ms           | 24.5±0.07ms         |    0.45 | solve.TimeSparseSystem.time_linsolve_Ab(30)                          |

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

Copy link
Contributor

@tjstienstra tjstienstra left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM

@moorepants moorepants merged commit 3b6df5e into sympy:master Apr 23, 2024
48 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants