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