UDF: Roots of 1D Polynomials of Degree N and Real or Complex Coeff


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.
 
Ad

Advertisements

M

monir

Hello;

1) I've just converted the two Fortran routines Zroots() and Laguer() to VBA.
The driver is the array Function Zroots (a, m, polish), and it is entered on
the w/s into a block of cells F11:G13:
{=Zroots(B11:C14,B9,TRUE)}
The function Zroots() returns #VALUE! to the 6 cells!!

2) Here's the general layout of the procedure: (I would be glad to post the
VBA code of the two short routines)

Option Base 1
Function Zroots(a, m, polish)
....... declaration & my code1
...................................................
j =
ad =
x(1) =
x(2) =
Call Laguer(ad, j, x(1), x(2), its)
....... my code 2
...................................................
roots(,) =
Zroots = roots
End Function

Sub Laguer(a, m, xx1, xx2, its)
....... declaration & my code3
If x(1) = x1(1) And x(2) = x1(2) Then xx1 = x(1): xx2 = x(2): Return
....... my code4
...................................................
MsgBox "Too many iterations in Sub Laguer. Very unusual"
Return
End Sub

3) The custom Function Zroots() is used in w/s cells and it includes actions
such as calling procedure Sub Laguer(). These actions don't necessarily
execute while XL is calculating/updating the w/s.

4) Could this be the cause of the problem forcing the array function
Zroots() to return the #VALUE! error value ??
Although I'm reasonably confident about the math involved, I'll continue to
re-check the code just in case.

5) Surprisingly, the Sub Laguer name doesn't appear under:
Tools::Macro::Macros::Macro name !! Why ??

Would appreciate your expert help.

Regards.


monir said:
(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.
 
S

SteveM

Hello;

1) I've just converted the two Fortran routines Zroots() and Laguer() to VBA.
The driver is the array Function Zroots (a, m, polish), and it is entered on
the w/s into a block of cells F11:G13:
{=Zroots(B11:C14,B9,TRUE)}
The function Zroots() returns #VALUE! to the 6 cells!!

2) Here's the general layout of the procedure: (I would be glad to post the
VBA code of the two short routines)

Option Base 1
Function Zroots(a, m, polish)
...... declaration & my code1
..................................................
j =
ad =
x(1) =
x(2) =
Call Laguer(ad, j, x(1), x(2), its)
...... my code 2
..................................................
roots(,) =
Zroots = roots
End Function

Sub Laguer(a, m, xx1, xx2, its)
...... declaration & my code3
If x(1) = x1(1) And x(2) = x1(2) Then xx1 = x(1): xx2 = x(2): Return
...... my code4
..................................................
MsgBox "Too many iterations in Sub Laguer. Very unusual"
Return
End Sub

3) The custom Function Zroots() is used in w/s cells and it includes actions
such as calling procedure Sub Laguer(). These actions don't necessarily
execute while XL is calculating/updating the w/s.

4) Could this be the cause of the problem forcing the array function
Zroots() to return the #VALUE! error value ??
Although I'm reasonably confident about the math involved, I'll continue to
re-check the code just in case.

5) Surprisingly, the Sub Laguer name doesn't appear under:
Tools::Macro::Macros::Macro name !! Why ??

Would appreciate your expert help.

Regards.

monir said:
(This is a cross-post.)

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 degreepoly
has (N+1) coefficients and N roots.
I apologize for the lengthy thread, but I thought a detailed descriptionis
warranted since a reliable UDF would be widely used .
Thank you kindly.

A VBA function only returns a single value. But you are trying to
return an array. Convert the function to a procedure and explicitly
assign each root to a cell.

SteveM
 
M

monir

Steve;

Function Zroots() is an array function returning values to multiple cells;
F11:G13 in this case.

Regards.


SteveM said:
Hello;

1) I've just converted the two Fortran routines Zroots() and Laguer() to VBA.
The driver is the array Function Zroots (a, m, polish), and it is entered on
the w/s into a block of cells F11:G13:
{=Zroots(B11:C14,B9,TRUE)}
The function Zroots() returns #VALUE! to the 6 cells!!

