Re: Partial Fraction Decomposition with imaginary coeff.
- To: mathgroup at yoda.physics.unc.edu
- Subject: Re: Partial Fraction Decomposition with imaginary coeff.
- From: Brian Evans <evans at gauss.eedsp.gatech.edu>
- Date: Sun, 8 Mar 1992 22:54:48 GMT
You just don't know how difficult Apart has made my life.
It is crucial to the implementation of inverse linear
transforms. Under 1.2, it would not handle polynomials
with rational coefficients expressed in decimal form.
Under 1.2 and 2.0, Apart does not break up terms like
1 / (x^2 + 1). About two years ago, I had to write a general
purpose routine to work around Apart's drawbacks. I call it
MyApart. It is embedded in the signal processing packages
for Mathematica and is used by the inverse z- and Laplace
transforms when Apart doesn't complete the decomposition.
MyApart is, of course, darn slow.
Brian L. Evans
Digital Signal Processing Laboratory
School of Electrical Engineering
Georgia Institute of Technology
Atlanta, GA 30332-0250
e-mail: evans at eedsp.gatech.edu
----------------------------- cut here ----------------------------------
(* :Title: Miscellaneous Routines *)
(* :Authors: Brian Evans, James McClellan *)
(* :Summary: To provide a robust partial fractions routine. *)
(* :Context: Misc` *)
(* :PackageVersion: 2.6 *)
(*
:Copyright: Copyright 1992 by Brian L. Evans
Georgia Tech Research Corporation
Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is
hereby granted, provided that the above copyright notice
appear in all copies and that both that copyright notice and
this permission notice appear in supporting documentation,
and that the name of the Georgia Tech Research Corporation,
Georgia Tech, or Georgia Institute of Technology not be used
in advertising or publicity pertaining to distribution of the
software without specific, written prior permission. Georgia
Tech makes no representations about the suitability of this
software for any purpose. It is provided "as is" without
express or implied warranty.
*)
(* :History: *)
(* :Keywords: partial fractions decomposition *)
(* :Source: *)
(* :Warning: *)
(* :Mathematica Version: 1.2 or 2.0 *)
(* :Limitation: *)
(*
:Functions: GetRoot
GetRootList
MyApart
*)
(* B E G I N P A C K A G E *)
BeginPackage[ "Misc`" ]
(* U S A G E I N F O R M A T I O N *)
GetRoot::usage =
"GetRoot[rule] extracts the value from an expression like \
{z -> 0.}, which is 0. in this case."
GetRootList::usage =
"GetRootList[p, x] returns a list of the approximate numerical roots \
of expression p, a function of x, with infinite roots removed. \
GetRootList[p, x, filter] applies filter to the list of roots \
returned by the Solve function (defaults to N)."
MyApart::usage =
"MyApart[ rational_polynomial, x ] decomposes the rational \
polynomial into a sum of fractions whose numerators are of \
the form (x + b)^n where b is a constant and n is an integer. \
MyApart[rational_polynomial, x, filter] specifies a filter to be \
placed on the output of the Solve command used to root the \
denominator: Identity for rational and N for real-valued roots. \
The default is Identity. \
In Mathematica 1.2, MyApart is about 25 times slower than Apart."
(* E N D U S A G E I N F O R M A T I O N *)
Begin[ "Private`" ]
(* GetRoot *)
GetRoot[{}] := {} (* no roots *)
GetRoot[rule_] := rule[[1]][[2]]
(* GetRootList *)
GetRootList[p_, x_, filter_:N] :=
Select[ Map[ GetRoot, filter[ Solve[ p == 0, x ] ] ],
(! MatchQ[#1, DirectedInfinity[___]])& ]
(* MyApart -- kludge around the way Apart does partial fractions *)
(* Root denominator and replace roots with symbols *)
MyApart[ratpoly_, x_, filter_:Identity] :=
Block [ {apart, denom, denomfactored, normfact, numer,
partfrac, rootlist, rootmult, rules},
numer = Numerator[ratpoly];
denom = Denominator[ratpoly];
normfact = Last[ CoefficientList[denom, x] ];
numer /= normfact;
denom /= normfact;
rootlist = Sort[ GetRootList[denom, x, filter] ];
{ denomfactored, rules } = multiplicityform[rootlist, x];
apart = Apart[numer / denomfactored, x];
partfrac = apart /. rules;
partfrac /. ( a_. / (b_ c_) :> a / ( Together[b] c ) /;
FreeQ[b, x] && ! FreeQ[c, x] ) ]
multiplicityform[ roots_, x_ ] :=
Block [ {count = 1, cur, denom = 1, i, last,
length, sublist = {}, sym = 1},
Clear[localvar]; (* localvar is global to package *)
length = Length[roots];
last = First[roots];
For [ i = 2, i <= length, i++,
cur = roots[[i]];
If [ SameQ[ cur, last ],
count++,
denom *= (x - localvar[sym])^count;
PrependTo[ sublist, localvar[sym] -> last ];
sym++;
count = 1 ];
last = cur ];
denom *= (x - localvar[sym])^count;
PrependTo[ sublist, localvar[sym] -> last ];
{ denom, sublist } ]
End[]
EndPackage[]
Protect[ GetRoot, GetRootList, MyApart ]
Null