M

#### monir

Hello;

I would very much appreciate your expert help in developing a UDF in XL VBA

to determine the real and imaginary roots of a polynomial of degree N and

real or complex coefficients.

1) This subject was already discussed using Goal Seek and Solver in the

thread:

"Goal Seek with Complex Numbers" located at:

http://www.mrexcel.com/forum/showthread.php?t=321367&goto=newpost

Goal Seek did not work with complex numbers, and I had limited success with

XL Solver. By reviewing the above thread, one would conclude that the Solver

procedure for this task had run its course and there would be no added value

in pursuing it any further. Hence the idea of developing a general and more

reliable UDF procedure.

Here're some of my thoughts on the subject and I would appreciate your

comments and suggestions.

2) One of the apparent difficulties in using the Solver procedure (item 1.

above) is that its success, if at all, crucially depends on having a good

first guess of the root, which by no means is readily available for a general

equation. Equally problematic is the fact that the user of the Solver

procedure has no control on hunting down the root even if it is accurately

bracketed. Meaning, there's no guarantee against having the (Solver and

similar) iterative procedure converges repeatedly to the same (nonmultiple)

root instead of separately to different roots.

3) But, which root-finding method one should implement in the UDF ?? In most

reliable algorithms based on Newton-Raphson method a user-supplied

root-bracketing is required, either by a subsidiary procedure or by providing

an array of guesses and let the procedure zooms in each root to within a

specified convergence criterion. To the best of my knowledge, almost all

Newton-based methods deal with real polynomial coefficients and determine

real roots.

4) I would suggest the Laguerre's method (code is provided in item 9 below).

It is guaranteed to converge to all types of roots in the - oo to + oo range,

handles polynomials with complex coefficients, and does not require an

initial guess or starting point or trial solutions. (Keep in mind, with

complex coefficients, complex roots may or may not occur in conjugate pairs.)

Sounds too good to be true ?? Well, it's actually true! and I successfully

used the method in the past.

5) The ref. NR, p. 264, 265 has the Laguer routine (~ 30 lines of code) and

its deriver routine Zroots (~ 30 lines). Both codes (provided in item 9

below) are in Fortran 77 and are relatively easy to follow and convert to VBA.

The UDF FUNCTION CubicEq() written by pgc01, MrExcel MVP, and posted

recently (in the thread quoted in item 1. above) could be used as a template

for the procedure. It is really that simple, specially for someone with

expertise in XL VBA.

6) The developed root-finding UDF VBA procedure would be applicable to any

one-dimensional equation of the form:

f(x) = sum [k=1 to k=N+1] A(k) x^(k-1)

where f(x) has only one independent variable "x", N is the degree, and

A(N+1) are the complex coefficients of the polynomial.

7) One should always try to get some idea of how the function behaves before

trying to find its roots. There's really no excuse for not to, since it is

one dimensional and the task can't be more simpler in XL.

The display (on the w/s) would identify where the function changes sign and

whether the change in sign is associated with a real root, a finite jump

discontinuity, or a singular point. Also, potential problems such as double

real roots (value of function is zero at its max or min), multiple real roots

in close proximity (oscillating function about the x-axis), and other

problems, could easily be identified from a simple display automatically

generated on the w/s by the UDF (just a thought).

8) One could validate the developed UDF using the analytical data for 3rd

and 4th degree equations. There're no general analytical solutions for

polynomials of degree higher than 4, though there're (very complicated)

solutions for particular families of polynomials of any degree. One can't,

however, justify the effort required in pursuing such solutions. One may

instead manufacture such solutions by multiplying a lower degree equation

(with known roots) by a known real or imaginary root.

9) Here're the two Fortran codes:

-------------------------------------------------------------------

SUBROUTINE Zroots (a,m,roots,polish)

! given the degree m and the complex coefficients a(1:m+1) of the polynomial:

! f(x) = a(1) + a(2) x^(1) + a(3) x^(2) + ... + a(m+1) x^(m)

