Skip to content
Snippets Groups Projects
Commit b7dd52c7 authored by nguyed99's avatar nguyed99
Browse files

Update UB4

parent 50e5d686
No related branches found
No related tags found
No related merge requests found
import numpy as np
import matplotlib.pyplot as plt
from tqdm import trange
# Ex1: Numerical integrators
def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
def implicit_euler(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwargs) -> np.ndarray:
"""
A-stable
:param f: function to be integrated
......@@ -14,10 +15,15 @@ def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarr
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
for i in trange(1, no_of_steps):
y_0 = y[i-1] + dt * f(y[i-1]) # predicted state obtained through EE
y_1 = y[i-1] + dt * f(y_0)
while abs(y_1 - y_0) > eps:
y_0 = y_1
y_1 = y[i-1] + dt * f(y_0)
y[i] = y_1
return y
......@@ -32,14 +38,13 @@ def explicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarr
y = np.zeros((no_of_steps, *y0.shape))
y[0] = y0
for i in range(1, no_of_steps):
for i in trange(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:
def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **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
......@@ -48,18 +53,24 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.nd
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)
for i in trange(1, no_of_steps):
y_0 = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1], **kwargs), **kwargs) # predicted state obtained through IM
y_1 = y[i-1] + dt * f( 1/2 * (y[i-1] + y_0) ,**kwargs)
while abs(y_1 - y_0) > eps:
y_0 = y_1
y_1 = y[i-1] + dt * f( 1/2 * (y[i-1] + y_0) ,**kwargs)
y[i] = y_1
return y
# Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
y0 = np.array([1])
y0 = np.array(1)
lamda = -1
t = [10, 100]
dt = [1.5, 1, 1e-1, 1e-2]
dt = [0.9, 1e-1, 1e-2]
f = lambda y: lamda * y
exact_solution = lambda t : np.exp(lamda * t) * y0
......@@ -84,14 +95,8 @@ exact_solution = lambda t : np.exp(lamda * t) * y0
#### 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 = 0.9 <-- Strange behaviour !
# IM decas as rapidly as ES, IE not as fast, EE too fast
### 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 !
......@@ -110,17 +115,14 @@ for i in range(len(dt)):
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 !
### t = 100, dt = 0.9 <-- Strange behaviour !
# Same as in the case t = 10, dt = 0.9
### t = 100, dt = 0.1 <-- Strange behaviour !
# Same as in the case t = 100, dt = 0.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'
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment