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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s