UB4.py

import numpy as np
import matplotlib.pyplot  as plt

# Ex1: Numerical integrators 

def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
    """
    A-stable
    :param f: function to be integrated
    :param y0: initial value
    :param t: time interval
    :param dt: time step
    """
    no_of_steps = int(t // dt)
    y = np.zeros((no_of_steps, *y0.shape))
    y[0] = y0
    
    for i in range(1, no_of_steps):
        y[i] = y[i-1] + f(y[i-1], **kwargs) * dt
        y[i] = y[i-1] + f(y[i], **kwargs) * dt

    return y

def explicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
    """
    :param f: function to be integrated
    :param y0: initial value
    :param t: time interval
    :param dt: time step
    """
    no_of_steps = int(t // dt)
    y = np.zeros((no_of_steps, *y0.shape))
    y[0] = y0

    for i in range(1, no_of_steps):
        y[i] = y[i-1] + f(y[i-1], **kwargs) * dt
    
    return y

def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
    """
    Not L-stable, doesn't decay properly - oscillation
    :param f: function to be integrated
    :param y0: initial value
    :param t: time interval
    :param dt: time step
    """
    no_of_steps = int(t // dt)
    y = np.zeros((no_of_steps, *y0.shape))
    y[0] = y0
    
    for i in range(1, no_of_steps):
        y[i] = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1], **kwargs), **kwargs)
        y[i] = y[i-1] + dt * f(1/2 *  (y[i-1] + y[i]), **kwargs)
    
    return y

# Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
y0 = np.array([1])
lamda = -1
t = [10, 100]
dt = [1.5, 1, 1e-1, 1e-2]
f = lambda y: lamda * y
exact_solution = lambda t : np.exp(lamda * t) * y0

# A-stability: the method should produce a numerical solution that does not grow in magnitude for any dt
# L-stability: stronger form of A-stability -> numerical solution decays (at least) as rapidly as the exact solution
# Expectation: 
## EE: A-stable (No), L-stable (No)
## IE: A-stable (Yes), L-stable (Yes)
## IM: A-stable (Yes), L-stable (No)

# for i in range(len(dt)):
#     plt.figure()
#     y_implicit_euler = implicit_euler(f, y0, t[0], dt[i])
#     plt.plot(implicit_euler(f, y0, t[0], dt[i]), '--', label = "implicit euler (IE)")
#     plt.plot(explicit_euler(f, y0, t[0], dt[i]), 'x', label = "explicit euler (EE)")
#     plt.plot(implicit_midpoint(f, y0, t[0], dt[i]), '>', label = "implicit midpoint (IM)")
#     time = np.arange(int(t[0] // dt[i]))
#     plt.plot(exact_solution(time), '3', label="exact solution (ES)")
#     plt.legend()
#     plt.title(f't={t[0]}, dt={dt[i]}')
#     plt.show()

#### For A-stability, dt < 1
#### Observation
### t = 10, dt = 1.5 <-- Strange behaviour !
# IE grows and EE decays!
# IM is robust

### t = 10, dt = 1 <-- Strange behaviour !
# IE gives a constant function (but it doesn't blow up as time progresses!).
# However, this means it's not A-stable because it does not appropriately damp the solution. 
# EE & IM decay faster than ES -> L-stable. EE gives a worse approximation in comparison to IM
### t = 10, dt = 1e-1  <-- Strange behaviour ! 
# They all behave similarly. They do not decay as rapidly as ES.
### t = 10, dt = 1e-2 <-- Strange behaviour !
# They all behave similarly. They do not decay as rapidly as ES.

for i in range(len(dt)):
    plt.figure()
    y_implicit_euler = implicit_euler(f, y0, t[1], dt[i])
    plt.plot(implicit_euler(f, y0, t[1], dt[i]), '--', label = "implicit euler (IE)")
    plt.plot(explicit_euler(f, y0, t[1], dt[i]), 'x', label = "explicit euler (EE)")
    plt.plot(implicit_midpoint(f, y0, t[1], dt[i]), '>', label = "implicit midpoint (IM)")
    time = np.arange(int(t[1] // dt[i]))
    plt.plot(exact_solution(time), '3', label="exact solution (ES)")
    plt.legend()
    plt.title(f't={t[1]}, dt={dt[i]}')
    plt.show()

#### Observation
### t = 100, dt = 1.5 <-- Strange behaviour !
# IE blows up at the end.
### t = 100, dt = 1 <-- Strange behaviour !
# Same as in the case t = 10, dt = 1
### t = 100, dt = 1e-1  <-- Strange behaviour ! 
# Same as in the case t = 100, dt = 1e-1
### t = 100, dt = 1e-2 <-- Strange behaviour !
# Same as in the case t = 100, dt = 1e-2



# Ex3: 
# def f(y: np.ndarray, k: float) -> np.ndarray:
#     assert y.shape[0] == 3, 'y has the wrong dimension. It should be 3'

#     A = np.array([[-1, 1, 0],
#                     [1, -1-k, k],
#                     [0, k, -k]])
    
#     return np.dot(A, y)

# dt = [1, 1e-1, 1e-2, 1e-3]
# t = [10,100]
# y0 = np.array([1,0,0])

# # short time
# fig = plt.figure()
# for dt_i in dt:
#     y = implicit_euler(f, y0, t[0], dt_i, k=1)
#     plt.plot(y[:,0], '--',label = 'State 1')
#     plt.plot(y[:,1], 'x', label = 'State 2')
#     plt.plot(y[:,2], '>', label = 'State 3')
#     plt.legend()
#     plt.title(f'{integrators[0]}, dt = {dt_i}')
#     plt.show()


# fig = plt.figure()
# for dt_i in dt:
#     y = explicit_euler(f, y0, t[0], dt_i, k=1)
#     plt.plot(y[:,0], '--', label = 'State 1')
#     plt.plot(y[:,1], 'x', label = 'State 2')
#     plt.plot(y[:,2], '>', label = 'State 3')
#     plt.legend()
#     plt.title(f'{integrators[1]}, dt = {dt_i}')
#     plt.show()


# fig = plt.figure()
# for dt_i in dt:
#     y = implicit_midpoint(f, y0, t[0], dt_i, k=1)
#     plt.plot(y[:,0], '--', label = 'State 1')
#     plt.plot(y[:,1], 'x', label = 'State 2')
#     plt.plot(y[:,2], '>', label = 'State 3')
#     plt.legend()
#     plt.title(f'{integrators[2]}, dt = {dt_i}')
#     plt.show()