CODEWITHSUNDEEP
pythonjupyterode
Savannah
Savannah

Reputation: 35

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

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!

Upvotes: 2

Views: 11168

Answers (2)

Martinus
Martinus

Reputation: 1

In your for-loop you need the range to start at 0, otherwise you are missing x[0]: for i in range[0,N-1]:

As a result I saw that in time the sinewave slowly goes down. This is what the 'actual' function does not do. It has to do with the size of td. If you make td smaller, the sinewave decends less and the error is less between 'actual' and 'forward'.

Upvotes: 0

Aziz
Aziz

Reputation: 20705

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.

Upvotes: 2

Related Questions