数理解析研究所講究録
第 643 巻 1988 年 203-225

203

A Series of Collocation Runge-Kutta Methods
Taketomo MITSUI and Hiroshi SUGIURA
三井斌友

杉浦

洋

Dept. of Information Eng., Fac. Eng., Nagoya Univ., Nagoya, JAPAN
名古屋大学

・

工学部・ 情報工学教室

1. lntroduction
In solving numerical initial value problem of ordinary differential equation(1.1)$\frac{dy}{dx}=f(x,y)(a<x<b)$ ,$y(a)=y_{0}$,the class of Runge-Kutta formulas (RK formulas, in short) is most important among many discrete variable methods. Especially the stiff problemof (1.1) requires highly sophisticated methods to solve, one of which is theimplicit RK formula.
$Ithas/$ the
form
$($

1.2

$a)$

$y_{n+1}=y_{n}+h \sum_{:=1}^{l}b:f(x_{n}+ch, Y)$

where$(1.2b)$$Y_{:}=y_{n}+h\sum_{j=1}^{l}a:jf(x_{n}+c_{j}h,Y_{j})$$(i=1, \ldots, s)$.

Here, the interval of integration, $[a, b]$ , is divided by the step-size
the step-points$x_{0}(=a)<x_{1}<\ldots<x_{n}<x_{n+1}<...$ $<x_{N}(=b)$

are given by(1.3)$x_{n}=x_{0}+nh$

1

.

$h$

so that

204

The real parameters $a_{ij},bc$ specify the method. Formula (1.2), which iscalled an s-stage one, is usually assumed to satisfy(1.4)$c_{i}=\sum_{j=1}^{l}a_{1j}$$(i=1, \ldots,s)$.

This condition can be motivated by the requirement that the approximate
solution $y_{n}(n=1,2, \ldots)$ has at least the order one for the exact value$y(x_{0}+nh)$ with respect to .
To specify the formula, i.e. to determine the parameters $ab$ and ,the collocation method seems to be most simple and powerful. Although
the exact definition will be given in the next section, the sense by the word”collocation” could be explained by considering the polynomial which interpolates the numerical solutionin (1.2b) at the collocation points
and$x_{n}+c:h(i=1, \ldots,s)$ . In the collocation method one will first choose thedinstinct abscissas, then the other parameters will be uniquelydetermined. Many known RK formulas, such as Gauss-Legendre, Radauand Lobatto types, belong to this class([2]).
The aim of the present paper is to analyse the process of the cho $os-$
ing the abscissas so as to increase the order of the approximate solution, to investigate the stability of the derived collocation methods, and toreconsider a method proposed by a group of physists ([1]), who called it
HIDM ( $=Higher$ order Implicit Difference Method).
$h$$c_{j}$$Y_{1}$$c_{1},$ $c_{2},$$\ldots,$$x_{n}$$c_{\iota}$$y_{n+1}$$-$

2. The Construction of CoIIocation RK Method
2.1 CoIIocation Method of Lower Order
As usual, the local discretization error $r(x., h)$ of the RK method (1.2)
is defined by the following:(2.1)$\tau(x_{n},h)$$=$$\frac{1}{h}\{y(x_{n+1})-\hat{y}_{n+1}\}$$=$$\frac{1}{h}\{y(x_{n+1})-y(x_{n})-h\sum_{:}b:f(ch,\hat{Y}_{i})\}$

where$\hat{Y}_{:}=y(x_{n})+h\sum_{j}:jch,\hat{Y}_{j})$2$(i=1, \ldots, s)$

.

,

205

Let us call the RK method to be consistent of order , if
positive integer such that
$p$

(2.2)

$\tau(x_{n}, h)=O(h^{p})$

as

$p$

is the largest$h\downarrow 0$

for any sufficiently differentiable solution $y(x)$ of (1.1). Butcher [4] has
introduced the simplifying conditions to easily determine the order fromvery complicated algebraic relations of the parameters.
Deflnition 1. An s-stage $RK$ formula is said to satisfy simplifying condition(2.3)$A(p)$

if $\tau(x_{n},h)=O(h^{p})$ ,(2.4)$B(p)$

if$\sum_{1=1}^{l}b_{i}c_{:}^{k-1}=\frac{1}{k}(k=1,2, \ldots,p)$(2.5)$C(p)$

if$\sum_{j=1}^{l}a:jc_{j}^{k-1}=\frac{1}{k}c_{:}^{k}(i=1, \ldots , s,k=1,2, \ldots ,p)$