! this driver routine successively calls Laguer routine and finds all m

complex roots.

! The logical variable polish should be input as:

! .true. if polishing (also by Laguerreâ€™s method) is desired,

! .false. if the roots will be subsequently polished by other means.

INTEGER m,MAXM

REAL EPS

COMPLEX a(m+1),roots(m)

LOGICAL polish

PARAMETER (EPS=1.e-6,MAXM=101)

INTEGER i,j,jj,its

COMPLEX ad(MAXM),x,b,c

do 11 j=1,m+1

ad(j)=a(j)

enddo 11

do 13 j=m,1,-1

x=cmplx(0.,0.)

call laguer(ad,j,x,its)

if(abs(aimag(x)).le.2.*EPS**2*abs(real(x))) x=cmplx(real(x),0.)

roots(j)=x

b=ad(j+1)

do 12 jj=j,1,-1

c=ad(jj)

ad(jj)=b

b=x*b+c

enddo 12

enddo 13

if (polish) then

do 14 j=1,m

call laguer(a,m,roots(j),its)

enddo 14

endif

do 16 j=2,m

x=roots(j)

do 15 i=j-1,1,-1

if(real(roots(i)).le.real(x)) goto 10

roots(i+1)=roots(i)

enddo 15

i=0

10 roots(i+1)=x

enddo 16

return

END

-------------------------------------------------------------------

SUBROUTINE Laguer (a,m,x,its)

! given the degree m and the complex coefficients a(1:m+1) of the polynomial:

! f(x) = sum [k=1 to k=m+1] a(k) x^(k-1)

! and given a complex value x, this routine improves x by Laguerreâ€™s method

until

! it converges, within the achievable roundoff limit, to a root of the given

polynomial

INTEGER m,its,MAXIT,MR,MT

REAL EPSS

COMPLEX a(m+1),x

PARAMETER (EPSS=2.e-7,MR=8,MT=10,MAXIT=MT*MR)

INTEGER iter, j

REAL abx,abp,abm,err,frac(MR)

COMPLEX dx,x1,b,d,f,g,h,sq,gp,gm,g2

SAVE frac

DATA frac /.5,.25,.75,.13,.38,.62,.88,1./

do 12 iter=1,MAXIT

its=iter

b=a(m+1)

err=abs(b)

d=cmplx(0.,0.)

f=cmplx(0.,0.)

abx=abs(x)

do 11 j=m,1,-1

f=x*f+d

d=x*d+b

b=x*b+a(j)

err=abs(b)+abx*err

enddo 11

err=EPSS*err

if(abs(b).le.err) then

return

else

g=d/b

g2=g*g

h=g2-2.*f/b

sq=sqrt((m-1)*(m*h-g2))

gp=g+sq

gm=g-sq

abp=abs(gp)

abm=abs(gm)

if(abp.lt.abm) gp=gm

if (max(abp,abm).gt.0.) then

dx=m/gp

else

dx=exp(cmplx(log(1.+abx),float(iter)))

endif

endif

x1=x-dx

if(x.eq.x1) return

if (mod(iter,MT).ne.0) then

x=x1

else

x=x-dx*frac(iter/MT)

endif

enddo 12

pause â€™too many iterations in Laguer. Very unusual'

return

END

10) Having Excel input/output complex values in two separate columns is a

practical option. As suggested in item 5 above, the referenced UDF array

function FUNCTION CubicEq() by pgc01 located at:

http://www.mrexcel.com/forum/showthread.php?t=321367&goto=newpost

is an excellent template as a starting point. The polynomial complex

coefficients (each in two cells) could be passed on as two dimensional array

and the root results returned as two dimensional array. The 2D arrays could

be dimensioned based on the degree of the polynomial since an Nth degree poly

has (N+1) coefficients and N roots.

I apologize for the lengthy thread, but I thought a detailed description is

warranted since a reliable UDF would be widely used .

Thank you kindly.