# Efficient Recursive Methods for Partial Fraction Expansion of General Rational Functions

Authors: Youneng Ma; Jinhua Yu; Yuanyuan Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Partial fraction expansion (pfe) is a classic technique used in many fields of pure or applied mathematics. The paper focuses on the pfe of general rational functions in both factorized and expanded form. Novel, simple, and recursive formulas for the computation of residues and residual polynomial coefficients are derived. The proposed pfe methods require only simple pure-algebraic operations in the whole computation process. They do not involve derivatives when tackling proper functions and require no polynomial division when dealing with improper functions. The methods are efficient and very easy to apply for both computer and manual calculation. Various numerical experiments confirm that the proposed methods can achieve quite desirable accuracy even for pfe of rational functions with multiple high-order poles or some tricky ill-conditioned poles.


---

## Body

## 1.  Introduction

Partial fraction expansion (pfe) has been a powerful tool widely used in the field of calculus, differential equations, control theory, and some other areas of pure or applied mathematics. It is also a classic topic studied by many scholars over the time. Though a sole solution is guaranteed, it is not an easy task to perform the computation of pfe effectively, especially when the functions to be expanded contain high-order or ill-conditioned poles. Many existing methods are difficult to apply and can lead to considerable errors. We aim to research this classic topic and propose several novel more efficient and simpler pfe methods for general rational functions in both factorized form and expanded form.SupposeR
(
s
)
=
Q
(
s
)
/
P
(
s
) is the rational function to be expanded in pfe. The polynomials, P
(
s
) and Q
(
s
), are denominator and numerator of R
(
s
), respectively. Generally, a polynomial (let us say Q
(
s
)), can be written in either expanded form Q
(
s
)
=
∑
i
a
i
s
i or factorized form Q
(
s
)
=
∏
i
(
s
-
z
i
)
n
i in terms of its zeros. Though theoretically equivalent, a method based on one of these forms may exhibit significant differences from the other in numerical accuracy and computational efficiency [1]. In pfe problems, as the poles of R
(
s
) have to be known, P
(
s
) is mostly written in factorized form. Though Q
(
s
) can be written in either factorized or expanded form, the latter seems to be more in favor in practice. Most current available pfe methods are designed for R
(
s
) with Q
(
s
) written in expanded form. However, one may still be motivated to design pfe methods for factorized R
(
s
). For one thing, many functions are originally written in factorized form. A pfe method for factorized R
(
s
) can deal with such functions directly and more efficiently, therefore, it is more preferable. For another, under many practical circumstances, even R
(
s
) is originally written in expanded form, the zeros and poles of R
(
s
) have to be obtained. In other words, R
(
s
) needs to be rewritten in factorized form. This can be often seen in the field of system analysis and design, where pfe methods are widely used in the derivation of Inverse Laplace transformation of transfer function. As the zeros and poles of the transfer function almost characterize the whole system, they are often obtained. The zero-pole plot of the transfer function is almost used as a routine in the analysis of signal and system. Considering that methods for factorized functions can be more easily performed, thus they can serve as an alternative scheme under such circumstances. Some other practical reasons for the development of pfe methods for factorized R
(
s
) are also mentioned in [2]. In a word, it is of significance to design methods for the pfe of R
(
s
) in both factorized and expanded form.The most well-knownpfe method is Heaviside’ cover-up formula. It is often introduced as a textbook method. It provides an elegant compact solution to pfe problems and often serves as a basis for other pfe methods. However, this classic method requires successive differentials when applying to functions with high-order poles, which gives rise to increasingly higher order polynomials. Evaluation of these high-order polynomials can eventually lead to large numerical errors [3]. Another standard pfe method is the method of undetermined coefficients. This method requires the construction and solution of a system of equations. It can also be very complicated and inconvenient when tackling functions with high-order poles. Besides the standard classic methods, many pfe methods are also proposed during the past years. Some of these methods tend to perform better than classical methods under certain conditions. However, many of them are only suitable for small-scale problems or some particular cases and are difficult to fit into computer program. Two desirable methods are specially designed for the pfe of functions in factorized form [3, 4]. In [4], Brugia developed an efficient noniterative method for pfe of functions with high-order poles by deriving the high-order derivative of the rational functions directly without passing its lower order ones. Compared with the classic methods and many other “tricky” methods, Brugia’s method is more applicable and efficient. It can be easily programed for computer computation. Linnér [3] simplified Brugia’s method by using operators and extended it to improper functions by Laurent expansion. In [5], a modification was also made to Brugia’s method, which leads to relative simpler implementation and a review of many nice “aged” methods can also be found in [5]. Most of existing methods assume that the rational function is written in expanded form. In [6], Karni proposed an algebraic procedure for pfe. Karni’s algorithm, unfortunately, is limited to just those cases where the degree of the numerator is no more than the number of poles. Karni’s method was further developed yet with much implementation costs in [7–9]. Pottle [10] proposed an iterative method for the digital computer use, but this method is subject to intolerable errors when applied to functions with high-order poles. Uraz and Nagy [11] proposed a method calculating the residues and coefficients utilizing matrix algebra. Besides the above relatively “aged” methods, many methods are also proposed in the recent years, such as the methods in [12–17]. Many of these methods are also only suitable for some particular cases and can become very complicated for large-scale problems. In [17], Ma et al. proposed an efficient pure algebraic method for pfe of general rational functions with high-order poles. It exhibits quite good performance and avoids long division when dealing with improper functions. The method does not involve derivative and can be coded conveniently. The calculation complexity of some pfe methods is discussed in [18, 19].In this paper, we develop efficient recursive algebraicpfe methods for rational functions in both factorized and expanded form. Simple, elegant, and pure-algebraic formulas are proposed to compute the pfe coefficients. The proposed methods can be used for pfe of both proper and improper rational functions with complex poles. They do not involve polynomial division or differentiation or the solution of a system of linear equations for the pfe of a general rational function. The remainder of the paper is organized as follows. In Section 2, we focus on the pfe of functions in factorized form. Brugia’s and Linnér’s methods are further developed and simplified for much easier implementation and extended to more general cases. We also develop novel simple algebraic methods to compute the coefficients of the residual polynomial of improper functions that avoids polynomial division. In Section 3, we focus on the pfe of rational function in expanded form. Two novel efficient methods are proposed. Several useful, compact formulas are derived for the computation of residues. Two efficient methods that avoid polynomial division are developed to compute the coefficients of the residual polynomial. In Section 4, numerical examples are provided to illustrate the usage and validity of the proposed methods, followed by the conclusion in Section 5.


## 2. Partial Faction Expansion of Rational Functions in Factorized Form

Let the factorized rational function be(1)
R
(
s
)
=
∏
j
=
1
N
(
s
-
z
j
)
n
j
∏
j
=
1
M
(
s
-
s
j
)
m
j
,
which is expanded as
(2)
R
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
c
)
h
+
∑
k
=
1
M
∑
L
=
1
m
k
c
k
L
(
s
-
s
k
)
L
.
Here c
k
L are the residues; z
j and s
j are the zeroes and poles of R
(
s
), respectively. The poles and zeroes are assumed distinct in the whole paper; n
j and m
j are the multiplicities of z
j and s
j, respectively; N is the number of zeros and M is the number of poles; e
h are the coefficients of the residual polynomial which will be zeros for a proper R
(
s
); s
c is a constant introduced to make the residual polynomial more general. s
c is often set as zero in other articles. From (1) and (2), we can observe that the degree of the numerator of R
(
s
) is V
=
∑
j
=
1
N
n
j; the degree of the denominator is K
=
∑
j
=
1
M
m
j; E
=
V
-
K. In a pfe problem, c
k
L, e
h are the unknown coefficients to be calculated. We first propose a simple recursive formula to obtain the residues and then provide two algebraic methods for obtaining coefficients of the residual polynomial for an improper function.


### 2.1. Computing Residues of Factorized Functions

Brugia [4] developed a classic pfe method for factorized rational functions. His method is more efficient and easier to apply for both manual and computer calculation compared with many other methods. Here we are to further develop and simplify Brugia’s method. To allow for context and illumination of differences and to make this paper self-contained, part of Brugia’s method is briefly introduced below. According to Heaviside’s formula, we have
(3)
c
k
L
=
R
k
(
K
)
|
s
=
s
k
K
!
,
where K
=
m
k
-
L and
(4)
R
k
(
s
)
=
∏
j
=
1
N
(
s
-
z
j
)
n
j
∏
j
=
1
,
j
≠
k
M
(
s
-
s
j
)
m
j
.
As shown in [4], the Kth derivatives of R
k
(
s
) are
(5)
R
k
(
K
)
=
(
R
k
D
k
)
(
K
-
1
)
=
∑
i
=
0
K
-
1
C
K
-
1
i
R
k
(
K
-
i
-
1
)
D
k
(
i
)
,
where K
>
0 and
(6)
D
k
(
i
)
=
(
-
1
)
i
i
!
(
∑
j
=
1
N
n
j
(
s
-
z
j
)
i
+
1
-
∑
j
=
1
j
≠
k
M
m
j
(
s
-
s
j
)
i
+
1
)
.
Here C
m
i
=
m
!
/
(
m
-
i
)
!
i
!, i
=
i
(
i
-
1
)
⋯
1, and 0
!
=
1. In the above equations, we denote R
k
(
s
) by R
k for simplicity of expression. Similar notations are also used for other functions in the following part. This is Brugia’s main result; more details can be found in [4]. Equations (5) and (6) can generate a linear system of equations and then Gramer’s method can used to solve the equations to obtain the derivatives of R
k
(
s
). It represents a very efficient pfe method and the computer implementation is easier to be accomplished compared with many other existing “tricky” methods. Here we propose a simpler method based on (3), (5), and (6). According to (3), we can obtain
(7)
R
k
(
i
)
|
s
=
s
k
=
i
!
c
k
(
m
k
-
i
)
.
Making some mathematical transformations on (6) yields
(8)
D
k
(
i
)
=
i
!
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
+
1
)
.
Here we mention that the modification on (6) to obtain (8) is relatively little, but in practice it can reduce computation costs and facilitate implementation to a considerable extent. Combining (5) and (3), we have
(9)
c
k
L
=
1
K
!
∑
i
=
0
K
-
1
C
K
-
1
i
R
k
(
K
-
i
-
1
)
D
k
(
i
)
|
s
=
s
k
.
Using (7) and (8) to substitute the derivatives of R
k and D
k in the summation of (9) yields
(10)
c
k
L
=
1
K
∑
i
=
0
K
-
1
c
k
(
m
k
-
(
K
-
i
-
1
)
)
×
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
+
1
)
|
s
=
s
k
.
Notice K
=
m
k
-
L. Equation (10) further yields(11a)
c
k
L
=
1
m
k
-
L
∑
i
=
1
m
k
-
L
c
k
(
L
+
i
)
λ
^
i
|
s
=
s
k
,
(
L
=
m
k
-
1
:
-
1
:
1
)
,
where
(11b)
λ
^
i
=
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
)
.Here we mention that in (11a) L
=
m
k
-
1
:
-
1
:
1 means L
=
m
k
-
1
,
m
k
-
2
,
m
k
-
3
,
…
,
2
,
1. This is a Matlab denotation which will also be used in the following part for the simplicity of expression. Obviously c
k
m
k
=
R
k
|
s
=
s
k. By using (11a) and (11b) and decreasing progressively the value of L from m
k
-
1 to 1, we can then derive other residues successively. Equations (11a) and (11b) can be implemented recursively in a straightforward way. It can be achieved without the aid of a comparatively complex table as done in [3]. It is much easier to apply for manual calculation than the Brugia’s method [4] and its modifications [3, 5]. Meanwhile, it can be programed more easily and reduce calculation costs.


### 2.2. Computing Residual Polynomial Coefficients of Factorized Functions

Generally, as practiced by most current articles, long division or polynomial division can be used to compute the coefficients of the residual polynomial, namely,e
h, in (2). However, as errors will propagate at each step of the polynomial division, this scheme is not quite suitable for large-scale problems. Moreover, for a given factorized rational function, the long division requires both the numerator and denominator be rewritten in expanded form, which can lead to additional errors and much computation costs. Here we provide two novel algebraic methods to compute e
h that avoid long division and the reformation of R
(
s
). The first method described in Section 2.2.1 is obtained by further developing Linnér’s method. The other method introduced in Section 2.2.2 is based on derivatives.


#### 2.2.1. Laurent Expansion Applied to Improper Rational Functions for Anys
c

