Home |
Search |
Today's Posts |
#1
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex Coeff
(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/showthr...7&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 Laguerres 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 Laguerres 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/showthr...7&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. |
#2
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex Coeff
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 procedu (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" wrote: (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/showthr...7&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 Laguerres 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 Laguerres 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/showthr...7&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. |
#3
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or ComplexCoeff
On Jun 10, 1:40 pm, monir wrote:
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 procedu (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" wrote: (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/showthr...7&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/showthr...7&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 |
#4
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
Steve;
Function Zroots() is an array function returning values to multiple cells; F11:G13 in this case. Regards. "SteveM" wrote: On Jun 10, 1:40 pm, monir wrote: 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 procedu (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" wrote: (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/showthr...7&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 Laguerres 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 Laguerres 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/showthr...7&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 |
#5
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex Coeff
"SteveM" wrote in message
news:e3e3f6a6-5758-4806-8c64- 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 |
#6
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
#7
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
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 ---------------- "monir" wrote: 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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
#8
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex Coeff
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. "Graeme Dennes" wrote: 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 ---------------- "monir" wrote: 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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
#9
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
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 --------------------------- "monir" wrote: 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. "Graeme Dennes" wrote: 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 ---------------- "monir" wrote: 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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
#10
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
Graeme;
I've just emailed you a copy of my demo w/b (XL 2007). Please acknowledge receipt. Regards. "Graeme" wrote: 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 --------------------------- "monir" wrote: 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. "Graeme Dennes" wrote: 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 ---------------- "monir" wrote: 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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
#11
Posted to microsoft.public.excel.programming
|
|||
|
|||
UDF: Roots of 1D Polynomials of Degree N and Real or Complex C
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 --------------------- "monir" wrote: Graeme; I've just emailed you a copy of my demo w/b (XL 2007). Please acknowledge receipt. Regards. "Graeme" wrote: 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 --------------------------- "monir" wrote: 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. "Graeme Dennes" wrote: 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 ---------------- "monir" wrote: 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. "Peter T" wrote: "SteveM" wrote in message news:e3e3f6a6-5758-4806-8c64- 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 |
Reply |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Forum | |||
Polynomials | Charts and Charting in Excel | |||
how make an angle to show in degree, minute,second instead of degree | Excel Programming | |||
Imputting an angle degree. Degree, minute, second | Excel Discussion (Misc queries) | |||
Imputting an angle degree. Degree, minute, second | Excel Worksheet Functions | |||
Convert decimal degree (lattitude/longitude) into Degree, | Excel Discussion (Misc queries) |