Characteristic Polynomials and Root Condition

6.2. Characteristic Polynomials and Root Condition#

It is possible to derive any number of multistep methods of various order, but the most important factor related to each method is its stability property.

Definition 6.1 (Stability)

Stability of a numerical method is its capability to keep errors at bay during computation. These errors are required to be dampened out eventually or remain bounded without contaminating the numerical solution significantly.

For simplicity, we study the stability by considering a general linear multistep method,

(6.1)#\[ \sum_{i=0}^{k}\alpha_i y_{j+i}~=~ h\sum_{i=0}^{k}\beta_i f_{j+i}, \qquad \text{with}~ \alpha_k =1 \]

and a test differential equation,

(6.2)#\[ \frac{dy}{dx}= f(x,y) = \lambda y, \quad y(x_0) \,=\,y_0 \]

where \(\lambda\) is a constant.

Note

The general solution to the test problem (6.2) is

\[ \begin{aligned} && y \,&=\, A\,e^{\lambda x}\,, & y(x_0) &= y_0 ~~\Rightarrow~~ A = y_0\,e^{-\lambda x_0} && \end{aligned} \]

thus, the exact solution is: \(\displaystyle{y_j ~=~ y_0\,e^{\lambda(x_j-x_0)}.}\)

According to the test problem (6.2), we can substitute \( f_j = \lambda y_j \) into equation (6.1) and get

(6.3)#\[ \sum_{i=0}^{k}\alpha_i y_{j+i}~=~ h\lambda\sum_{i=0}^{k}\beta_i y_{j+i}~, \]

which can be re-arranged as

(6.4)#\[ \sum_{i=0}^{k}\gamma_i y_{j+i}=0~, ~\text{with}~ \gamma_i =\alpha_i - h\lambda \beta_i \]

We can write equation (6.4) in the operator form as

\[ \L (\E) y_j = \left(\gamma_0 \E^0 + \gamma_1 \E^1 + \gamma_2 \E^2 + \cdots + \gamma_k \E^k \right) y_j =0~\, \]

and get its characteristic equation

(6.5)#\[ \L (z, h\lambda)= \gamma_0 + \gamma_1 z + \gamma_2 z^2 + \cdots + \gamma_k z^k = 0 \]

This is a polynomial equation of degree \(k\) and has \(k\) roots, \(\,z_1, z_2, \dots, z_k\,\), which may be real or complex. If these roots are distinct, the general solution of the difference equation (6.4) is

(6.6)#\[ \begin{aligned} y_j ~=~ A_1\, z_1^j + A_2\, z_2^j + \,\dots\, + A_k\,z_k^j \end{aligned} \]

where \(A_1\), \(A_2\), …, \(A_k\) are constant. See Section 2.4 for more details on difference and characteristic equations.

Question

What kind of conditions shall we impose to keep the general solution (6.6) bounded?

Definition 6.2 (Characteristic Polynomials)

For a general linear multistep method

\[ \sum_{i=0}^k \alpha_i y_{j+i} = h \sum_{i=0}^k f_{j+i}, \]

the expression in the form

(6.7)#\[ \L (z, h\lambda) = \sum_{i=0}^k \left(\alpha_i - h\lambda \beta_i\right) z^i \]

is called its characteristic polynomial.

Note

In some books, \(\pi(z, h\lambda)\) is used to represent the characteristic polynomial.

We call

(6.8)#\[ \rho(z) = \sum_{i=0}^{k} \alpha_i z^i \]

the first characteristic polynomial, and

(6.9)#\[ \sigma(z) = \sum_{i=0}^{k} \beta_i z^i \]

the second characteristic polynomial of the multistep method.

These polynomails given in equations (6.7), (6.8) and (6.9) have the following relation

(6.10)#\[ \L (z, h\lambda) = \rho(z) - h\lambda \sigma(z) \]

Example 6.2

Find the characteristic polynomial \(\L (z, h\lambda)\) for a 2-step method

\[ \alpha_0 y_{j-2} + \alpha_1 y_{j-1} + \alpha_2 y_{j} = h(\beta_0f_{j-2 } + \beta_1 f_{j-1} + \beta_2 f_{j}) \]

