# 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=u
du=-p*u + 2*p * u - 4*p*(u)^3 + p * cos(p*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)
xx_arr=push!(xx_arr,pair)
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.

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