Efficient Recursive Methods for Partial Fraction Expansion of General Rational Functions

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.


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 highorder 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.
Suppose () = ()/() is the rational function to be expanded inpfe.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 pfe 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 pfe methods are designed for () with () written in expanded form.However, one may still be motivated to design pfe methods for factorized ().For one thing, many functions are originally written in factorized form.A pfe 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 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 2 Journal of Applied Mathematics for factorized () are also mentioned in [2].In a word, it is of significance to design methods for the pfe of () in both factorized and expanded form.
The most well-known pfe method is Heaviside' coverup 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][8][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][13][14][15][16][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 algebraic pfe 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.

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  = ∑  =1   ; the degree of the denominator is  = ∑  =1   ;  =  − .In a pfe 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.[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

Computing Residues of Factorized Functions. Brugia
where  =   −  and As shown in [4], the Kth derivatives of   () are where  > 0 and Here    = !/( − )!!,  = ( − 1) ⋅ ⋅ ⋅ 1, and 0! = 1.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 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 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)  =   − 1 : −1 : 1 means  =   −1,   −2,   −3, . . ., 2, 1.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   − 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.

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.

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/, (1) transforms into where z = 1/(  −   ), s = 1/(  −   ), and  is a constant.Now write where R() does not include the factor .In (13), applying a Taylor expansion to R() at  = 0, we can obtain where  is an integer.Now expand () as where  2 () 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 R() 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  −   = 1/ in   (), we have where Here z = 1/(  −   ), s = 1/(  −   ), and By Taylor expanding R () at  = 0 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   (), R () and R() 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 D and R in (28) and then evaluating (28) at  = 0, 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   = 1 is obvious.Therefore, we just need to obtain the remaining coefficients  0 ,  1 ,  2 , . . .,  −1 .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  0 ,  1 ,  2 , . . .,  −1 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) ℎ = 1 :  means ℎ = 1, 2, 3, . . ., .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 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.

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;  0 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 pfe 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  ℎ .

Computing Residues of Functions in Expanded
Form.The method for computation of residues can be divided into two steps: first, calculate the residues of  0 (); second, calculate the residues of ().First, let us suppose  0 () can be expanded in pfe as Notice that  0 () 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   = 0 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  −1 () 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  0 ().It can be observed that the numerator of () can be rewritten as According to (52) and (61), we can denote By the above denotation, we have   = ().Suppose   can be expanded in pfe as Notice that  0 () = 1/().Thus according to (58), we easily obtain the residues of  0 If we have already obtained the residues of  −1 ,  (−1) , referring to [17], the residues of ( −  0 ) −1 can be obtained as  (−1)(+1) +(  − 0 )× (−1) .Therefore, according to (62b) we have Notice that   =   .As  0 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 + 1 >   or  +  −  >   , the residues  0(+1) and  0(+−) do not exist and we can simply set them as zero.Formula (65) can be proved by induction as below.
Setting  =  in (65) yields Equation (72) presents a direct relation between the residues of  0 () and the residues of ().It can be useful when calculating a particular residue of ().

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.

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 s = 1/(  −   ) and  is a constant.Let  =  + 1 in (73), we have From ( 73) and (75), we can observe that  ℎ =  (+1)(ℎ+1) , (ℎ = 0 :  − ) .
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).

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.2requires 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   ().(a) Calculate the residues of  0 (): (84) (c) Calculate the residual polynomial coefficients of   () using (79a) and (79b).(d) Calculate the residual polynomial coefficients using (78).

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.

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 ] . (87)
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 largescale 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 pfe, pfe(), together to validate the accuracy of the expansion results.And we also calculate the relative errors of pfe() at different values of  as |(() − pfe())/()| (() ̸ = 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 () (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)   = 8 (pole of (), method in Section 2.2.2(3) used); (b)   = 3.5 (zero of (), method in Section 2.2.2(2) used); (c)   = 0 (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. = [1000, 0, 0.1, 8, 9];  = [8,7,7,6,6];  = [−1, −2, 7, 11, 12, 990];  = [8,8,8,8,8,8].The degree of the numerator is 34.The degree of the denominator is 48.() contains three ill-conditioned high-order poles: two highorder 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, 10 50 ), 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.Relative error Figure 2: Relative errors of pfe with residual polynomial coefficients calculated using the methods in Sections 2.     Then referring to (78) or Table 1, we can obtain all the  ℎ as shown in Table 3.
Thus the desired expansion is  () = −32 + 12 − 4 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]. = [8,7,6,5,4,3,2,1];  = [1,2,3,4,5,6,7,8];  = [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. 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]. = [1000, 3, 2, 1, 1.1];  = [7,3,3,6,7];  = [0.1,0.2, . . ., 2.2].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 [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.

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 (),  = ∑  =1   the degree of denominator,  as number of poles,  = ∑  =1  2  , and  as the number of zeroes.Table 4 presents counts of long operations above onedegree 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 pfe algorithm (see Algorithm 1) using (11a) and (11b) to calculate residues and (30a) and (30b) to calculate residues to achieve the best performance.
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.

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 |(()−pfe())/()| 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

Journal of Applied Mathematics
Step Step 3.1.  = 1 Step 3. intolerable errors when functions with ill-conditioned highorder 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 pfe 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 Step 1. calculate residues of r 0 (s) using (59a) (59b).
Step   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 () can be more easily performed by manual calculation.
Finally,  0 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  0 and   of () are not zeros.Moreover, though a careful choice of  0 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,  0 and   can generally be simply set as zeros or other values at your convenience in practice.

Conclusions
In this paper, we developed efficient recursive methods for the pfe 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.

Method 1 .
The following procedures comprise Method 1.

Figure 1 :
Figure 1: Comparison of pfe 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.

Figure 3 :
Figure 3: Comparison of pfe with reference solution (a) and relative error of pfe (b) in Example 3.

Figure 4 :
Figure 4: Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example 5.

Figure 5 :
Figure 5: Comparison of expansion results of the proposed methods with reference (a) and relative errors of Method 1 (b) and Method 2 (c) in Example 6.

Algorithm 2 :
Algorithm for functions in expanded form.
ℎ 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,   =  0 (namely,   is a zero of ()) because firstly, this condition is most suitable as it requires the least computation and secondly,   =  0 = 0 in most practical cases.Notice that  > .Thus we can use a formula like (44) to calculate  ℎ when  0 is not equal to any pole of   ().When  0 is equal to a pole of   (), we just need to make a simple modification.  = 1 is obvious.Using the procedures in Section 2.2.2(2), we can finally have

Table 4 :
Counts of long operations above one-degree complexity for a proper function.

Table 5 :
Calculation time of large-scale examples.