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.

1. Introduction

Partial fraction expansion () 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 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 methods for general rational functions in both factorized form and expanded form.

Suppose is the rational function to be expanded in . The polynomials, and , are denominator and numerator of , respectively. Generally, a polynomial (let us say ), can be written in either expanded form or factorized form 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 problems, as the poles of have to be known, is mostly written in factorized form. Though can be written in either factorized or expanded form, the latter seems to be more in favor in practice. Most current available methods are designed for with written in expanded form. However, one may still be motivated to design methods for factorized . For one thing, many functions are originally written in factorized form. A method for factorized can deal with such functions directly and more efficiently, therefore, it is more preferable. For another, under many practical circumstances, even is originally written in expanded form, the zeros and poles of have to be obtained. In other words, needs to be rewritten in factorized form. This can be often seen in the field of system analysis and design, where 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 methods for factorized are also mentioned in [2]. In a word, it is of significance to design methods for the of in both factorized and expanded form.

The most well-known method is Heaviside’ cover-up formula. It is often introduced as a textbook method. It provides an elegant compact solution to problems and often serves as a basis for other 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 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 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 of functions in factorized form [3, 4]. In [4], Brugia developed an efficient noniterative method for 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 . 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 [79]. 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 [1217]. 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 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 methods is discussed in [18, 19].

In this paper, we develop efficient recursive algebraic methods for rational functions in both factorized and expanded form. Simple, elegant, and pure-algebraic formulas are proposed to compute the coefficients. The proposed methods can be used for 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 of a general rational function. The remainder of the paper is organized as follows. In Section 2, we focus on the 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 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 which is expanded as Here are the residues; and are the zeroes and poles of , respectively. The poles and zeroes are assumed distinct in the whole paper; and are the multiplicities of and , respectively; is the number of zeros and is the number of poles; are the coefficients of the residual polynomial which will be zeros for a proper ; is a constant introduced to make the residual polynomial more general. is often set as zero in other articles. From (1) and (2), we can observe that the degree of the numerator of is ; the degree of the denominator is ; . In a problem, , 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 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 where and As shown in [4], the th derivatives of are where and Here , , and . In the above equations, we denote by 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 . It represents a very efficient 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 Making some mathematical transformations on (6) yields 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 Using (7) and (8) to substitute the derivatives of and in the summation of (9) yields Notice . Equation (10) further yields where Here we mention that in (11a) means . This is a Matlab denotation which will also be used in the following part for the simplicity of expression. Obviously . By using (11a) and (11b) and decreasing progressively the value of from 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, , 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 that avoid long division and the reformation of . 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 Any

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 in (2) not be equal to zeros or poles of , which can cause inconvenience in practice. We will extend to any value. Linnér’s method is briefly introduced to allow for context, completeness, and illumination of differences. With the substitution , (1) transforms into where , , and is a constant. Now write where does not include the factor . In (13), applying a Taylor expansion to at , we can obtain where is an integer. Now expand as where is the proper fraction, which can be expanded as the second part on the right-side of (2). Comparing (14) and (15), we have This is Linnér’s main result. In the method, the derivatives of are calculated using equations like (5) and (6) and it also requires and , which can cause inconvenience in practice. In the following part, we first show how to extend to zeros and poles of . We then obtain much simpler recursive formula for the calculation of .

Suppose is the th root of , . We can rewrite where is as defined in (4). With the substitution in , we have where Here , , and By Taylor expanding at in (18), we can obtain From (17) and (21), we have Comparing (15) and (22), we can observe that Thus is extended to the th root of , . The same procedure can be used to extend to the th zero of , .

We can easily obtain by setting in (23) or by inspection: 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 , and are both of factorized form and that (3) and (23), which calculate and , respectively, are of the similar form. Therefore, the aforeproposed procedure in Section 2.1 for the calculation of is also suitable for the calculation of here. Take , for an instance. It can be proven that where From (23), we have According to (23) and (25), we have Using (26) and (27) to substitute the derivatives of and in (28) and then evaluating (28) at , we can obtain The same procedure can be used to compute when . In summary, can be calculated according to where 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 using (30b). Then can be immediately obtained using (30a).

