# Differential Equations

### Big Ideas
Since most differential equations cannot be solved exactly, we use *numerical methods* to find approximate solutions. The **finite difference method** is one such technique that approximates a differential equation with a system of algebraic equations. For a linear differential equation, this process results in a linear system of equations, $A \mathbf{y} = \mathbf{b}$, which we can solve to find the approximate solution.

## Second-Order Differential Equations
:::info
**Definition**
* An **ordinary differential equation (ODE)** is an equation that involves an unknown function, $y(t)$, and its derivatives.
* The **order** of an ODE is determined by the highest derivative present.
* In this section, we focus on **second-order linear ODEs**, which have the form: $$y'' + p(t)y' + q(t)y = r(t)$$
* To find a specific solution, we must impose **boundary conditions**, which are constraints on the solution or its derivatives at the domain's endpoints, such as $$y(t_0) = \alpha, \quad y(t_f) = \beta \quad \text{or} \quad y'(t_0) = \alpha \hspace{5mm} y(t_f) = \beta$$
:::

## Finite Difference Method
The [**finite difference method**](https://en.wikipedia.org/wiki/Finite_difference) uses formulas derived from [Taylor series]((https://en.wikipedia.org/wiki/Taylor_series)) to approximate the derivatives in a differential equation.

:::info
**Definition** Let $y(t)$ be a smooth function and consider the Taylor series
$$
\begin{align*}
y(t+h) &= y(t) + y'(t)h + \frac{y''(t)}{2}h^2 + \cdots \\
y(t-h) &= y(t) - y'(t)h + \frac{y''(t)}{2}h^2 + \cdots
\end{align*}
$$

* Truncate and rearrange the first series for the **forward difference formula**
$$
y'(t) \approx \frac{y(t+h) - y(t)}{h}
$$
* Truncate and rearrange the second series for the **backward difference formula**
$$
y'(t) \approx \frac{y(t) - y(t - h)}{h}
$$
* Subtract $y(t+h) - y(t-h)$ and truncate to get the (first order) **central difference formula**
$$
y'(t) \approx \frac{y(t+h) - y(t-h)}{2h}
$$
* Add $y(t+h) + y(t-h)$ and truncate to get the (second order) **central difference formula**
$$
y''(t) \approx \frac{y(t+h) -2y(t) + y(t-h)}{h^2}
$$
:::

:::info
**Definition.** The **finite difference method** uses formulas derived from Taylor series expansions to approximate derivatives.

The core of the method involves four steps for a second-order linear ODE with boundary conditions:

1. **Discretize the domain:** Divide the interval $[t_0,t_f]$ into $N+1$ equally spaced points, $t_k$, using a step size $h = \frac{t_f - t_0}{N+1}$.
2. **Approximate the solution:** Let $y_k \approx y(t_k)$ be the approximation of the solution at each point.
3. **Substitute finite difference formulas:** Replace the derivatives in the differential equation with their finite difference approximations. For example, the second-order central difference formula for $y''$ is: $$y''(t) \approx \frac{y(t+h) -2y(t) + y(t-h)}{h^2}$$
4. **Solve the linear system:** Arrange the resulting system of equations into the matrix form $A \mathbf{y} = \mathbf{b}$ and solve for the vector of approximate solutions, $\mathbf{y} = [y_1 \ \cdots \ y_N]^T$.
:::

This method is valuable for handling complex equations without analytical solutions. While it typically produces approximations, a smaller step size $h$ leads to higher accuracy, though at a greater computational cost.

:::warning
**Example**
Consider a second order linear ordinary differential equation with boundary conditions of the form
$$
y'' = r(t) \ \ , \ \ y(t_0) = \alpha \ , \ \ y(t_f) = \beta
$$

1. Choose $N$ and let $h = (t_f - t_0)/(N+1)$ and define $t_k = t_0 + kh$. Let $y_k$ denote an approximation of $y(t_k)$.
2. Note that the boundary conditions give us $y_0 = \alpha$ and $y_{N+1} = \beta$ and let $$\mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{bmatrix}$$
3. Let $r_k = r(t_k)$ and substitute the central difference formula at $t_k$ into the differential equation $$\frac{y_{k+1} -2y_k + y_{k-1}}{h^2} = r_k$$
4. Therefore we have $N$ equations and $N$ unknowns $y_k$ for $k=1,\dots,N$. Use the boundary conditions $y_0 = \alpha$ and $y_{N+1} = \beta$ and rearrange the equations in matrix form $A \mathbf{y} = \mathbf{b}$ where
$$
A =
\left[ \begin{array}{rrcrr}
-2 & 1 & & & \\
1 & -2 & 1 & & \\
& & \ddots & & \\
& & 1 & -2 & 1 \\
& & & 1 & -2
\end{array} \right]
\quad
\boldsymbol{b} =
\begin{bmatrix}
h^2 r_1 - \alpha \\ h^2 r_2 \\ \vdots \\ h^2 r_{N-1} \\ h^2 r_N - \beta
\end{bmatrix}
$$
:::

:::warning
**Example** Setup a linear system $A \boldsymbol{y} = \boldsymbol{b}$ for the equation with boundary conditions
$$
y'' = -2 \ \ , \ \ y(0) = 0 \ , \ \ y(1) = 0
$$
using step size $h=0.2$.

1. The step size $h$ corresponds to $N=4$ in our formulation, and $r(t) = -2$, $\alpha = \beta = 0$ therefore
$$
\left[ \begin{array}{rrrr}
-2 & 1 & & \\
1 & -2 & 1 & \\
& 1 & -2 & 1 \\
& & 1 & -2
\end{array} \right]
\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}
=\begin{bmatrix} -0.08 \\ -0.08 \\ -0.08 \\ -0.08 \end{bmatrix}
$$
2. Solve the system to find
$$
\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}
=\begin{bmatrix} 0.16 \\ 0.24 \\ 0.24 \\ 0.16 \end{bmatrix}
$$
3. The equation is very simple and we can solve exactly by integrating twice $$y(t) = -t^2 + C_1t + C_2$$ The boundary conditions imply $C_1 = 1$ and $C_2 = 0$ and therefore the exact solution is $$y(t) = t - t^2$$
:::