Brugia’s method is limited to proper functions. In [3], Linnér proposed an excellent procedure to compute the residual polynomial coefficient of improper functions by Laurent expansion. Here we will further develop Linnér’s method and obtain much simpler and more efficient calculation formulas. Linnér’s method demands that s
c in (2) not be equal to zeros or poles of R
(
s
), which can cause inconvenience in practice. We will extend s
c to any value. Linnér’s method is briefly introduced to allow for context, completeness, and illumination of differences. With the substitution s
-
s
c
=
1
/
t, (1) transforms into
(12)
R
(
t
)
=
t
-
E
W
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
,
where z
~
j
=
1
/
(
z
j
-
s
c
), s
~
j
=
1
/
(
s
j
-
s
c
), and W is a constant. Now write
(13)
R
(
t
)
=
t
-
E
R
~
(
t
)
,
where R
~
(
t
) does not include the factor t. In (13), applying a Taylor expansion to R
~
(
t
) at t
=
0, we can obtain
(14)
R
(
s
)
=
∑
h
=
0
H
R
~
(
h
)
(
0
)
h
!
(
s
-
s
c
)
E
-
h
+
o
[
(
s
-
s
c
)
E
-
H
]
,
where H is an integer. Now expand R
(
s
) as
(15)
R
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
c
)
h
+
R
2
(
s
)
,
where R
2
(
s
) is the proper fraction, which can be expanded as the second part on the right-side of (2). Comparing (14) and (15), we have
(16)
e
h
=
R
~
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
This is Linnér’s main result. In the method, the derivatives of R
~
(
t
) are calculated using equations like (5) and (6) and it also requires s
c
≠
z
j and s
c
≠
s
j, which can cause inconvenience in practice. In the following part, we first show how to extend s
c to zeros and poles of R
(
s
). We then obtain much simpler recursive formula for the calculation of e
h.Supposes
c is the kth root of R
(
s
), s
k. We can rewrite
(17)
R
(
s
)
=
(
s
-
s
k
)
-
m
k
R
k
(
s
)
,
where R
k
(
s
) is as defined in (4). With the substitution s
-
s
k
=
1
/
t in R
k
(
s
), we have
(18)
R
k
(
t
)
=
t
-
m
k
-
E
R
~
k
(
t
)
,
where
(19)
R
~
k
(
t
)
=
W
k
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
,
j
≠
k
M
(
t
-
s
~
j
)
m
j
.
Here z
~
j
=
1
/
(
z
j
-
s
k
), s
~
j
=
1
/
(
s
j
-
s
k
), and
(20)
W
k
=
∏
j
=
1
N
(
s
k
-
z
j
)
n
j
∏
j
=
1
,
j
≠
k
M
(
s
k
-
s
j
)
m
j
.
By Taylor expanding R
~
k
(
t
) at t
=
0 in (18), we can obtain
(21)
R
k
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
m
k
+
o
[
(
s
-
s
k
)
E
-
H
+
m
k
]
.
From (17) and (21), we have
(22)
R
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
o
[
(
s
-
s
k
)
E
-
H
]
.
Comparing (15) and (22), we can observe that
(23)
e
h
=
R
~
k
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
Thus s
c is extended to the kth root of R
(
s
), s
k. The same procedure can be used to extend s
c to the kth zero of R
(
s
), z
k.We can easily obtaine
E by setting h
=
E in (23) or by inspection:
(24)
e
E
=
1
.
In Section 2.1, we have simplified Brugia’s method and derived simple recursive formula for the calculation of the residues. Notice that, similar to R
k
(
s
), R
~
k
(
t
) and R
~
(
t
) are both of factorized form and that (3) and (23), which calculate c
i
j and e
h, respectively, are of the similar form. Therefore, the aforeproposed procedure in Section 2.1 for the calculation of c
i
j is also suitable for the calculation of e
h here. Take s
c
=
s
k, for an instance. It can be proven that
(25)
R
~
k
(
K
)
=
(
R
~
k
D
~
k
)
(
K
-
1
)
=
∑
i
=
0
K
-
1
C
K
-
1
i
R
~
k
(
K
-
i
-
1
)
D
~
k
(
i
)
,
where
(26)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
j
≠
k
M
m
j
(
s
~
j
-
t
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
~
j
-
t
)
i
+
1
)
.
From (23), we have
(27)
R
~
k
(
i
)
(
0
)
=
i
!
e
E
-
i
.
According to (23) and (25), we have
(28)
e
h
=
1
(
E
-
h
)
!
∑
i
=
0
E
-
h
-
1
C
E
-
h
-
1
i
R
~
k
(
E
-
h
-
i
-
1
)
D
~
k
(
i
)
|
t
=
0
.
Using (26) and (27) to substitute the derivatives of D
~
k and R
~
k in (28) and then evaluating (28) at t
=
0, we can obtain
(29)
e
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
k
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
k
)
i
)
,
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
E
-
1
:
-
1
:
0
)
.
The same procedure can be used to compute e
h when s
c
≠
s
k. In summary, e
h can be calculated according to(30a)
e
h
=
{
1
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
η
i
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(30b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
z
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
otherwise
)
.Equations (30a) and (30b) are such a simple and compact formula that its implementation is even very easy for hand calculation. It is also much more convenient for computer program and requires less computation costs than the original method in [3]. We just need first calculate η
i using (30b). Then e
h can be immediately obtained using (30a).


#### 2.2.2. Computinge
h through Derivatives

We provide another simple way to computee
h which also avoids the unfavorable polynomial division. That e
E
=
1 is obvious. Therefore, we just need to obtain the remaining coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1. Notice that as the residues c
k
L can be firstly obtained in Section 2.1, they can be seen as known quantities here. Three conditions are considered as below regarding the different values of s
c.(1)
s
c
≠
z
k
and
s
c
≠
s
k
. We first discuss the calculation of e
h when s
c
≠
z
k and s
c
≠
s
k. Here z
k and s
k are the kth zero and the kth pole, respectively. From (1) and (2), we observe that the coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1 can be obtained according to
(31)
e
h
=
1
h
!
(
R
(
h
)
+
T
(
h
)
)
|
s
=
s
c
,
where
(32)
T
(
s
)
=
-
∑
k
=
1
M
∑
L
=
1
m
k
c
k
L
(
s
-
s
k
)
L
.
Here R
=
R
(
s
) as defined in (1) and T
=
T
(
s
). The hth derivative of T
(
s
) can be easily obtained as
(33)
T
(
h
)
=
(
-
1
)
h
-
1
∑
k
=
1
M
∑
j
=
1
m
k
P
h
+
j
-
1
h
c
k
j
(
s
-
s
k
)
j
+
h
,
where P
m
n
=
m
!
/
(
m
-
n
)
!. As R
(
s
) is of factorized form, its derivatives can be derived using functions like (4) and (5). Thus the problem is solved. To further simplify the implementation, we introduce an auxiliary variable
(34)
e
~
h
=
1
h
!
R
(
h
)
|
s
=
s
c
.
From (31) and (34), we can easily observe that
(35)
e
h
=
e
~
h
+
1
h
!
T
(
h
)
|
s
=
s
c
.
Equation (34) has a similar formulation as (3); thus, e
~
h can be calculated and using the procedure we compute c
i
j in Section 2.1. Accordingly, we can finally obtain(36a)
e
~
h
=
{
R
(
s
c
)
(
h
=
0
)
,
1
h
∑
i
=
1
h
e
~
h
-
i
λ
i
(
h
=
1
:
E
)
,
where
(36b)
λ
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
.Here we mention that in (36a) h
=
1
:
E means h
=
1
,
2
,
3
,
…
,
E. This is a Matlab denotation which will also be used in the remaining part for simplicity of expression. We can first obtain e
~
h using (36a) and (36b); then, e
h can be easily obtained using (35). Notice the evaluation of (36b) requires s
c
≠
z
k and s
c
≠
s
k. Thus other two conditions remain to be discussed.(2)
s
c
=
z
k
. Here we discuss the condition when s
c is one of the zeros of R
(
s
) (Let us say z
k). Actually it will be seen that it is more preferable to let s
c
=
z
k. First we rewrite R
(
s
) as
(37)
R
(
s
)
=
(
s
-
z
k
)
n
k
G
k
(
s
)
,
where n
k is the multiplicity of z
k and
(38)
G
k
(
s
)
=
∏
j
=
1
,
j
≠
k
N
(
s
-
z
j
)
n
j
∏
j
=
1
M
(
s
-
s
j
)
m
j
.
Computing the hth derivative of both sides of (37) using Leibniz’s formula and then setting s as z
k yield
(39)
R
(
h
)
(
z
k
)
=
{
0
(
h
<
n
k
)
P
h
n
k
G
k
(
h
-
n
k
)
(
z
k
)
(
h
≥
n
k
)
.
As G
k
(
s
) is of factorized form, its derivatives can be computed recursively(40a)
G
k
(
h
-
n
k
)
=
{
G
k
(
h
=
n
k
)
,
∑
i
=
1
h
-
n
k
P
h
-
n
k
-
1
i
-
1
D
~
k
(
h
-
n
k
-
i
)
λ
~
i
(
h
>
n
k
)
,
where
(40b)
λ
~
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
,
(40c)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
+
1
)
.Combining (31) and (39), we have
(41)
e
h
=
{
T
(
h
)
(
z
k
)
h
!
(
h
=
0
:
n
k
-
1
)
,
G
k
(
h
-
n
k
)
(
z
k
)
(
h
-
n
k
)
!
+
T
(
h
)
(
z
k
)
h
!
(
h
=
n
k
:
E
)
.
As we did in Section 2.2.2(1), to further simply the implementation, we can also introduce an auxiliary variable e
~
h which satisfies
(42)
e
~
h
=
e
h
-
T
(
h
)
(
z
k
)
h
!
.
And then using the procedure that we proposed in Section 2.1 leads to
(43)
e
~
h
=
{
0
(
h
=
0
:
n
k
-
1
)
,
G
k
(
z
k
)
(
h
=
n
k
)
,
1
h
-
n
k
∑
i
=
1
h
e
~
h
-
i
λ
~
i
|
s
c
=
z
k
(
h
=
n
k
+
1
:
E
-
1
)
,
where λ
~
i is as defined in (40b). It can be easily observed that the larger n
k is, the less computation load the method will involve. Thus s
c can be chosen as z
k
=
max
⁡
{
z
1
,
z
2
,
…
,
z
N
}. Particularly, if n
k
≥
E, we have
(44)
e
h
=
T
(
h
)
(
z
k
)
h
!
.(3)
s
c
=
s
k
. Last, let us discuss the condition when s
c
=
s
k. Here s
k is the kth pole of R
(
s
). Multiplying both sides of (2) with (
s
-
s
k
)
m
k gives
(45)
R
k
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
k
)
h
+
m
k
-
(
s
-
s
k
)
m
k
T
(
s
)
,
where R
k
(
s
) and T
(
s
) are as defined in (4) and (32), respectively. According to (45), we have
(46)
e
h
=
R
k
(
h
+
m
k
)
(
s
)
+
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
(
h
+
m
k
)
!
|
s
=
s
k
.
Based on Leibniz’s formula, we can obtain
(47)
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
|
s
=
s
k
=
P
h
+
m
k
m
k
T
k
(
h
)
|
s
=
s
k
,
where
(48)
T
k
(
h
)
=
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
i
j
(
s
-
s
i
)
j
+
h
.
Let
(49)
e
~
h
=
R
k
(
h
+
m
k
)
(
s
)
|
s
=
s
k
(
h
+
m
k
)
!
.
Then according to (46), (47), and (49), we have
(50)
e
h
=
e
~
h
+
T
k
(
h
)
h
!
|
s
=
s
k
.
Notice that (49) has the similar formulation as (3) so that it can also be tackled using the procedure we proposed in Section 2.1. We can finally obtain
(51)
e
~
h
=
1
(
h
+
m
k
)
(
∑
i
=
1
h
e
~
(
h
-
i
)
λ
^
i
+
∑
i
=
h
+
1
m
k
+
h
c
k
(
i
-
h
)
λ
^
i
)
|
s
=
s
k
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
0
:
E
-
1
)
,
where λ
^
i are as defined in (11b) and c
k
(
i
-
h
) are residues related to the kth pole. It is worthy to mention that we do not have to calculate all the λ
^
i here if we stored them when calculating the residues in Section 2.1.Till this point, the method for the computation ofe
h is discussed for any s
c. Generally, it is more preferable to let s
c
=
z
k in practice as this condition requires the least computation costs.


## 2.1. Computing Residues of Factorized Functions

Brugia [4] developed a classic pfe method for factorized rational functions. His method is more efficient and easier to apply for both manual and computer calculation compared with many other methods. Here we are to further develop and simplify Brugia’s method. To allow for context and illumination of differences and to make this paper self-contained, part of Brugia’s method is briefly introduced below. According to Heaviside’s formula, we have
(3)
c
k
L
=
R
k
(
K
)
|
s
=
s
k
K
!
,
where K
=
m
k
-
L and
(4)
R
k
(
s
)
=
∏
j
=
1
N
(
s
-
z
j
)
n
j
∏
j
=
1
,
j
≠
k
M
(
s
-
s
j
)
m
j
.
As shown in [4], the Kth derivatives of R
k
(
s
) are
(5)
R
k
(
K
)
=
(
R
k
D
k
)
(
K
-
1
)
=
∑
i
=
0
K
-
1
C
K
-
1
i
R
k
(
K
-
i
-
1
)
D
k
(
i
)
,
where K
>
0 and
(6)
D
k
(
i
)
=
(
-
1
)
i
i
!
(
∑
j
=
1
N
n
j
(
s
-
z
j
)
i
+
1
-
∑
j
=
1
j
≠
k
M
m
j
(
s
-
s
j
)
i
+
1
)
.
Here C
m
i
=
m
!
/
(
m
-
i
)
!
i
!, i
=
i
(
i
-
1
)
⋯
1, and 0
!
=
1. In the above equations, we denote R
k
(
s
) by R
k for simplicity of expression. Similar notations are also used for other functions in the following part. This is Brugia’s main result; more details can be found in [4]. Equations (5) and (6) can generate a linear system of equations and then Gramer’s method can used to solve the equations to obtain the derivatives of R
k
(
s
). It represents a very efficient pfe method and the computer implementation is easier to be accomplished compared with many other existing “tricky” methods. Here we propose a simpler method based on (3), (5), and (6). According to (3), we can obtain
(7)
R
k
(
i
)
|
s
=
s
k
=
i
!
c
k
(
m
k
-
i
)
.
Making some mathematical transformations on (6) yields
(8)
D
k
(
i
)
=
i
!
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
+
1
)
.
Here we mention that the modification on (6) to obtain (8) is relatively little, but in practice it can reduce computation costs and facilitate implementation to a considerable extent. Combining (5) and (3), we have
(9)
c
k
L
=
1
K
!
∑
i
=
0
K
-
1
C
K
-
1
i
R
k
(
K
-
i
-
1
)
D
k
(
i
)
|
s
=
s
k
.
Using (7) and (8) to substitute the derivatives of R
k and D
k in the summation of (9) yields
(10)
c
k
L
=
1
K
∑
i
=
0
K
-
1
c
k
(
m
k
-
(
K
-
i
-
1
)
)
×
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
+
1
)
|
s
=
s
k
.
Notice K
=
m
k
-
L. Equation (10) further yields(11a)
c
k
L
=
1
m
k
-
L
∑
i
=
1
m
k
-
L
c
k
(
L
+
i
)
λ
^
i
|
s
=
s
k
,
(
L
=
m
k
-
1
:
-
1
:
1
)
,
where
(11b)
λ
^
i
=
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
)
i
)
.Here we mention that in (11a) L
=
m
k
-
1
:
-
1
:
1 means L
=
m
k
-
1
,
m
k
-
2
,
m
k
-
3
,
…
,
2
,
1. This is a Matlab denotation which will also be used in the following part for the simplicity of expression. Obviously c
k
m
k
=
R
k
|
s
=
s
k. By using (11a) and (11b) and decreasing progressively the value of L from m
k
-
1 to 1, we can then derive other residues successively. Equations (11a) and (11b) can be implemented recursively in a straightforward way. It can be achieved without the aid of a comparatively complex table as done in [3]. It is much easier to apply for manual calculation than the Brugia’s method [4] and its modifications [3, 5]. Meanwhile, it can be programed more easily and reduce calculation costs.


## 2.2. Computing Residual Polynomial Coefficients of Factorized Functions

Generally, as practiced by most current articles, long division or polynomial division can be used to compute the coefficients of the residual polynomial, namely,e
h, in (2). However, as errors will propagate at each step of the polynomial division, this scheme is not quite suitable for large-scale problems. Moreover, for a given factorized rational function, the long division requires both the numerator and denominator be rewritten in expanded form, which can lead to additional errors and much computation costs. Here we provide two novel algebraic methods to compute e
h that avoid long division and the reformation of R
(
s
). The first method described in Section 2.2.1 is obtained by further developing Linnér’s method. The other method introduced in Section 2.2.2 is based on derivatives.


### 2.2.1. Laurent Expansion Applied to Improper Rational Functions for Anys
c