2.2.2. Computing through Derivatives

We provide another simple way to compute which also avoids the unfavorable polynomial division. That is obvious. Therefore, we just need to obtain the remaining coefficients . Notice that as the residues 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 .

(1)   and . We first discuss the calculation of when and . Here and are the th zero and the th pole, respectively. From (1) and (2), we observe that the coefficients can be obtained according to where Here as defined in (1) and . The th derivative of can be easily obtained as where . As 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 From (31) and (34), we can easily observe that Equation (34) has a similar formulation as (3); thus, can be calculated and using the procedure we compute in Section 2.1. Accordingly, we can finally obtain where Here we mention that in (36a) means . This is a Matlab denotation which will also be used in the remaining part for simplicity of expression. We can first obtain using (36a) and (36b); then, can be easily obtained using (35). Notice the evaluation of (36b) requires and . Thus other two conditions remain to be discussed.

(2)  . Here we discuss the condition when is one of the zeros of (Let us say ). Actually it will be seen that it is more preferable to let . First we rewrite as where is the multiplicity of and Computing the th derivative of both sides of (37) using Leibniz’s formula and then setting as yield As is of factorized form, its derivatives can be computed recursively where Combining (31) and (39), we have As we did in Section 2.2.2(1), to further simply the implementation, we can also introduce an auxiliary variable which satisfies And then using the procedure that we proposed in Section 2.1 leads to where is as defined in (40b). It can be easily observed that the larger is, the less computation load the method will involve. Thus can be chosen as . Particularly, if , we have

(3)  . Last, let us discuss the condition when . Here is the th pole of . Multiplying both sides of (2) with gives where and are as defined in (4) and (32), respectively. According to (45), we have Based on Leibniz’s formula, we can obtain where Let Then according to (46), (47), and (49), we have 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 where are as defined in (11b) and are residues related to the th pole. It is worthy to mention that we do not have to calculate all the here if we stored them when calculating the residues in Section 2.1.

Till this point, the method for the computation of is discussed for any . Generally, it is more preferable to let 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 where , are the residues and coefficients of residual polynomial, respectively; is the number of poles; are the poles of ; are the multiplicities of ; are the polynomial coefficients of the numerator; and 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 , and the degree of the denominator is , . For a proper function, are zeros. It can be observed that can also be expressed as where Suppose can be expanded in as Inserting (55) into (53) and then compare (53) with (52), we can observe that the residues of can be obtained as and the residual polynomial coefficients can be calculated according to In the following part, we will first show how to compute and then to calculate .

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 of ; second, calculate the residues of . First, let us suppose can be expanded in as Notice that 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 in (11b). Referring to (11a) and (11b), the following formula is immediately obtained where

It has been proved in [17] that the residues of and those of satisfy Based on the above formula, all the residues of can be obtained. Then using (56), the residues of can be obtained. Here we develop novel simpler formulas to derive the residues of based on the residues of . It can be observed that the numerator of can be rewritten as According to (52) and (61), we can denoteBy the above denotation, we have . Suppose can be expanded in as Notice that . Thus according to (58), we easily obtain the residues of If we have already obtained the residues of , , referring to [17], the residues of can be obtained as . Therefore, according to (62b) we have Notice that . As can already be calculated using (59a) and (59b), we can finally obtain the residues of using formulas (64a) and (64b) successively by increasing the value of from 1 to . Equations (64a) and (64b) can be further developed and yield another beautiful formula: In (64a), (64b), and (65), if or , the residues and do not exist and we can simply set them as zero. Formula (65) can be proved by induction as below.

Proof of Formula (65). When , (65) is reduced to (64a). Thus it is true. Suppose (65) is true when ; namely, Then according to (64b), when , we have Applying (66) to (67) yields From (68), we have where Equation (70a) can further lead to From (69) and (70b), we can finally obtain Equation (71) is (65) with replaced by ; thus, by induction, it follows that (61) is true.

