Tag: polynomials

Complex Factors of Real and Complex Polynomials – part 2

Last week I posted a Maxima utility to factor polynomials using complex roots.  The idea was to use the built-in solver solve to find each root a with multiplicity m, and then construct a factor of the form (z-a)^m.  Typically, solve can fail to identify some roots that are returned by the more robust to_poly_solve.  My problem at that time was that I didn’t know how to access the multiplicities for the roots returned by to_poly_solve.

I’ve written a utility to determine multiplicities for the roots returned by to_poly_solve. With that, here’s an updated version of factorC that takes advantage of the more robust solver.

/* factor a polynomial using all its complex roots */
factorC(_f,_z):=block(
[s,n,m,fp,j],
fp:1,
ss:to_poly_solve(_f,_z),
s:create_list(args(ss)[k][1],k,1,length(ss)),
m:tpsmult1(_f,_z,ss),
/*These lines were needed before I figured out how
 how to handle multiplicities with to_poly_solve
s:solve(_f,_z),
m:multiplicities,*/
n:length(s),
for j:1 thru n do 
 if lhs(s[j])#0 
 then fp:fp*(_z-(rhs(s[j])))^m[j],
 fp:fp*divide(_f,fp)[1],
fp
);
Advertisement

Roots, Multiplicities, and to_poly_solve

A few days ago I posted some code to factor and perform partial fraction decomposition using all available complex roots.  The approach is to find roots a — taking note of the multiplicity m of each —  and then construct factors of the form (z-a)^m.

As I’ve noted before, the built-in command solve can fail to find solutions.  The more robust Maxima solver is to_poly_solve.    That said, after a call to to_poly_solve, I don’t know how to get my hands on the multiplicities.  So, I’ve written a little utility to compute and return a list of multiplicities for the roots returned by to_poly_solve.  Here’s an example, followed by the code.  Notice that because to_poly_solve returns solutions using the %union construction, extracting the jth root requires something like root:rhs(part(s,j,1))

tpsmult1

/* tpsmult1 assumes results of to_poly_solve applied to _f, 
 are passed in the argument s */
tpsmult1(_f,_z,s):=block(
[d,root,nroots,mult,m1],
nroots:length(s),
mult:[],
for j:1 thru nroots do (
root:rhs(part(s,j,1)),
m1:rootmult(_f,_z,root),
mult: endcons(m1,mult)
),
mult
);

/* find the multiplicity of a root _a for the polynomial _f */
rootmult(_f,_z,_a):=block(
[d,val,m,mm],
m:1,
d:divide(_f,(_z-_a)),
val:subst(_z=_a,part(d,1)),
if (val=0) then (
 mm:rootmult(part(d,1),_z,_a),
 m:m+mm
 ),
m
);

Partial Fraction Decomposition Using Complex Factors

A few days ago I wrote about my solution for factoring a polynomial into factors corresponding to complex roots,  factorC.

One of the uses for that utility is to perform a partial fractions decomposition of rational functions.  Here’s a complex partial fraction function that uses factorC and the built-in Maxima function partfrac.

partfracC

partfracC(_f,_z):=block(
[d,fd],
d:denom(_f),
fd:factorC(d,_z),
partfrac(1/fd,_z)
);

Complex Factors of Real and Complex Polynomials

I’m getting ready for my fall Complex Variables class.  I noticed that the built-in Maxima function residue doesn’t reliably do the right thing.  My goal is to make some improvements to  Maxima residue calculations in Maxima over the course of the next month.

As I started to look at some test cases, I realized I didn’t know how to factor a polynomial into complex factors.  In the simplest case:

realfactor

but I wanted to see

manual_complexfactor

Maybe someone will find this and let me know that there’s a simple way to make that happen using existing Maxima commands.  Until then, I’ve written a little utility  to identify the roots (both real and complex) of a polynomial and return a factorization.  Here’s an example, and the code.  Notice that it only works as well as the root finder solve.  I tried to upgrade to the more robust to_poly_solve,  but I don’t yet know how to handle multiplicities in that case.

factorCexample

factorC(_f,_z):=block(
[s,n,m,fp,j],
fp:1,
/* This commented code was meant to use the
more robust solver to_poly_solve, but 
I couldn't understand how to handle multiplicities
ss:args(to_poly_solve(_f,_z)),
s:create_list(ss[k][1],k,1,length(ss)),*/
s:solve(_f,_z),
m:multiplicities,
n:length(s),
for j:1 thru n do 
 if lhs(s[j])#0 
 then fp:fp*(_z-(rhs(s[j])))^m[j],
 fp:fp*divide(_f,fp)[1],
fp
);