Problem in solving differential equation

[Oliver321]

Hello everyone!
I was studying chaotic systems and therefore made some computer simulations in python. I simulated the driven damped anhatmonic oscillator.

The problem I am facing is with solving the differential equation for t=0s-200s. I used numpy.linspace(0,200,timesteps) for generate a time scale. I did this for three different timesteps: 10^5,10^6 and 10^7. For every timestep i get different results, and they jump between one another. If I increase it further, it changes again (see pictures).

Now I was wondering what was going on here! I thought the accuracy of the solution should be independent of the timesteps? If not, I can conclude: As I increase the timesteps, the result gets more accurate and if I increase it to high, there are going to be numerical errors. But how could I tell which solution is the right one? Shouldn’t 10^6 steps be enough?

Thanks for every answer!

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig

F=63
xv0=[4.346,0]
                                   
                           
res=10000000                                
tmax=200                                
t=np.linspace(0,tmax,res)
t2=np.linspace(0,250,res)
w=2                                        

def ddp(xv,t):                            
    x=xv[0]      
    v=xv[1]        
    b=0.3
    x2=-x-x**3-b*v+F*np.cos(w*t)
    return v,x2
x=odeint(ddp,xv0,np.linspace(0,50*2*np.pi/w,1000))    
xv0=[x[-1,0],x[-1,1]]

x=odeint(ddp,xv0,t)    
x2=odeint(ddp,xv0,t2)                                
n=int((res/tmax*2*np.pi/w))                          
peak,_=sig.find_peaks(-1*x[:,0],4)
x_s=np.take(x[:,0],peak)
v_s=np.take(x[:,1],peak)
t_s=np.take(t,peak)

fig, axs = plt.subplots(ncols=2,figsize=(12,4.8))
axs[0].plot(x2[:,0],x2[:,1],linewidth=0.5)
axs[0].plot(x_s,v_s,'.r',markersize='1')
axs[0].set(xlabel="x",ylabel="v")
axs[1].plot(x[:,0],x[:,1],linewidth=0.5)
axs[1].plot(x_s,v_s,'.r',markersize='3')
axs[1].set_xlim([-6.5,-4])
axs[1].set(xlabel="x",ylabel="v")
plt.ylabel('v')
#plt.savefig('phase_plot',dpi=300)
plt.show()

 

[S.G. Janssens]

Now I was wondering what was going on here! I thought the accuracy of the solution should be independent of the timesteps? If not, I can conclude: As I increase the timesteps, the result gets more accurate and if I increase it to high, there are going to be numerical errors. But how could I tell which solution is the right one? Shouldn’t 10^6 steps be enough?

Long time integrations of chaotic systems are very difficult, precisely because they are chaotic: For few time steps (large step size), the algorithm is usually not accurate enough. Increasing the number of time steps leads to more floating point calculations, so more round-off errors, and the effect of these is amplified by the chaotic nature of the system ("sensitive dependence on initial conditions").

For certain systems, you can use specialized integrations, e.g. symplectic integrators for chaotic Hamiltonian systems, but your system is dissipative. You should in any case compare different integrators, vary the step sizes (as you already did), and see if you can discover consistent behavior among integrators, for step sizes that are somewhere between "not too lange" and "not too small".

 

[Oliver321]

Thank you very much, this is helping me a lot! Unfortunately I am very inexperienced in this topic.
At one hand, the step size 10^6 would be optimal from a theoretical standpoint (the Lorenz section would (if plotted in the right way) represent a u-sequence). On the other hand result of stepsize 10^7 would be more constant over a long range of the stepsizes. But maybe this is to large.
Could you give me a clue what you think would be a optimal range?

Edit: I now tried it with different solvers, and it seems that the plot from stepsize 10^7 is most accurate!

 

[valenumr]

The optimal time step would be related to the time constant of you system. Systems with multiple time constants that have a large separation are called "rigid", iirc, and don't behave well with numerical integration. You can try rk4, or multi step methods such as adams-moulton or adams-bashforth.