(2.6)

$D(q)$

if$\sum_{:=1}^{l}bc^{\ell-1}a=\frac{1}{p}b_{j}$$(2.7)/$$E(p,q)$

if

(

$1$

一

$c_{j}^{\ell}$,)($j=1,$$\ldots$,, 馬 $l=1$ , 名...

$\sum_{i,j=1}^{l}bc^{\ell-1}ac_{j}^{k-1}=\frac{1}{k(\ell+k)}(\ell=1, \ldots,p, k=1, .

$q$

),

.., q)$

Definition 2 (Burrage,3). If an s-stage $RK$ formula satisfies the conditions $B(s)$ and $C(s)$ , then it is called a collocation method.
The implication of the collocation RK method can be considred as follows.
For the sake of notational convenience, we introduce a scaling for the
independent variable by $x=x_{0}+th$ . Then we will consider the differential equation and its approximate solution in the term of the variable. Moreover, for the RK methods, it suffices to consider on the interval$t\in[0,1]$ . Hereafter, within the present section, we will restrict ourselves on
this situation.
If the RK method (1.2) is regarded as a method of numerical quadrature,
we may apply the method to a differential equation depending only on ,
$t$

$t$

(2.8)

$\frac{dy}{dt}=hf(t)(0<t<1),$3$y(0)=0$.

.

:$\backslash$’$20^{t}6$

Then it reduces the RK method to the quarature rule(2.9)

$y_{1}=h \sum_{:=1}^{l}b:f(c:)$,

where the right-hand side is the approximation of the integral$y(1)=h \int_{0}^{1}f(t)dt$.

The condition $B(p)$ implies that (2.9) is exact if is merely a polynomialof of degree at most $p-1$ . This means the quadrature (2.9) is of order .
From (2.8), we have
$f$

$t$

$p$

(2.10)$Y_{:}=h\sum_{j=1}^{\iota}a;jf(c_{j})$for each , which should be considered as an approximation of the integral$i$$y( c:)=h\int_{0}^{e}:f(t)dt$.

Let $P(t)$ be an interpolating polynomial of degree at most . Then thederivative of $P(t)$ is integrated exactly at every internal point . Thus we
obtain from (2.10)
$s$

$c:$

$Y_{:}-y_{0}=P(c_{i})-P(0)=h\int_{0}^{\iota}:P’(t)dt=h\sum_{j=1}^{l}a:jP^{t}(c_{j})$.

By virtue of (1.2b), the equation$f(c_{j})=P’(c_{j})(j=1, \ldots, s)$holds, which means that the condition $C(s)$ implies the $\dot{u}\iota terpolancy$ of thederivative $P’(t)$ at every internal point $c_{i}(i=1, \ldots,s)$ .
In conclusion, the collocation method gives the exact values at $t=1$
for the solution, at $t=c_{i}$ for its derivative if it is a polynomial of degree at

most

$s$

.

Further, a stronger statement holds. Here Butcher made an essential
contribution in showing the following
1

4

207

Theorem 1. Let an s-stage $RK$ formula have distinct abscissas
and nonzero weights. Then the following implications hold:$c_{1},$$b_{1},$$c_{l}$(2.11)(2.12)(2.13)(2.14)(2.15)(2.16)(2.17)$A(p)$

今$A(p+q)$,$E(p, q)$$B(p+q)$ and $C(q)$
$p$

$s$

十

