DiffEq result for a dynamical system differs from wolframalpha

I want to solve this differential equation:

(d^2 y(t))/(dt^2) = -0.1×(dy(t))/(dt) + 0.5 y(t) - 2 y(t)^3 + 2 cos(2.4 t)

How I tried to solve it using DiffEq:

using DifferentialEquations
using Plots

function duffing!(du,u,p,t)
    du[1]=u[2]
    du[2]=-p[1]*u[2] + 2*p[2] * u[1] - 4*p[3]*(u[1])^3 + p[4] * cos(p[5]*t)
end

function Duffing(x0,xx0,F0,ω,γ,max_k=5000,m=1,a=0.25,b=0.5)
    u0=[x0,xx0]
    tspan=(0.0,max_k*2*pi/ω)
    p=[γ,a,b,F0,ω]
    prob=ODEProblem(duffing!,u0,tspan,p)
    return solve(prob,DPRKN6(),saveat=pi/ω)
end

function plotter(x0,xx0,F0,ω,γ,max_k)
    closeall()
    sol=Duffing(x0,xx0,F0,ω,γ,max_k)
    x_arr=[];xx_arr=[];
    period=2*pi/ω
    for k=0:max_k
        pair=sol(k*period)
        x_arr=push!(x_arr,pair[1])
        xx_arr=push!(xx_arr,pair[2])
    end
    scatter(x_arr,xx_arr,markersize=0.8,label="")
    #savefig("img.png")
end



However, what I get looks like a big clump rather than what wolframalpha gives out. It doesn’t look like my system is chaotic but rather goes towards a fixed point.

EDIT:What the Book says I should get for the same conditions as the figure above

1 Like

We have an example of the stroboscopic map of the duffing system in the DynamicalSystems.jl documentation, and it shows the chaotic attractor without a problem: https://juliadynamics.github.io/DynamicalSystems.jl/dev/chaos/orbitdiagram/#Stroboscopic-Map-1

And the way we define the dynamical rule is as follows:

(trajectory is just a wrapper around solve that explicitly gives saveat)

my equations are slightly different than yours, because I use the out of place form (recommended for small systems). I also use slightly different parameters, but that shouldn’t matter. Perhaps check your code for typos? Also, does DPRKN6() allow interpolation? (since the syntaxx sol(some_time) interpolates)

1 Like

Yes

1 Like

It appears the problem is solved if I change the constants. Perhaps the book is wrong (Thijssen Comp. Phys.).

2 Likes