Explicit (Forward) and Implicit (Backward) Euler Methods in Python

Question

I am trying to implement both the explicit and implicit Euler methods to approximate a solution for the following ODE:

dx/dt = -kx, where k = cos(2 pi t), and x(0) = 1

Euler's methods use finite differencing to approximate a derivative:

dx/dt = (x(t+dt) - x(t)) / dt

The forward method explicitly calculates x(t+dt) based on a previous solution x(t):

x(t+dt) = x(t) + f(x,t)dt

The backwards method is implicit, and finds the solution x(t+dt) by solving an equation involving the current state of the system x(t) and the later one x(t+dt):

x(t) = x(t+dt) - f(x,(t+dt))dt

My code for approximating a solution to dx/dt = -kx, x(0) = 1 and plotting it alongside the actual solution is given below:

### Import Necessary Packages
import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_format = 'retina'
plt.rcParams['figure.figsize'] = (10.0, 6.0)

### Defining Basic Data
t0 = 0                                # initial t
tf = 4*np.pi                          # final t       
N = 1000                              # factor affecting time step size
dt = (4*np.pi)/N                      # time step size (as a factor of N)
t = np.linspace(t0,tf,N)              # defining a vector of t values from t0 to tf  
x0 = 1                                # initial x
x = np.zeros([N])                     # initializing array for x values
f = lambda x,t: -np.cos(2*np.pi*t)*x  # defining f(x,t) on RHS on ODE

### Define a Function for Euler's Forward Method ###
def ForwardEuler(f,x0,t):
    x[0] = x0
    for i in range(1,N-1):
        x[i+1] = x[i] + (f(x[i],t[i]))*dt
    return x

# Plot Solution
forward = ForwardEuler(f,x0,t)
actual = 1/np.exp((1/(2*np.pi))*np.sin(2*np.pi*t))
plt.plot(t,actual,'r-',t,forward,'b-')
plt.legend(['Actual','Backward Euler'])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title("Solution to $x'=-kx$,   $x$(0)=1")
plt.grid(True);plt.show()

My question lies in how to adapt the for-loop section of the code to display the backward Euler method instead of forward Euler method. I am having trouble with this since the equations require you to know x[i+1] in order to solve x[i+1].

I believe the backwards for-loop would be what is given below, but I am unsure:

def BackwardEuler(f,x0,t):
    x[0] = x0
    for i in range(1,N-1):
        x[i] = x[i+1] - (f(x[i+1],t[i+1]))*dt
    return x

I have found very few resources online and am at a complete loss. Any help on this would be appreciated! Thank you!

Aziz0000_83605 · Answer

Usually, for Backward Euler and Trapezoidal Rule, you write the expression as a equation (or a system of equations), then solve it (i.e. find its zeros). The zeros represent the value of x[i+1].

For example, for Backward Euler, the system is:

x[i+1] = x[i] + (f(x[i+1],t[i+1]))*dt

Which you can rewrite as:

x[i+1] - x[i] - dt*f(x[i+1], t[i+1]) = 0

The values x[i] and t[i+1] are known. The only unknown is x[i+1]. You can solve this system numerically (using something like fsolve), and the solution would be your x[i+1]. It is possible, of course, that you get more than one solution. You have to select the one that fits your problem (i.e. x cannot be an imaginary number, or x cannot be negative, etc.)

The same technique can be applied for Trapezoidal Rule, with the system being:

x[i+1] - x[i] - (f(x[i],t[i]) + f(x[i+1],t[i+1]))*(dt/2) = 0

PS: Check out Computational Science StackExchange. It is more suitable for question related to numerical and computational methods.

Explicit (Forward) and Implicit (Backward) Euler Methods in Python

1 Answers