diff --git a/UB4/UB4.py b/UB4/UB4.py index a31f962baf1f002d6ca1d110c83097550a9d30c3..7e087a7fbe08f77ab0d7c175160e44ff6df4a034 100644 --- a/UB4/UB4.py +++ b/UB4/UB4.py @@ -67,6 +67,36 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwar return y +def IE(y0, t, dt, lamda = -1): + 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] = 1 / (1 - dt*lamda) * y[i-1] + + return y + +def EE(y0, t, dt, lamda = -1): + 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 * lamda*y[i-1] + + return y + +def IM(y0, t, dt, lamda = -1): + 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] = (1 + lamda * dt/2) /(1 - lamda * dt/2) * y[i-1] + + return y + # Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1 y0 = np.array(1) @@ -83,17 +113,19 @@ exact_solution = lambda t : np.exp(lamda * t) * y0 ## 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 i in range(len(dt)): + plt.figure() + # 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)") + plt.plot(IE(y0, t[1], dt[i]), '--', label = "implicit euler (IE)") + plt.plot(EE(y0, t[1], dt[i]), 'x', label = "explicit euler (EE)") + plt.plot(IM(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[0]}, dt={dt[i]}') + plt.show() #### For A-stability, dt < 1 #### Observation @@ -254,84 +286,84 @@ y0 = np.array([0.5,0.2,0.3]) # p_dot = - partial H / partial q = -q # q_dot = partial H / partial p = p -p_0 = np.array([0]) -q_0 = np.array([1]) -t = 100 -dt = 0.001 - -def HO_IE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): - """ - Integrator = implicit Euler - p_dot = - partial H / partial q = -q - q_dot = partial H / partial p = p - """ - no_of_steps = int(t // dt) - p = np.zeros((no_of_steps, *p_0.shape)) - q = np.zeros((no_of_steps, *q_0.shape)) - p[0] = p_0 - q[0] = q_0 - - for i in trange(1, no_of_steps): - p[i] = 1/(1 + dt**2) * (p[i-1] - dt * q[i-1]) - q[i] = 1/(1 + dt**2) * (q[i-1] + dt * p[i-1]) - - return p, q - -def HO_EE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): - """ - Integrator = explicit Euler - p_dot = - partial H / partial q = -q - q_dot = partial H / partial p = p - """ - no_of_steps = int(t // dt) - p = np.zeros((no_of_steps, *p_0.shape)) - q = np.zeros((no_of_steps, *q_0.shape)) - p[0] = p_0 - q[0] = q_0 - - for i in trange(1, no_of_steps): - p[i] = p[i-1] + dt * -q[i-1] - q[i] = q[i-1] + dt * p[i-1] +# p_0 = np.array([0]) +# q_0 = np.array([1]) +# t = 100 +# dt = 0.001 + +# def HO_IE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): +# """ +# Integrator = implicit Euler +# p_dot = - partial H / partial q = -q +# q_dot = partial H / partial p = p +# """ +# no_of_steps = int(t // dt) +# p = np.zeros((no_of_steps, *p_0.shape)) +# q = np.zeros((no_of_steps, *q_0.shape)) +# p[0] = p_0 +# q[0] = q_0 + +# for i in trange(1, no_of_steps): +# p[i] = 1/(1 + dt**2) * (p[i-1] - dt * q[i-1]) +# q[i] = 1/(1 + dt**2) * (q[i-1] + dt * p[i-1]) + +# return p, q + +# def HO_EE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): +# """ +# Integrator = explicit Euler +# p_dot = - partial H / partial q = -q +# q_dot = partial H / partial p = p +# """ +# no_of_steps = int(t // dt) +# p = np.zeros((no_of_steps, *p_0.shape)) +# q = np.zeros((no_of_steps, *q_0.shape)) +# p[0] = p_0 +# q[0] = q_0 + +# for i in trange(1, no_of_steps): +# p[i] = p[i-1] + dt * -q[i-1] +# q[i] = q[i-1] + dt * p[i-1] - return p, q - -def HO_IM(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): - """ - Integrator = implicit midpoint - p_dot = - partial H / partial q = -q - q_dot = partial H / partial p = p - """ - no_of_steps = int(t // dt) - p = np.zeros((no_of_steps, *p_0.shape)) - q = np.zeros((no_of_steps, *q_0.shape)) - p[0] = p_0 - q[0] = q_0 - - for i in trange(1, no_of_steps): - q[i] = (q[i-1] + dt*(p[i-1] - dt/4 * q[i-1])) / (1 + dt**2 / 4) - p[i] = (p[i-1] - dt*(q[i-1] - dt/4 * p[i-1])) / (1 + dt**2 / 4) - - return p, q - -p, q = HO_IE(p_0, q_0, t, dt) -plt.figure() -plt.plot(p, q) -plt.title("Implicit Euler") -plt.show() -### Observation -> Could see that phase space volume increases - -p, q = HO_EE(p_0, q_0, t, dt) -plt.figure() -plt.plot(p, q) -plt.title("Explicit Euler") -plt.show() - -p, q = HO_IM(p_0, q_0, t, dt) -plt.figure() -plt.plot(p, q) -plt.title("Implicit midpoint") -plt.show() -### Observation -> Conserve phase space volume better than the other 2 methods +# return p, q + +# def HO_IM(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float): +# """ +# Integrator = implicit midpoint +# p_dot = - partial H / partial q = -q +# q_dot = partial H / partial p = p +# """ +# no_of_steps = int(t // dt) +# p = np.zeros((no_of_steps, *p_0.shape)) +# q = np.zeros((no_of_steps, *q_0.shape)) +# p[0] = p_0 +# q[0] = q_0 + +# for i in trange(1, no_of_steps): +# q[i] = (q[i-1] + dt*(p[i-1] - dt/4 * q[i-1])) / (1 + dt**2 / 4) +# p[i] = (p[i-1] - dt*(q[i-1] - dt/4 * p[i-1])) / (1 + dt**2 / 4) + +# return p, q + +# p, q = HO_IE(p_0, q_0, t, dt) +# plt.figure() +# plt.plot(p, q) +# plt.title("Implicit Euler") +# plt.show() +# ### Observation -> Could see that phase space volume increases + +# p, q = HO_EE(p_0, q_0, t, dt) +# plt.figure() +# plt.plot(p, q) +# plt.title("Explicit Euler") +# plt.show() + +# p, q = HO_IM(p_0, q_0, t, dt) +# plt.figure() +# plt.plot(p, q) +# plt.title("Implicit midpoint") +# plt.show() +# ### Observation -> Conserve phase space volume better than the other 2 methods