juliadifferentialequations.jl

Error when attempting to solve Ordinary Differential equation using Symplectic solver in Julia

using LinearAlgebra 
using PyPlot
using Random
using DifferentialEquations

 N = 100
 σ, a, J, K = 10.0, 1.0, 1.0,-0.1;
 ω = zeros(N);  

 function heaviside(x::Real)
 return x >= 0.0 ? 1.0 : 0.0
 end 

 function rhs(du,u,p,t)

  u1 = @view u[1:N]
  du1 = @view du[1:N]
  u2 = @view u[N+1:2*N]
  du2 = @view du[N+1:2*N]
 
  σ,a,J,K=p

  for i in 1:N
    du1[i]=(1/N)*sum(j->( (u1[j]-u1[i])*(1+J*cos(u2[j]-u2[i]))-sign(u1[j]-u1[i]) ),1:N)
    
    du2[i]=ω[i]+ (K/N)* sum(j->( sin(u2[j]-u2[i])*(1-(u1[j]-u1[i])^2/σ^2)*heaviside(σ - 
 abs((u1[j]-u1[i])))),1:N)
end

return du 

end

tspan = (0.0, 200)  # Assuming you want to integrate from 0 to T
dts=0.1 #savetime
p=[σ,a,J,K]
Random.seed!(123)
u0= [(rand() * 20.0 - 10.0, rand() * (2.0 * π) - π) for j in 1:N]
#println(u0)
u0=vcat([x[1] for x in u0], [x[2] for x in u0]);

 prob = ODEProblem(rhs, u0, tspan,p);

 @time sol = solve(prob,KahanLi8(),dt=0.1,saveat=dts);

While I solve this ODEproblem using Tsit5() solver, I am able to integrate it and it gives me result. However, when I try to solve it with Symplectic Solver Like KahanLi8(), the code shows errors. From the error it seems that the problem is maybe in defining the function which is not compatible with the KahanLi8() solver, but I am unable to figure it out in any way. Any help will be appreciated.

The error is as follows

type Array has no field x

Solution

  • The symplectic integrators that you are trying require a partitioned ODE, as documented here:

    https://docs.sciml.ai/DiffEqDocs/stable/solvers/dynamical_solve/

    If you want to solve a non-partitioned ODE with these methods, you need to use a different solver that is designed to be symplectic for non-split ODEs. One such solver in the SciML ecosystem is https://github.com/SciML/IRKGaussLegendre.jl. Thus you could try changing

    @time sol = solve(prob,KahanLi8(),dt=0.1,saveat=dts);

    to:

    using IRKGaussLegendre
@time sol = solve(prob,IRKGL16(),dt=0.1,saveat=dts);