2) Here's the general layout of the procedure: (I would be glad to post the
VBA code of the two short routines)

Option Base 1
Function Zroots(a, m, polish)
...... declaration & my code1
..................................................
j =
ad =
x(1) =
x(2) =
Call Laguer(ad, j, x(1), x(2), its)
...... my code 2
..................................................
roots(,) =
Zroots = roots
End Function

Sub Laguer(a, m, xx1, xx2, its)
...... declaration & my code3
If x(1) = x1(1) And x(2) = x1(2) Then xx1 = x(1): xx2 = x(2): Return
...... my code4
..................................................
MsgBox "Too many iterations in Sub Laguer. Very unusual"
Return
End Sub

3) The custom Function Zroots() is used in w/s cells and it includes actions
such as calling procedure Sub Laguer(). These actions don't necessarily
execute while XL is calculating/updating the w/s.

4) Could this be the cause of the problem forcing the array function
Zroots() to return the #VALUE! error value ??
Although I'm reasonably confident about the math involved, I'll continue to
re-check the code just in case.

5) Surprisingly, the Sub Laguer name doesn't appear under:
Tools::Macro::Macros::Macro name !! Why ??

Would appreciate your expert help.

Regards.

monir said:
(This is a cross-post.)

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.

A VBA function only returns a single value. But you are trying to
return an array. Convert the function to a procedure and explicitly
assign each root to a cell.

SteveM
 
P

Peter T

SteveM said:
A VBA function only returns a single value. But you are trying to
return an array.

Indeed a UDF can return a array to cells

Select a similar size array of cells, with the formula in the input bar
Array enter with Ctrl-shift-Enter

Regards,
Peter T
 
M

monir

Hello;

1) I've ironed out ALL the problems with the VBA procedure .... at last!
In addition to accessing the VBA atpvbaen.xls library (as per Dana DeLouis'
suggestion), I've also replaced the input array argument in the array
Function Zroots2() from a 2D array to two 1D arrays.

2) The VBA procedure now works perfectly and as desired. I've successfully
tested the procedure on numerous 1D polynomials for up to degree 20 with real
and complex coefficients. ALL real and complex roots are determined
correctly with almost no effort from the user. Just enter the deg of the
poly and the poly coeffs on a w/s.
(The max poly deg of 20 is only for ease formatting of the w/s. The VBA
works fine for any 1D poly of any deg.)

If there's interest, please let me know and I'd be glad to post the VBA
procedure with a demo.

Thank you all for your help in resolving the issue.

Kind regards.
 
Ad

Advertisements

G

Graeme Dennes

I would be pleased to see the code for the VBA procedure you've developed,
along with the demo. I'm sure others would too.

Graeme

----------------
 
M

monir

Graeme;

I'd be glad to post a copy of my self-contained demo w/b with the VBA
procedure.

1) I'm not sure if I could post the w/b as an Attachment in this MS DG. Do
you know how ??
If this posting feature is not available, please mail me your work email
address.
2) Are you interested in accessing the VBA code or in the use of the
procedure ??
3) Which version of Excel are you running ??
I've a demo version for XL 2007 and another for XL 2003 (and possibly
earlier). There're slight differences in the developed VBA code.

Regards.
 
G

Graeme

Hi monir,

1. I don't know how to attach files, or if it can be done, and I could not
find the answer in the Help.
2. Re the code: I'm interested mainly in the VBA code from an algorithmic
perspective, but having the procedure to run some tests and demos would be
most useful (private research purposes only).
3. Excel 2007.
4. my email address is gdennes_at_hotmail_dot_com Let me know if this is
not clear!!

Thank you for offering to make it available.

Graeme

---------------------------
 
M

monir

Graeme;

I've just emailed you a copy of my demo w/b (XL 2007).
Please acknowledge receipt.

Regards.
 
Ad

Advertisements

G

Graeme

monir,

I have received the file. Unfortunately, I won't be able to spend any time
with it for the next eight weeks (away from home).

Thank you for your generous offer to share it.

Graeme

---------------------
 
Ad

Advertisements


Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Top