logo

Stability and Root Conditions of Consistent Multistep Methods 📂Numerical Analysis

Stability and Root Conditions of Consistent Multistep Methods

Theorem

If a multistep method is consistent, the method is stable     \iff the method satisfies the root condition.

Explanation

When creating node points by cutting the closed interval [x0,b][x_{0} , b] at intervals of hh, let’s call it x0x1xN(h)1xN(h)bx_{0} \le x_{1} \le \cdots \le x_{N(h) -1} \le x_{N(h) } \le b. Here, N(h)N(h) represents the index of the last node point that changes according to hh.

Consider z0,,zpz_{0} , \cdots , z_{p}, a very slight change to the original given initial value y0,,ypy_{0} , \cdots , y_{p}. That the method is stable means for sufficiently small positive numbers hh and ϵ\epsilon, if it is said that max0npynznϵ\displaystyle \max_{0 \le n \le p} | y_{ n} - z_{n} | \le \epsilon, then there exists some constant CC that satisfies max0nN(h)ynzncϵ\displaystyle \max_{0 \le n \le N(h) } | y_{ n} - z_{n} | \le c \epsilon, independent of hh. If the initial small change can become a big change as hh is adjusted, it is said that the method is not stable.

Proof 1

Consistent Multistep Method: For the initial value problem {y=f(x,y)(y(x0),,y(xp))=(Y0,,Yp)\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ) , \cdots , y(x_{p}) ) = (Y_{0}, \cdots , Y_{p} ) \end{cases}, a multistep method yn+1=j=0pajynj+hj=1pbjf(xnj,ynj) y_{n+1} = \sum_{j=0}^{p} a_{j} y_{n-j} + h \sum_{j = -1}^{p} b_{j} f (x_{n-j} , y_{n-j} ) satisfies the following. {j=0paj=1j=0pjaj+j=1pbj=1 \begin{cases} \displaystyle \sum_{j = 0}^{p} a_{j} = 1 \\ \displaystyle - \sum_{j = 0}^{p} j a_{j} + \sum_{j = -1}^{p} b_{j} = 1 \end{cases}

Root Condition: Let’s say ρ(r)=rp+1j=0pajrpj\displaystyle \rho ( r) = r^{p+1} - \sum_{j=0}^{p} a_{j} r^{p-j} for a multistep method that is consistent. When the roots r0,,rpr_{0} , \cdots , r_{p} of equation ρ(r)=0\rho (r ) = 0 satisfy the following conditions, the given multistep method is said to satisfy the Root Condition.

  • (i): rj1| r_{j} | \le 1
  • (ii): rj=1    ρ(rj)0|r_{j}| = 1 \implies \rho ‘(r_{j}) \ne 0

()(\Rightarrow)

Assume that the method is stable but does not satisfy the root condition.

This assumption will be shown to be contradicted by the initial value problem {y=0y(0)=0\begin{cases} y ' = 0 \\ y(0) = 0 \end{cases}.

The numerical solution to the given problem is trivially yn=0y_{n} = 0 for all n0n \ge 0.

When considering z0,,zpz_{0} , \cdots , z_{p}, a very slight change given to y0,,ypy_{0} , \cdots , y_{p} zn+1=a0zn+a1zn1++ap1znp1+apznp z_{n+1} = a_{0} z_{n} + a_{1} z_{n-1} + \cdots + a_{p-1} z_{n-p-1} + a_{p} z_{n -p} To obtain the characteristic solution, let’s say zn:=rnz_{n} := r^{n} rn+1=a0rn+a1rn1++ap1rnp1+aprnp r^{n+1} = a_{0} r^{n} + a_{1} r^{n-1} + \cdots + a_{p-1} r^{n-p-1} + a_{p} r^{n -p} Dividing both sides by rnpr^{n-p} rp+1=a0rp+a1rn1++ap1r1+ap r^{p+1} = a_{0} r^{p} + a_{1} r^{n-1} + \cdots + a_{p-1} r^{1} + a_{p} This p+1p+1th degree equation has r1,r2,,rpr_{1}, r_{2} , \cdots , r_{p} roots in addition to r0=1r_{0}=1. Therefore, the general solution is zn=c0r0n+c1r1n++cprpn z_{n} = c_{0} r_{0}^{n} + c_{1} r_{1}^{n} + \cdots + c_{p} r_{p}^{n} for some c0,,cpc_{0} , \cdots , c_{p}.

Consider the following two cases for some 0jp0 \le j \le p.

  • Case 1.
    If condition (i) is not satisfied and for j j, it is rj>1|r_{j}| > 1, then z0=0++cjrj0++0=ϵz_{0} = 0 + \cdots + c_{j} r_{j}^{0} + \cdots + 0 = \epsilon, that is, it’s said z0=cjrj0=ϵ z_{0} = c_{j} r_{j}^{0} = \epsilon . In other words, since z0=cj(rj)0    zn=cj(rj)nz_{0} = c_{j} ( r_{j} ) ^{0} \implies z_{ n} = c_{j} ( r_{j} ) ^{n} z0=ϵ,z1=ϵrj,,zp=ϵrjp z_{0} = \epsilon, z_{1} = \epsilon r_{j } , \cdots , z_{p} = \epsilon r_{j}^{p} and max0npynzn=ϵrjp \max_{0 \le n \le p } | y_{n} - z_{n} | = \epsilon | r_{j} |^{p} Applying the method to [x0,b][x_{0} , b] maxx0xnbynzn=ϵrjN(h) \max_{x_{0} \le x_{n} \le b } | y_{n} - z_{n} | = \epsilon |r_{j}|^{N(h)} But when h0h \to 0, N(h)N(h) \to \infty and since rj>1|r_{j}| > 1, maxx0xnbynzn=cϵrjp \max_{x_{0} \le x_{n} \le b } | y_{n} - z_{n} | = c \cdot \epsilon |r_{j}|^{p} there cannot exist any C>0C>0 that satisfies. Therefore, the method becomes unstable.
  • Case 2.
    If condition (ii) is not satisfied and for jj, it is rj=1    ρ(rj)=0|r_{j}| = 1 \implies \rho’(r_{j}) = 0, this means if rj=1|r_{j}| = 1, then rjr_{j} is a repeated root of the characteristic equation, and the multiplicity of rjr_{j} must be at least 22. In this case, the general solution is represented as a linear combination including linearly independent particular solutions rjn,nrjn,,nν1rjnr_{j}^{n} , n r_{j}^{n} , \cdots , n^{ \nu - 1} r_{j}^{n}. Mathematically, zn=cj0rjn+cj1nrjn++cjν1nν1rjn z_{n} = c_{j_{0}} r_{j}^{n} + c_{j_{1}} n r_{j}^{n} + \dots + c_{j_{\nu-1}} n^{\nu-1} r_{j}^{n} Here, by setting ϵ:=max0kν1cjk\displaystyle \epsilon := \max_{ 0 \le k \le \nu-1} | c_{j_{k}} |, since rj=1| r_{j} | = 1 z0ϵ,z12ϵ,,zpϵ(1+p++pν1)=ϵpν1p1 | z_{0 } | \le \epsilon , | z_{1 } | \le 2 \epsilon , \cdots , \displaystyle | z_{ p} | \le \epsilon ( 1 + p + \cdots + p^{\nu -1} ) = \epsilon {{ p^{\nu} -1 } \over { p - 1}} As it’s a multistep method, for p2p \ge 2 max0npynzn=ϵpν1p1 \max_{0 \le n \le p } | y_{n} - z_{n} | = \epsilon {{ p^{\nu} -1 } \over { p - 1}} Applying the method to [x0,b][x_{0} , b] maxx0xnbynzn=ϵN(h)ν1N(h)1 \max_{x_{0} \le x_{n} \le b } | y_{n} - z_{n} | = \epsilon {{ N(h)^{\nu} -1 } \over { N(h) - 1}} But when h0h \to 0, N(h)N(h) \to \infty, therefore maxx0xnbynzn=cϵpν1p1 \max_{x_{0} \le x_{n} \le b } | y_{n} - z_{n} | = c \cdot \epsilon {{ p^{\nu} -1 } \over { p - 1}} there cannot exist any C>0C>0 that satisfies. Therefore, the method becomes unstable.

()(\Leftarrow) The original proof is too difficult and contains many leaps2.

If it is considered en:=ynzne_{n} : = y_{n} - z_{n}, it can be thought of as solving the non-homogeneous linear differential equation y=z+ey = z + e.

Conversely, z=yez = y - e becomes a homogeneous linear differential equation. yn+1=j=0pajynj+hj=1pbjf(xnj,ynj) \begin{align} \displaystyle y_{n+1} = \sum_{j=0}^{p} a_{j} y_{n-j} + h \sum_{j = -1}^{p} b_{j} f (x_{n-j} , y_{n-j} ) \end{align}

zn+1=j=0pajznj+hj=1pbjf(xnj,znj) \begin{align} \displaystyle z_{n+1} = \sum_{j=0}^{p} a_{j} z_{n-j} + h \sum_{j = -1}^{p} b_{j} f (x_{n-j} , z_{n-j} ) \end{align} If (1)(1) is subtracted from (2)(2), en+1=i=0najenj+hj=1pbj[f(xnj,ynj)f(xnj,znj)] e_{n+1} = \sum_{i=0}^{n} a_{j} e_{n-j} + h \sum_{j= -1}^{p} b_{j} \left[ f(x_{n-j} , y_{n- j }) - f(x_{n-j} , z_{n- j }) \right] is obtained. If it is said that yn=d0s0n+d1s1n++dpspny_{n} = d_{0} s_{0}^{n} + d_{1} s_{1}^{n} + \cdots + d_{p} s_{p}^{n}, en=g0(r0+s0)n+g1(r1+s1)n++gp(rp+sp)n e_{n} = g_{0} ( r_{0} + s_{0} )^n + g_{1} ( r_{1} + s_{1} )^n + \cdots + g_{p} ( r_{p} + s_{p} )^n Let’s consider max0npynzn=ϵ\displaystyle \max_{0 \le n \le p} | y_{ n} - z_{n} | = \epsilon for this. The given method en=g0(r0+s0)n+g1(r1+s1)n++gp(rp+sp)n e_{n} = g_{0} ( r_{0} + s_{0} )^n + g_{1} ( r_{1} + s_{1} )^n + \cdots + g_{p} ( r_{p} + s_{p} )^n also satisfies the root condition for (r0+s0),,(rp+sp)( r_{0} + s_{0} ) , \cdots , ( r_{p} + s_{p} ), en=g0r0+s0n+g1r1+s1n++gprp+spnpmax1ipgimax1ipri+sinpϵ1 \begin{align*} \displaystyle | e_{n} | =& | g_{0} | | r_{0} + s_{0} |^n + | g_{1} | | r_{1} + s_{1} | ^n + \cdots + | g_{p} | | r_{p} + s_{p} |^n \\ \le & p \max_{1 \le i \le p} | g_{i}| \max_{1 \le i \le p} | r_{i} + s_{i} |^n \\ \le & p \cdot \epsilon \cdot 1 \end{align*} Therefore, maxx0xnben=pϵ\displaystyle \max_{ x_{0} \le x_{n} \le b} | e_{n} | = p \epsilon holds, and the method is stable.


  1. Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p398~401. ↩︎

  2. Isaacson. (2012). ANALYSIS OF NUMERICAL METHODS: p405~417. https://www.researchgate.net/file.PostFileLoader.html?id=56c583ac5e9d97577f8b458e&assetKey=AS:330399868833792@1455784875579 ↩︎