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