Setting in (65) yields Equation (72) presents a direct relation between the residues of and the residues of . It can be useful when calculating a particular residue of .

3.2. Computing Residual Polynomial Coefficients of Functions in Expanded Form

We will propose novel methods to calculate the residual polynomial coefficients, namely, 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 of through Laurent Expansion

As shown in (53), is a sum of , if we can obtain the residual polynomial coefficients of , , can be then obtained using (57). Obviously, will contribute to the of only when . Thus is assumed to be no less than in this part. The problem is then reduced to how to calculate the of (). Notice all the 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 where Here and is a constant. Let in (73), we have From (73) and (75), we can observe that Notice that are the residual polynomial coefficients of . It is clear that if we have already obtained the coefficients of , we can immediately obtain the coefficients of by (76). It hereby follows that if we have already obtained the coefficients of , , we can successively obtain . From (76), the following formula can be obtained: Equation (77) presents the relationship between the () and , which can be seen more clearly in Table 1. Here we mention that we calculate the with decreasing from to , while the method in [17] calculates with increasing from to . It will be seen that the method proposed here is simpler and more convenient.

According to (57) and (77), we have From the above analysis, the whole problem is reduced to how to obtain . As is of factorized form, can be calculated using the methods we proposed in Section 2.2.1 through Laurent expansion. Thus, we have where In summary, the residual polynomial coefficients can be obtained by first using (79a) and (79b) and then using (78).

3.2.2. Computation of through Derivatives

As we discussed in Section 3.2.1, we just have to calculate the residual polynomial coefficients of , and then we can obtain . As is of factorized form, it is also possible to use the method we proposed in Section 2.2.2 to calculate . In Section 2.2.2, three conditions are considered. Here we only discuss the most favorable condition, (namely, is a zero of ) because firstly, this condition is most suitable as it requires the least computation and secondly, in most practical cases. Notice that . Thus we can use a formula like (44) to calculate when is not equal to any pole of . When is equal to a pole of , we just need to make a simple modification. is obvious. Using the procedures in Section 2.2.2(2), we can finally have where

3.3. Recommendations on the Choices of Procedures

In Sections 3.1 and 3.2, several procedures are proposed for the computation of both and . As the calculation of can have a relation with that of , 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 . As we can see from (80) and (81), the method in Section 3.2.2 requires that we first obtain the residues of . Therefore, we use (60) (rather than (64a) and (64b)) to obtain the residues because the residues of 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 . It is more preferable to use (79a) and (79b) to calculate as it is not based on the residues of .

Method 1. The following procedures comprise Method 1.(a) Calculate the residues of : where .(b) Calculate the residues of using and then obtain the residues of using .(c) Calculate the residues of using (80) and (81).(d) Calculate the residual polynomial coefficients using (78).

Method 2. The following procedures comprise Method 2.(a) Calculate the residues of the same way as the first step of Method 1.(b) Calculate the residues of . Here notice that . Consider (c) Calculate the residual polynomial coefficients of 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, and are set as zeros, if not mentioned particularly. For simplicity of expression, for large-scale functions of , we denote the coefficients of as follows: array of poles, ; array of the multiplicities of poles, ; array of zeros, ; array of the multiplicities of zeros, ; polynomial coefficients of the numerator of in expanded form, . 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 Section 2. 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: Obviously, the desired expansion is given by where

4.1.1. Calculation of Residues

The following residues can be calculated in a direct way: Using (11a) and (11b), we can calculate the remaining residues recursively at . Notice that −1 is the second pole of . First, calculate using (11b) Then using (11a), Thus we have

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 be −1, the second pole of . Referring to (30a) and (30b), we have Thus, we have

(b) Method through Derivatives. Here we use the method proposed in Section 2.2.2 to compute . Obviously, the second zero of , 1, has the largest multiplicity of 3. Thus it is chosen as here. As , we can use (44) to derive the residual polynomial coefficients Thus we have The complete expansion is now given by

