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 = 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}).$$