title: Applying Modified Euler to solve a Pendulum ODE in Python
body: <p>So I am trying to implement some numerical methods into python and I am having some issues where all of my functions output more or less the same thing as the regular euler method. I assume this is because I am messing up in some way when I am implementing the method into code.</p>
  
  <p>My pendulum is defined as this:</p>
  
  <pre><code>def func(y,t):
      ### Simplified the Function to remove friction since it cancelled out
  
      x,v = y[:3],y[3:6]
      grav = np.array([0.,  0., -9.8 ])
      lambd = (grav.dot(x)+v.dot(v))/x.dot(x)
  
      return np.concatenate([v, grav - lambd*x] )
  
  def dF_matrix(y):
  
      n=y.size
      dF=np.zeros((6,6))
  
      xp=np.array([y[1],y[2],y[3]])[np.newaxis]
  
      mass=1.
      F1=0.
      F2=0.
      F3=-mass*9.8
      F=np.array([F1,F2,F3])[np.newaxis]
  
      phix=2.*y[0]
      phiy=2.*y[4]
      phiz=2.*y[5]
      G=np.array([phix,phiy,phiz])[np.newaxis]
  
      H=2.*np.eye(3)
      lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)
  
      dF[0,3]=1
      dF[1,4]=1
      dF[2,5]=1
      dF[3,0]=(y[0]*F1+2*lambd)/mass
      dF[3,1]=(y[0]*F2)/mass
      dF[3,2]=(y[0]*F3)/mass
      dF[3,3]=phix*y[1]
      dF[3,4]=phix*y[2]
      dF[3,5]=phix*y[3]
      dF[4,0]=(y[4]*F1)/mass
      dF[4,1]=(y[4]*F2+2*lambd)/mass
      dF[4,2]=(y[4]*F3)/mass
      dF[4,3]=phiy*y[1]
      dF[4,4]=phiy*y[2]
      dF[4,5]=phiy*y[3]
      dF[5,0]=(y[5]*F1)/mass
      dF[5,1]=(y[5]*F2)/mass
      dF[5,2]=(y[5]*F3+2*lambd)/mass
      dF[5,3]=phiz*y[1]
      dF[5,4]=phiz*y[2]
      dF[5,5]=phiz*y[3]
  
      return dF
  </code></pre>
  
  <p>The functions that I have made to integrate the ODE function are as follows (with help from others in previous a thread):</p>
  
  <pre><code>from scipy.integrate import odeint
  from mpl_toolkits import mplot3d
  import matplotlib.pyplot as plt
  </code></pre>
  
  <h3>Forward Euler Method</h3>
  
  <pre><code>def forward_Euler(function, y_matrix, time):
      y = np.zeros((np.size(time), np.size(y_matrix)))
      y[0, :] = y_matrix
  
      for i in range(len(time) - 1):
          dt = time[i + 1] - time[i]
          y[i + 1, :] = y[i, :] + np.asarray(function(y[i, :], time[i])) * dt
  
      return y
  
  </code></pre>
  
  <h3>Modified Euler Method</h3>
  
  <p>ERROR STARTS HERE </p>
  
  <p>The error I am getting is:
   RuntimeWarning: invalid value encountered in double_scalars
    lambd = (grav.dot(x)+v.dot(v))/x.dot(x)</p>
  
  <pre><code>def modified_Euler(function, y_matrix, time):
      y = np.zeros((np.size(time), np.size(y_matrix)))  # creates the matrix that we will fill
      y[0, :] = y_matrix # sets the initial values of the matrix
  
      for i in range(len(time) - 1):  # apply the Euler
          dt = time[i + 1] - time[i]  # Step size
          k1 = np.asarray(function(y[i, :], time[i])*dt)
          k2 = np.asarray(function(y[i] + k1, time[i+1])*dt)
          y[i + 1, :] = y[i, :] + .5 * (k1 + k2)
  
      return y
  </code></pre>
  
  <h3>Adams-Bashforth  2nd order</h3>
  
  <pre><code>def Adams_Bash_2nd(function, y_matrix, time):
      y = np.zeros((np.size(time), np.size(y_matrix)))
      y[0, :] = y_matrix
  
      dt = time[1] - time[0]
      f_0 = function(y[0], time[0])
      y[1] = y[0] + dt * f_0
      y[1] = y[0] + 0.5*dt * (f_0+function(y[1], time[1]))
  
      for i in range(len(time) - 1):
          dt = time[i + 1] - time[i]
          f_1 = function(y[i, :], time[i])
          f_2 = function(f_1-1, time[i-1])
          y[i + 1] = y[i] + 0.5 * dt * (3 * f_1 - f_2)
  
      return y
  </code></pre>
  
  <h3>Adams Bashforth Moulton Method</h3>
  
  <pre><code>
  def Adams_Moulton(function, y_matrix, time):
      y = np.zeros((np.size(time), np.size(y_matrix)))
      y[0, :] = y_matrix
  
  ### predictor formula
  
      for i in range(len(time) - 1):
          dt = time[i + 1] - time[i]
          f_1 = function(y[i, :], time[i])
          f_2 = function(f_1-1, time[i-1])
          y[i + 1, :] = y[i, :] + dt * f_1 + ((dt**2)/2) * f_2
  
  ### Corrector formula
  
      for i in range(len(time) - 1):
          dt = time[i + 1] - time[i]
          k_1 = 9 * (function(y[i, :], time[i+1]))
          k_2 = 19 * (function(y[i, :], time[i]))
          k_3 = 5 * (function(y[i, :], time[i-1]))
          k_4 = (function(y[i, :], time[i-2]))
          y[i + 1, :] = y[i] + (dt/24) * (k_1 + k_2 - k_3 + k_4)
  
      return y
  </code></pre>
  
  <h3>RK4 step to use in next function</h3>
  
  <pre><code>def RK4_step(f,y,t,dt, N=1):
      dt /= N;
      for k in range(N):
          k1=f(y,t)*dt; k2=f(y+k1/2,t+dt/2)*dt; k3=f(y+k2/2,t+dt/2)*dt; k4=f(y+k3,t+dt)*dt;
          y, t = y+(k1+2*(k2+k3)+k4)/6, t+dt
      return y
  
  </code></pre>
  
  <h3>Adams-Bashforth Moulton Method 4th order</h3>
  
  <pre><code>def Adams_Moulton_4th(function, y_matrix, time):
      y = np.zeros((np.size(time), np.size(y_matrix)))
      y[0] = y_matrix
          ### bootstrap steps with 4th order one-step method
      dt = time[4] - time[0]
      y[4] = RK4_step(function, y[0], time[0], dt, N=4)
      y[5] = RK4_step(function, y[4], time[4], dt, N=4)
      y[1] = RK4_step(function, y[5], time[5], dt, N=4)
  
      f_m2 = function(y[0], time[0])
      f_m1 = function(y[4], time[4])
      f_0 = function(y[5], time[5])
      f_1 = function(y[1], time[1])
      for i in range(3, len(time) - 1):
          ### predictor formula 4th order [ 55/24, -59/24, 37/24, -3/8 ]
          f_m3, f_m2, f_m1, f_0 = f_m2, f_m1, f_0, f_1
          y[i + 1] = y[i] + (dt / 24) * (55 * f_0 - 59 * f_m1 + 37 * f_m2 - 9 * f_m3)
          f_1 = function(y[i + 1], time[i + 1])
          ### Corrector formula 4th order [ 3/8, 19/24, -5/24, 1/24 ]
          y[i + 1] = y[i] + (dt / 24) * (9 * f_1 + 19 * f_0 - 5 * f_m1 + f_m2)
          f_1 = function(y[i + 1], time[i + 1])
      return y
  </code></pre>
  
  <p>I decided to program the way I am testing the functions into a with a function eliminating a good amount of lines from the previous iteration</p>
  
  <pre><code># initial condition
  y0 = np.array([0.0, 1.0, 0.0, 0.8, 0.0, 1.2])
  
  
  def test_function(test_function):
      print(test_function.__name__ + "...")
      nt = 2500
      time_start = process_time()
      # time points
      t = np.linspace(0, 25, nt)
      # solve ODE
      y1 = test_function(func, y0, t)
      time_elapsed = (process_time() - time_start)
      print('elapsed time', time_elapsed)
      # compute residual:
      r = y1[:, 0] ** 2 + y1[:, 1] ** 2 + y1[:, 2] ** 2 - 1
      rmax1 = np.max(np.abs(r))
      print('error', rmax1)
  
      fig = plt.figure()
      ax = plt.axes(projection='3d')
      ax.plot3D(y1[:, 0], y1[:, 1], y1[:, 2], 'gray')
  
      plt.show()
  
  
  test_function(odeint)
  test_function(forward_Euler)
  test_function(modified_Euler)
  test_function(Adams_Bash_2nd)
  test_function(Adams_Moulton)
  test_function(Adams_Moulton_4th)
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/62CHb.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/62CHb.png" alt="odeint outputs"></a>
  <a href="https://i.stack.imgur.com/qyHY0.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qyHY0.png" alt="forward_euler outputs"></a>
  <a href="https://i.stack.imgur.com/BLC47.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/BLC47.png" alt="modified_euler outputs"></a>
  <a href="https://i.stack.imgur.com/P94YL.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/P94YL.png" alt="Adams_Bash_2nd outputs"></a>
  <a href="https://i.stack.imgur.com/LFbZE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/LFbZE.png" alt="Adams_Moulton outputs"></a>
  <a href="https://i.stack.imgur.com/vgI5P.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vgI5P.png" alt="Adams_Moulton_4th outputs"></a></p>