M
monir
(This is a cross-post.)
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.
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.