But Why Does This Matter

The reason that we want to reframe ODEs in this way is because of the following fact:

For all constant-linear ODEs, we can express the ODE as a matrix such that applying it to any point in the vector space would map any point to a valid point on the curve defined by the ODE

Looking at the equations above, these matrices (\(T\)) are:

1. Simple Harmonic Motion: \[T = \begin{bmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & -\omega^2 & 0 & 0\\ \end{bmatrix} \]

2. Radioactive Decay: \[T = \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & -\lambda & 0\\ \end{bmatrix}\]

3. RC Circuit Equation: \[T = \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & \frac{-1}{RC} & 0\\ \end{bmatrix}\]

4. Damped Harmonic Oscillator: \[T = \begin{bmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & -\omega^2 & -2\gamma & 0\\ \end{bmatrix}\]

5. Heat Equation (One-Dimensional): \[T = \begin{bmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & \frac{1}{\alpha} & 0\\ \end{bmatrix}\]

6. Exponential Growth or Decay: \[T = \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & k & 0\\ \end{bmatrix}\]

Using these to fit ODEs

Now that we can express the ODEs in the form of a matrix, we can implement these matriexies in the ODE solver package to make the solution fit the ode. It’s important here to note that I’ve diverted from my old definitions of \(Y\) here, where the first element of the vector is \(y(x)\).

To make a step in the approximation we use the following equation:

\[ \begin{bmatrix} 1 \\ x+h \\ y(x+h)\\ y'(x+h)\\ y''(x+h)\\ ...\\ y^{n}(x+h)\\ \end{bmatrix} = S \cdot \begin{bmatrix} 1 \\ x\\ y(x)\\ y'(x)\\ y''(x)\\ ...\\ y^{n}(x)\\ \end{bmatrix} \epsilon \]

Where \(S\) is: \[ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & ... & 0 \\ h & 1 & 0 & 0 & 0 & ... & 0 \\ 0 & 0 & 1 & \frac{h}{1!} & \frac{h^2}{2!} & ... & \frac{h^n}{n!}\\ 0 & 0 & 0 & 1 & \frac{h}{1!} & ... & \frac{h^{n-1}}{(n-1)!}\\ 0 & 0 & 0 & 0 & 1 & ... & \frac{h^{n-2}}{(n-2)!}\\ ... & ... & ... & ... & ... & ...\\ 0 & 0 & 0 & 0 & 0 & ... & 1\\ \end{bmatrix}\]

When making this step, the error in the approximation will move the point away from the plane that contains all valid solutions to the ODE, and therefore we will have to snap it back using one of the transformation matrices (\(T\)).

Implementing this method in our python library:

def expanded_euler(dims, h):
    step_matrix = np.zeros((dims, dims))
    for i in range(dims):
        for j in range(i, dims):
            # Is 1, and h at j-i =0, 1 respectively
            step_matrix[i, j] = h ** (j - i) / math.factorial(j - i)
    expanded_matrix = add_x_and_1(step_matrix, h)
    return expanded_matrix


def add_x_and_1(original_matrix, h):
    new_size = len(original_matrix) + 2
    new_matrix = np.zeros((new_size, new_size), dtype=original_matrix.dtype)

    # Set the 2x2 top left matrix
    new_matrix[0:2, 0:2] = [[1, 0], [h, 1]]

    # Copy the original matrix to the bottom right of the new matrix.
    new_matrix[2:, 2:] = original_matrix
    return new_matrix


def linear(y, step_matrix_generator, transformation_matrix, steps=10, h=0.1):
    dims = len(y) - 2
    step_matrix = transformation_matrix @ step_matrix_generator(dims, h)
    output_list = []

    y_n = y.copy()
    i = 0
    while i < steps:
        y_n = step_matrix @ y_n
        output_list.append(y_n)
        i += 1

Bind this machinery together, and you get a tool capable of tackling the initial example of \(y' = 2x\) passing through the point (0,0):

init_y = [1,0,0,0] #[1,x,y,y']
transformation_matrix = np.array([
   [ 1,0,0,0 ],
   [ 0,1,0,0 ],
   [ 0,0,1,0 ],
   [ 0,2,0,0 ]
])
solution = linear(
    init_y,
    expanded_euler,
    transformation_matrix,
    steps=100, h=0.01)