# Instability of perpetual manifolds of coupled oscillators, problem relevant to numerical integration algorithms

**URL:** <https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [July 31, 2023, 6:15pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328 "2023-07-31T18:15:25Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Prof\_FG](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prof_fg/32/51969_2.png) [@Prof\_FG](https://discourse.julialang.org/u/Prof_FG)\
**Post date:** [July 31, 2023, 6:15pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/1 "2023-07-31T18:15:26Z")

</div>

Dear all,  
I downloaded Julia to try solving numerically a system of nonlinear  
differential equations, which is based on my article about perpetual  
points/manifolds, in Journal of Mechanical Science [https://doi.org/10.1177/0954406220934833](https://doi.org/10.1177/0954406220934833).  
In the 2 degrees of freedom nonlinear example there are three Perpetual  
Manifolds, restricting the study to the simple one, which is the 2nd,  
the numerical simulation due to instability cannot capture rigid body  
motion with displacements described by two straight lines (same slope)  
with constant space difference, and of course constant equal velocities  
(the perpetual manifolds are defined by zero accelerations).  
The third example is even more difficult to obtain the numerical solution,  
and I skip it for the time being, I guess dealing with the 2nd example  
the third one might be solved too.  
Regarding Julia solutions, I tried  
to solve the problem with high accuracy (attached script and the plot)  
but didn’t work properly. I have to highlight two more things, a) using  
semi analytically the finite difference written in a paper the solution  
of straight lines can be easily verified, b) none of the software could  
capture totally this solution Matlab, Mathematica, Maple, Scilab,  
Octave.  
My motivation for this work, I am working for the development  
of the Perpetual Mechanics Theory, targeting a wide audience too, that  
is based on these solutions and without commercial software verifying  
them my work in this particular direction is quite uneasy…  
Please I would appreciate it if it is possible to confirm my solution and, feel  
free to contact me for any relevant question, and I would be glad to  
explain the problem in mathematical terms too, accompanied with its  
significance in mechanics.  
Thank you in advance for the time and effort spent on this matter.  
with my best regards and wishes  
Prof. Fotios Georgiades/Georgiadis

Julia code

using DifferentialEquations  
m1=1000.00;  
m2=1000.00;  
k1=1.00e6;  
k2=-5.00e5;  
c=427.21;  
y0=1.00;  
x0=y0+(-k1/k2)^0.5;  
dx0=1.00;  
dy0=1.00;  
params=(m1, m2, k1, k2, c);  
function dof2(dx, x, params, t)  
m1, m2, k1, k2, c=params  
dx[1] = x[3]  
dx[2] = x[4]  
dx[3] = (-k1\*(x[1]-x[2])-k2\*(x[1]-x[2])^3-c\*(x[3]-x[4]))/m1  
dx[4] = (k1\*(x[1]-x[2])+k2\*(x[1]-x[2])^3+c\*(x[3]-x[4]))/m2  
end  
xtot0=[x0, y0, dx0, dy0]  
tspn=[0.00, 1.00]  
prb=ODEProblem(dof2, xtot0, tspn, params)  
sln = solve(prb, abstol = 1e-16, reltol = 1e-16)  
using Plots  
plot(sln)  
savefig(“plotall.svg”)

 ![plots_all](https://global.discourse-cdn.com/julialang/original/3X/e/5/e569aa5a4a7bf27690feda64c2aa3d7dcf0340fe.jpeg)

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [July 31, 2023, 7:42pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/2 "2023-07-31T19:42:02Z")

</div>

Based on the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems) I tried `RadauIIA5`

```julia
sln = solve(prb, RadauIIA5(), abstol = 1e-16, reltol = 1e-16);

```

And got this solution  
 ![plotall](https://global.discourse-cdn.com/julialang/original/3X/c/7/c7fa524f95889f63adde7e12108ddf103629261f.png)

It is probably also worth trying Arbitrary precision floats as well.

---

<div class="post-metadata">

**Author:** ![Prof\_FG](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prof_fg/32/51969_2.png) [@Prof\_FG](https://discourse.julialang.org/u/Prof_FG)\
**Post date:** [July 31, 2023, 8:19pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/3 "2023-07-31T20:19:13Z")

</div>

Please check by changing the limits e.g. up to 10 or 20. In another software use of three different limits led to three different solutions with one of them being correct.

Στις Δευ, 31 Ιουλ 2023, 10:47 μ.μ. ο χρήστης Russel Howe via Julia Programming Language \<[notifications@julialang.discoursemail.com](mailto:notifications@julialang.discoursemail.com)\> έγραψε:

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 1, 2023, 3:30am UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/4 "2023-08-01T03:30:08Z")

</div>

Just using a solver for stiff systems seems fine here? There’s nothing else to it from what I can see.

---

<div class="post-metadata">

**Author:** ![Prof\_FG](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prof_fg/32/51969_2.png) [@Prof\_FG](https://discourse.julialang.org/u/Prof_FG)\
**Post date:** [August 1, 2023, 8:53am UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/5 "2023-08-01T08:53:54Z")

</div>

Dears Russel and Chris,

First of all thanks for your time and effort in finding a solution. Regarding the solution itself rather obvious that is a well solved problem. Regarding the analysis even treating the problem as a stiff one in other software (as I mentioned) the solution couldn’t be captured. I did a search for the Randau 17 and I found an article,  
[A 17th-order Radau IIA method for package RADAU. Applications in mechanical systems - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S0898122109007925)  
that even in the introduction (attached with yellow highlighting of the points) mentioning that there is a certain class of implicit Runge-Kutta algorithms (RADAU, STRIKE, DIRK, SDIRK) that are not only dealing with stiff problems but they are rather stable and this is the case in here.

 ![IMPLICIT_RUNGE_KUTTA](https://global.discourse-cdn.com/julialang/original/3X/d/f/df61c0e43d9315f8c8ec9ba1419b1510f3951ce7.jpeg)

Thanks for all, with my best wishes and regards

Fotis

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 1, 2023, 2:21pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/6 "2023-08-01T14:21:07Z")

</div>

> [@Prof\_FG](#):
>
> that even in the introduction (attached with yellow highlighting of the points) mentioning that there is a certain class of implicit Runge-Kutta algorithms (RADAU, STRIKE, DIRK, SDIRK) that are not only dealing with stiff problems but they are rather stable and this is the case in here.

Yes, these are in DifferentialEquations.jl. Radau5 is just one choice. See [ODE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)

---

<div class="post-metadata">

**Author:** ![Prof\_FG](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prof_fg/32/51969_2.png) [@Prof\_FG](https://discourse.julialang.org/u/Prof_FG)\
**Post date:** [August 1, 2023, 7:10pm UTC](https://discourse.julialang.org/t/instability-of-perpetual-manifolds-of-coupled-oscillators-problem-relevant-to-numerical-integration-algorithms/102328/7 "2023-08-01T19:10:24Z")

</div>

Thanks a lot for all.

Fotis