Applying Modified Euler to solve a Pendulum ODE in Python

Question

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.

My pendulum is defined as this:

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

The functions that I have made to integrate the ODE function are as follows (with help from others in previous a thread):

from scipy.integrate import odeint
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt

Forward Euler Method

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

Modified Euler Method

ERROR STARTS HERE

The error I am getting is: RuntimeWarning: invalid value encountered in double_scalars lambd = (grav.dot(x)+v.dot(v))/x.dot(x)

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

Adams-Bashforth 2nd order

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

Adams Bashforth Moulton Method

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

RK4 step to use in next function

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

Adams-Bashforth Moulton Method 4th order

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

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

# 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)

Answer 1

The modified Euler method Does Not access points outside the step i -> i+1, there is no i-1 (note that in your source document the step, in the python code, not the formulas, is i-1 -> i with the loops starting at an appropriately increased index). It simply is (as you can find everywhere the mod. Euler or Heun method is discussed)

k1 = f(y[i]   , t[i  ])*dt;
k2 = f(y[i]+k1, t[i+1])*dt;
y[i+1] = y[i] + 0.5*(k1+k2);

In contrast, the Adams-Bashford method of order 2 and Adams-Moulton methods of order greater 2 Do access points from before the step i -> i+1, formally one has in AB2

y[i+1] = y[i] + 0.5*dt * (3*f[i] - f[i-1])

For a first implementation it would make sense to declare the f array the same way as the y array to implement this formula verbatim. It can be more economical to only keep a short array of f values that is shifted or rotated to give access to the last few f values.

Note that you need to initialize y[1] and f[1] with some other method of similar or higher order. Or if you want to have a "pure" run of the method, you need to initialize y[-1] and f[-1] and further back so that y[1] can be computed with the method formula.