# Ordinary differential equations

Ordinary linear differential equation is a equation in form of:

$$\sum_n^M a_n \frac{d^n}{dx^n} y(x) =0$$

where $M$ is the order of differential equation. There are few methods that can be used to solve ODE. In this text I will show how to solve the differential equation analytically with characteristic equation and how to solve it with the numerical methods.

## The analytical method

Let assume that we want to solve a first order ODE:

$$a \frac{dy(x)}{dx} + b y(x)=0. \hspace{1cm} (1)$$

We can “guess” the solution by choosing a function in form of $y(x) = e^{rx}$ where $r$ is yet unknown. After plugging our trial solution into the equation we get

$$a r e^{rx} + be^{rx}=0 .$$

If we assume that for all $x$ our guess function $y(x) \ne 0$ then we can divide the above equation by $y(x) = e^{rx}$ to get:

$$a r + b=0$$

which has one solution

$$r = -\frac{b}{a}.$$

The above results can be obtained in different way. Instead of plugging the test function we divide (1) by the $a$ and move the second term to the right hand side

$$\frac{dy(x)}{dx}=-\frac{b}{a}y(x)$$

Next we “multiply” by $dx$, divide by $y(x)$ and integrate over $dx$. Again the we take the assumption that $y(x) \ne 0$.

$$\int \frac{dy(x)}{y(x)}dx=-\int\frac{b}{a} dx .$$

Both integrals are trivial and the result is

$$ln~y(x)=-\frac{b}{a}x + C$$

where $C$ is the integration constant. Next step is to take the left and right hand side to the exponent

$$e^{ln~y(x)}=e^{-\frac{b}{a}x + C}$$

to get

$$y(x)=e^{-\frac{b}{a} + C}=e^{-\frac{b}{a}x}e^C$$

which is of course

$$y(x)=ce^{-\frac{b}{a}x}$$

where $c=e^C$ is a constant.

## The numerical method

In order to solve (1) we need to approximate somehow the derivative and the function $y(x)$ in equation (1). First we need to expresses this function in terms of discrete points in space instead of continuous function. We do that by defining the points $x_i$ separated by distance of $\Delta x$, where $i$ numerates the the points on the number line. The idea is presented below. The $x_{i+1} = x_{i}+\Delta x$.

We approximate
$$\frac{dy(x)}{dx} \approx \frac{y(x_i+\Delta x)-y(x_i)}{\Delta x}.$$
We can put the above approximation to the (1) to get
$$\frac{y(x_i+\Delta x)-y(x_i)}{\Delta x}=-\frac{b}{a}y(x_i).\hspace{1cm} (2)$$
This can bo solved to get
$$y(x_i+\Delta x)=y(x_i)-\Delta x\frac{b}{a}y(x_i)=\left(1-\Delta x\frac{b}{a}\right)y(x_i) \hspace{1cm} (3)$$

The above formula is known as a Forward Euler Method (link to Wikipedia article).
To solve this we must know the initial value of function $y(x)$. Lets assume that we have $N_x$ points, $y(x_0=0)=1$ and $a=b=1$. Using formula (3) we can get the following table:

 x y 0 1 0.1 0.9 0.2 0.81 0.3 0.729 0.4 0.6561 0.5 0.59049 0.6 0.531441 0.7 0.4782969 0.8 0.43046721 0.9 0.387420489 1 0.3486784401 1.1 0.3138105961 1.2 0.2824295365 1.3 0.2541865828 1.4 0.2287679245 1.5 0.2058911321 1.6 0.1853020189 1.7 0.166771817 1.8 0.1500946353 1.9 0.1350851718 2 0.1215766546 2.1 0.1094189891 2.2 0.0984770902 2.3 0.0886293812 2.4 0.0797664431 2.5 0.0717897988

This table generates a graph that is solution to this equation:

The red line represents the curve $y(x)=e^{-x}$ which is the exact solution to this problem.
It is worth noting that we used a first order formula to solve this equation and its accuracy is proportional to the $\Delta x$.

Below you can find the Python code for the example above:

from numpy import *
import scipy as sp
import matplotlib
import matplotlib.pyplot as plt

y0 = 1.
xmin = 0.
xmax = 3.
npoints = 16

dx = (xmax-xmin)/float(npoints)

x = sp.linspace(xmin,xmax,npoints)
sol = sp.zeros(npoints)
sol[0] = y0

for i in xrange(1,npoints):
sol[i] = sol[i-1]-dx*sol[i-1]

sol_exact = exp(-x)

plt.figure()
plt.plot(x, sol, 'ok')
plt.plot(x, sol_exact, '--k')
plt.ylabel("forward")
plt.show()

The problem can be solved in different way by changing the right hand side of (2)

$$\frac{y(x_i+\Delta x)-y(x_i)}{\Delta x}=-\frac{b}{a}y(x_{i+1}).$$