# Numerical stability and ODEs

### General solution

Let’s assume we are trying to solve the following differential equation:

$\displaystyle U\frac{dC}{dx} = D\frac{d^2 C}{dx^2} - K \cdot C$.

As it’s a simple, linear ODE, it can be solved by using as ansatz

$\displaystyle C(x) = e^{\lambda\,x}$.

If we replace $C(x)$ by it, we get

$\displaystyle U\,\lambda\,e^{\lambda\,x} = D\,\lambda^2\,e^{\lambda\,x} - K\,e^{\lambda x}$.

As $e^{\lambda\,x}$ is nonzero, lambda must be a root of the characteristic polynomial,

$D\,\lambda^2 - U\,\lambda - K = 0$.

Using the well-known “quadratic formula”,

$\displaystyle \lambda = \frac{U \pm \sqrt{U^2 + 4\,D\,K}}{2D}$.

If we call the two values of $\lambda_1$ and $\lambda_2$, assuming they are different, the general solution will be a linear combination:

$\displaystyle C(x) = A_1\,e^{\lambda_1\,x} + A_2\,e^{\lambda_2\,x}$.

### Specific example

If we take the following specific values for the parameters,

$\displaystyle U = 36\,{\rm km \cdot h^{-1}}$

$\displaystyle D = 0.14$ (assuming appropriate units)

$\displaystyle K = 10$ (assuming appropriate units),

we get the following values for $\lambda$:

$\displaystyle \lambda = \frac{36 \pm \sqrt{36^2 + 4 \cdot 0.14 \cdot 10 }}{2 \cdot 0.14}$

$\displaystyle = \frac{36 \pm \sqrt{1300}}{0.28}$.

$\displaystyle \lambda_1 = 257.34$

$\displaystyle \lambda_2 = -0.19825$

Now we have to use the given boundary conditions:

$\displaystyle C(0) = 1000000$

$\displaystyle C'(0) = -1000$

By replacing:

$\displaystyle A_1 + A_2 = 1000000$

$\displaystyle \lambda_1 A_1 + \lambda_2 A_2 = -1000$

Substituting:

$\displaystyle \lambda_1 A_1 + \lambda_2 (1000000 - A_1) = -1000$

$\displaystyle (\lambda_1 - \lambda_2) A_1 + 1000000 \lambda_2 = -1000$

$\displaystyle A_1 = \frac{-1000 - 1000000 \lambda_2}{\lambda_1 - \lambda_2}$

$\displaystyle A_1 = \frac{-1000 - 1000000 (-0.19825)}{257.34 + 0.19825}$

$\displaystyle A_1 = 765.9$

$\displaystyle A_2 = 999234$

Then the full solution will be:

$\displaystyle C(x) = 765.9\,e^{257.34\,x} + 999234\,e^{-0.19825\,x}$.

This solution will quickly diverge for values of $x$ not much smaller than one.

Even if the initial conditions were finely tuned for getting $A_1 = 0$, any small numerical error would give a quickly diverging solution.

Advertisements