Basic syntax

For example, to solve the initial value problem $$ \begin{align} u''(t) &= -u \\ u(0) &= 0 \\ u'(0) &= 1 \end{align} $$ we can express it as a first-order system $$ \begin{align} \vec{u}' &= \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \vec{u} \\ \vec{u} &= \begin{bmatrix}0 \\ 1 \end{bmatrix} \end{align} $$ then solve the system with a time-integrator like RK4 like so:

import numpy as np
from odeiter import TimeDomain, RK4
u0 = np.array([0.0, 1.0])

def rhs(t, u):
    return np.array([[0, 1], [-1, 0]]) @ u

time = TimeDomain(0, 0.1, 4)
solver = RK4()
for u in solver.solution_generator(u0, rhs, time):
    # so something with the solution at this time
    print(u)

# OUTPUT
# [0. 1.]
# [0.09983333 0.99500417]
# [0.19866917 0.9800666 ]
# [0.29551996 0.95533654]
# [0.38941803 0.9210611 ]

Interactive Plotting

This shows some of the flexibly by solving the system above, plotting the solution and plotting the error.

The plotting occurs as the system is solved, rather than after it is solved and requires pyqt5 and matplotlib to function properly.

"""
Solve the second order equation
x'' = -x
x(0) = 0
x'(0) = 1

which has the solution
x(t) = sin(t)
"""
import matplotlib.pyplot as plt
import numpy as np
from odeiter import TimeDomain_Start_Stop_MaxSpacing, RK4

x0 = 0
y0 = 1
A = np.array(
    [
        [0, 1],
        [-1, 0],
    ]
)
u0 = np.array([x0, y0])


def rhs(t, u):
    return A @ u


def exact(t):  # exact solution for testing
    return np.sin(t)


t0, tf = 0, 2 * np.pi
max_time_step = 1e-2
# Create a TimeDomain object
time = TimeDomain_Start_Stop_MaxSpacing(t0, tf, max_time_step)
# Choose a solver
solver = RK4()
plt.ion()  # requires pyqt5 for interactive plotting
for t, u in zip(time.array, solver.solution_generator(u0, rhs, time)):
    # do whatever you want with the solution
    x, y = u
    plt.plot(t, x, "k.")
    plt.pause(1e-3)
plt.show(block=True)