Brugia’s method is limited to proper functions. In [3], Linnér proposed an excellent procedure to compute the residual polynomial coefficient of improper functions by Laurent expansion. Here we will further develop Linnér’s method and obtain much simpler and more efficient calculation formulas. Linnér’s method demands that s
c in (2) not be equal to zeros or poles of R
(
s
), which can cause inconvenience in practice. We will extend s
c to any value. Linnér’s method is briefly introduced to allow for context, completeness, and illumination of differences. With the substitution s
-
s
c
=
1
/
t, (1) transforms into
(12)
R
(
t
)
=
t
-
E
W
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
,
where z
~
j
=
1
/
(
z
j
-
s
c
), s
~
j
=
1
/
(
s
j
-
s
c
), and W is a constant. Now write
(13)
R
(
t
)
=
t
-
E
R
~
(
t
)
,
where R
~
(
t
) does not include the factor t. In (13), applying a Taylor expansion to R
~
(
t
) at t
=
0, we can obtain
(14)
R
(
s
)
=
∑
h
=
0
H
R
~
(
h
)
(
0
)
h
!
(
s
-
s
c
)
E
-
h
+
o
[
(
s
-
s
c
)
E
-
H
]
,
where H is an integer. Now expand R
(
s
) as
(15)
R
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
c
)
h
+
R
2
(
s
)
,
where R
2
(
s
) is the proper fraction, which can be expanded as the second part on the right-side of (2). Comparing (14) and (15), we have
(16)
e
h
=
R
~
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
This is Linnér’s main result. In the method, the derivatives of R
~
(
t
) are calculated using equations like (5) and (6) and it also requires s
c
≠
z
j and s
c
≠
s
j, which can cause inconvenience in practice. In the following part, we first show how to extend s
c to zeros and poles of R
(
s
). We then obtain much simpler recursive formula for the calculation of e
h.Supposes
c is the kth root of R
(
s
), s
k. We can rewrite
(17)
R
(
s
)
=
(
s
-
s
k
)
-
m
k
R
k
(
s
)
,
where R
k
(
s
) is as defined in (4). With the substitution s
-
s
k
=
1
/
t in R
k
(
s
), we have
(18)
R
k
(
t
)
=
t
-
m
k
-
E
R
~
k
(
t
)
,
where
(19)
R
~
k
(
t
)
=
W
k
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
,
j
≠
k
M
(
t
-
s
~
j
)
m
j
.
Here z
~
j
=
1
/
(
z
j
-
s
k
), s
~
j
=
1
/
(
s
j
-
s
k
), and
(20)
W
k
=
∏
j
=
1
N
(
s
k
-
z
j
)
n
j
∏
j
=
1
,
j
≠
k
M
(
s
k
-
s
j
)
m
j
.
By Taylor expanding R
~
k
(
t
) at t
=
0 in (18), we can obtain
(21)
R
k
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
m
k
+
o
[
(
s
-
s
k
)
E
-
H
+
m
k
]
.
From (17) and (21), we have
(22)
R
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
o
[
(
s
-
s
k
)
E
-
H
]
.
Comparing (15) and (22), we can observe that
(23)
e
h
=
R
~
k
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
Thus s
c is extended to the kth root of R
(
s
), s
k. The same procedure can be used to extend s
c to the kth zero of R
(
s
), z
k.We can easily obtaine
E by setting h
=
E in (23) or by inspection:
(24)
e
E
=
1
.
In Section 2.1, we have simplified Brugia’s method and derived simple recursive formula for the calculation of the residues. Notice that, similar to R
k
(
s
), R
~
k
(
t
) and R
~
(
t
) are both of factorized form and that (3) and (23), which calculate c
i
j and e
h, respectively, are of the similar form. Therefore, the aforeproposed procedure in Section 2.1 for the calculation of c
i
j is also suitable for the calculation of e
h here. Take s
c
=
s
k, for an instance. It can be proven that
(25)
R
~
k
(
K
)
=
(
R
~
k
D
~
k
)
(
K
-
1
)
=
∑
i
=
0
K
-
1
C
K
-
1
i
R
~
k
(
K
-
i
-
1
)
D
~
k
(
i
)
,
where
(26)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
j
≠
k
M
m
j
(
s
~
j
-
t
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
~
j
-
t
)
i
+
1
)
.
From (23), we have
(27)
R
~
k
(
i
)
(
0
)
=
i
!
e
E
-
i
.
According to (23) and (25), we have
(28)
e
h
=
1
(
E
-
h
)
!
∑
i
=
0
E
-
h
-
1
C
E
-
h
-
1
i
R
~
k
(
E
-
h
-
i
-
1
)
D
~
k
(
i
)
|
t
=
0
.
Using (26) and (27) to substitute the derivatives of D
~
k and R
~
k in (28) and then evaluating (28) at t
=
0, we can obtain
(29)
e
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
k
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
k
)
i
)
,
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
E
-
1
:
-
1
:
0
)
.
The same procedure can be used to compute e
h when s
c
≠
s
k. In summary, e
h can be calculated according to(30a)
e
h
=
{
1
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
η
i
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(30b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
z
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
otherwise
)
.Equations (30a) and (30b) are such a simple and compact formula that its implementation is even very easy for hand calculation. It is also much more convenient for computer program and requires less computation costs than the original method in [3]. We just need first calculate η
i using (30b). Then e
h can be immediately obtained using (30a).


### 2.2.2. Computinge
h through Derivatives

We provide another simple way to computee
h which also avoids the unfavorable polynomial division. That e
E
=
1 is obvious. Therefore, we just need to obtain the remaining coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1. Notice that as the residues c
k
L can be firstly obtained in Section 2.1, they can be seen as known quantities here. Three conditions are considered as below regarding the different values of s
c.(1)
s
c
≠
z
k
and
s
c
≠
s
k
. We first discuss the calculation of e
h when s
c
≠
z
k and s
c
≠
s
k. Here z
k and s
k are the kth zero and the kth pole, respectively. From (1) and (2), we observe that the coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1 can be obtained according to
(31)
e
h
=
1
h
!
(
R
(
h
)
+
T
(
h
)
)
|
s
=
s
c
,
where
(32)
T
(
s
)
=
-
∑
k
=
1
M
∑
L
=
1
m
k
c
k
L
(
s
-
s
k
)
L
.
Here R
=
R
(
s
) as defined in (1) and T
=
T
(
s
). The hth derivative of T
(
s
) can be easily obtained as
(33)
T
(
h
)
=
(
-
1
)
h
-
1
∑
k
=
1
M
∑
j
=
1
m
k
P
h
+
j
-
1
h
c
k
j
(
s
-
s
k
)
j
+
h
,
where P
m
n
=
m
!
/
(
m
-
n
)
!. As R
(
s
) is of factorized form, its derivatives can be derived using functions like (4) and (5). Thus the problem is solved. To further simplify the implementation, we introduce an auxiliary variable
(34)
e
~
h
=
1
h
!
R
(
h
)
|
s
=
s
c
.
From (31) and (34), we can easily observe that
(35)
e
h
=
e
~
h
+
1
h
!
T
(
h
)
|
s
=
s
c
.
Equation (34) has a similar formulation as (3); thus, e
~
h can be calculated and using the procedure we compute c
i
j in Section 2.1. Accordingly, we can finally obtain(36a)
e
~
h
=
{
R
(
s
c
)
(
h
=
0
)
,
1
h
∑
i
=
1
h
e
~
h
-
i
λ
i
(
h
=
1
:
E
)
,
where
(36b)
λ
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
.Here we mention that in (36a) h
=
1
:
E means h
=
1
,
2
,
3
,
…
,
E. This is a Matlab denotation which will also be used in the remaining part for simplicity of expression. We can first obtain e
~
h using (36a) and (36b); then, e
h can be easily obtained using (35). Notice the evaluation of (36b) requires s
c
≠
z
k and s
c
≠
s
k. Thus other two conditions remain to be discussed.(2)
s
c
=
z
k
. Here we discuss the condition when s
c is one of the zeros of R
(
s
) (Let us say z
k). Actually it will be seen that it is more preferable to let s
c
=
z
k. First we rewrite R
(
s
) as
(37)
R
(
s
)
=
(
s
-
z
k
)
n
k
G
k
(
s
)
,
where n
k is the multiplicity of z
k and
(38)
G
k
(
s
)
=
∏
j
=
1
,
j
≠
k
N
(
s
-
z
j
)
n
j
∏
j
=
1
M
(
s
-
s
j
)
m
j
.
Computing the hth derivative of both sides of (37) using Leibniz’s formula and then setting s as z
k yield
(39)
R
(
h
)
(
z
k
)
=
{
0
(
h
<
n
k
)
P
h
n
k
G
k
(
h
-
n
k
)
(
z
k
)
(
h
≥
n
k
)
.
As G
k
(
s
) is of factorized form, its derivatives can be computed recursively(40a)
G
k
(
h
-
n
k
)
=
{
G
k
(
h
=
n
k
)
,
∑
i
=
1
h
-
n
k
P
h
-
n
k
-
1
i
-
1
D
~
k
(
h
-
n
k
-
i
)
λ
~
i
(
h
>
n
k
)
,
where
(40b)
λ
~
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
,
(40c)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
+
1
)
.Combining (31) and (39), we have
(41)
e
h
=
{
T
(
h
)
(
z
k
)
h
!
(
h
=
0
:
n
k
-
1
)
,
G
k
(
h
-
n
k
)
(
z
k
)
(
h
-
n
k
)
!
+
T
(
h
)
(
z
k
)
h
!
(
h
=
n
k
:
E
)
.
As we did in Section 2.2.2(1), to further simply the implementation, we can also introduce an auxiliary variable e
~
h which satisfies
(42)
e
~
h
=
e
h
-
T
(
h
)
(
z
k
)
h
!
.
And then using the procedure that we proposed in Section 2.1 leads to
(43)
e
~
h
=
{
0
(
h
=
0
:
n
k
-
1
)
,
G
k
(
z
k
)
(
h
=
n
k
)
,
1
h
-
n
k
∑
i
=
1
h
e
~
h
-
i
λ
~
i
|
s
c
=
z
k
(
h
=
n
k
+
1
:
E
-
1
)
,
where λ
~
i is as defined in (40b). It can be easily observed that the larger n
k is, the less computation load the method will involve. Thus s
c can be chosen as z
k
=
max
⁡
{
z
1
,
z
2
,
…
,
z
N
}. Particularly, if n
k
≥
E, we have
(44)
e
h
=
T
(
h
)
(
z
k
)
h
!
.(3)
s
c
=
s
k
. Last, let us discuss the condition when s
c
=
s
k. Here s
k is the kth pole of R
(
s
). Multiplying both sides of (2) with (
s
-
s
k
)
m
k gives
(45)
R
k
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
k
)
h
+
m
k
-
(
s
-
s
k
)
m
k
T
(
s
)
,
where R
k
(
s
) and T
(
s
) are as defined in (4) and (32), respectively. According to (45), we have
(46)
e
h
=
R
k
(
h
+
m
k
)
(
s
)
+
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
(
h
+
m
k
)
!
|
s
=
s
k
.
Based on Leibniz’s formula, we can obtain
(47)
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
|
s
=
s
k
=
P
h
+
m
k
m
k
T
k
(
h
)
|
s
=
s
k
,
where
(48)
T
k
(
h
)
=
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
i
j
(
s
-
s
i
)
j
+
h
.
Let
(49)
e
~
h
=
R
k
(
h
+
m
k
)
(
s
)
|
s
=
s
k
(
h
+
m
k
)
!
.
Then according to (46), (47), and (49), we have
(50)
e
h
=
e
~
h
+
T
k
(
h
)
h
!
|
s
=
s
k
.
Notice that (49) has the similar formulation as (3) so that it can also be tackled using the procedure we proposed in Section 2.1. We can finally obtain
(51)
e
~
h
=
1
(
h
+
m
k
)
(
∑
i
=
1
h
e
~
(
h
-
i
)
λ
^
i
+
∑
i
=
h
+
1
m
k
+
h
c
k
(
i
-
h
)
λ
^
i
)
|
s
=
s
k
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
0
:
E
-
1
)
,
where λ
^
i are as defined in (11b) and c
k
(
i
-
h
) are residues related to the kth pole. It is worthy to mention that we do not have to calculate all the λ
^
i here if we stored them when calculating the residues in Section 2.1.Till this point, the method for the computation ofe
h is discussed for any s
c. Generally, it is more preferable to let s
c
=
z
k in practice as this condition requires the least computation costs.


## 2.2.1. Laurent Expansion Applied to Improper Rational Functions for Anys
c

Brugia’s method is limited to proper functions. In [3], Linnér proposed an excellent procedure to compute the residual polynomial coefficient of improper functions by Laurent expansion. Here we will further develop Linnér’s method and obtain much simpler and more efficient calculation formulas. Linnér’s method demands that s
c in (2) not be equal to zeros or poles of R
(
s
), which can cause inconvenience in practice. We will extend s
c to any value. Linnér’s method is briefly introduced to allow for context, completeness, and illumination of differences. With the substitution s
-
s
c
=
1
/
t, (1) transforms into
(12)
R
(
t
)
=
t
-
E
W
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
,
where z
~
j
=
1
/
(
z
j
-
s
c
), s
~
j
=
1
/
(
s
j
-
s
c
), and W is a constant. Now write
(13)
R
(
t
)
=
t
-
E
R
~
(
t
)
,
where R
~
(
t
) does not include the factor t. In (13), applying a Taylor expansion to R
~
(
t
) at t
=
0, we can obtain
(14)
R
(
s
)
=
∑
h
=
0
H
R
~
(
h
)
(
0
)
h
!
(
s
-
s
c
)
E
-
h
+
o
[
(
s
-
s
c
)
E
-
H
]
,
where H is an integer. Now expand R
(
s
) as
(15)
R
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
c
)
h
+
R
2
(
s
)
,
where R
2
(
s
) is the proper fraction, which can be expanded as the second part on the right-side of (2). Comparing (14) and (15), we have
(16)
e
h
=
R
~
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
This is Linnér’s main result. In the method, the derivatives of R
~
(
t
) are calculated using equations like (5) and (6) and it also requires s
c
≠
z
j and s
c
≠
s
j, which can cause inconvenience in practice. In the following part, we first show how to extend s
c to zeros and poles of R
(
s
). We then obtain much simpler recursive formula for the calculation of e
h.Supposes
c is the kth root of R
(
s
), s
k. We can rewrite
(17)
R
(
s
)
=
(
s
-
s
k
)
-
m
k
R
k
(
s
)
,
where R
k
(
s
) is as defined in (4). With the substitution s
-
s
k
=
1
/
t in R
k
(
s
), we have
(18)
R
k
(
t
)
=
t
-
m
k
-
E
R
~
k
(
t
)
,
where
(19)
R
~
k
(
t
)
=
W
k
∏
j
=
1
N
(
t
-
z
~
j
)
n
i
∏
j
=
1
,
j
≠
k
M
(
t
-
s
~
j
)
m
j
.
Here z
~
j
=
1
/
(
z
j
-
s
k
), s
~
j
=
1
/
(
s
j
-
s
k
), and
(20)
W
k
=
∏
j
=
1
N
(
s
k
-
z
j
)
n
j
∏
j
=
1
,
j
≠
k
M
(
s
k
-
s
j
)
m
j
.
By Taylor expanding R
~
k
(
t
) at t
=
0 in (18), we can obtain
(21)
R
k
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
m
k
+
o
[
(
s
-
s
k
)
E
-
H
+
m
k
]
.
From (17) and (21), we have
(22)
R
(
s
)
=
∑
h
=
0
H
R
~
k
(
h
)
(
0
)
h
!
(
s
-
s
k
)
E
-
h
+
o
[
(
s
-
s
k
)
E
-
H
]
.
Comparing (15) and (22), we can observe that
(23)
e
h
=
R
~
k
(
E
-
h
)
(
0
)
(
E
-
h
)
!
.
Thus s
c is extended to the kth root of R
(
s
), s
k. The same procedure can be used to extend s
c to the kth zero of R
(
s
), z
k.We can easily obtaine
E by setting h
=
E in (23) or by inspection:
(24)
e
E
=
1
.
In Section 2.1, we have simplified Brugia’s method and derived simple recursive formula for the calculation of the residues. Notice that, similar to R
k
(
s
), R
~
k
(
t
) and R
~
(
t
) are both of factorized form and that (3) and (23), which calculate c
i
j and e
h, respectively, are of the similar form. Therefore, the aforeproposed procedure in Section 2.1 for the calculation of c
i
j is also suitable for the calculation of e
h here. Take s
c
=
s
k, for an instance. It can be proven that
(25)
R
~
k
(
K
)
=
(
R
~
k
D
~
k
)
(
K
-
1
)
=
∑
i
=
0
K
-
1
C
K
-
1
i
R
~
k
(
K
-
i
-
1
)
D
~
k
(
i
)
,
where
(26)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
j
≠
k
M
m
j
(
s
~
j
-
t
)
i
+
1
-
∑
j
=
1
N
n
j
(
z
~
j
-
t
)
i
+
1
)
.
From (23), we have
(27)
R
~
k
(
i
)
(
0
)
=
i
!
e
E
-
i
.
According to (23) and (25), we have
(28)
e
h
=
1
(
E
-
h
)
!
∑
i
=
0
E
-
h
-
1
C
E
-
h
-
1
i
R
~
k
(
E
-
h
-
i
-
1
)
D
~
k
(
i
)
|
t
=
0
.
Using (26) and (27) to substitute the derivatives of D
~
k and R
~
k in (28) and then evaluating (28) at t
=
0, we can obtain
(29)
e
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
k
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
k
)
i
)
,
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
E
-
1
:
-
1
:
0
)
.
The same procedure can be used to compute e
h when s
c
≠
s
k. In summary, e
h can be calculated according to(30a)
e
h
=
{
1
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
η
i
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(30b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
z
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
otherwise
)
.Equations (30a) and (30b) are such a simple and compact formula that its implementation is even very easy for hand calculation. It is also much more convenient for computer program and requires less computation costs than the original method in [3]. We just need first calculate η
i using (30b). Then e
h can be immediately obtained using (30a).


## 2.2.2. Computinge
h through Derivatives

We provide another simple way to computee
h which also avoids the unfavorable polynomial division. That e
E
=
1 is obvious. Therefore, we just need to obtain the remaining coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1. Notice that as the residues c
k
L can be firstly obtained in Section 2.1, they can be seen as known quantities here. Three conditions are considered as below regarding the different values of s
c.(1)
s
c
≠
z
k
and
s
c
≠
s
k
. We first discuss the calculation of e
h when s
c
≠
z
k and s
c
≠
s
k. Here z
k and s
k are the kth zero and the kth pole, respectively. From (1) and (2), we observe that the coefficients e
0
,
e
1
,
e
2
,
…
,
e
E
-
1 can be obtained according to
(31)
e
h
=
1
h
!
(
R
(
h
)
+
T
(
h
)
)
|
s
=
s
c
,
where
(32)
T
(
s
)
=
-
∑
k
=
1
M
∑
L
=
1
m
k
c
k
L
(
s
-
s
k
)
L
.
Here R
=
R
(
s
) as defined in (1) and T
=
T
(
s
). The hth derivative of T
(
s
) can be easily obtained as
(33)
T
(
h
)
=
(
-
1
)
h
-
1
∑
k
=
1
M
∑
j
=
1
m
k
P
h
+
j
-
1
h
c
k
j
(
s
-
s
k
)
j
+
h
,
where P
m
n
=
m
!
/
(
m
-
n
)
!. As R
(
s
) is of factorized form, its derivatives can be derived using functions like (4) and (5). Thus the problem is solved. To further simplify the implementation, we introduce an auxiliary variable
(34)
e
~
h
=
1
h
!
R
(
h
)
|
s
=
s
c
.
From (31) and (34), we can easily observe that
(35)
e
h
=
e
~
h
+
1
h
!
T
(
h
)
|
s
=
s
c
.
Equation (34) has a similar formulation as (3); thus, e
~
h can be calculated and using the procedure we compute c
i
j in Section 2.1. Accordingly, we can finally obtain(36a)
e
~
h
=
{
R
(
s
c
)
(
h
=
0
)
,
1
h
∑
i
=
1
h
e
~
h
-
i
λ
i
(
h
=
1
:
E
)
,
where
(36b)
λ
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
.Here we mention that in (36a) h
=
1
:
E means h
=
1
,
2
,
3
,
…
,
E. This is a Matlab denotation which will also be used in the remaining part for simplicity of expression. We can first obtain e
~
h using (36a) and (36b); then, e
h can be easily obtained using (35). Notice the evaluation of (36b) requires s
c
≠
z
k and s
c
≠
s
k. Thus other two conditions remain to be discussed.(2)
s
c
=
z
k
. Here we discuss the condition when s
c is one of the zeros of R
(
s
) (Let us say z
k). Actually it will be seen that it is more preferable to let s
c
=
z
k. First we rewrite R
(
s
) as
(37)
R
(
s
)
=
(
s
-
z
k
)
n
k
G
k
(
s
)
,
where n
k is the multiplicity of z
k and
(38)
G
k
(
s
)
=
∏
j
=
1
,
j
≠
k
N
(
s
-
z
j
)
n
j
∏
j
=
1
M
(
s
-
s
j
)
m
j
.
Computing the hth derivative of both sides of (37) using Leibniz’s formula and then setting s as z
k yield
(39)
R
(
h
)
(
z
k
)
=
{
0
(
h
<
n
k
)
P
h
n
k
G
k
(
h
-
n
k
)
(
z
k
)
(
h
≥
n
k
)
.
As G
k
(
s
) is of factorized form, its derivatives can be computed recursively(40a)
G
k
(
h
-
n
k
)
=
{
G
k
(
h
=
n
k
)
,
∑
i
=
1
h
-
n
k
P
h
-
n
k
-
1
i
-
1
D
~
k
(
h
-
n
k
-
i
)
λ
~
i
(
h
>
n
k
)
,
where
(40b)
λ
~
i
=
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
,
(40c)
D
~
k
(
i
)
=
i
!
(
∑
j
=
1
M
m
j
(
s
j
-
s
)
i
+
1
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
)
i
+
1
)
.Combining (31) and (39), we have
(41)
e
h
=
{
T
(
h
)
(
z
k
)
h
!
(
h
=
0
:
n
k
-
1
)
,
G
k
(
h
-
n
k
)
(
z
k
)
(
h
-
n
k
)
!
+
T
(
h
)
(
z
k
)
h
!
(
h
=
n
k
:
E
)
.
As we did in Section 2.2.2(1), to further simply the implementation, we can also introduce an auxiliary variable e
~
h which satisfies
(42)
e
~
h
=
e
h
-
T
(
h
)
(
z
k
)
h
!
.
And then using the procedure that we proposed in Section 2.1 leads to
(43)
e
~
h
=
{
0
(
h
=
0
:
n
k
-
1
)
,
G
k
(
z
k
)
(
h
=
n
k
)
,
1
h
-
n
k
∑
i
=
1
h
e
~
h
-
i
λ
~
i
|
s
c
=
z
k
(
h
=
n
k
+
1
:
E
-
1
)
,
where λ
~
i is as defined in (40b). It can be easily observed that the larger n
k is, the less computation load the method will involve. Thus s
c can be chosen as z
k
=
max
⁡
{
z
1
,
z
2
,
…
,
z
N
}. Particularly, if n
k
≥
E, we have
(44)
e
h
=
T
(
h
)
(
z
k
)
h
!
.(3)
s
c
=
s
k
. Last, let us discuss the condition when s
c
=
s
k. Here s
k is the kth pole of R
(
s
). Multiplying both sides of (2) with (
s
-
s
k
)
m
k gives
(45)
R
k
(
s
)
=
∑
h
=
0
E
e
h
(
s
-
s
k
)
h
+
m
k
-
(
s
-
s
k
)
m
k
T
(
s
)
,
where R
k
(
s
) and T
(
s
) are as defined in (4) and (32), respectively. According to (45), we have
(46)
e
h
=
R
k
(
h
+
m
k
)
(
s
)
+
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
(
h
+
m
k
)
!
|
s
=
s
k
.
Based on Leibniz’s formula, we can obtain
(47)
[
(
s
-
s
k
)
m
k
T
(
s
)
]
(
h
+
m
k
)
|
s
=
s
k
=
P
h
+
m
k
m
k
T
k
(
h
)
|
s
=
s
k
,
where
(48)
T
k
(
h
)
=
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
i
j
(
s
-
s
i
)
j
+
h
.
Let
(49)
e
~
h
=
R
k
(
h
+
m
k
)
(
s
)
|
s
=
s
k
(
h
+
m
k
)
!
.
Then according to (46), (47), and (49), we have
(50)
e
h
=
e
~
h
+
T
k
(
h
)
h
!
|
s
=
s
k
.
Notice that (49) has the similar formulation as (3) so that it can also be tackled using the procedure we proposed in Section 2.1. We can finally obtain
(51)
e
~
h
=
1
(
h
+
m
k
)
(
∑
i
=
1
h
e
~
(
h
-
i
)
λ
^
i
+
∑
i
=
h
+
1
m
k
+
h
c
k
(
i
-
h
)
λ
^
i
)
|
s
=
s
k
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
(
h
=
0
:
E
-
1
)
,
where λ
^
i are as defined in (11b) and c
k
(
i
-
h
) are residues related to the kth pole. It is worthy to mention that we do not have to calculate all the λ
^
i here if we stored them when calculating the residues in Section 2.1.Till this point, the method for the computation ofe
h is discussed for any s
c. Generally, it is more preferable to let s
c
=
z
k in practice as this condition requires the least computation costs.


## 3. Partial Fraction Expansion of Rational Functions in Expanded Form

Let the rational function in expanded form be(52)
R
(
s
)
=
Q
(
s
)
P
(
s
)
=
∑
k
=
0
N
a
k
(
s
-
s
0
)
k
∏
i
=
1
M
(
s
-
s
i
)
m
i
=
∑
h
=
0
E
e
h
(
s
-
s
c
)
h
+
∑
i
=
1
M
∑
j
=
1
m
i
c
i
j
(
s
-
s
i
)
j
,
where c
i
j, e
h are the residues and coefficients of residual polynomial, respectively; M is the number of poles; s
i are the poles of R
(
s
); m
i are the multiplicities of s
i; a
k are the polynomial coefficients of the numerator; s
0 and s
c are constants which are often zeros in practice and they are introduced to make the polynomials more general. Obviously, the degree of the numerator is N, and the degree of the denominator is K
=
∑
m
j, E
=
N
-
K. For a proper function, e
h are zeros. It can be observed that R
(
s
) can also be expressed as
(53)
R
(
s
)
=
∑
k
=
0
N
a
k
r
k
(
s
)
,
where
(54)
r
k
(
s
)
=
(
s
-
s
0
)
k
∏
i
=
1
M
(
s
-
s
i
)
m
i
.
Suppose r
k
(
s
) can be expanded in pfe as
(55)
r
k
(
s
)
=
∑
h
=
0
k
-
K
e
k
h
(
s
-
s
c
)
h
+
∑
i
=
1
M
∑
j
=
1
m
i
c
k
i
j
(
s
-
s
i
)
j
.
Inserting (55) into (53) and then compare (53) with (52), we can observe that the residues of R
(
s
) can be obtained as
(56)
c
i
j
=
∑
k
=
0
N
a
k
c
k
i
j
,
and the residual polynomial coefficients can be calculated according to
(57)
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
(
K
+
h
+
i
)
h
,
(
h
=
0
:
E
)
.
In the following part, we will first show how to compute c
i
j and then to calculate e
h.


### 3.1. Computing Residues of Functions in Expanded Form

The method for computation of residues can be divided into two steps: first, calculate the residues ofr
0
(
s
); second, calculate the residues of R
(
s
). First, let us suppose r
0
(
s
) can be expanded in pfe as
(58)
r
0
(
s
)
=
∑
i
=
1
M
∑
j
=
1
m
i
c
0
i
j
(
s
-
s
i
)
j
.
Notice that r
0
(
s
) is a function of factorized form. Thus it can be expanded using the method we proposed in Section 2.1. We just need to set the numerator as 1 in (1); namely, all the n
i
=
0 in (11b). Referring to (11a) and (11b), the following formula is immediately obtained(59a)
c
0
i
j
=
1
m
i
-
j
∑
k
=
1
m
i
-
j
c
0
i
(
j
+
k
)
λ
k
|
s
=
s
i
,
(
j
=
m
i
-
1
:
-
1
:
1
)
,
where
(59b)
λ
k
=
∑
L
=
1
L
≠
i
M
m
L
(
s
L
-
s
)
k
.It has been proved in [17] that the residues of r
k
(
s
) and those of r
k
-
1
(
s
) satisfy
(60)
c
k
i
j
=
c
(
k
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
c
(
k
-
1
)
i
j
.
Based on the above formula, all the residues of r
k
(
s
) can be obtained. Then using (56), the residues of R
(
s
) can be obtained. Here we develop novel simpler formulas to derive the residues of R
(
s
) based on the residues of r
0
(
s
). It can be observed that the numerator of R
(
s
) can be rewritten as
(61)
Q
(
s
)
=
a
0
+
(
s
-
s
0
)
(
(
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
)
a
1
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
a
2
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
a
N
-
2
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
)
)
.
According to (52) and (61), we can denote(62a)
f
0
=
a
N
P
(
s
)
,
(62b)
f
L
=
a
N
-
L
P
(
s
)
+
(
s
-
s
0
)
f
L
-
1
,
(
L
=
1
:
N
)
.By the above denotation, we have f
N
=
R
(
s
). Suppose f
L can be expanded in pfe as
(63)
f
L
=
∑
k
=
0
L
-
Κ
e
L
k
(
s
-
s
c
)
k
+
∑
i
=
1
M
∑
j
=
1
m
i
d
L
i
j
(
s
-
s
i
)
j
.
Notice that r
0
(
s
)
=
1
/
P
(
s
). Thus according to (58), we easily obtain the residues of f
0
(64a)
d
0
i
j
=
a
N
×
c
0
i
j
.
If we have already obtained the residues of f
L
-
1, d
(
L
-
1
)
i
j, referring to [17], the residues of (
s
-
s
0
)
f
L
-
1 can be obtained as d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j. Therefore, according to (62b) we have
(64b)
d
L
i
j
=
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
,
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
(
L
=
1
:
N
)
.Notice that c
i
j
=
d
N
i
j. As c
0
i
j can already be calculated using (59a) and (59b), we can finally obtain the residues of R
(
s
) using formulas (64a) and (64b) successively by increasing the value of L from 1 to N. Equations (64a) and (64b) can be further developed and yield another beautiful formula:
(65)
d
L
i
j
=
∑
m
=
0
L
a
N
-
(
L
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
In (64a), (64b), and (65), if j
+
1
>
m
i or j
+
m
-
k
>
m
i, the residues c
0
i
(
j
+
1
) and c
0
i
(
j
+
m
-
k
) do not exist and we can simply set them as zero. Formula (65) can be proved by induction as below.Proof of Formula (65).
WhenL
=
0, (65) is reduced to (64a). Thus it is true. Suppose (65) is true when L
=
n; namely,
(66)
d
n
i
j
=
∑
m
=
0
n
a
N
-
(
n
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
Then according to (64b), when L
=
n
+
1, we have
(67)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
d
n
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
n
i
j
.
Applying (66) to (67) yields
(68)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
∑
m
=
0
n
a
N
-
(
n
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
1
+
m
-
k
)
+
(
s
i
-
s
0
)
×
∑
m
=
0
n
a
N
-
(
n
-
m
)
×
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
From (68), we have
(69)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
Λ
,
where(70a)
Λ
=
∑
m
=
0
n
a
N
-
(
n
-
m
)
×
{
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
1
+
m
-
k
)
a
a
a
A
a
a
+
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
+
1
c
0
i
(
j
+
m
-
k
)
}
.
Equation (70a) can further lead to
(70b)
Λ
=
∑
m
=
1
n
+
1
a
N
-
(
n
+
1
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.From (69) and (70b), we can finally obtain
(71)
d
(
n
+
1
)
i
j
=
∑
m
=
0
n
+
1
a
N
-
(
n
+
1
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
Equation (71) is (65) with L replaced by n
+
1; thus, by induction, it follows that (61) is true.SettingL
=
N in (65) yields
(72)
c
i
j
=
∑
m
=
0
N
a
m
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
[
j
+
m
-
k
]
.
Equation (72) presents a direct relation between the residues of r
0
(
s
) and the residues of R
(
s
). It can be useful when calculating a particular residue of R
(
s
).


### 3.2. Computing Residual Polynomial Coefficients of Functions in Expanded Form

We will propose novel methods to calculate the residual polynomial coefficients, namely,e
h in (52). Most existing methods calculate such coefficients using polynomial division. An algebraic procedure to calculate such coefficients that avoids polynomial division can be found in [17]. Here we propose two novel simpler methods which also avoid polynomial division. We first propose a method through Laurent expansion and then propose a method through derivatives.


#### 3.2.1. Computation ofe
h through Laurent Expansion

As shown in (53), R
(
s
) is a sum of r
k
(
s
), if we can obtain the residual polynomial coefficients of r
k
(
s
), e
k
h, e
h can be then obtained using (57). Obviously, r
k
(
s
) will contribute to the e
h of R
(
s
) only when k
≥
K
=
∑
m
j. Thus k is assumed to be no less than K in this part. The problem is then reduced to how to calculate the e
k
h of r
k
(
s
) (k
=
K
:
N). Notice all the r
k
(
s
) are in factorized form; thus, they can be expanded using the method we proposed in Section 2. However, it is not advisable to do so, for this can lead to cumbersome implementation. Actually, simpler methods can be obtained. Referring to Section 2.2.1, we can prove that
(73)
e
k
h
=
r
~
(
k
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
-
K
)
,
where
(74)
r
~
(
t
)
=
W
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
.
Here s
~
j
=
1
/
(
s
j
-
s
c
) and W is a constant. Let k
=
k
+
1 in (73), we have
(75)
e
(
k
+
1
)
h
=
r
~
(
k
+
1
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
+
1
-
K
)
.
From (73) and (75), we can observe that
(76)
e
k
h
=
e
(
k
+
1
)
(
h
+
1
)
,
(
h
=
0
:
k
-
K
)
.
Notice that e
(
k
+
1
)
(
h
+
1
) are the residual polynomial coefficients of r
k
+
1
(
s
). It is clear that if we have already obtained the coefficients of r
k
+
1
(
s
), we can immediately obtain the coefficients of r
k
(
s
) by (76). It hereby follows that if we have already obtained the coefficients of r
N
(
s
), e
N
h, we can successively obtain e
(
N
-
1
)
h
,
e
(
N
-
2
)
h
,
…
,
e
K
h. From (76), the following formula can be obtained:
(77)
e
(
K
+
i
)
h
=
e
N
(
E
-
i
+
h
)
,
(
i
=
0
:
E
,
h
=
0
:
i
)
.
Equation (77) presents the relationship between the e
k
h (k
<
N) and e
N
h, which can be seen more clearly in Table 1. Here we mention that we calculate the e
k
h with k decreasing from N to K, while the method in [17] calculates e
k
h with k increasing from K to N. It will be seen that the method proposed here is simpler and more convenient.Table 1
The relationship between thee
k
h (k
<
N) and e
N
h.
e
k
h
h
0
1
2
⋯
E
-
2
E
-
1
E
k
N
e
N
0
e
N
1
e
N
2
⋯
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
N
-
1
e
N
1
e
N
2
e
N
3
⋯
e
N
(
E
-
1
)
e
N
E
N
-
2
e
N
2
e
N
2
e
N
4
⋯
e
N
E
⋯
⋯
⋯
⋯
⋯
K
+
2
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
K
+
1
e
N
(
E
-
1
)
e
N
E
K
e
N
EAccording to (57) and (77), we have
(78)
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
N
(
E
-
i
)
,
(
h
=
0
:
E
)
.
From the above analysis, the whole problem is reduced to how to obtain e
N
h. As r
N
(
s
) is of factorized form, e
N
h can be calculated using the methods we proposed in Section 2.2.1 through Laurent expansion. Thus, we have(79a)
e
N
h
=
{
1
,
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
N
(
h
+
i
)
η
i
,
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(79b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
,
(
s
c
=
s
0
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
otherwise
)
.In summary, the residual polynomial coefficients can be obtained by first using (79a) and (79b) and then using (78).


#### 3.2.2. Computation ofe
h through Derivatives

As we discussed in Section3.2.1, we just have to calculate the residual polynomial coefficients of r
N
(
s
), e
N
h and then we can obtain e
h. As r
N
(
s
) is of factorized form, it is also possible to use the method we proposed in Section 2.2.2 to calculate e
N
h. In Section 2.2.2, three conditions are considered. Here we only discuss the most favorable condition, s
c
=
s
0 (namely, s
c is a zero of r
N
(
s
)) because firstly, this condition is most suitable as it requires the least computation and secondly, s
c
=
s
0
=
0 in most practical cases. Notice that N
>
K. Thus we can use a formula like (44) to calculate e
N
h when s
0 is not equal to any pole of r
N
(
s
). When s
0 is equal to a pole of r
N
(
s
), we just need to make a simple modification. e
N
E
=
1 is obvious. Using the procedures in Section 2.2.2(2), we can finally have
(80)
e
N
h
=
T
N
(
h
)
(
s
)
h
!
|
s
=
s
c
,
(
h
=
1
:
E
-
1
)
,
where
(81)
T
N
(
h
)
(
s
)
=
{
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
s
=
s
k
)
,
(
-
1
)
h
-
1
∑
i
=
1
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
otherwise
)
.


### 3.3. Recommendations on the Choices of Procedures

In Sections3.1 and 3.2, several procedures are proposed for the computation of both e
h and c
i
j. As the calculation of e
h can have a relation with that of c
i
j, we recommend two methods comprised of several aforeproposed procedures for improper rational functions. We only suggest two typical methods which we consider to be preferable and suitable as below.Method 1 uses (80) and (81) to calculate e
h. As we can see from (80) and (81), the method in Section 3.2.2 requires that we first obtain the residues of r
N
(
s
). Therefore, we use (60) (rather than (64a) and (64b)) to obtain the residues because the residues of r
N
(
s
) are obtained in the calculation process. Method 2 uses (59a), (59b), (64a), and (64b) to obtain the residues, (79a) and (79b), and (78) to calculate e
h. It is more preferable to use (79a) and (79b) to calculate e
N
h as it is not based on the residues of r
N
(
s
).Method 1. The following procedures comprise Method 1.(a)
Calculate the residues ofr
0
(
s
):
(82)
c
0
i
j
=
{
1
∏
k
=
1
,
k
≠
i
M
(
s
i
-
s
k
)
m
k
,
(
i
=
1
:
M
;
j
=
m
i
)
1
m
i
-
j
∑
k
=
1
m
i
-
j
c
0
i
(
j
+
k
)
λ
k
|
s
=
s
i
,
(
i
=
1
:
M
;
j
=
m
i
-
1
:
-
1
:
1
)
,
whereλ
k
=
∑
L
=
1
,
L
≠
i
M
(
m
L
/
(
s
L
-
s
)
k
).
(b)
Calculate the residues ofr
k
(
s
) using
(83)
c
k
i
j
=
c
(
k
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
c
(
k
-
1
)
i
j
,
aa
(
k
=
1
:
N
;
i
=
1
:
M
;
j
=
1
:
m
i
)
,
and then obtain the residues ofR
(
s
) using c
i
j
=
∑
k
=
0
N
a
k
c
k
i
j.
(c)
Calculate the residues ofr
N
(
s
) using (80) and (81).
(d)
Calculate the residual polynomial coefficients using (78).Method 2. The following procedures comprise Method 2.(a)
Calculate the residues ofr
0
(
s
) the same way as the first step of Method 1.
(b)
Calculate the residues ofR
(
s
). Here notice that c
i
j
=
d
N
i
j. Consider
(84)
d
L
i
j
=
{
a
N
×
c
0
i
j
,
(
L
=
0
)
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
,
(
L
=
1
:
N
)
.
(c)
Calculate the residual polynomial coefficients ofr
N
(
s
) using  (79a) and (79b).
(d)
Calculate the residual polynomial coefficients using (78).


## 3.1. Computing Residues of Functions in Expanded Form

The method for computation of residues can be divided into two steps: first, calculate the residues ofr
0
(
s
); second, calculate the residues of R
(
s
). First, let us suppose r
0
(
s
) can be expanded in pfe as
(58)
r
0
(
s
)
=
∑
i
=
1
M
∑
j
=
1
m
i
c
0
i
j
(
s
-
s
i
)
j
.
Notice that r
0
(
s
) is a function of factorized form. Thus it can be expanded using the method we proposed in Section 2.1. We just need to set the numerator as 1 in (1); namely, all the n
i
=
0 in (11b). Referring to (11a) and (11b), the following formula is immediately obtained(59a)
c
0
i
j
=
1
m
i
-
j
∑
k
=
1
m
i
-
j
c
0
i
(
j
+
k
)
λ
k
|
s
=
s
i
,
(
j
=
m
i
-
1
:
-
1
:
1
)
,
where
(59b)
λ
k
=
∑
L
=
1
L
≠
i
M
m
L
(
s
L
-
s
)
k
.It has been proved in [17] that the residues of r
k
(
s
) and those of r
k
-
1
(
s
) satisfy
(60)
c
k
i
j
=
c
(
k
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
c
(
k
-
1
)
i
j
.
Based on the above formula, all the residues of r
k
(
s
) can be obtained. Then using (56), the residues of R
(
s
) can be obtained. Here we develop novel simpler formulas to derive the residues of R
(
s
) based on the residues of r
0
(
s
). It can be observed that the numerator of R
(
s
) can be rewritten as
(61)
Q
(
s
)
=
a
0
+
(
s
-
s
0
)
(
(
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
)
a
1
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
a
2
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
a
N
-
2
+
(
s
-
s
0
)
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
×
(
a
N
-
1
+
a
N
(
s
-
s
0
)
)
⋯
)
)
)
.
According to (52) and (61), we can denote(62a)
f
0
=
a
N
P
(
s
)
,
(62b)
f
L
=
a
N
-
L
P
(
s
)
+
(
s
-
s
0
)
f
L
-
1
,
(
L
=
1
:
N
)
.By the above denotation, we have f
N
=
R
(
s
). Suppose f
L can be expanded in pfe as
(63)
f
L
=
∑
k
=
0
L
-
Κ
e
L
k
(
s
-
s
c
)
k
+
∑
i
=
1
M
∑
j
=
1
m
i
d
L
i
j
(
s
-
s
i
)
j
.
Notice that r
0
(
s
)
=
1
/
P
(
s
). Thus according to (58), we easily obtain the residues of f
0
(64a)
d
0
i
j
=
a
N
×
c
0
i
j
.
If we have already obtained the residues of f
L
-
1, d
(
L
-
1
)
i
j, referring to [17], the residues of (
s
-
s
0
)
f
L
-
1 can be obtained as d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j. Therefore, according to (62b) we have
(64b)
d
L
i
j
=
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
,
aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
(
L
=
1
:
N
)
.Notice that c
i
j
=
d
N
i
j. As c
0
i
j can already be calculated using (59a) and (59b), we can finally obtain the residues of R
(
s
) using formulas (64a) and (64b) successively by increasing the value of L from 1 to N. Equations (64a) and (64b) can be further developed and yield another beautiful formula:
(65)
d
L
i
j
=
∑
m
=
0
L
a
N
-
(
L
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
In (64a), (64b), and (65), if j
+
1
>
m
i or j
+
m
-
k
>
m
i, the residues c
0
i
(
j
+
1
) and c
0
i
(
j
+
m
-
k
) do not exist and we can simply set them as zero. Formula (65) can be proved by induction as below.Proof of Formula (65).
WhenL
=
0, (65) is reduced to (64a). Thus it is true. Suppose (65) is true when L
=
n; namely,
(66)
d
n
i
j
=
∑
m
=
0
n
a
N
-
(
n
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
Then according to (64b), when L
=
n
+
1, we have
(67)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
d
n
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
n
i
j
.
Applying (66) to (67) yields
(68)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
∑
m
=
0
n
a
N
-
(
n
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
1
+
m
-
k
)
+
(
s
i
-
s
0
)
×
∑
m
=
0
n
a
N
-
(
n
-
m
)
×
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
From (68), we have
(69)
d
(
n
+
1
)
i
j
=
a
N
-
(
n
+
1
)
×
c
0
i
j
+
Λ
,
where(70a)
Λ
=
∑
m
=
0
n
a
N
-
(
n
-
m
)
×
{
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
1
+
m
-
k
)
a
a
a
A
a
a
+
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
+
1
c
0
i
(
j
+
m
-
k
)
}
.
Equation (70a) can further lead to
(70b)
Λ
=
∑
m
=
1
n
+
1
a
N
-
(
n
+
1
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.From (69) and (70b), we can finally obtain
(71)
d
(
n
+
1
)
i
j
=
∑
m
=
0
n
+
1
a
N
-
(
n
+
1
-
m
)
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
(
j
+
m
-
k
)
.
Equation (71) is (65) with L replaced by n
+
1; thus, by induction, it follows that (61) is true.SettingL
=
N in (65) yields
(72)
c
i
j
=
∑
m
=
0
N
a
m
∑
k
=
0
m
C
m
k
(
s
i
-
s
0
)
k
c
0
i
[
j
+
m
-
k
]
.
Equation (72) presents a direct relation between the residues of r
0
(
s
) and the residues of R
(
s
). It can be useful when calculating a particular residue of R
(
s
).


## 3.2. Computing Residual Polynomial Coefficients of Functions in Expanded Form

We will propose novel methods to calculate the residual polynomial coefficients, namely,e
h in (52). Most existing methods calculate such coefficients using polynomial division. An algebraic procedure to calculate such coefficients that avoids polynomial division can be found in [17]. Here we propose two novel simpler methods which also avoid polynomial division. We first propose a method through Laurent expansion and then propose a method through derivatives.


### 3.2.1. Computation ofe
h through Laurent Expansion

As shown in (53), R
(
s
) is a sum of r
k
(
s
), if we can obtain the residual polynomial coefficients of r
k
(
s
), e
k
h, e
h can be then obtained using (57). Obviously, r
k
(
s
) will contribute to the e
h of R
(
s
) only when k
≥
K
=
∑
m
j. Thus k is assumed to be no less than K in this part. The problem is then reduced to how to calculate the e
k
h of r
k
(
s
) (k
=
K
:
N). Notice all the r
k
(
s
) are in factorized form; thus, they can be expanded using the method we proposed in Section 2. However, it is not advisable to do so, for this can lead to cumbersome implementation. Actually, simpler methods can be obtained. Referring to Section 2.2.1, we can prove that
(73)
e
k
h
=
r
~
(
k
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
-
K
)
,
where
(74)
r
~
(
t
)
=
W
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
.
Here s
~
j
=
1
/
(
s
j
-
s
c
) and W is a constant. Let k
=
k
+
1 in (73), we have
(75)
e
(
k
+
1
)
h
=
r
~
(
k
+
1
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
+
1
-
K
)
.
From (73) and (75), we can observe that
(76)
e
k
h
=
e
(
k
+
1
)
(
h
+
1
)
,
(
h
=
0
:
k
-
K
)
.
Notice that e
(
k
+
1
)
(
h
+
1
) are the residual polynomial coefficients of r
k
+
1
(
s
). It is clear that if we have already obtained the coefficients of r
k
+
1
(
s
), we can immediately obtain the coefficients of r
k
(
s
) by (76). It hereby follows that if we have already obtained the coefficients of r
N
(
s
), e
N
h, we can successively obtain e
(
N
-
1
)
h
,
e
(
N
-
2
)
h
,
…
,
e
K
h. From (76), the following formula can be obtained:
(77)
e
(
K
+
i
)
h
=
e
N
(
E
-
i
+
h
)
,
(
i
=
0
:
E
,
h
=
0
:
i
)
.
Equation (77) presents the relationship between the e
k
h (k
<
N) and e
N
h, which can be seen more clearly in Table 1. Here we mention that we calculate the e
k
h with k decreasing from N to K, while the method in [17] calculates e
k
h with k increasing from K to N. It will be seen that the method proposed here is simpler and more convenient.Table 1
The relationship between thee
k
h (k
<
N) and e
N
h.
e
k
h
h
0
1
2
⋯
E
-
2
E
-
1
E
k
N
e
N
0
e
N
1
e
N
2
⋯
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
N
-
1
e
N
1
e
N
2
e
N
3
⋯
e
N
(
E
-
1
)
e
N
E
N
-
2
e
N
2
e
N
2
e
N
4
⋯
e
N
E
⋯
⋯
⋯
⋯
⋯
K
+
2
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
K
+
1
e
N
(
E
-
1
)
e
N
E
K
e
N
EAccording to (57) and (77), we have
(78)
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
N
(
E
-
i
)
,
(
h
=
0
:
E
)
.
From the above analysis, the whole problem is reduced to how to obtain e
N
h. As r
N
(
s
) is of factorized form, e
N
h can be calculated using the methods we proposed in Section 2.2.1 through Laurent expansion. Thus, we have(79a)
e
N
h
=
{
1
,
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
N
(
h
+
i
)
η
i
,
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(79b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
,
(
s
c
=
s
0
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
otherwise
)
.In summary, the residual polynomial coefficients can be obtained by first using (79a) and (79b) and then using (78).


### 3.2.2. Computation ofe
h through Derivatives

As we discussed in Section3.2.1, we just have to calculate the residual polynomial coefficients of r
N
(
s
), e
N
h and then we can obtain e
h. As r
N
(
s
) is of factorized form, it is also possible to use the method we proposed in Section 2.2.2 to calculate e
N
h. In Section 2.2.2, three conditions are considered. Here we only discuss the most favorable condition, s
c
=
s
0 (namely, s
c is a zero of r
N
(
s
)) because firstly, this condition is most suitable as it requires the least computation and secondly, s
c
=
s
0
=
0 in most practical cases. Notice that N
>
K. Thus we can use a formula like (44) to calculate e
N
h when s
0 is not equal to any pole of r
N
(
s
). When s
0 is equal to a pole of r
N
(
s
), we just need to make a simple modification. e
N
E
=
1 is obvious. Using the procedures in Section 2.2.2(2), we can finally have
(80)
e
N
h
=
T
N
(
h
)
(
s
)
h
!
|
s
=
s
c
,
(
h
=
1
:
E
-
1
)
,
where
(81)
T
N
(
h
)
(
s
)
=
{
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
s
=
s
k
)
,
(
-
1
)
h
-
1
∑
i
=
1
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
otherwise
)
.


## 3.2.1. Computation ofe
h through Laurent Expansion

As shown in (53), R
(
s
) is a sum of r
k
(
s
), if we can obtain the residual polynomial coefficients of r
k
(
s
), e
k
h, e
h can be then obtained using (57). Obviously, r
k
(
s
) will contribute to the e
h of R
(
s
) only when k
≥
K
=
∑
m
j. Thus k is assumed to be no less than K in this part. The problem is then reduced to how to calculate the e
k
h of r
k
(
s
) (k
=
K
:
N). Notice all the r
k
(
s
) are in factorized form; thus, they can be expanded using the method we proposed in Section 2. However, it is not advisable to do so, for this can lead to cumbersome implementation. Actually, simpler methods can be obtained. Referring to Section 2.2.1, we can prove that
(73)
e
k
h
=
r
~
(
k
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
-
K
)
,
where
(74)
r
~
(
t
)
=
W
∏
j
=
1
M
(
t
-
s
~
j
)
m
j
.
Here s
~
j
=
1
/
(
s
j
-
s
c
) and W is a constant. Let k
=
k
+
1 in (73), we have
(75)
e
(
k
+
1
)
h
=
r
~
(
k
+
1
-
K
-
h
)
(
0
)
h
!
,
(
h
=
0
:
k
+
1
-
K
)
.
From (73) and (75), we can observe that
(76)
e
k
h
=
e
(
k
+
1
)
(
h
+
1
)
,
(
h
=
0
:
k
-
K
)
.
Notice that e
(
k
+
1
)
(
h
+
1
) are the residual polynomial coefficients of r
k
+
1
(
s
). It is clear that if we have already obtained the coefficients of r
k
+
1
(
s
), we can immediately obtain the coefficients of r
k
(
s
) by (76). It hereby follows that if we have already obtained the coefficients of r
N
(
s
), e
N
h, we can successively obtain e
(
N
-
1
)
h
,
e
(
N
-
2
)
h
,
…
,
e
K
h. From (76), the following formula can be obtained:
(77)
e
(
K
+
i
)
h
=
e
N
(
E
-
i
+
h
)
,
(
i
=
0
:
E
,
h
=
0
:
i
)
.
Equation (77) presents the relationship between the e
k
h (k
<
N) and e
N
h, which can be seen more clearly in Table 1. Here we mention that we calculate the e
k
h with k decreasing from N to K, while the method in [17] calculates e
k
h with k increasing from K to N. It will be seen that the method proposed here is simpler and more convenient.Table 1
The relationship between thee
k
h (k
<
N) and e
N
h.
e
k
h
h
0
1
2
⋯
E
-
2
E
-
1
E
k
N
e
N
0
e
N
1
e
N
2
⋯
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
N
-
1
e
N
1
e
N
2
e
N
3
⋯
e
N
(
E
-
1
)
e
N
E
N
-
2
e
N
2
e
N
2
e
N
4
⋯
e
N
E
⋯
⋯
⋯
⋯
⋯
K
+
2
e
N
(
E
-
2
)
e
N
(
E
-
1
)
e
N
E
K
+
1
e
N
(
E
-
1
)
e
N
E
K
e
N
EAccording to (57) and (77), we have
(78)
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
N
(
E
-
i
)
,
(
h
=
0
:
E
)
.
From the above analysis, the whole problem is reduced to how to obtain e
N
h. As r
N
(
s
) is of factorized form, e
N
h can be calculated using the methods we proposed in Section 2.2.1 through Laurent expansion. Thus, we have(79a)
e
N
h
=
{
1
,
(
h
=
E
)
,
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
N
(
h
+
i
)
η
i
,
(
h
=
E
-
1
:
-
1
:
0
)
,
where
(79b)
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
,
(
s
c
=
s
0
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
,
(
otherwise
)
.In summary, the residual polynomial coefficients can be obtained by first using (79a) and (79b) and then using (78).


## 3.2.2. Computation ofe
h through Derivatives

As we discussed in Section3.2.1, we just have to calculate the residual polynomial coefficients of r
N
(
s
), e
N
h and then we can obtain e
h. As r
N
(
s
) is of factorized form, it is also possible to use the method we proposed in Section 2.2.2 to calculate e
N
h. In Section 2.2.2, three conditions are considered. Here we only discuss the most favorable condition, s
c
=
s
0 (namely, s
c is a zero of r
N
(
s
)) because firstly, this condition is most suitable as it requires the least computation and secondly, s
c
=
s
0
=
0 in most practical cases. Notice that N
>
K. Thus we can use a formula like (44) to calculate e
N
h when s
0 is not equal to any pole of r
N
(
s
). When s
0 is equal to a pole of r
N
(
s
), we just need to make a simple modification. e
N
E
=
1 is obvious. Using the procedures in Section 2.2.2(2), we can finally have
(80)
e
N
h
=
T
N
(
h
)
(
s
)
h
!
|
s
=
s
c
,
(
h
=
1
:
E
-
1
)
,
where
(81)
T
N
(
h
)
(
s
)
=
{
(
-
1
)
h
-
1
∑
i
=
1
i
≠
k
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
s
=
s
k
)
,
(
-
1
)
h
-
1
∑
i
=
1
M
∑
j
=
1
m
i
P
h
+
j
-
1
h
c
N
i
j
(
s
-
s
i
)
j
+
h
,
(
otherwise
)
.


## 3.3. Recommendations on the Choices of Procedures

In Sections3.1 and 3.2, several procedures are proposed for the computation of both e
h and c
i
j. As the calculation of e
h can have a relation with that of c
i
j, we recommend two methods comprised of several aforeproposed procedures for improper rational functions. We only suggest two typical methods which we consider to be preferable and suitable as below.Method 1 uses (80) and (81) to calculate e
h. As we can see from (80) and (81), the method in Section 3.2.2 requires that we first obtain the residues of r
N
(
s
). Therefore, we use (60) (rather than (64a) and (64b)) to obtain the residues because the residues of r
N
(
s
) are obtained in the calculation process. Method 2 uses (59a), (59b), (64a), and (64b) to obtain the residues, (79a) and (79b), and (78) to calculate e
h. It is more preferable to use (79a) and (79b) to calculate e
N
h as it is not based on the residues of r
N
(
s
).Method 1. The following procedures comprise Method 1.(a)
Calculate the residues ofr
0
(
s
):
(82)
c
0
i
j
=
{
1
∏
k
=
1
,
k
≠
i
M
(
s
i
-
s
k
)
m
k
,
(
i
=
1
:
M
;
j
=
m
i
)
1
m
i
-
j
∑
k
=
1
m
i
-
j
c
0
i
(
j
+
k
)
λ
k
|
s
=
s
i
,
(
i
=
1
:
M
;
j
=
m
i
-
1
:
-
1
:
1
)
,
whereλ
k
=
∑
L
=
1
,
L
≠
i
M
(
m
L
/
(
s
L
-
s
)
k
).
(b)
Calculate the residues ofr
k
(
s
) using
(83)
c
k
i
j
=
c
(
k
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
c
(
k
-
1
)
i
j
,
aa
(
k
=
1
:
N
;
i
=
1
:
M
;
j
=
1
:
m
i
)
,
and then obtain the residues ofR
(
s
) using c
i
j
=
∑
k
=
0
N
a
k
c
k
i
j.
(c)
Calculate the residues ofr
N
(
s
) using (80) and (81).
(d)
Calculate the residual polynomial coefficients using (78).Method 2. The following procedures comprise Method 2.(a)
Calculate the residues ofr
0
(
s
) the same way as the first step of Method 1.
(b)
Calculate the residues ofR
(
s
). Here notice that c
i
j
=
d
N
i
j. Consider
(84)
d
L
i
j
=
{
a
N
×
c
0
i
j
,
(
L
=
0
)
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
,
(
L
=
1
:
N
)
.
(c)
Calculate the residual polynomial coefficients ofr
N
(
s
) using  (79a) and (79b).
(d)
Calculate the residual polynomial coefficients using (78).


## 4. Examples and Discussions

We provide six examples to illustrate the usage and validity of the proposed methods. We use similar practices as used in [17] to validate the efficiency of the proposed methods for large-scale problems. Matlab codes that perform the proposed methods can be easily obtained. The calculation of the numerical examples is performed by Matlab (2010b) on a PC of 32-bit word length. In those examples, s
c and s
0 are set as zeros, if not mentioned particularly. For simplicity of expression, for large-scale functions of R
(
s
), we denote the coefficients of R
(
s
) as follows: array of poles, S
=
[
s
1
,
s
2
,
…
,
s
M
]; array of the multiplicities of poles, m
=
[
m
1
,
m
2
,
…
,
m
M
]; array of zeros, Z
=
[
z
1
,
z
2
,
…
,
z
M
]; array of the multiplicities of zeros, n
=
[
n
1
,
n
2
,
…
,
n
M
]; polynomial coefficients of the numerator of R
(
s
) in expanded form, A
=
[
a
0
,
a
1
,
…
,
a
N
]. The reader can refer to (1) and (52) to see clearly the meaning of those coefficients.


### 4.1. Usage and Validation of Methods for Factorized Functions

Example 1.
The following simple rational function is provided as an instance to illustrate the usage of the method proposed in Section2. The example is also used by Linnér in [3]. One may refer to [3] to see the differences between the proposed method and Linnér’ method:
(85)
R
(
s
)
=
(
s
+
3
)
(
s
-
1
)
3
(
s
-
2
)
3
(
s
-
3
)
s
(
s
+
1
)
4
(
s
+
2
)
.
Obviously, the desired expansion is given by
(86)
R
(
s
)
=
y
(
s
)
-
T
(
s
)
,
where
(87)
y
(
s
)
=
e
0
+
e
1
(
s
-
s
c
)
+
e
2
(
s
-
s
c
)
2
,
T
(
s
)
=
-
[
c
1
s
+
c
21
(
s
+
1
)
+
c
22
(
s
+
1
)
2
a
a
a
a
a
+
c
23
(
s
+
1
)
3
+
c
24
(
s
+
1
)
4
+
c
3
s
+
2
]
.


#### 4.1.1. Calculation of Residues

The following residues can be calculated in a direct way:(88)
c
1
=
s
R
(
s
)
|
s
=
0
=
-
36
,
c
24
=
(
s
+
1
)
4
R
(
s
)
|
s
=
-
1
=
1728
,
c
3
=
(
s
+
2
)
R
(
s
)
|
s
=
-
2
=
4320
.
Using (11a) and (11b), we can calculate the remaining residues recursively at s
=
-
1. Notice that −1 is the second pole of R
(
s
). First, calculate λ
^
i using (11b)
(89)
λ
^
1
=
(
1
-
2
+
1
+
1
0
+
1
)
-
(
1
-
3
+
1
+
3
1
+
1
+
3
2
+
1
+
1
3
+
1
)
=
-
9
4
,
λ
^
2
=
(
1
(
-
1
)
2
+
1
1
2
)
-
(
1
(
-
2
)
2
+
3
2
2
+
3
3
2
+
1
4
2
)
=
29
48
,
λ
^
3
=
(
1
(
-
1
)
3
+
1
1
3
)
-
(
1
(
-
2
)
3
+
3
2
3
+
3
3
3
+
1
4
3
)
=
-
217
576
.
Then using (11a),
(90)
c
23
=
c
24
λ
^
1
=
-
3888
;
c
22
=
1
2
(
c
23
λ
^
1
+
c
24
λ
^
2
)
=
4896
;
c
21
=
1
3
(
c
22
λ
^
1
+
c
23
λ
^
2
+
c
24
λ
^
3
)
=
-
4672
.
Thus we have
(91)
T
(
s
)
=
-
[
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
a
a
a
a
a
+
1728
(
s
+
1
)
4
+
4320
s
+
2
]
.


#### 4.1.2. Computing Coefficients of the Residual Polynomial

(a) Method Based on Laurent Expansion. Here we use the method proposed in Section 2.2.1. Let s
c be −1, the second pole of R
(
s
). Referring to (30a) and (30b), we have
(92)
η
1
=
(
1
·
(
-
2
+
1
)
+
1
·
(
1
+
0
)
)
-
(
1
·
(
-
3
+
1
)
+
3
·
(
1
+
1
)
aaaaaa
+
3
·
(
2
+
1
)
+
1
·
(
3
+
1
)
)
=
-
17
,
η
2
=
(
1
+
1
)
-
(
(
-
2
)
2
+
3
·
2
2
+
3
·
3
2
+
4
2
)
=
-
57
,
e
2
=
1
,
e
1
=
e
2
η
1
=
-
17
,
e
0
=
1
2
(
e
1
η
1
+
e
2
η
2
)
=
116
.
Thus, we have
(93)
y
(
s
)
=
116
-
17
(
s
+
1
)
+
(
s
+
1
)
2
.(b) Method through Derivatives. Here we use the method proposed in Section 2.2.2 to compute e
h. Obviously, the second zero of R
(
s
), 1, has the largest multiplicity of 3. Thus it is chosen as s
c here. As 2
=
E
<
3
=
n
2, we can use (44) to derive the residual polynomial coefficients
(94)
e
0
=
T
(
1
)
=
86
,
e
1
=
T
(
1
)
(
1
)
=
-
13
.
Thus we have
(95)
y
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
.
The complete expansion is now given by
(96)
R
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
+
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
+
1728
(
s
+
1
)
4
+
4320
s
+
2
.Example 2.
We consider here a large-scale factorized improperR
(
s
) with S
=
[
1
,
2
,
…
,
9
,
10
]; m
=
[
10
,
10
,
…
,
10
,
10
]; Z
=
[
0.5
,
1.5
,
…
,
8.5
,
9.5
]; n
=
[
11
,
11
,
…
,
11
,
11
]. In this case, the degree of the numerator is 100. The degree of the denominator reaches 110. It is quite a large-scale problem. Here, one may see the significance of developing a method which can deal with factorized R
(
s
) directly. If we use a method designed for rational functions in expanded form, the reformulation of the R
(
s
) into expanded form can lead to huge computation costs. It can also introduce considerable additional error. As done in [17], in the large-scale experiments, we also do not display the expansion coefficients because they are too many and too tedious to display and references of the coefficients are very difficult to find. Instead, we plot the figure of R
(
s
) and that of its pfe, pfe
(
s
), together to validate the accuracy of the expansion results. And we also calculate the relative errors of pfe
(
s
) at different values of s as |
(
R
(
s
)
-
pfe
(
s
)
)
/
R
(
s
)
|  (
R
(
s
)
≠
0
). The results are shown in Figures 1 and 2. The figures of the functions are plotted within the interval [
25,100
]. The residues are calculated using the methods proposed in Section 2.1. In Figure 1, the residual polynomial coefficients are calculated by the method we proposed in Section 2.2.1. Notice that, in this figure, the curve of pfe, pfe is in perfect agreement with that of R
(
s
) (the reference solution), demonstrating the high accuracy of the expansion results. And the relative error is nearly negligible. Figure 2 presents the relative errors when residual polynomial coefficients are calculated using the method proposed in Section 2.2.2. Three experiments regarding the values of s
c are considered: (a) s
c
=
8 (pole of R
(
s
), method in Section 2.2.2(3) used); (b) s
c
=
3.5 (zero of R
(
s
), method in Section 2.2.2(2) used); (c) s
c
=
0 (neither pole nor zero of R
(
s
), method in Section 2.2.2(1) used). As we can see, the relative errors are also quite small, demonstrating the good performance of those proposed methods. From Figures 1 and 2, one may notice that the method through Laurent expansion seems to perform better than the method through derivatives to some extent. This is because the former calculates the residues without the usage of the residues.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 2. Residual polynomial coefficients are calculated using the method proposed in Section 2.2.1.
(a)
(b)Relative errors ofpfe with residual polynomial coefficients calculated using the methods in Sections 2.2.2(3)  (a), 2.2.2(2)  (b), and 2.2.2(1)  (c) in Example 2.
(a)
(b)
(c)Example 3.
In this example, we consider a large-scale factorized improper rational function with ill-conditioned poles.S
=
[
1000
,
0
,
0.1
,
8
,
9
]; m
=
[
8
,
7
,
7
,
6
,
6
]; Z
=
[
-
1
,
-
2
,
7
,
11
,
12
,
990
]; n
=
[
8
,
8
,
8
,
8
,
8
,
8
]. The degree of the numerator is 34. The degree of the denominator is 48. R
(
s
) contains three ill-conditioned high-order poles: two high-order poles (0 and 0.1), very close to each other, and an 8th pole (1000) much larger than the other poles. This problem can be rather tricky and the ill-conditioned poles may often lead to intolerable large errors. In this example, the residual polynomial coefficients are calculated using the method proposed in Section 2.2.1. We find the method proposed in Section 2.2.2 is not quite suitable here as some of the residues are extremely large (nearly, 1050), the manipulation of which will inevitably introduce large errors. As can be seen in Figure 3, our method can provide quite satisfactory results. The figure of pfe function is in perfect agreement with the reference solution. And the nearly negligible relative error well demonstrates the method’s good performance in tackling functions with ill-conditioned poles.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 3.
(a)
(b)


### 4.2. Usage and Validation of Methods for Functions in Expanded Form

Example 4.
The following simple rational function is provided to illustrate the usage of the method we proposed for thepfe of functions in expanded form proposed in Section 3:
(97)
R
(
s
)
=
8
+
3
s
+
s
2
+
4
s
3
+
5
s
4
+
s
6
+
2
s
7
+
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
.
It can be estimated that
(98)
R
(
s
)
=
e
0
+
e
1
s
+
e
2
s
2
+
e
2
s
3
+
c
11
s
+
c
21
(
s
+
1
)
+
c
22
(
s
+
1
)
2
+
c
31
(
s
+
2
)
+
c
32
(
s
+
2
)
2
.
Notice the polynomial coefficients of the numerator are A
=
[
8
,
3
,
1
,
4
,
5
,
0
,
1
,
2
,
1
]. We will use the two methods to solve this problem. The procedures of the methods are described in Section 3.3.Method 1. (a) Use (59a) and (59b) to expand r
0
(
s
)
(99)
r
0
(
s
)
=
1
s
(
s
+
1
)
2
(
s
+
2
)
2
=
0.25
s
+
1
(
s
+
1
)
+
-
1
(
s
+
1
)
2
+
-
1.25
(
s
+
2
)
+
-
0.5
(
s
+
2
)
2
.(b) Use (60) to calculate the residues of r
k
(
s
), c
k
i
jbased on the residues of r
0
(
s
). The results are shown in Table 2.Table 2
The residues ofr
k
(
s
), c
k
i
j in Example 4.
k
c
k
11
c
k
21
c
k
22
c
k
31
c
k
32
1
0
-
2
1
2
1
2
0
3
-
1
-
3
-
2
3
0
-
4
1
4
4
4
0
5
-
1
-
4
-
8
5
0
-
6
1
0
16
6
0
7
-
1
16
-
32
7
0
-
8
1
-
64
64
8
0
9
-
1
192
-
128Then use (56) to obtain the residues of R
(
s
). 								(100)
c
11
=
2
,
c
21
=
14
,
c
22
=
-
7
,
c
31
=
69
,
c
32
=
-
59
.(c) From Table2, we can find the residues of r
8
(
s
), c
8
i
j. 								(101)
r
8
(
s
)
=
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
=
e
80
+
e
81
s
-
e
82
s
2
+
s
3
-
T
8
(
s
)
,
where
(102)
T
8
(
s
)
=
-
[
9
(
s
+
1
)
+
-
1
(
s
+
1
)
2
+
192
(
s
+
2
)
+
-
128
(
s
+
1
)
2
]
.
Using (80) and (81), we have
(103)
e
80
=
T
8
(
0
)
=
-
72
,
e
81
=
T
8
(
1
)
(
0
)
=
23
,
e
82
=
T
8
(
2
)
(
0
)
=
-
4
,
e
83
=
1
.
Then referring to (78) or Table 1, we can obtain all the e
k
h as shown in Table 3.Table 3
The residual polynomial coefficients ofr
k
(
s
), e
k
h in Example 4.
e
k
h
h
0
1
2
3
k
8
-
72
23
-
6
1
7
-
23
-
6
1
6
-
6
1
5
1(d) Using (78), we have e
0
=
-
32, e
1
=
12, e
2
=
-
4, e
3
=
1.Method 2. (a) Use (59a) and (59b) to expand r
0
(
s
).(b) Using (64a) and (64b) by increasing the value of L from 0 to N, we can then obtain the residues of R
(
s
) as
(104)
c
11
=
2
,
c
21
=
14
,
c
22
=
-
7
,
c
31
=
69
,
c
32
=
-
59
.(c) Obtain the residual polynomial coefficients ofr
N
(
s
) (here, k
=
N
=
8), e
N
h using (79a) an (79b)
(105)
r
8
(
s
)
=
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
=
s
7
(
s
+
1
)
2
(
s
+
2
)
2
=
-
72
+
23
s
-
6
s
2
+
s
3
+
T
8
(
s
)
,
where T
8
(
s
) is proper residual fraction. Then using (78) or Table 1, we can obtain all the e
k
h as shown in Table 3.(d) Finally, using (78), we have e
0
=
-
32, e
1
=
12, e
2
=
-
4, e
3
=
1.Thus the desired expansion is(106)
R
(
s
)
=
-
32
+
12
s
-
4
s
2
+
s
3
+
2
s
+
14
(
s
+
1
)
+
-
7
(
s
+
1
)
2
+
69
(
s
+
2
)
+
-
59
(
s
+
2
)
2
.Example 5.
In this example, we consider a large-scale rational function in expanded form. The example has also been used in [17]. Thus one can make a comparison between the performance of the proposed methods and the method in [17]. S
=
[
8
,
7
,
6
,
5
,
4
,
3
,
2
,
1
]; m
=
[
1
,
2
,
3
,
4
,
5
,
6
,
7
,
8
]; A
=
[
0.1
,
0.2
,
…
,
5.9
,
6
]. In this case, the degree of the denominator is 36, and the degree of the numerator is 59. The two methods proposed in Section 3 are validated. The results are shown in Figure 4. As shown in the figure, the methods we proposed can provide quite good results for such a large-scale problem. The relative error is desirable for both the methods, demonstrating its good performance. From Figures 4(b) and 4(c), we can see that Method 2 performs better than Method 1.Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example5.
(a)
(b)
(c)Example 6.
In this example, we considera proper rational function in expanded form with ill-conditioned poles. This example is also used in [17]. S
=
[
1000
,
3
,
2
,
1
,
1.1
]; m
=
[
7
,
3
,
3
,
6
,
7
]; A
=
[
0.1
,
0.2
,
…
,
2.2
]. Thus R
(
s
) is a function containing three ill-conditioned high-order poles: two close high-order poles (1 and 1.1) and a 7th pole (1000) much larger than the other poles. Again we validate the two methods proposed in Section 3. As shown in Figure 5, our methods can provide quite satisfactory results and the relative errors are small enough. While the relative error is comparable to the method in [17], we find that the proposed methods are more robust. As mentioned or shown in [8, 9, 17], the input order of poles can have significant influence on the accuracy of the expansion results especially when the rational functions to be expanded contain ill-conditioned poles. However, the accuracy of our method tends to be less affected by the order of the poles. In this Example, we can rearrange the sequence of poles S as [
1
,
1.1
,
2
,
3
,
1000
] or other sequences and achieve comparable accuracy as the sequence used this example, while the method in [17] does not perform quite well under all these conditions.Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example6.
(a)
(b)
(c)


### 4.3. Computational Efficiency

To provide an insight into the calculation efficiency of the foregoing methods, we count the number of long operations (multiplication and division) of the proposed methods. Suppose that the zeroes and poles are known and the functions to be expanded are proper. We only calculate the main part of operations. Operations of constant andO
(
n
) complexity are not included. Let V be the degree of numerator of R
(
s
), K
=
∑
i
=
1
M
m
i the degree of denominator, M as number of poles, H
=
∑
i
=
1
M
m
i
2, and N as the number of zeroes. Table 4 presents counts of long operations above one-degree complexity for a proper function.Table 4
Counts of long operations above one-degree complexity for a proper function.
Counts of long operations above one-degree complexity
Method for factorized functions (Section2.2.1)
1
/
2
(
H
M
-
M
K
)
+
M
N
+
M
2
Method for expanded functions
Method 1:1
/
2
(
H
M
-
M
K
)
+
M
2
+
2
V
K Method 2:    1
/
2
(
H
M
-
M
K
)
+
M
2
+
V
KThe calculation time of the examples can also provide an insight into the efficiency of the methods. The calculation time of the examples is displayed in Table5. Overall, we can easily notice that the computing time for each method is very small, which demonstrates their good computational efficiency. For factorized Examples (Examples 2 and 3), we only provide the results when the residual polynomial coefficients are calculated by the method we proposed in Section 2.2.1 in Table 5. If the residual polynomial coefficients are calculated using the methods in Section 2.2.2, the calculation time is between 0.06 s and 0.07 s for Example 2. In Example 3, methods in Section 2.2.2 perform poorly and are not applicable. Obviously method in Section 2.2.1 is more computational efficient and accurate than the methods in Section 2.2.2. Thus for factorized function, we recommend pfe algorithm (see Algorithm 1) using (11a) and (11b) to calculate residues and (30a) and (30b) to calculate residues to achieve the best performance.Table 5
Calculation time of large-scale examples.
Example2
Example3
Example5
Example6
Calculation time (s)
0.0446
0.0384
0.0752 (Method 1)0.0172 (Method 2)
0.0160 (Method 1)0.0155 (Method 2)Algorithm 1:Algorithm for factorized functions.
Step  1.
calculate
c
k
m
k
directly
Fork
=
1
:
M
c
k
m
k
=
∏
j
=
1
N
(
s
k
-
z
j
)
n
j
∏
j
=
1
j
≠
k
M
(
s
k
-
s
j
)
m
j.
End
Step  2.
calculate
c
k
L
(
L
=
m
k
-
1
:
-
1
:
1
)
using (11a) and (11b)
Fork
=
1
:
M
Fori
=
1
:
m
k
-
1
λ
^
i
=
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
k
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
k
)
i
)
End
ForL
=
m
k
-
1
:
-
1
:
1
c
k
L
=
0
Fori
=
1
:
m
k
-
L
c
k
L
=
c
k
L
+
c
k
(
L
+
i
)
λ
^
i
End
c
k
L
=
c
k
L
(
m
k
-
L
)
End
End
Step  3.
Calculate
e
h
using (30a) and (30b)
Step  3.1.
e
E = 1
Step  3.2. calculate
η
i  using (30b)
Fori
=
1
:
E
-
1
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
z
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
otherwise
)
.
End
Step  3.3. Calculate
e
h  (
h
=
E
-
1
:
-
1
:
1
)  using (30a)
Forh
=
E
-
1
:
-
1
:
1
e
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
η
i
EndFor functions in expanded form, both methods perform well for the proper function (Example6). The calculation time is less than 0.02 s for such a large-scale tricky problem, which proves their good calculation efficiency. In Example 5 which involves an improper function, though both perform well, Method 2 is obviously more efficient as it uses less than 1/4 of the time of Method 1. Meanwhile, from the expansion results, we can also find Method 2 can produce more accurate results than Method 1 for improper functions. Thus it is more preferable in practical use. The algorithm for Method 2 is provided in Algorithm 2.Algorithm 2:Algorithm for functions in expanded form.
Step  1.
calculate residues of
r
0
(
s
)
using (59a) (59b).
Step  1.1.  calculate
c
0
i
m
i  directly
For  i
=
1
:
M
c
0
i
m
i
=
1
∏
j
=
1
j
≠
i
M
(
s
i
-
s
j
)
m
j.
End
Step  1.2. calculate
c
0
i
j  (
j
=
m
i
-
1
:
-
1
:
1
)  using
For i
=
1
:
M
Fork
=
1
:
m
i
-
1
λ
^
i
=
∑
L
=
1
L
≠
i
M
m
L
(
s
L
-
s
i
)
k
End
Forj
=
m
i
-
1
:
-
1
:
1
c
0
i
j
=
0
Fork
=
1
:
m
i
-
j
c
0
i
j
=
c
0
i
j
+
c
0
i
(
j
+
k
)
λ
^
k
End
c
0
i
j
=
c
0
i
j
(
m
i
-
j
)
End
End
Step 2.
calculate the residues of
R
(
s
)
using (64a) and (64b)
d
0
i
j
=
a
N
×
c
0
i
j. (
i
=
1
:
M
,
j
=
1
:
m
i
)
ForL
=
1
:
N
Fori
=
1
:
M
Forj
=
1
:
m
i
d
L
i
j
=
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
End
End
End
Step 3.
calculate the residual polynomial coefficients of
r
N
(
s
)
,
e
N
h
using (79a) (79b).
Step  3.1. e
N
E = 1
Step  3.2. calculate
η
i  using (79b)
Fori
=
1
:
E
-
1
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
(
s
c
=
s
0
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
(
otherwise
)
.
End
Step  3.3. Calculate
e
N
h  (
h
=
E
-
1
:
-
1
:
1
)  using (79a)
Forh
=
E
-
1
:
-
1
:
1
e
N
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
N
(
h
+
i
)
η
i
End
Step 4.
calculate the residual polynomial coefficients of
R
(
s
)
using
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
N
(
E
-
i
)
(
h
=
0
:
E
).


### 4.4. Discussions

We would like to discuss some other issues about the usage and numerical validation of the proposed methods. Firstly, in the experiments, we define relative error as|
(
R
(
s
)
-
pfe
(
s
)
)
/
R
(
s
)
| to estimate quantitatively the accuracy of the expansion results. This criterion is valid because the original function and its pfe are supposed to be theoretically equivalent. There is no doubt that the more coherent their figures are the more accurate the expansion results are. However, incoherence between the figures does not always indicate low-accuracy results. The errors calculated by this definition are not the exact errors of the expansion methods. The actual errors can be much smaller than those shown in the figures, especially when ill-conditioned poles are involved. This is because that in floating point arithmetic, evaluation of partial fraction functions itself (even when the expansion result is 100% correct) can lead to considerably large, even intolerable errors when functions with ill-conditioned high-order poles [20]. Hence, it does not necessarily indicate inaccurate expansion results even when the relative errors are comparatively large at some values of s. On the other hand, small relative errors can guarantee good expansion results. When using the defined relative error to validate a pfe method, we suggest that s should not have the values that make R
(
s
) have an extreme small (or extreme large) absolute value. Otherwise the relative errors will inevitably grow large due to the properties of floating point arithmetic of the computer.Secondly, we provide several methods to calculate the residual polynomial coefficients of improperR
(
s
) in Section 2.2. In terms of accuracy, those methods perform almost equally for small- or median-scale problems and some large-scale problems. But when R
(
s
) contains ill-conditioned poles and when the numerator of R
(
s
) has much larger degree than its denominator, the method proposed in Section 2.2.1 performs better than the methods proposed Section 2.2.2, as it does not use the residues as a basis.Thirdly, though we proposed methods forR
(
s
) in both factorized and expanded form, for a given factorized or polynomial function, it is not necessary to transform the function into another form and then perform pfe. Whether in factorized or expanded form, the expansion results using the proposed methods are comparable. Meanwhile, we mention that the methods for factorized R
(
s
) can be more easily performed by manual calculation.Finally,s
0 and s
c in (2) and (52) are set as zeros in almost all the existing articles. In our paper, they are introduced to generalize a polynomial. They can facilitate implementation and reduce calculation in many practical situations when s
0 and s
c of R
(
s
) are not zeros. Moreover, though a careful choice of s
0 and s
c may improve the accuracy to some extent in certain cases, it will not lead to a significant improvement of accuracy in most cases. Thus, s
0 and s
c can generally be simply set as zeros or other values at your convenience in practice.


## 4.1. Usage and Validation of Methods for Factorized Functions

Example 1.
The following simple rational function is provided as an instance to illustrate the usage of the method proposed in Section2. The example is also used by Linnér in [3]. One may refer to [3] to see the differences between the proposed method and Linnér’ method:
(85)
R
(
s
)
=
(
s
+
3
)
(
s
-
1
)
3
(
s
-
2
)
3
(
s
-
3
)
s
(
s
+
1
)
4
(
s
+
2
)
.
Obviously, the desired expansion is given by
(86)
R
(
s
)
=
y
(
s
)
-
T
(
s
)
,
where
(87)
y
(
s
)
=
e
0
+
e
1
(
s
-
s
c
)
+
e
2
(
s
-
s
c
)
2
,
T
(
s
)
=
-
[
c
1
s
+
c
21
(
s
+
1
)
+
c
22
(
s
+
1
)
2
a
a
a
a
a
+
c
23
(
s
+
1
)
3
+
c
24
(
s
+
1
)
4
+
c
3
s
+
2
]
.


### 4.1.1. Calculation of Residues

The following residues can be calculated in a direct way:(88)
c
1
=
s
R
(
s
)
|
s
=
0
=
-
36
,
c
24
=
(
s
+
1
)
4
R
(
s
)
|
s
=
-
1
=
1728
,
c
3
=
(
s
+
2
)
R
(
s
)
|
s
=
-
2
=
4320
.
Using (11a) and (11b), we can calculate the remaining residues recursively at s
=
-
1. Notice that −1 is the second pole of R
(
s
). First, calculate λ
^
i using (11b)
(89)
λ
^
1
=
(
1
-
2
+
1
+
1
0
+
1
)
-
(
1
-
3
+
1
+
3
1
+
1
+
3
2
+
1
+
1
3
+
1
)
=
-
9
4
,
λ
^
2
=
(
1
(
-
1
)
2
+
1
1
2
)
-
(
1
(
-
2
)
2
+
3
2
2
+
3
3
2
+
1
4
2
)
=
29
48
,
λ
^
3
=
(
1
(
-
1
)
3
+
1
1
3
)
-
(
1
(
-
2
)
3
+
3
2
3
+
3
3
3
+
1
4
3
)
=
-
217
576
.
Then using (11a),
(90)
c
23
=
c
24
λ
^
1
=
-
3888
;
c
22
=
1
2
(
c
23
λ
^
1
+
c
24
λ
^
2
)
=
4896
;
c
21
=
1
3
(
c
22
λ
^
1
+
c
23
λ
^
2
+
c
24
λ
^
3
)
=
-
4672
.
Thus we have
(91)
T
(
s
)
=
-
[
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
a
a
a
a
a
+
1728
(
s
+
1
)
4
+
4320
s
+
2
]
.


### 4.1.2. Computing Coefficients of the Residual Polynomial

(a) Method Based on Laurent Expansion. Here we use the method proposed in Section 2.2.1. Let s
c be −1, the second pole of R
(
s
). Referring to (30a) and (30b), we have
(92)
η
1
=
(
1
·
(
-
2
+
1
)
+
1
·
(
1
+
0
)
)
-
(
1
·
(
-
3
+
1
)
+
3
·
(
1
+
1
)
aaaaaa
+
3
·
(
2
+
1
)
+
1
·
(
3
+
1
)
)
=
-
17
,
η
2
=
(
1
+
1
)
-
(
(
-
2
)
2
+
3
·
2
2
+
3
·
3
2
+
4
2
)
=
-
57
,
e
2
=
1
,
e
1
=
e
2
η
1
=
-
17
,
e
0
=
1
2
(
e
1
η
1
+
e
2
η
2
)
=
116
.
Thus, we have
(93)
y
(
s
)
=
116
-
17
(
s
+
1
)
+
(
s
+
1
)
2
.(b) Method through Derivatives. Here we use the method proposed in Section 2.2.2 to compute e
h. Obviously, the second zero of R
(
s
), 1, has the largest multiplicity of 3. Thus it is chosen as s
c here. As 2
=
E
<
3
=
n
2, we can use (44) to derive the residual polynomial coefficients
(94)
e
0
=
T
(
1
)
=
86
,
e
1
=
T
(
1
)
(
1
)
=
-
13
.
Thus we have
(95)
y
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
.
The complete expansion is now given by
(96)
R
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
+
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
+
1728
(
s
+
1
)
4
+
4320
s
+
2
.Example 2.
We consider here a large-scale factorized improperR
(
s
) with S
=
[
1
,
2
,
…
,
9
,
10
]; m
=
[
10
,
10
,
…
,
10
,
10
]; Z
=
[
0.5
,
1.5
,
…
,
8.5
,
9.5
]; n
=
[
11
,
11
,
…
,
11
,
11
]. In this case, the degree of the numerator is 100. The degree of the denominator reaches 110. It is quite a large-scale problem. Here, one may see the significance of developing a method which can deal with factorized R
(
s
) directly. If we use a method designed for rational functions in expanded form, the reformulation of the R
(
s
) into expanded form can lead to huge computation costs. It can also introduce considerable additional error. As done in [17], in the large-scale experiments, we also do not display the expansion coefficients because they are too many and too tedious to display and references of the coefficients are very difficult to find. Instead, we plot the figure of R
(
s
) and that of its pfe, pfe
(
s
), together to validate the accuracy of the expansion results. And we also calculate the relative errors of pfe
(
s
) at different values of s as |
(
R
(
s
)
-
pfe
(
s
)
)
/
R
(
s
)
|  (
R
(
s
)
≠
0
). The results are shown in Figures 1 and 2. The figures of the functions are plotted within the interval [
25,100
]. The residues are calculated using the methods proposed in Section 2.1. In Figure 1, the residual polynomial coefficients are calculated by the method we proposed in Section 2.2.1. Notice that, in this figure, the curve of pfe, pfe is in perfect agreement with that of R
(
s
) (the reference solution), demonstrating the high accuracy of the expansion results. And the relative error is nearly negligible. Figure 2 presents the relative errors when residual polynomial coefficients are calculated using the method proposed in Section 2.2.2. Three experiments regarding the values of s
c are considered: (a) s
c
=
8 (pole of R
(
s
), method in Section 2.2.2(3) used); (b) s
c
=
3.5 (zero of R
(
s
), method in Section 2.2.2(2) used); (c) s
c
=
0 (neither pole nor zero of R
(
s
), method in Section 2.2.2(1) used). As we can see, the relative errors are also quite small, demonstrating the good performance of those proposed methods. From Figures 1 and 2, one may notice that the method through Laurent expansion seems to perform better than the method through derivatives to some extent. This is because the former calculates the residues without the usage of the residues.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 2. Residual polynomial coefficients are calculated using the method proposed in Section 2.2.1.
(a)
(b)Relative errors ofpfe with residual polynomial coefficients calculated using the methods in Sections 2.2.2(3)  (a), 2.2.2(2)  (b), and 2.2.2(1)  (c) in Example 2.
(a)
(b)
(c)Example 3.
In this example, we consider a large-scale factorized improper rational function with ill-conditioned poles.S
=
[
1000
,
0
,
0.1
,
8
,
9
]; m
=
[
8
,
7
,
7
,
6
,
6
]; Z
=
[
-
1
,
-
2
,
7
,
11
,
12
,
990
]; n
=
[
8
,
8
,
8
,
8
,
8
,
8
]. The degree of the numerator is 34. The degree of the denominator is 48. R
(
s
) contains three ill-conditioned high-order poles: two high-order poles (0 and 0.1), very close to each other, and an 8th pole (1000) much larger than the other poles. This problem can be rather tricky and the ill-conditioned poles may often lead to intolerable large errors. In this example, the residual polynomial coefficients are calculated using the method proposed in Section 2.2.1. We find the method proposed in Section 2.2.2 is not quite suitable here as some of the residues are extremely large (nearly, 1050), the manipulation of which will inevitably introduce large errors. As can be seen in Figure 3, our method can provide quite satisfactory results. The figure of pfe function is in perfect agreement with the reference solution. And the nearly negligible relative error well demonstrates the method’s good performance in tackling functions with ill-conditioned poles.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 3.
(a)
(b)


## 4.1.1. Calculation of Residues

The following residues can be calculated in a direct way:(88)
c
1
=
s
R
(
s
)
|
s
=
0
=
-
36
,
c
24
=
(
s
+
1
)
4
R
(
s
)
|
s
=
-
1
=
1728
,
c
3
=
(
s
+
2
)
R
(
s
)
|
s
=
-
2
=
4320
.
Using (11a) and (11b), we can calculate the remaining residues recursively at s
=
-
1. Notice that −1 is the second pole of R
(
s
). First, calculate λ
^
i using (11b)
(89)
λ
^
1
=
(
1
-
2
+
1
+
1
0
+
1
)
-
(
1
-
3
+
1
+
3
1
+
1
+
3
2
+
1
+
1
3
+
1
)
=
-
9
4
,
λ
^
2
=
(
1
(
-
1
)
2
+
1
1
2
)
-
(
1
(
-
2
)
2
+
3
2
2
+
3
3
2
+
1
4
2
)
=
29
48
,
λ
^
3
=
(
1
(
-
1
)
3
+
1
1
3
)
-
(
1
(
-
2
)
3
+
3
2
3
+
3
3
3
+
1
4
3
)
=
-
217
576
.
Then using (11a),
(90)
c
23
=
c
24
λ
^
1
=
-
3888
;
c
22
=
1
2
(
c
23
λ
^
1
+
c
24
λ
^
2
)
=
4896
;
c
21
=
1
3
(
c
22
λ
^
1
+
c
23
λ
^
2
+
c
24
λ
^
3
)
=
-
4672
.
Thus we have
(91)
T
(
s
)
=
-
[
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
a
a
a
a
a
+
1728
(
s
+
1
)
4
+
4320
s
+
2
]
.


## 4.1.2. Computing Coefficients of the Residual Polynomial

(a) Method Based on Laurent Expansion. Here we use the method proposed in Section 2.2.1. Let s
c be −1, the second pole of R
(
s
). Referring to (30a) and (30b), we have
(92)
η
1
=
(
1
·
(
-
2
+
1
)
+
1
·
(
1
+
0
)
)
-
(
1
·
(
-
3
+
1
)
+
3
·
(
1
+
1
)
aaaaaa
+
3
·
(
2
+
1
)
+
1
·
(
3
+
1
)
)
=
-
17
,
η
2
=
(
1
+
1
)
-
(
(
-
2
)
2
+
3
·
2
2
+
3
·
3
2
+
4
2
)
=
-
57
,
e
2
=
1
,
e
1
=
e
2
η
1
=
-
17
,
e
0
=
1
2
(
e
1
η
1
+
e
2
η
2
)
=
116
.
Thus, we have
(93)
y
(
s
)
=
116
-
17
(
s
+
1
)
+
(
s
+
1
)
2
.(b) Method through Derivatives. Here we use the method proposed in Section 2.2.2 to compute e
h. Obviously, the second zero of R
(
s
), 1, has the largest multiplicity of 3. Thus it is chosen as s
c here. As 2
=
E
<
3
=
n
2, we can use (44) to derive the residual polynomial coefficients
(94)
e
0
=
T
(
1
)
=
86
,
e
1
=
T
(
1
)
(
1
)
=
-
13
.
Thus we have
(95)
y
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
.
The complete expansion is now given by
(96)
R
(
s
)
=
86
-
13
(
s
-
1
)
+
(
s
-
1
)
2
+
-
36
s
+
-
4672
(
s
+
1
)
+
4896
(
s
+
1
)
2
+
-
3888
(
s
+
1
)
3
+
1728
(
s
+
1
)
4
+
4320
s
+
2
.Example 2.
We consider here a large-scale factorized improperR
(
s
) with S
=
[
1
,
2
,
…
,
9
,
10
]; m
=
[
10
,
10
,
…
,
10
,
10
]; Z
=
[
0.5
,
1.5
,
…
,
8.5
,
9.5
]; n
=
[
11
,
11
,
…
,
11
,
11
]. In this case, the degree of the numerator is 100. The degree of the denominator reaches 110. It is quite a large-scale problem. Here, one may see the significance of developing a method which can deal with factorized R
(
s
) directly. If we use a method designed for rational functions in expanded form, the reformulation of the R
(
s
) into expanded form can lead to huge computation costs. It can also introduce considerable additional error. As done in [17], in the large-scale experiments, we also do not display the expansion coefficients because they are too many and too tedious to display and references of the coefficients are very difficult to find. Instead, we plot the figure of R
(
s
) and that of its pfe, pfe
(
s
), together to validate the accuracy of the expansion results. And we also calculate the relative errors of pfe
(
s
) at different values of s as |
(
R
(
s
)
-
pfe
(
s
)
)
/
R
(
s
)
|  (
R
(
s
)
≠
0
). The results are shown in Figures 1 and 2. The figures of the functions are plotted within the interval [
25,100
]. The residues are calculated using the methods proposed in Section 2.1. In Figure 1, the residual polynomial coefficients are calculated by the method we proposed in Section 2.2.1. Notice that, in this figure, the curve of pfe, pfe is in perfect agreement with that of R
(
s
) (the reference solution), demonstrating the high accuracy of the expansion results. And the relative error is nearly negligible. Figure 2 presents the relative errors when residual polynomial coefficients are calculated using the method proposed in Section 2.2.2. Three experiments regarding the values of s
c are considered: (a) s
c
=
8 (pole of R
(
s
), method in Section 2.2.2(3) used); (b) s
c
=
3.5 (zero of R
(
s
), method in Section 2.2.2(2) used); (c) s
c
=
0 (neither pole nor zero of R
(
s
), method in Section 2.2.2(1) used). As we can see, the relative errors are also quite small, demonstrating the good performance of those proposed methods. From Figures 1 and 2, one may notice that the method through Laurent expansion seems to perform better than the method through derivatives to some extent. This is because the former calculates the residues without the usage of the residues.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 2. Residual polynomial coefficients are calculated using the method proposed in Section 2.2.1.
(a)
(b)Relative errors ofpfe with residual polynomial coefficients calculated using the methods in Sections 2.2.2(3)  (a), 2.2.2(2)  (b), and 2.2.2(1)  (c) in Example 2.
(a)
(b)
(c)Example 3.
In this example, we consider a large-scale factorized improper rational function with ill-conditioned poles.S
=
[
1000
,
0
,
0.1
,
8
,
9
]; m
=
[
8
,
7
,
7
,
6
,
6
]; Z
=
[
-
1
,
-
2
,
7
,
11
,
12
,
990
]; n
=
[
8
,
8
,
8
,
8
,
8
,
8
]. The degree of the numerator is 34. The degree of the denominator is 48. R
(
s
) contains three ill-conditioned high-order poles: two high-order poles (0 and 0.1), very close to each other, and an 8th pole (1000) much larger than the other poles. This problem can be rather tricky and the ill-conditioned poles may often lead to intolerable large errors. In this example, the residual polynomial coefficients are calculated using the method proposed in Section 2.2.1. We find the method proposed in Section 2.2.2 is not quite suitable here as some of the residues are extremely large (nearly, 1050), the manipulation of which will inevitably introduce large errors. As can be seen in Figure 3, our method can provide quite satisfactory results. The figure of pfe function is in perfect agreement with the reference solution. And the nearly negligible relative error well demonstrates the method’s good performance in tackling functions with ill-conditioned poles.Comparison ofpfe with reference solution (a) and relative error of pfe (b) in Example 3.
(a)
(b)


## 4.2. Usage and Validation of Methods for Functions in Expanded Form

Example 4.
The following simple rational function is provided to illustrate the usage of the method we proposed for thepfe of functions in expanded form proposed in Section 3:
(97)
R
(
s
)
=
8
+
3
s
+
s
2
+
4
s
3
+
5
s
4
+
s
6
+
2
s
7
+
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
.
It can be estimated that
(98)
R
(
s
)
=
e
0
+
e
1
s
+
e
2
s
2
+
e
2
s
3
+
c
11
s
+
c
21
(
s
+
1
)
+
c
22
(
s
+
1
)
2
+
c
31
(
s
+
2
)
+
c
32
(
s
+
2
)
2
.
Notice the polynomial coefficients of the numerator are A
=
[
8
,
3
,
1
,
4
,
5
,
0
,
1
,
2
,
1
]. We will use the two methods to solve this problem. The procedures of the methods are described in Section 3.3.Method 1. (a) Use (59a) and (59b) to expand r
0
(
s
)
(99)
r
0
(
s
)
=
1
s
(
s
+
1
)
2
(
s
+
2
)
2
=
0.25
s
+
1
(
s
+
1
)
+
-
1
(
s
+
1
)
2
+
-
1.25
(
s
+
2
)
+
-
0.5
(
s
+
2
)
2
.(b) Use (60) to calculate the residues of r
k
(
s
), c
k
i
jbased on the residues of r
0
(
s
). The results are shown in Table 2.Table 2
The residues ofr
k
(
s
), c
k
i
j in Example 4.
k
c
k
11
c
k
21
c
k
22
c
k
31
c
k
32
1
0
-
2
1
2
1
2
0
3
-
1
-
3
-
2
3
0
-
4
1
4
4
4
0
5
-
1
-
4
-
8
5
0
-
6
1
0
16
6
0
7
-
1
16
-
32
7
0
-
8
1
-
64
64
8
0
9
-
1
192
-
128Then use (56) to obtain the residues of R
(
s
). 								(100)
c
11
=
2
,
c
21
=
14
,
c
22
=
-
7
,
c
31
=
69
,
c
32
=
-
59
.(c) From Table2, we can find the residues of r
8
(
s
), c
8
i
j. 								(101)
r
8
(
s
)
=
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
=
e
80
+
e
81
s
-
e
82
s
2
+
s
3
-
T
8
(
s
)
,
where
(102)
T
8
(
s
)
=
-
[
9
(
s
+
1
)
+
-
1
(
s
+
1
)
2
+
192
(
s
+
2
)
+
-
128
(
s
+
1
)
2
]
.
Using (80) and (81), we have
(103)
e
80
=
T
8
(
0
)
=
-
72
,
e
81
=
T
8
(
1
)
(
0
)
=
23
,
e
82
=
T
8
(
2
)
(
0
)
=
-
4
,
e
83
=
1
.
Then referring to (78) or Table 1, we can obtain all the e
k
h as shown in Table 3.Table 3
The residual polynomial coefficients ofr
k
(
s
), e
k
h in Example 4.
e
k
h
h
0
1
2
3
k
8
-
72
23
-
6
1
7
-
23
-
6
1
6
-
6
1
5
1(d) Using (78), we have e
0
=
-
32, e
1
=
12, e
2
=
-
4, e
3
=
1.Method 2. (a) Use (59a) and (59b) to expand r
0
(
s
).(b) Using (64a) and (64b) by increasing the value of L from 0 to N, we can then obtain the residues of R
(
s
) as
(104)
c
11
=
2
,
c
21
=
14
,
c
22
=
-
7
,
c
31
=
69
,
c
32
=
-
59
.(c) Obtain the residual polynomial coefficients ofr
N
(
s
) (here, k
=
N
=
8), e
N
h using (79a) an (79b)
(105)
r
8
(
s
)
=
s
8
s
(
s
+
1
)
2
(
s
+
2
)
2
=
s
7
(
s
+
1
)
2
(
s
+
2
)
2
=
-
72
+
23
s
-
6
s
2
+
s
3
+
T
8
(
s
)
,
where T
8
(
s
) is proper residual fraction. Then using (78) or Table 1, we can obtain all the e
k
h as shown in Table 3.(d) Finally, using (78), we have e
0
=
-
32, e
1
=
12, e
2
=
-
4, e
3
=
1.Thus the desired expansion is(106)
R
(
s
)
=
-
32
+
12
s
-
4
s
2
+
s
3
+
2
s
+
14
(
s
+
1
)
+
-
7
(
s
+
1
)
2
+
69
(
s
+
2
)
+
-
59
(
s
+
2
)
2
.Example 5.
In this example, we consider a large-scale rational function in expanded form. The example has also been used in [17]. Thus one can make a comparison between the performance of the proposed methods and the method in [17]. S
=
[
8
,
7
,
6
,
5
,
4
,
3
,
2
,
1
]; m
=
[
1
,
2
,
3
,
4
,
5
,
6
,
7
,
8
]; A
=
[
0.1
,
0.2
,
…
,
5.9
,
6
]. In this case, the degree of the denominator is 36, and the degree of the numerator is 59. The two methods proposed in Section 3 are validated. The results are shown in Figure 4. As shown in the figure, the methods we proposed can provide quite good results for such a large-scale problem. The relative error is desirable for both the methods, demonstrating its good performance. From Figures 4(b) and 4(c), we can see that Method 2 performs better than Method 1.Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example5.
(a)
(b)
(c)Example 6.
In this example, we considera proper rational function in expanded form with ill-conditioned poles. This example is also used in [17]. S
=
[
1000
,
3
,
2
,
1
,
1.1
]; m
=
[
7
,
3
,
3
,
6
,
7
]; A
=
[
0.1
,
0.2
,
…
,
2.2
]. Thus R
(
s
) is a function containing three ill-conditioned high-order poles: two close high-order poles (1 and 1.1) and a 7th pole (1000) much larger than the other poles. Again we validate the two methods proposed in Section 3. As shown in Figure 5, our methods can provide quite satisfactory results and the relative errors are small enough. While the relative error is comparable to the method in [17], we find that the proposed methods are more robust. As mentioned or shown in [8, 9, 17], the input order of poles can have significant influence on the accuracy of the expansion results especially when the rational functions to be expanded contain ill-conditioned poles. However, the accuracy of our method tends to be less affected by the order of the poles. In this Example, we can rearrange the sequence of poles S as [
1
,
1.1
,
2
,
3
,
1000
] or other sequences and achieve comparable accuracy as the sequence used this example, while the method in [17] does not perform quite well under all these conditions.Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example6.
(a)
(b)
(c)


## 4.3. Computational Efficiency

To provide an insight into the calculation efficiency of the foregoing methods, we count the number of long operations (multiplication and division) of the proposed methods. Suppose that the zeroes and poles are known and the functions to be expanded are proper. We only calculate the main part of operations. Operations of constant andO
(
n
) complexity are not included. Let V be the degree of numerator of R
(
s
), K
=
∑
i
=
1
M
m
i the degree of denominator, M as number of poles, H
=
∑
i
=
1
M
m
i
2, and N as the number of zeroes. Table 4 presents counts of long operations above one-degree complexity for a proper function.Table 4
Counts of long operations above one-degree complexity for a proper function.
Counts of long operations above one-degree complexity
Method for factorized functions (Section2.2.1)
1
/
2
(
H
M
-
M
K
)
+
M
N
+
M
2
Method for expanded functions
Method 1:1
/
2
(
H
M
-
M
K
)
+
M
2
+
2
V
K Method 2:    1
/
2
(
H
M
-
M
K
)
+
M
2
+
V
KThe calculation time of the examples can also provide an insight into the efficiency of the methods. The calculation time of the examples is displayed in Table5. Overall, we can easily notice that the computing time for each method is very small, which demonstrates their good computational efficiency. For factorized Examples (Examples 2 and 3), we only provide the results when the residual polynomial coefficients are calculated by the method we proposed in Section 2.2.1 in Table 5. If the residual polynomial coefficients are calculated using the methods in Section 2.2.2, the calculation time is between 0.06 s and 0.07 s for Example 2. In Example 3, methods in Section 2.2.2 perform poorly and are not applicable. Obviously method in Section 2.2.1 is more computational efficient and accurate than the methods in Section 2.2.2. Thus for factorized function, we recommend pfe algorithm (see Algorithm 1) using (11a) and (11b) to calculate residues and (30a) and (30b) to calculate residues to achieve the best performance.Table 5
Calculation time of large-scale examples.
Example2
Example3
Example5
Example6
Calculation time (s)
0.0446
0.0384
0.0752 (Method 1)0.0172 (Method 2)
0.0160 (Method 1)0.0155 (Method 2)Algorithm 1:Algorithm for factorized functions.
Step  1.
calculate
c
k
m
k
directly
Fork
=
1
:
M
c
k
m
k
=
∏
j
=
1
N
(
s
k
-
z
j
)
n
j
∏
j
=
1
j
≠
k
M
(
s
k
-
s
j
)
m
j.
End
Step  2.
calculate
c
k
L
(
L
=
m
k
-
1
:
-
1
:
1
)
using (11a) and (11b)
Fork
=
1
:
M
Fori
=
1
:
m
k
-
1
λ
^
i
=
(
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
k
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
k
)
i
)
End
ForL
=
m
k
-
1
:
-
1
:
1
c
k
L
=
0
Fori
=
1
:
m
k
-
L
c
k
L
=
c
k
L
+
c
k
(
L
+
i
)
λ
^
i
End
c
k
L
=
c
k
L
(
m
k
-
L
)
End
End
Step  3.
Calculate
e
h
using (30a) and (30b)
Step  3.1.
e
E = 1
Step  3.2. calculate
η
i  using (30b)
Fori
=
1
:
E
-
1
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
j
≠
k
N
n
j
(
z
j
-
s
c
)
i
(
s
c
=
z
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
∑
j
=
1
N
n
j
(
z
j
-
s
c
)
i
(
otherwise
)
.
End
Step  3.3. Calculate
e
h  (
h
=
E
-
1
:
-
1
:
1
)  using (30a)
Forh
=
E
-
1
:
-
1
:
1
e
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
(
h
+
i
)
η
i
EndFor functions in expanded form, both methods perform well for the proper function (Example6). The calculation time is less than 0.02 s for such a large-scale tricky problem, which proves their good calculation efficiency. In Example 5 which involves an improper function, though both perform well, Method 2 is obviously more efficient as it uses less than 1/4 of the time of Method 1. Meanwhile, from the expansion results, we can also find Method 2 can produce more accurate results than Method 1 for improper functions. Thus it is more preferable in practical use. The algorithm for Method 2 is provided in Algorithm 2.Algorithm 2:Algorithm for functions in expanded form.
Step  1.
calculate residues of
r
0
(
s
)
using (59a) (59b).
Step  1.1.  calculate
c
0
i
m
i  directly
For  i
=
1
:
M
c
0
i
m
i
=
1
∏
j
=
1
j
≠
i
M
(
s
i
-
s
j
)
m
j.
End
Step  1.2. calculate
c
0
i
j  (
j
=
m
i
-
1
:
-
1
:
1
)  using
For i
=
1
:
M
Fork
=
1
:
m
i
-
1
λ
^
i
=
∑
L
=
1
L
≠
i
M
m
L
(
s
L
-
s
i
)
k
End
Forj
=
m
i
-
1
:
-
1
:
1
c
0
i
j
=
0
Fork
=
1
:
m
i
-
j
c
0
i
j
=
c
0
i
j
+
c
0
i
(
j
+
k
)
λ
^
k
End
c
0
i
j
=
c
0
i
j
(
m
i
-
j
)
End
End
Step 2.
calculate the residues of
R
(
s
)
using (64a) and (64b)
d
0
i
j
=
a
N
×
c
0
i
j. (
i
=
1
:
M
,
j
=
1
:
m
i
)
ForL
=
1
:
N
Fori
=
1
:
M
Forj
=
1
:
m
i
d
L
i
j
=
a
N
-
L
×
c
0
i
j
+
d
(
L
-
1
)
i
(
j
+
1
)
+
(
s
i
-
s
0
)
×
d
(
L
-
1
)
i
j
End
End
End
Step 3.
calculate the residual polynomial coefficients of
r
N
(
s
)
,
e
N
h
using (79a) (79b).
Step  3.1. e
N
E = 1
Step  3.2. calculate
η
i  using (79b)
Fori
=
1
:
E
-
1
η
i
=
{
∑
j
=
1
j
≠
k
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
(
s
c
=
s
k
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
(
s
c
=
s
0
)
,
∑
j
=
1
M
m
j
(
s
j
-
s
c
)
i
-
N
(
z
j
-
s
c
)
N
(
otherwise
)
.
End
Step  3.3. Calculate
e
N
h  (
h
=
E
-
1
:
-
1
:
1
)  using (79a)
Forh
=
E
-
1
:
-
1
:
1
e
N
h
=
1
(
E
-
h
)
∑
i
=
1
E
-
h
e
N
(
h
+
i
)
η
i
End
Step 4.
calculate the residual polynomial coefficients of
R
(
s
)
using
e
h
=
∑
i
=
0
E
-
h
a
K
+
h
+
i
e
N
(
E
-
i
)
(
h
=
0
:
E
).


## 4.4. Discussions

We would like to discuss some other issues about the usage and numerical validation of the proposed methods. Firstly, in the experiments, we define relative error as|
(
R
(
s
)
-
pfe
(
s
)
)
/
R
(
s
)
| to estimate quantitatively the accuracy of the expansion results. This criterion is valid because the original function and its pfe are supposed to be theoretically equivalent. There is no doubt that the more coherent their figures are the more accurate the expansion results are. However, incoherence between the figures does not always indicate low-accuracy results. The errors calculated by this definition are not the exact errors of the expansion methods. The actual errors can be much smaller than those shown in the figures, especially when ill-conditioned poles are involved. This is because that in floating point arithmetic, evaluation of partial fraction functions itself (even when the expansion result is 100% correct) can lead to considerably large, even intolerable errors when functions with ill-conditioned high-order poles [20]. Hence, it does not necessarily indicate inaccurate expansion results even when the relative errors are comparatively large at some values of s. On the other hand, small relative errors can guarantee good expansion results. When using the defined relative error to validate a pfe method, we suggest that s should not have the values that make R
(
s
) have an extreme small (or extreme large) absolute value. Otherwise the relative errors will inevitably grow large due to the properties of floating point arithmetic of the computer.Secondly, we provide several methods to calculate the residual polynomial coefficients of improperR
(
s
) in Section 2.2. In terms of accuracy, those methods perform almost equally for small- or median-scale problems and some large-scale problems. But when R
(
s
) contains ill-conditioned poles and when the numerator of R
(
s
) has much larger degree than its denominator, the method proposed in Section 2.2.1 performs better than the methods proposed Section 2.2.2, as it does not use the residues as a basis.Thirdly, though we proposed methods forR
(
s
) in both factorized and expanded form, for a given factorized or polynomial function, it is not necessary to transform the function into another form and then perform pfe. Whether in factorized or expanded form, the expansion results using the proposed methods are comparable. Meanwhile, we mention that the methods for factorized R
(
s
) can be more easily performed by manual calculation.Finally,s
0 and s
c in (2) and (52) are set as zeros in almost all the existing articles. In our paper, they are introduced to generalize a polynomial. They can facilitate implementation and reduce calculation in many practical situations when s
0 and s
c of R
(
s
) are not zeros. Moreover, though a careful choice of s
0 and s
c may improve the accuracy to some extent in certain cases, it will not lead to a significant improvement of accuracy in most cases. Thus, s
0 and s
c can generally be simply set as zeros or other values at your convenience in practice.


## 5. Conclusions

In this paper, we developed efficient recursive methods for thepfe of general rational functions in both factorized and expanded form. Simple, elegant, recursive formulas that describe the relation of the residues and the coefficients of the residual polynomial are obtained. These methods tend to be simpler and more applicable than many existing methods for pfe of rational functions with multiple high-order poles. They can be easily programmed for computer use with desirable efficiency and accuracy. They are also very stable whose accuracy is less affected by the input-order of poles. The methods are also very suitable for manual calculation.



---

*Source: 895036-2014-10-12.xml*