:::warning
**Example** Setup a linear system $A \mathbf{y} = \mathbf{b}$ for the equation with boundary conditions $$y'' = \cos(t), \quad y(0) = 0, \quad y(2\pi) = 1$$ using 7 equally spaced points from $t_0 = 0$ to $t_f = 2\pi$.

1. The value $N=5$ corresponds to 7 equally spaced points in our formulation with step size
$$
h = \frac{t_f - t_0}{N + 1} = \frac{2\pi - 0}{5 + 1} = \frac{\pi}{3}
$$
2. We have $r(t) = \cos(t)$ and $\alpha = 0$ and $\beta = 1$. Note that $r_k = \cos(k\pi/3)$ therefore
$$
\left[ \begin{array}{rrrrr}
-2 & 1 & & & \\
1 & -2 & 1 & & \\
& 1 & -2 & 1 & \\
& & 1 & -2 & 1 \\
& & & 1 & -2
\end{array} \right]
\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{bmatrix}
=
\begin{bmatrix} \pi^2/18 \\ -\pi^2/18 \\ -\pi^2/9 \\ -\pi^2/18 \\ \pi^2/18 - 1 \end{bmatrix}
$$
3. Use `scipy.linalg.solve` to compute the solution. The equation is elementary and we can solve exactly by integrating twice $$y(t) = 1 - \cos(t) + \frac{t}{2 \pi}$$
4. Plot the exact solution together with our approximation
:::

:::success

1. Determine whether the statement is **True** or **False**.
* The finite difference method applied to a linear second order differential equation with boundary conditions will compute exact values of the solution $y_k = y(t_k)$ if the step size $h$ is chosen to be small enough.
* The finite difference method applied to a linear second order differential equation with boundary conditions will never compute exact values of the solution $y_k = y(t_k)$ for any differential equation and step size $h$.
2. Setup a linear system $A \mathbf{y} = \mathbf{b}$ to approximate the solution of the equation with boundary conditions $$y'' = 2^t, \quad y'(0) = 1, \quad y(1) = 0$$ using step size $h=0.25$. Use the forward difference formula and the boundary condition $y'(0)=1$ to approximate the boundary value $y_0$.
3. Explain why it is not possible to derive a linear system $A \mathbf{y} = \mathbf{b}$ for the equation $$y'' = \cos(y), \quad y(0) = 0, \quad y(1) = \frac{\pi}{4}$$ by applying finite difference formulas.
4. Setup a linear system $A \mathbf{y} = \mathbf{b}$ for the finite difference approximation of $$y'' + ty = 0, \quad y(1) = 1, \quad y(3) = -1$$ using 5 equally spaced points from $t_0 = 1$ to $t_4 = 3$.
5. Setup the linear system $A \mathbf{y} = \mathbf{b}$ corresponding to the finite difference method applied to the equation $$y'' + y' = t^2, \quad y(-1) = y(1) = 0$$ using 9 equally spaced points on the domain $[-1,1]$.
:::