HW2 solution

import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

def lv(u, p, t):

    x, y = u

    a, b, c, d = p

    dx = a * x - b * x * y

    dy = -c * y + d * x * y

    return np.array([dx, dy])

# ODE solver

def solve_euler(f, u0, tspan, h, p):

    ts = np.arange(tspan[0], tspan[1] + h, h)

    us = np.zeros((len(ts), len(u0)))

    us[0, :] = np.asarray(u0)

    for i in range(1, len(ts)):

        u = us[i-1,:]

        t = ts[i-1]

        nextu = u + h * f(u, p, t)

        us[i,:] = nextu

    return {"t": ts, "y": np.asarray(us)}

h = 0.001

u0 = [2.0, 1.0]

ts = (0.0, 30.0)

p = (2/3, 4/3, 1, 1)

sol = solve_euler(lv, u0, ts, h, p)

plt.figure(figsize=(10, 7))

plt.plot(sol["t"], sol["y"][:, 0], label="x")

plt.plot(sol["t"], sol["y"][:, 1], label="y")

plt.legend(loc="best")

# Comparison between time steps (h) 

u0 = [2.0, 1.0]

ts = (0.0, 30.0)

p = (2/3, 4/3, 1, 1)

sols = [solve_euler(lv, u0, ts, h, p) for h in (0.001, 0.005, 0.01, 0.05, 0.1)]

plt.figure(figsize=(15, 10))

for sol in sols:

    plt.plot(sol["t"], sol["y"][:, 0], label="x")

    plt.plot(sol["t"], sol["y"][:, 1], label="y")

def solve_rk4(f, u0, tspan, h, p):

    ts = np.arange(tspan[0], tspan[1] + h, h)

    us = np.zeros((len(ts), len(u0)))

    us[0, :] = np.asarray(u0)

    for i in range(1, len(ts)):

        u = us[i-1,:]

        t = ts[i-1]

        k1 = h * f(u, p, t)

        k2 = h * f(u + 0.5 * k1, p, t + 0.5 * h)

        k3 = h * f(u + 0.5 * k2, p, t + 0.5 * h)

        k4 = h * f(u + k3, p, t + h)

        nextu = u + (k1 + 2*k2 + 2*k3+ k4) / 6

        us[i,:] = nextu

    return {"t": ts, "y": np.asarray(us)}

# Comparison between FE and RK4 (h(RK4) = 0.01, h(FE) = 0.001)

u0 = [2.0, 1.0]

ts = (0.0, 30.0)

p = (2/3, 4/3, 1, 1)

solrk = solve_rk4(lv, u0, ts, 0.01, p)

solfe = solve_euler(lv, u0, ts, 0.001, p)

plt.figure(figsize=(10, 7))

plt.plot(solrk["t"], solrk["y"][:, 0], label="x (RK)")

plt.plot(solrk["t"], solrk["y"][:, 1], label="y (RK)")

plt.plot(solfe["t"], solfe["y"][:, 0], '--', label="x (FE)")

plt.plot(solfe["t"], solfe["y"][:, 1], '--', label="y (FE)")

plt.legend(loc="best")

# Comparison between FE and RK4 (h = 0.01)

u0 = [2.0, 1.0]

ts = (0.0, 30.0)

p = (2/3, 4/3, 1, 1)

solrk = solve_rk4(lv, u0, ts, 0.01, p)

solfe = solve_euler(lv, u0, ts, 0.01, p)

plt.figure(figsize=(10, 7))

plt.plot(solrk["t"], solrk["y"][:, 0], label="x (RK)")

plt.plot(solrk["t"], solrk["y"][:, 1], label="y (RK)")

plt.plot(solfe["t"], solfe["y"][:, 0], '--', label="x (FE)")

plt.plot(solfe["t"], solfe["y"][:, 1], '--', label="y (FE)")

plt.legend(loc="best")

# I use solve_ivp() instead of odeint() so that the solver could determine the time step itself

from scipy.integrate import solve_ivp

sol = solve_ivp(lambda t, y: lv(y, p, t), ts, u0, dense_output=True, method="RK45")

sol

# To get smooth curves,  use interpolants

t = np.linspace(ts[0], ts[1], num=500)

y = sol.sol(t).transpose()

plt.figure(figsize=(15, 10))

plt.plot(t, y)

plt.plot(solrk["t"], solrk["y"], '--')

plt.plot(solfe["t"], solfe["y"])

using Pkg

pkg"up; add DifferentialEquations Plots"

function lv(u, p, t)

    x, y = u

    a, b, c, d = p

    dx = a * x - b * x * y

    dy = -c * y + d * x * y

    return [dx, dy]

end

# ODE solver

function solve_euler(f, u0, tspan, h, p)

    ts = tspan[1]:h:tspan[end]

    us = zeros(length(ts), length(u0))

    us[1, :] .= u0

    for i in 2:length(ts)

        us[i, :] .= us[i-1, :] .+ h .* f(us[i-1, :], p, ts[i-1])

    end

    return (t = ts, u = us)

end

h = 0.001

u0 = [2.0, 1.0]

ts = (0.0, 30.0)

p = (2/3, 4/3, 1, 1)

ts, us = solve_euler(lv, u0, ts, h, p)

using Plots

plot(ts, us)

function solve_rk4(f, u0, tspan, h, p)

    ts = tspan[1]:h:tspan[end]

    us = zeros(length(ts), length(u0))

    us[1, :] .= u0

    for i in 2:length(ts)

        u = us[i-1, :]

        t = ts[i-1]

        k1 = h * f(u, p, t)

        k2 = h * f(u + 0.5 * k1, p, t + 0.5 * h)

        k3 = h * f(u + 0.5 * k2, p, t + 0.5 * h)

        k4 = h * f(u + k3, p, t + h)

        us[i,:] = u + (k1 + 2* k2 + 2*k3 + k4) / 6

    end

    return (t = ts, u = us)

end

solrk = solve_rk4(lv, u0, ts, 0.01, p)

plot(solrk.t, solrk.u)

solfe = solve_euler(lv, u0, ts, 0.01, p)

plot(solrk.t, solrk.u, labels=["x (RK4)" "y (RK4)"])

plot!(solfe.t, solfe.u, labels=["x (Euler)" "y (Euler)"], legend=:topleft)