Francais | English | Espanõl

Stiff equation

From Wikipedia, the free encyclopedia

Jump to: navigation, search

In mathematics, a stiff equation is a differential equation for which certain numerical methods for solving the equation are numerically unstable, unless the step size is taken to be extremely small. It has proved difficult to formulate a precise definition of stiffness, but the main idea is that the equation includes some terms that can lead to rapid variation in the solution.

Contents

[edit] Motivating example

Consider the initial value problem

<math>\,y'(t)=-15y(t),\quad t \ge 0, y(0)=1.</math>

The exact solution is <math>y(t)=e^{-15t}</math>, and clearly <math>y(t)</math> → <math>0</math> as <math>t</math> → ∞. We want the numerical solutions to exhibit the same behaviour.

Attempting to solve it using Euler's method, <math>y_{n+1}=y_n+hf(t_n,y_n)=y_n+h(-15y_n)</math> with a step size <math>h=1/4</math>

<math>t_0=0</math> <math>y_0=1</math>
<math>t_1=0.25</math> <math>y_1=1+1/4(-15)=-2.75</math>
<math>t_2=0.5</math> <math>y_2=-2.75+1/4(-15)(-2.75)=7.5625</math>
<math>t_3=0.75</math> <math>y_3=7.5625+1/4(-15)(7.5625)=-20.7969.</math>

It is obvious that for this choice of step size the desired limiting behaviour is not present. If we reattempt, halving the step size to give <math>h=1/8</math>:

<math>t_0=0</math> <math>y_0=1</math>
<math>t_1=0.125</math> <math>y_1=1+1/8(-15)=-0.875</math>
<math>t_2=0.25</math> <math>y_2=-0.875+1/8(-15)(-0.875)=0.765625</math>
<math>t_3=0.375</math> <math>y_3=0.765625+1/8(-15)(0.765625)=-0.669922</math>
<math>t_4=0.5</math> <math>y_4=-0.669922+1/8(-15)(-0.669922)=0.586182</math>
<math>t_5=0.625</math> <math>y_5=0.586182+1/8(-15)(0.586182)=-0.512909</math>
<math>t_6=0.75</math> <math>y_6=-0.512909+1/8(-15)(-0.512909)=0.448795.</math>

The oscillatory behaviour is still present.

Still using the same problem, we reattempt with the trapezoidal method, that is, the two-stage Adams–Moulton multistep method,

<math>y_{n+1}=y_n+{1\over 2}h\left(f(t_n,y_n)+f(t_{n+1},y_{n+1})\right)</math>

for this problem,

<math>y_{n+1}=y_n+{1\over 2}h\left((-15y_n)+(-15y_{n+1})\right)=y_n-{15\over 2}h(y_n+y_{n+1})</math>
<math>y_{n+1}=y_n-{15\over 2}hy_n-{15\over 2}hy_{n+1}</math>
<math>y_{n+1}+{15\over 2}hy_{n+1}=y_n-{15\over 2}hy_n</math>
<math>\left(1+{15h\over 2}\right)y_{n+1}=\left(1-{15h\over 2}\right)y_n</math>
<math>y_{n+1}={1-{15h\over 2}\over 1+{15h\over 2}}y_n</math>

and for <math>h=1/8</math>,

<math>y_{n+1}=1/31 y_n</math>
<math>t_0=0</math> <math>y_0=1</math>
<math>t_1=0.125</math> <math>y_1=1/31=0.0322581</math>
<math>t_2=0.25</math> <math>y_2=1/961=0.00104058</math>
<math>t_3=0.375</math> <math>y_3=1/29791=0.0000335672.</math>

The trapezoidal method clearly gives us the desired results. What differentiates methods like the trapezoidal method from methods like the Euler method?

[edit] Test equation and A-stability

Certain types of problems can be characterized as stiff:

problems of the form
<math>y' = ky + f(t),\,</math>
where <math>k</math>∈C and Re(<math>k</math>) is large and negative.
systems of the form
<math>\mathbf{y}' = K\mathbf{y}+\mathbf{f}(t),</math>
where <math>K</math> is a square matrix having at least one eigenvalue <math>l</math>∈C and Re(<math>l</math>) is large and negative.
systems of the form
<math>\mathbf{y}' = \mathbf{f}(\mathbf{y}, t),</math>
with the Jacobian of f having at least one eigenvalue <math>m</math>∈C and Re(<math>m</math>) is large and negative.

For problems of the first form, it is clear that we can, if we attempt to solve a problem of the form <math>y' = ky</math> with some method, determine the behaviour of a method. For a one-step method, applying the method to the test equation gives <math>y_{n+1}=\phi(hk)y_n</math>, and by induction, <math>y_n=\phi(hk)^ny_0</math>. Now the required behaviour is that <math>y_n</math> → 0 as <math>n</math> → ∞, so we require <math>|\phi(hk)|<1</math>.

