 def euler_cromer(state, rhs, dt):
   
     mid_state = state + rhs(state)*dt 
     mid_derivs = rhs(mid_state)       
     
     next_state = np.array([mid_state[0], state[1] + mid_derivs[1]*dt])
     
     return next_state

模拟数据

 w = 2
 period = 2*np.pi/w
 dt = period/200  
 T = 800*period  
 N = round(T/dt)
 
 print('The number of time steps is {}.'.format( N ))
 print('The time increment is {}'.format( dt ))
 
 t = np.linspace(0, T, N)
 
 x0 = 2    
 v0 = 0    
 
 num_sol = np.zeros([N,2])
 num_sol[0,0] = x0
 num_sol[0,1] = v0
 
 for i in range(N-1):
     num_sol[i+1] = euler_cromer(num_sol[i], springmass, dt)

The number of time steps is 160000. The time increment is 0.015707963267948967

首先，得到解析解。然后，您选择绘制振荡运动的前几个周期：数值和解析。

 x_an = x0*np.cos(w * t)

 iend = 800 
 fig = plt.figure(figsize=(6,4))
 plt.plot(t[:iend], num_sol[:iend, 0], linewidth=2, linestyle='--', label='Numerical solution')
 plt.plot(t[:iend], x_an[:iend], linewidth=1, linestyle='-', label='Analytical solution')
 plt.xlabel('Time [s]')
 plt.ylabel('$x$ [m]')
 plt.title('Spring-mass system, with Euler-Cromer method.\n');

该图显示，欧拉-克罗默不存在振幅增大的问题。从这个意义上讲，你应该对此感到满意。但是，如果你绘制一段较长模拟的末尾，你就会发现它确实开始偏离解析解。

 istart = 400
 
 fig = plt.figure(figsize=(6,4))
 
 plt.plot(t[-istart:], num_sol[-istart:, 0], linewidth=2, linestyle='--', label='Numerical solution')
 plt.plot(t[-istart:], x_an[-istart:], linewidth=1, linestyle='-', label='Analytical solution')
 plt.xlabel('Time [s]')
 plt.ylabel('$x$ [m]')
 plt.title('Spring-mass system, with Euler-Cromer method. \n');

观察一段很长的运行中的最后几次振荡，即使时间增量很小，也会发现轻微的相位差。因此，尽管欧拉-克罗默方法解决了欧拉方法的一个大问题，但它仍然存在一些错误。它仍然是一阶方法！

PreviousC++(Python)肥皂泡沫普拉托边界膜曲面模型算法 NextPython湍流隐式模型耗散粘性方程和大涡流模拟

Last updated 4 months ago

def euler_cromer(state, rhs, dt): mid_state = state + rhs(state)*dt mid_derivs = rhs(mid_state) next_state = np.array([mid_state[0], state[1] + mid_derivs[1]*dt]) return next_state

w = 2 period = 2*np.pi/w dt = period/200 T = 800*period N = round(T/dt) print('The number of time steps is {}.'.format( N )) print('The time increment is {}'.format( dt )) t = np.linspace(0, T, N) x0 = 2 v0 = 0 num_sol = np.zeros([N,2]) num_sol[0,0] = x0 num_sol[0,1] = v0 for i in range(N-1): num_sol[i+1] = euler_cromer(num_sol[i], springmass, dt)

iend = 800 fig = plt.figure(figsize=(6,4)) plt.plot(t[:iend], num_sol[:iend, 0], linewidth=2, linestyle='--', label='Numerical solution') plt.plot(t[:iend], x_an[:iend], linewidth=1, linestyle='-', label='Analytical solution') plt.xlabel('Time [s]') plt.ylabel('$x$ [m]') plt.title('Spring-mass system, with Euler-Cromer method.\n');

istart = 400 fig = plt.figure(figsize=(6,4)) plt.plot(t[-istart:], num_sol[-istart:, 0], linewidth=2, linestyle='--', label='Numerical solution') plt.plot(t[-istart:], x_an[-istart:], linewidth=1, linestyle='-', label='Analytical solution') plt.xlabel('Time [s]') plt.ylabel('$x$ [m]') plt.title('Spring-mass system, with Euler-Cromer method. \n');