CODEWITHSUNDEEP
pythonmathnumerical-methodsode
killer croc
killer croc

Reputation: 70

Implementing the Backwards Euler method in python to solve a pendulum

I am trying to set up an implicit solver to a pendulum F and dF are defined as:

def func(y,t):
    ### Simplified the Function to remove friction since it canceled 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[3],y[4],y[5]])[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[1]
    phiz=2.*y[2]
    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[3]
    dF[3,4]=phix*y[4]
    dF[3,5]=phix*y[5]
    dF[4,0]=(y[1]*F1)/mass
    dF[4,1]=(y[1]*F2+2*lambd)/mass
    dF[4,2]=(y[1]*F3)/mass
    dF[4,3]=phiy*y[3]
    dF[4,4]=phiy*y[4]
    dF[4,5]=phiy*y[5]
    dF[5,0]=(y[2]*F1)/mass
    dF[5,1]=(y[2]*F2)/mass
    dF[5,2]=(y[2]*F3+2*lambd)/mass
    dF[5,3]=phiz*y[3]
    dF[5,4]=phiz*y[4]
    dF[5,5]=phiz*y[5]

    return dF

I am trying to set up an implicit solver that uses the function and the dF in order to guess values in order to create a more accurate solver for my pendulum function.

the error that I get is the following

Traceback (most recent call last): line 206, in test_function(backwards_Euler) line 186, in test_function y1 = test_function(func, y0, t) line 166, in backwards_Euler F = np.asarray(y[i] + dt * function(zold, time[i+1])-zold) line 13, in func lambd = (grav.dot(x)+v.dot(v))/x.dot(x)

ValueError: shapes (3,6) and (3,6) not aligned: 6 (dim 1) != 3 (dim 0)

I think it is because one of the calculations I am doing is diverging to infinity but I cannot figure out the issue that might be causing this.

Here is my current backward Euler function:

### Backwards Euler Method

def backwards_Euler(function, y_matrix, time):
    y = np.zeros((np.size(time), np.size(y_matrix)))
    y[0, :] = y_matrix

    dt = time[1] - time[0]

    for i in range(len(time-1)):
        err = 1
        zold = y[i] + dt*function(y[i],time[i])  ## guess with forward euler
        I = 0

        while err > 10**(-10) and I < 5:
            F = y[i] + dt * function(zold, time[i+1])-zold  ## Here is where my error occurs
            dF = dt*dF_matrix(y[i+1])-1
            znew = zold - F/dF
            zold = znew
            I+=1
        y[i+1]=znew
    return y

Here is how I am testing the function although currently, I have no output

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

Upvotes: 1

Views: 6483

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

No, the error is that the matrix shapes in some matrix operations don't match. The immediate cause is that x.dot(x) works when x is a simple array, but not for two-dimensional arrays. Then the rules of matrix multiplication have to be strictly obeyed.

However, there is no reason that there is a matrix in place of x, you do not call the function in a "vectorized" fashion. The true reason is that python is not Matlab, the division in

znew = zold - F/dF

does not do what you think it does. Instead it does elementwise division, repeating the smaller object to fill the shape of the larger, which also happens with the subtraction, so that znew in the end has the same shape as the matrix dF, with the observed results in the next ODE function evaluation.

You need to explicitly invoke a linear solver, the standard one in in np.linalg (there are others there and in scipy.sparse.linalg using other matrix factorizations) is simply called solve(A,b) and computes the solution of A*x=b.

znew = zold - np.linalg.solve(dF, F)

Implicit Euler gives a diverging solution, the length of the pendulum increases rapidly. Applying these methods to the similar implicit trapezoidal method, which is also Adams-Moulton 2nd order, gives the code

def Adams_Moulton_2nd(function, y_init, time):
    # solve F(z)=0 with simplified Newton, where
    # F    = y + 0.5*dt*(f(y,t)+f(z,t+dt)) - z
    # is the defect of the implicit trapezium method. Using the derivative
    # dF_z = 0.5*dt*df_z(z,t+dt) - I
    # the non-linear solver step is solution of the linear system
    # dF_z*dz = -F
    y = np.zeros((np.size(time), np.size(y_init)))
    y[0, :] = y_init
    Id = np.eye(len(y_init))
    dt = time[1] - time[0]
    dt4 = dt**4

    for i in range(len(time)-1):
        f_0 = function(y[i],time[i])
        norm_0 = sum(f_0**2)**0.5 + 1e-6
        z = y[i] + dt*f_0  ## guess with forward euler
        dF = 0.5*dt*dF_matrix(z, t[i+1]) - Id
        for I in range(5):
            F = y[i] + 0.5*dt * (f_0+function(z, time[i+1])) - z
            dz = - np.linalg.solve(dF,F)
            if max(abs(dz)) < norm_0*dt4: break
            z += dz
        y[i+1]=z
    return y### Adams_Moulton_2nd

Upvotes: 1

Related Questions