Example 2. We consider here a large-scale factorized improper with ; ; ; . 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 directly. If we use a method designed for rational functions in expanded form, the reformulation of the 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 and that of its , , together to validate the accuracy of the expansion results. And we also calculate the relative errors of at different values of as   . The results are shown in Figures 1 and 2. The figures of the functions are plotted within the interval . 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 , is in perfect agreement with that of (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 are considered: (a) (pole of , method in Section 2.2.2(3) used); (b) (zero of , method in Section 2.2.2(2) used); (c) (neither pole nor zero of , 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.

Example 3. In this example, we consider a large-scale factorized improper rational function with ill-conditioned poles. ; ; ; . The degree of the numerator is 34. The degree of the denominator is 48. 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 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.

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 the of functions in expanded form proposed in Section 3: It can be estimated that Notice the polynomial coefficients of the numerator are . 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

(b) Use (60) to calculate the residues of , based on the residues of . The results are shown in Table 2.

Then use (56) to obtain the residues of .

(c) From Table 2, we can find the residues of , . where Using (80) and (81), we have Then referring to (78) or Table 1, we can obtain all the as shown in Table 3.

(d) Using (78), we have , , , .

Method 2. (a) Use (59a) and (59b) to expand .

(b) Using (64a) and (64b) by increasing the value of from 0 to , we can then obtain the residues of as

(c) Obtain the residual polynomial coefficients of (here, ), using (79a) an (79b) where is proper residual fraction. Then using (78) or Table 1, we can obtain all the as shown in Table 3.

(d) Finally, using (78), we have , , , .

Thus the desired expansion is

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]. ; ; . 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.

Example 6. In this example, we consider a proper rational function in expanded form with ill-conditioned poles. This example is also used in [17]. ; ; . Thus 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 as 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.

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 and complexity are not included. Let be the degree of numerator of , the degree of denominator, as number of poles, , and as the number of zeroes. Table 4 presents counts of long operations above one-degree complexity for a proper function.

The 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 Table 5. 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 algorithm (see Algorithm 1) using (11a) and (11b) to calculate residues and (30a) and (30b) to calculate residues to achieve the best performance.

Step  1. calculate directly
For
.
End
Step  2. calculate    using (11a) and (11b)
For
 For
  
 End
 For
  
  For
   
  End
  
 End
End
Step  3. Calculate   using (30a) and (30b)
Step  3.1. = 1
Step  3.2. calculate using (30b)
 For
  
 End
Step  3.3. Calculate   using (30a)
 For
   
 End

For functions in expanded form, both methods perform well for the proper function (Example 6). 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.

Step  1. calculate residues of using (59a) (59b).
Step  1.1. calculate directly
For
.
End
Step  1.2. calculate   using
For
  For
   
  End
  For
   
   For
    
   End
   
  End
End
Step 2. calculate the residues of using (64a) and (64b)
.
For
   For
    For
    
    End
   End
End
Step 3. calculate the residual polynomial coefficients of , using (79a) (79b).
Step  3.1. = 1
Step  3.2. calculate using (79b)
 For
   
 End
Step  3.3. Calculate   using (79a)
 For
    
End
Step 4. calculate the residual polynomial coefficients of using .

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 to estimate quantitatively the accuracy of the expansion results. This criterion is valid because the original function and its 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 . On the other hand, small relative errors can guarantee good expansion results. When using the defined relative error to validate a method, we suggest that should not have the values that make 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 improper 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 contains ill-conditioned poles and when the numerator of 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 for 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 . Whether in factorized or expanded form, the expansion results using the proposed methods are comparable. Meanwhile, we mention that the methods for factorized can be more easily performed by manual calculation.

Finally, and 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 and of are not zeros. Moreover, though a careful choice of and may improve the accuracy to some extent in certain cases, it will not lead to a significant improvement of accuracy in most cases. Thus, and 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 the 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 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.

Appendix

See Algorithms 1 and 2.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by The National Natural Science Foundation of China (81101049, 61271071), Shanghai Pujiang Talent Program (12PJ1401200), Doctoral Fund of Ministry of Education (20110071 120019), and Fudan University ASIC and System State Key Laboratory Project (11MS007).