ODE_methods (1).ipynb

Kernel: Python 3

!pip install torchdiffeq
!pip install git+https://github.com/rtqichen/torchdiffeq
import numpy as np
import matplotlib.pyplot as plt

Collecting torchdiffeq
  Downloading torchdiffeq-0.2.5-py3-none-any.whl.metadata (440 bytes)
Requirement already satisfied: torch>=1.5.0 in /usr/local/lib/python3.10/dist-packages (from torchdiffeq) (2.5.1+cu121)
Requirement already satisfied: scipy>=1.4.0 in /usr/local/lib/python3.10/dist-packages (from torchdiffeq) (1.13.1)
Requirement already satisfied: numpy<2.3,>=1.22.4 in /usr/local/lib/python3.10/dist-packages (from scipy>=1.4.0->torchdiffeq) (1.26.4)
Requirement already satisfied: filelock in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (3.16.1)
Requirement already satisfied: typing-extensions>=4.8.0 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (4.12.2)
Requirement already satisfied: networkx in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (3.4.2)
Requirement already satisfied: jinja2 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (3.1.4)
Requirement already satisfied: fsspec in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (2024.10.0)
Requirement already satisfied: sympy==1.13.1 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq) (1.13.1)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from sympy==1.13.1->torch>=1.5.0->torchdiffeq) (1.3.0)
Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from jinja2->torch>=1.5.0->torchdiffeq) (3.0.2)
Downloading torchdiffeq-0.2.5-py3-none-any.whl (32 kB)
Installing collected packages: torchdiffeq
Successfully installed torchdiffeq-0.2.5
Collecting git+https://github.com/rtqichen/torchdiffeq
  Cloning https://github.com/rtqichen/torchdiffeq to /tmp/pip-req-build-haygst08
  Running command git clone --filter=blob:none --quiet https://github.com/rtqichen/torchdiffeq /tmp/pip-req-build-haygst08
  Resolved https://github.com/rtqichen/torchdiffeq to commit a88aac53cae738addee44251288ce5be9a018af3
  Preparing metadata (setup.py) ... done
Requirement already satisfied: torch>=1.5.0 in /usr/local/lib/python3.10/dist-packages (from torchdiffeq==0.2.5) (2.5.1+cu121)
Requirement already satisfied: scipy>=1.4.0 in /usr/local/lib/python3.10/dist-packages (from torchdiffeq==0.2.5) (1.13.1)
Requirement already satisfied: numpy<2.3,>=1.22.4 in /usr/local/lib/python3.10/dist-packages (from scipy>=1.4.0->torchdiffeq==0.2.5) (1.26.4)
Requirement already satisfied: filelock in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (3.16.1)
Requirement already satisfied: typing-extensions>=4.8.0 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (4.12.2)
Requirement already satisfied: networkx in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (3.4.2)
Requirement already satisfied: jinja2 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (3.1.4)
Requirement already satisfied: fsspec in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (2024.10.0)
Requirement already satisfied: sympy==1.13.1 in /usr/local/lib/python3.10/dist-packages (from torch>=1.5.0->torchdiffeq==0.2.5) (1.13.1)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from sympy==1.13.1->torch>=1.5.0->torchdiffeq==0.2.5) (1.3.0)
Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from jinja2->torch>=1.5.0->torchdiffeq==0.2.5) (3.0.2)

Implementations of Euler's and RK4 methods

Euler's method is one of the oldest and simplest algorithms. The core idea is approximating the solution function step by step using tangents lines.

$\frac{dy}{dt} = f(t, y)$ and $y(t_0) = y_0$, then $y_t = y_0 + f(t_0, y_0)(t-t_0)$

---

Runge Kutta 4 method

The Runge-Kutta 4th Order method is a numerical technique for solving ordinary differential equations (ODEs). It provides a balance between accuracy and computational efficiency by approximating the solution at each time step using four intermediate evaluations of the derivative. The 4th order Runge-Kutta method approximates the curve of the function f between two points by a 4th degree polynomial. It requires 4 slope estimates at each step. Geometrically, we can visualize the curve between two points as part of a 4th degree polynomial. The 4th order Runge-Kutta method finds 4 tangents to this curve at evenly spaced points, and uses the average of their slopes as an approximation of the entire curve between the two points.

Here we are going to write function for each method and then utilize it initial value problem solver function. Next, set a table and graph for showing how far the hidden point of t from the exact point:

def solveIVP(f, tspan, y0, h, solver):

    # Initialise t and y arrays
    t = np.arange(tspan[0], tspan[1] + h, h)
    y = np.zeros((len(t),len(y0)))
    y[0,:] = y0

    # Loop through steps and calculate single step solver solution
    for n in range(len(t) - 1):
        y[n+1,:] = solver(f, t[n], y[n,:], h)

    return t, y


def euler(f, t, y, h):
    return y + h * f(t, y)

def rk4(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
    k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
    k4 = f(t + h, y + h * k3)
    return y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

# Define the ODE function
def f(t, y):
    return t * y

def exact(t):
  return np.exp(t**2/2)


# Define IVP parameters
tspan = [0, 1]  # boundaries of the t domain
y0 = [1]        # initial values
h = 0.2        # step length

# Calculate the solution to the IVP
t, yEuler = solveIVP(f, tspan, y0, h, euler)
t, yRK4 = solveIVP(f, tspan, y0, h, rk4)

# Print table of solution values
print("|  t   |     yEuler     |    RK4    |   Exact   |")
print("|------|-----------|-----------|-----------|")
for n in range(len(t)):
    print(f"| {t[n]:4.2f} | {yEuler[n,0]:9.6f} | {yRK4[n,0]:9.6f} | {exact(t[n]):9.6} |" )

#calculate exact solution
tExact=np.linspace(tspan[0], tspan[1], 200)
yExact=exact(tExact)
# Plot solution
fig, ax = plt.subplots()
plt.plot(tExact, yExact, "k-", label="Exact Solution")
plt.plot(t, yEuler, "bo-", label="Euler")
plt.plot(t, yRK4, "r-o", label="RK4 method")
plt.xlabel("$t$", fontsize=12)
plt.ylabel("$y$", fontsize=12)
plt.legend()
plt.show()

|  t   |     y     |    RK4    |   Exact   |
|------|-----------|-----------|-----------|
| 0.00 |  1.000000 |  1.000000 |       1.0 |
| 0.20 |  1.000000 |  1.020201 |    1.0202 |
| 0.40 |  1.040000 |  1.083287 |   1.08329 |
| 0.60 |  1.123200 |  1.197217 |   1.19722 |
| 0.80 |  1.257984 |  1.377126 |   1.37713 |
| 1.00 |  1.459261 |  1.648717 |   1.64872 |

# Define ODE function and the exact solution
def f(t,y):
    return t * y


def exact(t):
    return np.exp(t ** 2 / 2)


# Define IVP parameters
tspan = [0, 1]
y0 = [1]
hvals = [ 0.2, 0.1, 0.05, 0.025 ]  # step length values

# Loop through h values and calculate errors
errors = []
errorsRk4 = []
print(f"Exact solution: y(1) = {exact(1):0.6f}\n")
print("|   h   |  Euler   |  Error   |   RK4   |  Error   |")
print("|:-----:|:--------:|:--------:|:-------:|:--------:|")
for h in hvals:
    t, y = solveIVP(f, tspan, y0, h, euler)
    t, yrk4 = solveIVP(f, tspan, y0, h, rk4)
    errors.append(abs(y[-1,0] - exact(t[-1])))
    errorsRk4.append(abs(yrk4[-1,0] - exact(t[-1])))
    print(f"| {h:0.3f} | {y[-1,0]:0.6f} | {errors[-1]:0.2e} | {yrk4[-1,0]:0.6f} | {errorsRk4[-1]:0.2e} |")

# Plot errors
fig, ax = plt.subplots()
plt.plot(hvals, errors, "bo-")
plt.plot(hvals, errorsRk4, "ro-")
plt.xlabel("$h$", fontsize=12)
plt.ylabel("$E(h)$", fontsize=12)
plt.show()

Exact solution: y(1) = 1.648721

|   h   |  Euler   |  Error   |   RK4   |  Error   |
|:-----:|:--------:|:--------:|:-------:|:--------:|
| 0.200 | 1.459261 | 1.89e-01 | 1.648717 | 4.59e-06 |
| 0.100 | 1.547110 | 1.02e-01 | 1.648721 | 2.64e-07 |
| 0.050 | 1.595942 | 5.28e-02 | 1.648721 | 1.55e-08 |
| 0.025 | 1.621801 | 2.69e-02 | 1.648721 | 9.33e-10 |