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

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
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,
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
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
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

Based on the documentation I tried RadauIIA5

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

And got this solution

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

1 Like

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 <> έγραψε:

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

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
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.

Thanks for all, with my best wishes and regards


Yes, these are in DifferentialEquations.jl. Radau5 is just one choice. See ODE Solvers · DifferentialEquations.jl

Thanks a lot for all.