Solution

  • Method 1:

    Substituting \(\,f = \lambda\,y\,\) into the 2-step method formula, we have

    \[ \alpha_0 y_{j-2} + \alpha_1 y_{j-1} + \alpha_2 y_{j} = h\lambda(\beta_0 y_{j-2} + \beta_1 y_{j-1} + \beta_2 y_{j}) \]

    and re-arrange it as

    \[ (\alpha_0-h\lambda \beta_0) y_{j-2} + (\alpha_1-h\lambda \beta_1) y_{j-1} + (\alpha_2-h\lambda \beta_2) y_{j} = 0 \]

    Writing the equation in shift operator form

    \[ \left[ (\alpha_0-h\lambda \beta_0) \E^0 + (\alpha_1-h\lambda \beta_1) \E^1 + (\alpha_2-h\lambda \beta_2) \E^2 \right] y_{j-2} =0 \]

    so the characteristic polynomial is

    \[\begin{split} \begin{aligned} \L (z, h\lambda) &=\sum_{i=0}^{2} (\alpha_i - h\lambda \beta_i) z^i\\ & = (\alpha_0-h\lambda \beta_0) z^0 + (\alpha_1-h\lambda \beta_1) z^1 + (\alpha_2-h\lambda \beta_2) z^2 \end{aligned} \end{split}\]

  • Method 2:

    Using (6.7) given in Definition 6.2, we can quickly work out the characteristic polynomial for this 2-step method

    \[\begin{split} \begin{aligned} \L (z, h\lambda) &=\sum_{i=0}^{2} (\alpha_i - h\lambda \beta_i) z^i\\ & = (\alpha_0-h\lambda \beta_0) z^0 + (\alpha_1-h\lambda \beta_1) z^1 + (\alpha_2-h\lambda \beta_2) z^2 \end{aligned} \end{split}\]

Definition 6.3 (Root Condition)

Let \(\,z_1,z_2,\dots,z_k\,\) denote the roots of the first characteristic polynomial \(\,\rho(z)\,\) associated with a linear multistep method. If \(\,|z_i|\leq1\,\), for each \(\,i=1,2,...,k\,\), and all roots with absolute value \(1\) are simple roots (i.e. not repeated), then the difference method is said to satisfy the root condition.

Figure made with TikZ

Fig. 6.1 Unit circle in the complex plane \(\C\).

Theorem 6.1 (Stability)

A linear multistep method of order \(\,\geq1\,\) is stable if the roots of its first characteristic polynomial \(\,\rho(z)\,\) satisfies the root condition, i.e. the roots of \(\,\rho(z)\,\) lie inside or on the unit circle, and the roots on the unit circle are distinct (or simple) - this is called the root condition.

Definition 6.4 (Strong Stability)

A method is strongly stable if the roots of \(\rho(z) = 0\) are inside the unit circle except for only one root, call principal root, with the magnitude of one(i.e. \(|z| = 1\)).

Definition 6.5 (Weak Stability)

A method is weakly stable if it is stable (i.e. satisfies the root condition) but has more than one root on the unit circle.

Example 6.3

Find the roots associated with the following characteristic polynomials:

  1. \(\rho(z)= z^3 + z \)

    Solution (click to show)

    Let \(z^3 + z ~=~ 0\)
    \(\therefore\quad z(z^2 + 1) ~=~ 0\)
    \(\therefore\quad z_1 ~=~ 0, ~~ z_2 ~=~ i, ~~z_3 ~= -i\)

    We can also use Matlab to find the roots of the polynomial

    p=[1 0 1 0];     %vector defining the coefficients of the polynomial
roots(p)         %solve the polynomial equation

    Output

    ans =

    0.0000 + 0.0000i
    0.0000 + 1.0000i
    0.0000 - 1.0000i

    \(z_1\), \(z_2\), and \(z_3\) are shown on the unit circle. Because there are more than one root on the unit circle, the method is weakly stable.

    Figure made with TikZ

  2. \(\rho(z)=z^4 + z\)

    Solution (click to show)

    Let \(z^4 + z ~=~ 0\)

    \[\begin{split} \begin{aligned} z^4 + z & = z(z^3 + 1) \\ & = z (z+1) (z^2-z+1) \\ \end{aligned} \end{split}\]

    \(\therefore\quad z_1 ~=~ 0\), \(z_2 = -1\),  \(z_3 = \dfrac{1}{2} + i\dfrac{\sqrt{3}}{2}\),  \(z_4 = \dfrac{1}{2} - i\dfrac{\sqrt{3}}{2}\).

    We can also use Matlab to find the roots of the polynomial

    syms z
solve(z^4+z)

    Output

    ans =
                -1
                0
1/2 - (3^(1/2)*1i)/2
(3^(1/2)*1i)/2 + 1/2

    Because there are more than one root on the unit circle, the method is weakly stable.

    Figure made with TikZ