,
$E(p,q)$ ,$\Rightarrow$) and $D(p)$十 ) 伽 $dE(s, q)$$B(p+s)$ and $E(p,s)$$B(p),$ $C(q)$ and $D(r)$
(
$B$ (

$B$

,

$\ldots,$$B(p)$$\Rightarrow$$\ldots$$b_{\iota}$$\Rightarrow$$q$,$C(q)$ ,$D(p)$ ,$E(p, q)$$q$$\Rightarrow$$A(p)$$\Rightarrow$$(p \leq\min(q+r+1,2q+2))$.

From (2.17), the collocation RK method of s-stage is consistent of order at
least for general initial value problem (1.1) if it has distinct abscissas andnonzero weights. This is easily verified by putting $p=q=s$ and $r=0$ in(2.17). Moreover, it is not difficult to show$s$beTheorem 2 (Dekker-Verwer,2) . Let distinct abscissas
given. Then the simplifying conditions $B(s)$ and $C(s)$ uniquely define an$RK$ method with nonzero weights.
$s$$c_{1},$$\ldots,$$c_{\iota}$

Hence, we can focus on the determination of the abscissas for the collocation
RK method. Looking back to the interpolating property at the internal
point $t=c_{1}$ , we have to consider a polynomial on $[0,1]$ of degree whichinterpolates the solution $y(t)$ . Usually it can be constructed on the fixed. Its Lagran$ge$ form is given by
in terpolat $ing$ points$s$$\xi_{0},$$\xi_{1},$$\xi_{l}$$\ldots,$(2.18)$L_{\iota}(t)= \sum_{j=0}^{l}y(\xi_{j})\ell_{j}(t)$where(2.19)$l_{j}(t)= \prod_{k\neq j}\frac{t-\xi_{k}}{\xi_{j}-\xi_{k}}$.

on $[0,1]$ , the absc issas
Theorem 3. For fixed distinct points.. ., of the s-stage collocation method satisfy the equation$\xi_{0},$$\xi_{\iota}$$\ldots,$$c_{\iota}$(2.20)$\omega_{l}’(c:)=0$5$c_{1}$

,

208

where
$\omega_{\iota}(t)=\prod_{j}(t-\xi_{j})=(t-\xi_{0})(t-\xi_{1})\cdots(t-\xi_{\iota})$.

Proof will be given by considering the well-known residual formula of the
Lagrange interpolating form.
Most conceivable choice for the interpolating points
is the New ton-Cotes $type$ equidistant distribution$\xi_{j}=\frac{j}{s}$$(j=0,1, \ldots,s)$$\varphi_{l}^{1}(t)=0$$\xi_{\iota}$$\ldots,$

on $[0,1]$.

In this case, the deternining equation for the abscissas
by the algebraic equation(2.21)$\xi_{0},$$c_{1},$$\ldots,c_{\iota}$

is given,where(2.22)$\varphi,(t)=\prod_{j=0}^{l}(t-\frac{j}{s})=t(t-\frac{1}{s})\cdots(t-\frac{s-1}{s})(t-1)$.

It is obvious that the equation (2.21) has distinct roots on $(0,1)$ , and they
are situated as $(i-1)/s<c_{i}<i/s(i=1, \ldots,s)$ if they are labelled as$s$$c_{1}<c_{2}<\cdots<c,$.

Since the Newton-Cotes type formula of s-stage is a collocation one, it
is consistent of order at least . Here we arrive at a question what is theactual order of consistency of the Newton-Cotes type. Before attacking thisquestion, we will state the following Lemma.
$s$Lemma 1 Let$\phi_{p,q}(t)$be a polynomial of degree $(p+2q-1)$ defined by$\phi_{p,q}(t)=t^{q}(t-1)^{q}\prod_{j=1}^{p-1}(t-\frac{j}{p})$with $p,q\in N$ . Then we have$\int_{0}^{1}\phi_{p.q}(t)dt\{\begin{array}{l}=0forevenp\neq 0foroddp\end{array}$

6

209

Further, let$\phi_{p,q}^{*}(t)$be a polynomial of degre $e(p+2q)$ given by$\phi_{p,q}^{r}(t)=t\phi_{p,q}(t)$,$\int_{0}^{1}\phi_{p,q}^{*}(t)dt\neq 0$.

then we have

Proof will be carried out through an elementary integral calculus.
Theorem 4. The s-stage Newton-Cotes type formula is consistent of order $s+1$ for odd , of order $s+2$ for even .
$s$

$s$

Proof. By virtue of (2.17), it suffices to prove that the condition $B(s+1)$

and $B(s+2)$ are satisfied for the odd case and the even case, respectively.
The condtion $B(s+1)$ means that the Newton-Cotes formula is exact
for a polynomial $f(t)$ of degree at most as the numerical quadraute rulefor (2.8). So, let $f(t)$ be a polynomial of degree . Then, in the division of$f(t)$ by, the quotient and remainder are polynomials of degree andof degree less than , respectively:
$s$

$s$

$0$

$\varphi_{\iota}’(t)$$s$$f(t)=q\varphi_{l}’(t)+r_{\iota-1}(t)$.

Thus$y(1)$$=$$h \int_{0}^{1}\{q\varphi’,(t)+r_{\iota-1}(t)\}dt$$=$$qh[ \varphi_{\iota}(1)-\varphi_{\iota}(0)]+h\int_{0}^{1}r_{\iota-1}(t)dt$$=$$h \int_{0}^{1}r_{\iota-1}(t)dt$.

On the other hand, the identity$\sum_{i=1}^{l}b_{i}f(c_{i})$

$=$

$=$

:

$\sum_{:}$$\sum_{:}b_{i}r_{e-1}(c;)$

holds. The polynomial
is of degree less than andis equal to
.
Thus
$s$

$r_{\iota-1}$$\sum b_{i}r_{\iota-1}(c_{i})$$\int_{0^{1}}r_{\iota-1}(t)dt$7$c_{1},$ $\ldots,c_{\iota}$

are distinct.

210$\sim$

Next, let $f(t)$ be a polynomial of degree $s+1$ . Then we have$f(t)=q(t)\varphi’.(t)+r_{\iota-1}(t)$,where $q(t)=q_{1}t+q_{0}$ andis a polynomial of degree less than
to the above, we can deduce as$r_{-1}$$y(1)$$=$$h \int_{0}^{1}\{q(t)\varphi_{l}^{t}(t)+r_{\iota-1}(t)\}dt$$=$$h[q(1) \varphi_{\iota}(1)-q(0)\varphi.(0)]-q_{1}h\int_{0}^{1}\varphi_{\iota}(t)dt$$=$$-q_{1}h \int_{0}^{1}\varphi_{\iota}(t)dt+h\int_{0}^{1}r_{-1}(t)dt$$s$. Similar

十 $h \int_{0}^{1}r_{\iota-1}(t)dt$

because.
By virtue of Lemma 1 in the case of $p=s$ and $q=1$ , we have$\varphi_{l}(1)=\varphi_{\iota}(0)=0$$\int_{0}^{1}\varphi,(t)dt\{\begin{array}{l}=0forevenp\neq 0foroddp\end{array}$This implies that the condition $B(s+2)$ is satisfied for even while itcannot be satisfied for odd . But if $g(t)$ is of degree 2, i.e. if $f(t)$ is ofdegree $s+2$ , Lemma 1 states
$s$

$s$

$\int_{0}^{1}t\varphi_{l}(t)dt\neq 0$

Hence $B(s+3)$ never be satisfied for even

$s$

.

.口

For some $s’ s$ , the graph ofcan be found in Fig. 1. Theyandshow the position of the interpolating points
and the abscissasof the Newton-Cotes type collocation method.
$\varphi_{l}^{t}(t)$$\varphi_{\iota}(t)$$\xi_{0},$$\xi_{\iota}$$\ldots,$$c_{1},$$\ldots,$$c_{l}$A physists group ([1]) has presented a discrete variable method for ODE’s
which can be found to be equivalent to the Newton-Cotes type collocation
method. They called it HIDM, which sounds slightly unsuitable. (SeeAppendix.)

2.2 A Series of Methods of the Same Stage Number8$21_{1}^{s,}$To increase the order of consistency for the s-stage collocation method,
we will make a gradual change of the distribution of the interpolating points
which causes the change of the distribution of the collocation
points.
Letbe the polynomial of degree $(s-k)$ given by$\xi_{0},$$\xi_{\epsilon}$$\ldots,$$c_{1},$$\ldots,c_{\iota}$$\psi_{k}(t)$(2.23)$\psi_{k}(t)=\{\begin{array}{l}\prod_{j=1}^{s-k-1}(t--\underline{\dot{L}}_{\overline{k}})fork=0,\ldots,s-21fork=s-1\end{array}$Further, define the polynomial(2.24)$\varphi_{\iota,k}(t)$of degree $s+1$ by$\varphi_{\iota.k}(t)=\frac{d^{k}}{dt^{k}}\{t^{k+1}(t-1)^{k+1}\psi_{k}(t)\}$$(k=0,1, \ldots, s-1)$.

Note that the identity$\{\rho_{\iota,0}(t)=\varphi_{\iota}(t)$holds. Since $t=0$ and $t=1$ are both $(k+1)$ -ple root of $t^{k+1}(t-1)^{k+1}\psi_{k}(t)$ ,are. Obviously $s+1dist\dot{u}lct$ real roots ofsituated as$\varphi_{e,k}(0)=\varphi_{\iota,k}(1)=0$$\varphi_{\iota,k}(t)=0$$\xi_{0}(=0)<\xi_{1}<\cdots<\xi_{\iota-1}<\xi_{l}(=1)$.

Thus, $\varphi_{*,k}’(t)=0$ can be the determining equation for the abscissas of thecollocation method.

Theorem 5. The s-stage collocation formula determined by the equation(2.25)$\varphi_{\iota,k}^{\iota}(t)=0$

is consistent of order $s+k+1$ if $s-k$ is odd, and of order $s+k+2$ if $s-k$
is even.

Proof is similar to that of Th.4. Hence it is omitted here.
From Theorems 4 and 5, we can show the table of increasing the orderof consistency in this class. Refer to Table 1.

9

212

Table 1. The order of consistency for $(s,k)$ .

Theorem 6. The s-stage collocation method determ ined by the equation$\varphi_{\iota.\iota-1}^{t}(t)=0$

is nothing but the s-stage Gauss-Legendre method.
This theorem can be proved through the Rodrigues’ formula for the
Legendre polynomials.
In Fig.2, one can see the change of the zeros for the polynomial$\frac{d^{\ell}}{dt^{\ell}}${t ん+l $(t-1)^{k+1}\psi_{k}(t)\}$ $(p=0,1,$$\ldots$$k+1)$

for some pairs of $(s,k)$ .
On closing the present section, we have established a way to generate aseries of s-stage collocation method starting from the Newton-Cotes type
( $(s+1)- st$ or $(s+2)-nd$ order) up to the Gauss-Legendre formula ( $2s$ -thorder), increasing the order of consistency two by two.

3. A-Stability of the Collocation RK Method
3.1 Derivation of the Stability Factor
For the consideration of the stability of RK method, the matrix andvector notations make the process simple and elegant. Let and be thematrix and the vector given by$A$$A=(a:j)$$(1 \leq i,j\leq s)$

10

and $b=(b_{1}, \ldots,b,)^{T}$ ,

$b$

213

respectively. Defining the vector all of whose components are equal tounity, the stability factor $R(z)$ of the RK method (1.2) is given by
$e$

(3.1)

,

$R(z)= \frac{\det(I-zA+zeb^{T})}{\det(I-zA)}$

which can be found in many literatures (e.g. $Dekker-Verwer[2]$ ).
and are uniquely determined
In the s-stage collocation method,
providedby the set of the abscissasthat they are distinct (Th.2).
Thus, let us define the polynomial
$A$

$b$

$c_{1},\ldots,c_{\iota}$(3.2)$\pi_{\iota}(t)=\prod_{j=1}^{l}(t-c_{j})=(t-c_{1})(t-c_{2})\cdots(t-c_{\iota})$,which is the derivative ofin Th.3 divided by its leading coefficient tomake it monic. Thus we obtain$\omega_{\iota}(t)$(3.3)$\pi_{\iota}(t)=t^{*}+\sum_{j=0}^{s-1}p_{j}t^{j}$,$(j=0, \ldots,s-1)$ can be represented by thewhere the coefficients.
fundamental symmetric functions of
We introduce the Vandermonde matrix$p_{j}$$c_{1},\ldots,c_{\iota}$(3.4)$V=\{\begin{array}{lll}1 c_{1} c_{1}^{s-1}1 c_{2} c_{2^{-1}}^{l}\vdots \vdots \vdots 1 c_{l} c_{l}^{s-1}\end{array}\}$and the diagonal matrices(3.5)$C=diag[c_{1},c_{2}, \cdots,c_{\iota}]$

,

(3.6)

$S=diag[1,$.

$\frac{1}{2},$$\cdots,$$\frac{1}{s}]$Then, the simplifying conditions $B(s)$ and $C(s)$ are equivalent to the matrix
identities given by the following (Dekker-Verwer[2]):(3.7)$B(s)$:$b^{T}V=e^{T}S$

11

,

214

(3.8)

$AV=CVS$.

:$C(s)$The direct calculation leads to the equation$CVS=[\dot{d}_{:}/j]$.

On the other hand, the identity$[1,c_{i},$ $\cdots,c_{:}^{\iota-1}]\{\begin{array}{lllllll}$=$$[c:,$$\frac{c_{:}^{2}}{2}$)\end{array}\}$$\frac{c_{:}^{\epsilon-1}}{s-1}-\frac{\sum_{j=0}^{\epsilon-1}p_{j}c_{:}^{j}}{s}]$holds. By virtue of (3.3), the right-hand side of the above identity is equal
to the i-th row vector of the matrix product $CVS$ . Hence, we can establish
the equation$0 \frac{1}{2}$(3.9)$V^{-1}CVS=[01$.. ..

$-p_{\iota-1}^{l-2^{S}}-P-P..../_{/^{s}}-p^{0}/_{/}^{s_{S}}:^{1}\ovalbox{\tt\small REJECT}$$\frac{01}{\iota-1}$

Let $W$ be the diagonal matrix given by(3.10)$W=diag[1,1,2],3!,$12$\cdots,$$(s-1)!$ ].

.

@i 5

From (3.8),(3.9) and (3.10), we obtain$WV^{-1}CVSW^{-1}=WV^{-1}AVW^{-1}=\{\begin{array}{llllll} \end{array}\}$.

The right-hand side matrix is the companion matrix for the polynomial(3.11)$d(z)= \frac{p_{0}}{s!}+\frac{p_{1}}{s!}z+\cdots+\frac{(s-1)!p_{\epsilon-1}}{s!}z^{\iota-1}+z^{l}=\frac{1}{s!}\sum_{j=0}^{l}j!p_{j}z^{j}$.

That is, the identity(3.12)$d(z)=\det(zI-WV^{-1}AVW^{-1})$

holds. Hence, we arrive at the following lemma.
Lemma 2. The denominator polynomial
is given by(3.13)$\Delta(z)$$\Delta(z)=z^{\iota}d(\frac{1}{z})$of the stability factor $R(z)$.

Proof. By calculation, we have$\Delta(z)$$=$$\det(I-zA)=z^{\iota}\det(\frac{1}{z}I-A)$$=$$z^{\iota} \det(\frac{1}{z}I-WV^{-1}AVW^{-1})=z^{\iota}d(\frac{1}{z}).\square$

Next, we will consider the numerator of $R(z)$ . The numerator $N(z)$is given by $N(z)=\det(I-zA+zeb^{T})$ . From (3.7) and (3.8), we obtain$eb^{T}=ee^{T}SV^{-1}$ and $A=CVSV^{-1}$ , which implies(3.14)$A-eb^{T}=(CV-ee^{T})SV^{-1}$

13

.

216

The $(i,j)$ -th component of the matrix $CV-ee^{T}$ is given by. Then,putting $\gamma:=c:-1(i=1, \ldots,s)$ , we introduce the Vandermonde matrix
and the diagonal matrix by$c_{:}^{j}-1$$U$$\Gamma$(3.15)$U=\{\begin{array}{llll}1 \gamma_{1} \cdots \gamma i^{-1}1 \gamma_{2} \cdots \gamma_{2^{-1}}\vdots \vdots \vdots 1 \gamma_{l} \cdots \gamma_{l}^{l-l}\end{array}\}$(3.16)$\Gamma=diag[\gamma_{1},\gamma_{2}, \cdots,\gamma,]$,.

Taking the identity$\not\simeq_{:}=(c_{i}-1)^{j}=\sum_{\llcorner-0}^{j}(-1)^{j-\ell(\begin{array}{l}jp\end{array})c_{:}^{\ell}}$

into the consideration, the two matrices
In fact, we obtain(3.17)

$V$

and

$U$

are. in the linear relation.

$U=V\{\begin{array}{lll}11-1 -1 (-1)_{\iota+}_{\iota_{2}^{-1}}^{\iota_{1}}\{}(-1)^{\iota+2}(-1)_{-1}^{\iota+1}1-2 3 |1-3 \cdots \cdots|1 | -(_{-2}^{\iota-1})O 1\end{array}\}$On the other hand,$(CV-ee^{T})S=\{\begin{array}{lll}c_{l}-1 (c_{1}^{2}-1)/2 (c_{1}^{t}-1)/s\vdots \vdots \vdotsc_{l}-1 (c_{l}^{2}-1)/2 (c_{l}^{l}-1)/s\end{array}\}$

14

.

.

217

Thus, from (3.17) we have$(CV^{-}-e\ell^{T})SV^{-1}U$$=$$\{\begin{array}{lll}\gamma_{1} \gamma_{l}^{2}/2 \gamma i/s\vdots \vdots \vdots\gamma_{l} \gamma_{l}^{2}/2 \gamma_{l}^{l}/s\end{array}\}$$\Gamma US$$=$,

which, with (3.14), implies the identity$(A-eb^{T})U=rus$(3.18).

Hence, similar to the denominator case, we can establish the following, let us define the polynomial
result. By the shift of the variable inby$t$$\pi_{\iota}(t)$$p_{s}(t)$(3.20)$\rho_{\epsilon}(t)$$=\pi_{\iota}(t+1)$$=$$t^{\iota}+ \sum_{j=0}^{l-1}r_{j}t^{j}$,where $r_{j}(j=0, \ldots,s-1)$ is the coefficient.Further, introducing the polynomial(3.21)$e(z)= \frac{r_{0}}{s!}+\frac{r_{1}}{s!}z+\cdots+\frac{(s-1)!r_{\iota-1}}{s!}z^{\iota-1}+z^{\iota}=\frac{1}{s!}\sum_{j=0}^{l}j!r_{j}z^{j}$.

we have the following
Lemma 3. The numerator polynomial $N(z)$ of the stability factor $R(z)$is given by(3.22)$N(z)=z^{\iota}e( \frac{1}{z})$.

Due to the above Lemmas, $R(z)$ is represented by(3.23)$R(z)= \frac{N(z)}{\Delta(z)}$,and $N(z)$ are given by (3.13) and (3.22), respectively. The
where
magnitude
principle implies the necessary and sufficient condimaximum
tions $(S1)$ and $(S2)$ for the A-stability of the collocation RK method as
follows:$\Delta(z)$

15

218

$z$

.

$(S1)$All the roots of$\Delta(z)$belong to the right half-plane of the complex$(S2)|N(iy)|\leq|\Delta(iy)|$ for every real $y$. (The imaginary unit is denoted

by )
$i.$

, is symmetricDeflnition 3 If the $distr;bution$ of the abscissas$RK$with respect to 1/2, the colloc ation
method of these abscissas is called
as symmetric.
$c_{1},c_{2},$ $\ldots,$$c$symmetric, it is A-stable iff all
Theorem 7 If the collocation method
the zeros ofare in the right half-plane.
$|s$$\Delta(z)$3.2 A-Stability of the Newton-Cotes Type Methods
In the case of the Newton-Cotes type method, the collocation points are
given by (2.21). The polynomialis then given by$\pi_{\iota}(t)$$\pi_{\iota}(t)=\frac{1}{s+1}\varphi’,(t)$whereis defined by (2.22). Utilizing Lemma 2, we can calculate thefor various number of $s’ s$ . By Theorem 7,
denominator polynomials
whether all of them are situated
we are to investigate the zeros ofin the right half-plane. In other words, the question is whether all of the
The symbolic and algebraiczeros of $\Delta(-z)$ are in the left half-planecomputation by computer gives the polynomials $\Delta(-z)$ as in Table 2. As
can be seen, all the coefficiets of one polynonial have the same sign withthe leading coefficient. Thus they satisfy the necessary condition for therequired root distribution. But, as is well known, the exact solution is only
possible by $e.g$ . the Routh-Hu rwitz crit eria.
Given a polynomial $f(z)$ as$\varphi_{\iota}(t)$$\Delta(z)$$\Delta(z)$$C^{-}J$$f(z)=f_{0}z^{n}+f_{1}z^{n-1}+\cdots+f_{n-1}z+f_{n}$,where $f_{0}>0$ , then the Routh-Hurwitz criteria first defines the following
n-square matrix .
$Q$

16

a19

$Q=[f_{0}^{\iota}f0_{:}0_{:}f_{1}^{2}f_{0}^{\}ff:...\cdot$$f_{2}^{6}f^{\}f_{4}f:.\cdot.\cdot$...

$f_{2n-4}^{2n-1}f_{2n-\}^{2n-2}ff::.]$.

Here $f_{j}=0$ for $j>n$ . The necessary and sufficient condition for all theis that all the principal minors of are positive.
roots of $f(z)$ being in
We have calculated the matrix
and its principal minors for every$\Delta(-z)$ , whose results are shown in Table 3. This leads us to the following.
$C^{-}$$Q$$Q$

Theorem 8 The Newton-Cotes type collocation methods whose stage num$ber$ is less than 9 are A-stable.

Table

$e$

. $A(-z)$ for the Newton-Cotes type formulas

$z)$

$=$

$z^{2}+6z+12$$\Delta_{\}(-z)$$=$$z^{\}+11z^{2}+54z+108$$\Delta_{4}(-z)$$=$$3z^{4}+50z^{3}+420z^{2}+1920z+3840$$\Delta_{6}(-z)$$=$$12z^{\epsilon}+274z^{4}+3375z^{\}+25500z^{2}+112500z+225000$$\Delta_{6}(-z)$$=$$5z^{6}+147z^{6}+2436z^{4}+26460z^{\}+189000z^{2}+816480z+1632960$$\Delta_{7}(-z)$$=$$30z^{7}+1089z^{6}+22981z^{\epsilon}+331681z^{4}+3361400z^{\}+23193660z^{2}$$\Delta_{8}(-z)$$=$$105z^{8}+4566z^{7}+118124z^{6}+2153088z^{6}+28734720z^{4}$$\Delta_{2}($ 一$+98825160z+197650320$$+278691840z^{\}+1878589440z^{2}+7927234560z+15854469120$$\Delta_{9}(-z)$$=$$560z^{9}+28516z^{8}+879525z^{7}+19539360z^{6}+327229875z^{6}$$+4151341530z^{4}+39060913500z^{\}+258918055200z^{2}$$+1084777369200z+2169554738400$

17

220

Table S. Principal minors of the Newton-Cotes type formulas

Appendix is omitted because of the limit of space. Hence interested
reader will be given the full manuscript on his request by the authors.

References[1] Abe,K., Ishida,A., Watanabe,T., Kanada, Y. &Nishikawa,K., HIDM
–New numerical method for differntial equations, Kakuyugo-kenkyu,
57(1987), 85-95 (in Japanese).
[2] Dekker,K. &Verwer, J.G., Stability of Runge-Kutta Methods for Stiff
Nonlinear Differntial Equations, North-Holland, 1984.
[3] Burrage, K., High order algebraically stable Runge-Kutta methods,
BIT,18(1978), 373-383.

[4] Butcher, J.C., Implicit Runge-Kutta processes,
18(1964), 50-64.

18Math. Comp.,

221
Fig.1 The graphs of$\varphi_{*}(t)$and$\varphi_{s}’(t)$.

Fig.l.l The case of $s=3$ .

$\theta^{(XI}$$\rho’\iota\cross)$

Fig.1.2 The case of $s=4$ .

$\phi t\cross)$$\phi’1X1$$\nearrow$?

222
Fig.1.3 The case of $s=5$ .

$\ell tXl$$\uparrow’\iota\cross)$

Fig.1.4 The case of $s=6$ .

$\rho(\cross)$$\phi’1\cross)$$-\lambda_{O}$223
Fig.2 The zeros of the polynomials$\neg^{\ell}d^{d_{t}}\{t^{k+1}(t-1)^{k+1}\psi_{k}(t)\}$.

Fig.2.1 The case of $k=0(P=0,1)$ .

$\phi$

1XI$\backslash \backslash \backslash _{\backslash }$$\phi’1\cross)$$\vdash^{\backslash }A\backslash \backslash \wedgearrow-||\ovalbox{\tt\small REJECT}_{\vee\vee}^{\wedge}\wedge-|\backslash \backslash \backslash |$$\iota$$\backslash \backslash$$-\lambda..\backslash \backslash _{\backslash }$$\backslash \backslash _{\iota}$$\backslash$$\backslash \backslash$Fig.2.2 The case of$k\cdot\overline{\backslash ^{-}}1(p=0,1,2)$.

$\backslash$$\phi 1X)$$\ell’1XI$$\frac{|\wedge\wedge\wedge|\wedge\wedge\wedge|}{|\vee\vee\vee|\vee\vee\vee|}$冫1

磨 $Zg$Fig.2.3 The case of $k=2(\ell=0,1,2,3)$ .

$t^{(X1}$$\rho^{r}t\cross)$

Fig.2.4 The case of $k=3(\ell=0,1, \ldots,4)$ .

$\phi 1X)$$\uparrow^{r}tX1$$\frac{|\wedge\wedge\wedge|\wedge\wedge\wedge|}{|arrow\vee\vee|Rarrow|}$$-\rangle$

冫225
Fig.2.5 The case of $k=4(P=0,1, \ldots,5)$ .

$1=1$$\phi$

1X1$\ell’1XI$$|^{\vee\vee}|\vee||_{\wedge\wedge}||\ovalbox{\tt\small REJECT}_{\hat{\vee}}^{\wedge}rightarrow-$

Fig.2.6 The case of $k=5(\ell=0,1, \ldots, 6)$ .

$\theta^{(X1}$$\phi’$

1X1$\frac{|\wedge\wedge\wedge|\wedge\wedge\wedge|}{|\vee\vee\vee|\vee\vee\vee|}$冫$\cdot 3$