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.

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”)