We say the set {<math>z</math> ∈ C | <math>|\phi(z)|<1</math> } with <math>z=hk</math> is the region of absolute stability. If the region of absolute stability contains the set { <math>z</math> ∈ C | Re(<math>z</math>) < 0 }, that is, the left half plane, the method is said to be A-stable. A-stable methods do not exhibit the instability problems as described in the motivating example.

[edit] Example

Consider both the Euler and trapezoidal methods above. Applying the test equation:

<math>y_{n+1} = y_n + hf(t_n, y_n) = y_n + h(ky_n) = y_n + hky_n = (1+hk)y_n</math>
<math>y_n = (1+hk)^ny_0</math>
<math>\phi(z) = 1+z.</math>

The region of absolute stability for this method is { <math>z</math> ∈ C | |1+z| < 1 } which is a circle, thus the Euler method is not A-stable.

However, for the trapezoidal method:

<math>y_{n+1}=y_n+{1\over 2}h\left(f(t_n,y_n)+f(t_{n+1},y_{n+1})\right),</math>

applying the test equation

<math>y_{n+1}=y_n+{1\over 2}h\left(ky_n+ky_{n+1})\right)</math>
<math>y_{n+1}-{1\over 2}hky_{n+1}=y_n+{1\over 2}hky_n</math>
<math>\left( 1-{1\over 2}hk \right)y_{n+1}=\left(1+{1\over 2}hk \right) y_n</math>
<math>y_{n+1}={1+{1\over 2}hk \over 1-{1\over 2}hk}y_n</math>
<math>y_n=\left( {1-{1\over 2}hk \over 1-{1\over 2}hk} \right) ^ny_0</math>
<math>\phi(z)={1+{1\over 2}z \over 1-{1\over 2}z}.</math>

The region of absolute stability then is

<math>\left\{ z \in \Bbb{C} \left|\ \left| {1+{1\over 2}z \over 1-{1\over 2}z} \right| < 1 \right.\right\},</math>

which is indeed identical to the left half complex plane.

[edit] Multistep methods

For a multistep method of the form

<math>y_{n+1}=\sum_{i=0}^s a_i y_{n-i}+h\sum_{j=-1}^s b_j f(t_{n-j},y_{n-j})</math>

the test equation is applied,

<math>y_{n+1}=\sum_{i=0}^s a_i y_{n-i}+h\sum_{j=-1}^s b_jky_{n-j}=\sum_{i=0}^s a_i y_{n-i}+hk\sum_{j=-1}^s b_jy_{n-j}</math>

which can be simplified to

<math>(1-b_{-1}hk)y_{n+1}-\sum_{j=0}^s (a_j+b_jhk)y_{n-j}=0.</math>

Making the substitutions <math>z=hk</math> and <math>y_n=w^n</math> and simplifying further,

<math>w^{r+1}-\sum_{i=0}^s a_iw^{s-j} - z\sum_{j=-1}^s b_jw^{r-j}=0.</math>

But the first sum corresponds to the characteristic equation for the method. Define

<math>\Phi(z, w) = w^{r+1}-\sum_{i=0}^s a_iw^{s-j} - z\sum_{j=-1}^s b_jw^{r-j}.</math>

It is a fact that if the roots of <math>\Phi(z,w)</math> in <math>w</math> have modulus less than 1 for all roots, then the method has the required limiting behaviour.

The region of absolute stability for a multistep method of the above form is then

<math>\left\{ z \in \Bbb{C} \left|\ |w| < 1, \mathrm{for\ all}\ w\ \mathrm{such\ that\ } \Phi(z,w)=0\right.\right\}.</math>

Again, if this set contains the left half complex plane, the multi-step method is said to be A-stable.

[edit] Example

Let us determine the region of absolute stability for the two step Adams-Bashforth method

<math>y_{n+1} = y_n + h \left( {3\over 2}f(t_n, y_n)-{1\over 2}f(t_{n-1}, y_{n-1}) \right) .</math>

Applying the test equation

<math>y_{n+1} = y_n + h \left( {3 \over 2}ky_n - {1\over 2}ky_{n-1} \right) </math>
<math>=y_n + {3\over 2}hky_n - {1\over 2}hky_{n-1}</math>
<math>y_{n+1} - (1+{3\over 2}hk)y_n + {1\over 2}hky_{n-1}=0</math>
<math>w^2-(1+{3\over 2}z)w+{z\over 2}=0</math>

which has roots

<math>w={1\over 2}(1 + {3\over 2}z \pm\sqrt{1 + z + {9\over 4}z^2}</math>

thus the region of absolute stability is

<math>\left\{ z \in \Bbb{C} \left|\ \left|{1\over 2}(1 + {3\over 2}z \pm\sqrt{1 + z + {9\over 4}z^2}\right| < 1 \right.\right\}.</math>
Retrieved from "http://en.wikivisual.com/index.php/Stiff_equation"
Personal tools