- #1
Oliver321
- 59
- 5
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!
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!
Python code:
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()
Attachments
Last edited by a moderator: