Complex Variables: Practical Applications by D. James Benton

Author:D. James Benton [Benton, D. James]
Language: eng
Format: azw3
Publisher: UNKNOWN
Published: 2019-01-16T00:00:00+00:00

Traub-Jenkins Method

Real and complex roots of a polynomial having complex coefficients are found with the Traub-Jenkins Method.[22] This algorithm is similar to Bairstow's in that a Newton iteration is used to solve for the roots. Besides the complex coefficients, the biggest difference between the two is that the Traub-Jenkins method yields a single complex root with each step. Multiple initialization attempts are also used to handle the convergence problem.

The residual (i.e., complex error associated with each provisional root) is calculated and iteratively minimized. The polynomial is also scaled in order to better estimate the root. This approach typically finds the roots in order of increasing magnitude, but this is not always the case, so don't count on it. As is often the case, convergence is roughly quadratic, making this a fast algorithm.

The algorithm passes various modified forms of the ever-shrinking polynomial to several functions through global variables. In the earliest implementations, these arrays were static; however, in the code provided here (traubj.c), the arrays are dynamically allocated, eliminating any constraing on the degree. The core algorithm is listed below:

int TraubJenkins(double*a,int ndeg,double*z)

{

BOOL conv;

int i,icnt1,icnt2,ii,inx,inxi,j,n2,nn2;

double bnd,cosr,sinr,xx,xxx,yy,zi,zr;

if(ndeg<1)

return(__LINE__);

tj.are=DBL_EPSILON;

tj.rmre=2.*M_SQRT2*DBL_EPSILON;

xx=M_SQRT1_2;

yy=-xx;

cosr=-0.06975647;

sinr=0.9975641;

tj.nn=ndeg+1;

if(fabs(a[0])>FLT_MIN||fabs(a[1])>FLT_MIN)

{

while(1)

{

nn2=tj.nn+tj.nn;

if(fabs(a[nn2-2])>FLT_MIN||fabs(a[nn2-1])>FLT_MIN)

break;

inx=ndeg-tj.nn+2;

inxi=inx+ndeg;

z[inxi-1]=z[inx-1]=0.;

tj.nn--;

if(tj.nn==1)

return(0);

}

for(i=1;i<=tj.nn;i++)

{

ii=i+i;

tj.pr[i-1]=a[ii-2];

tj.pi[i-1]=a[ii-1];

tj.shr[i-1]=modulus(tj.pr[i-1],tj.pi[i-1]);

}

bnd=ScaleFactor(tj.nn,tj.shr);

if(bnd!=1.)

{

for(i=0;i<tj.nn;i++)

{

tj.pr[i]*=bnd;

tj.pi[i]*=bnd;

}

}

while(tj.nn>2)

{

for(i=0;i<tj.nn;i++)

tj.shr[i]=modulus(tj.pr[i],tj.pi[i]);

bnd=LowerBound(tj.nn,tj.shr,tj.shi);

for(icnt1=1;icnt1<=2;icnt1++)

{

dPdH(5);

for(icnt2=1;icnt2<=9;icnt2++)

{

xxx=cosr*xx-sinr*yy;

yy=sinr*xx+cosr*yy;

xx=xxx;

tj.sr=bnd*xx;

tj.si=bnd*yy;

PolyZero2(10*icnt2,&zr,&zi,&conv);

if(conv)

goto next;

}

}

return(__LINE__);

next:

inx=ndeg+2-tj.nn;

inxi=inx+ndeg;

z[inx-1]=zr;

z[inxi-1]=zi;

tj.nn--;

for(i=0;i<tj.nn;i++)

{

tj.pr[i]=tj.qpr[i];

tj.pi[i]=tj.qpi[i];

}

}

ComplexDivide(-tj.pr[1],-tj.pi[1],tj.pr[0],tj.pi[0],&z[ndeg-1],&z[ndeg+ndeg-1]);

for(i=1;i<=ndeg;i++)

tj.pi[i-1]=z[ndeg+i-1];

n2=ndeg+ndeg;

j=ndeg;

for(i=1;i<=ndeg;i++)

{

z[n2-2]=z[j-1];

z[n2-1]=tj.pi[j-1];

n2-=2;

j--;

}

return(0);

}

return(__LINE__);

}

This code can be found in the folder \examples\Traub-Jenkins. Typical output is as follows:

a=(1+2i),(3+4i),(5+6i),(7+8i)

z=(-0.120851+1.7188i),(-1.75701+0.207073i),(-0.322141-1.52588i)

a=(-9.959+8.467i),(-3.666-3.501i),(9.169+5.724i),(1.478-0.643i)

z=(-0.0832328+0.111289i),(0.762881+0.306844i),(-0.719835-0.803841i)

a=(-3.039-5.537i),(-4.295-1.856i),(-6.72+6.827i),(-0.039-9.509i)

z=(0.410517-0.484651i),(-1.60023-0.370736i),(0.604929+1.31012i)

a=(-7.005+1.942i),(-5.173-4.564i),(2.39+4.604i),(-6.098-9.847i)

z=(0.15504+0.936837i),(0.694713-0.680799i),(-1.36779-1.05119i)

a=(-9.708+2.382i),(7.421+8.716i),(9.718+9.895i),(-4.553-8.275i)

z=(0.517487+0.0547907i),(-1.0417-0.329409i),(1.03744+1.29836i)

a=(4.771+1.538i),(-8.131+9.912i),(-4.334-3.702i),(7.035-0.106i)

z=(0.650561+0.184551i),(-0.760413-0.427239i),(1.04699-2.13696i)

a=(-1.298-6.19i),(1.321+0.332i),(7.673-5.336i),(5.141-2.289i)

z=(-0.497039-0.170566i),(-0.36336+1.09371i),(0.954641-1.11679i)

a=(-1.748-3.132i),(-4.454-2.357i),(2.661+2.756i),(-9.964+2.859i)

z=(0.971482+0.806244i),(-0.0407683-0.989109i),(-2.10971+0.94695i)

a=(-1.277-0.259i),(-2.472-9.222i),(2.316-6.965i),(-7.811-8.158i)

z=(0.0362865+0.807379i),(-0.837316-1.44827i),(-2.46508-5.91829i)

a=(-9.712+0.105i),(-0.96-1.058i),(9.264-7.353i),(-2.555-6.196i)

z=(-0.0692952+0.47559i),(-1.13555+0.0243696i),(1.10719-0.609953i)

Download

Copyright Disclaimer:
This site does not store any files on its server. We only index and link to content provided by other sites. Please contact the content providers to delete copyright contents if any and email us, we'll remove relevant links or contents immediately.