BoundsError in FOChaos.jl

[finmod]

Trying t run the series of examples of Chaos,

result = solve(prob, h, GL()) generates BoundsError in every example.

[ErikQQY]

Hello there! Sorry for the errors, could you provide your FOChaos.jl example and the version of FractionalDiffEq.jl you are using?

[finmod]

I managed to solve the problem with the following code for all examples in FOChaos.jl with FractionalDiffEq.jl v0.3.0

# Fractional Order Chaotic examples come from 
#=
@inproceedings{Petr2011FractionalOrderNS,
title={Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation},
author={Ivo Petr{\'a}{\vs}},
year={2011}
}
=#
#

#=======Fractional Order Liu System=======#
using FractionalDiffEq, Plots

function Liu!(du, u, p, t)
    a, b, c, e, n, m = p
    du[1] = -a*u[1]-e*u[2]^2
    du[2] = b*u[2]-n*u[1]*u[3]
    du[3] = -c*u[3]+m*u[1]*u[2]
end

α = [0.95, 0.95, 0.95]
x0 = [0.2, 0, 0.5]
h = 0.01; tspan = (0, 100);
p = [1.0, 2.5, 5.0, 1.0, 4, 4]
prob = FODESystem(Liu!, α, x0, tspan, p)
result = solve(prob, h, NonLinearAlg())


plot(result, vars=(1,2,3), title="Fractional Order Liu System")

This only allows me to test NonLinearAlg() and not GL() or others that would be the interesting part of these examples.

[ErikQQY]

It is very stange🤔, I tested Liu example on my computer and it works fine.

[finmod]

You are right. It runs. It is just that GL() is 50 x slower than the other two solvers.
Still, on running the complete set, I have a bug in the Van der Pol example. Can you find it?

#=======Fractional Order Van der Pol Oscillator=======#
using FractionalDiffEq, Plots

function VanderPol!(du, u, p, t)
    ϵ = p
    du[1] = u[2]
    du[2] = -u[1]-ϵ*(u[1]^2-1)*u[2]
end

α = [1.2, 0.8]
x0 = [0.2, -0.2]
h = 0.01; tspan = (0, 60);
p = [1.0]
prob = FODESystem(VanderPol!, α, x0, tspan, p)
result = solve(prob, h, NonLinearAlg())
#result = solve(prob, h, GL())
#result = solve(prob, h, PIEX())

#plot(result, vars=(1,2), title="Fractional Order Van der Pol Oscillator")

[ErikQQY]

Indeed, the GL solver is much solver than the other solvers.

There is an MethodError here, this error is caused by p=[1.0], since there is ϵ = p in the definition of VanderPol!, we should set p=1.0 instead p=[1.0] to avoid this error.

[finmod]

Spot on! Now, the PIEX() solver gives a BoundsError: attempt to access 2-element Vector {Float 64} at index [[1],2]

#=======Fractional Order Van der Pol Oscillator=======#
using FractionalDiffEq, Plots

function VanderPol!(du, u, p, t)
    ϵ = p
    du[1] = u[2]
    du[2] = -u[1]-ϵ*(u[1]^2-1)*u[2]
end

α = [1.2, 0.8]
x0 = [0.2, -0.2]
h = 0.01; tspan = (0, 60);
p = 1.0
prob = FODESystem(VanderPol!, α, x0, tspan, p)
#result = solve(prob, h, NonLinearAlg())
#result = solve(prob, h, GL())
result = solve(prob, h, PIEX())

plot(result, vars=(1,2), title="Fractional Order Van der Pol Oscillator")

[ErikQQY]

Yes, it is a real bug here🤔, when we pass one parameter in the problem definition, it will cause BoundsError. Sorry for the errors, I will find some time to solve them. You can try to pass an array of parameters to the problem definition for now to avoid the bug.

[ErikQQY]

Ha! Just figure it out! I think the PIEX, PECE and FLMM solvers don't support high-order FDEs(here one of the orders in the problem is 1.2, bigger than 1), I think I should add that in the codebase

[finmod]

Correct. If the order 1.2 is changed to 0.9, it runs.
Thank you for the list of solvers where this happens: PIEX, PECE, FLMM. Any others?

[ErikQQY]

I think other solvers should work fine, I will immediately patch a release to solve this problem.

[ErikQQY]

Just patched a new release to address this problem: JuliaRegistries/General#80729

Thanks a lot for spotting this! @finmod