Müller's Method (find a complex root of an analytic function)

The three functions below will attempt to find a complex root to an analytic function. While Igor's ` FindRoots ` operation will provide complex roots to polynomials with real coefficients, it does not work for other functions with complex roots. The code is light on user-interface features and has not been extensively tested. Please email suggestions or comments.

Müller's method uses three initial evaluations of a function but does not require the derivative of the function. Essentially, the method fits a parabola through the last three points and uses the appropriate intersection with zero as its next value. Unlike other methods, you can start from a real guess and still find a complex root. The code is based on Numerical Recipes, 2nd ed. See also the Wikipedia entry.

Example:

``` •print sin(pi/3) 0.866025 •print /C Muller(cmplx(1,0.1),cmplx(1,0.2),cmplx(1,0.3)) // 1st root of z^3+1=0 (<pre><code class="language-igor"> FindRoots ``` would work on this example) (0.5,0.866025) •print /C Muller(cmplx(1,-0.1),cmplx(1,-0.2),cmplx(1,-0.3)) // 2nd root of z^3+1=0 (0.5,-0.866025) •print /C Muller(cmplx(-0.7,0),cmplx(-0.8,0),cmplx(-0.9,0)) // 3rd root of z^3+1=0 (the obvious one) (-1,0) •print /C Muller(cmplx(0,-0.6),cmplx(0,-0.7),cmplx(0,-0.8)) // root of sqrt(z+1+i)+1 (this is something ` FindRoots ` can't do) (3,-1)

Function/C Muller(x0,x1,x2)
Variable/c x0,x1,x2                                                    // initial values for root iteration
Variable tol=1e-9,nmax=100,k=2                                // tol = tolerance; nmax = max iterations (adjust to taste)
Make /c/d/o/n=(nmax+3) xx=0,f=0                 // temporary waves [keep when debugging]
xx[0]=x0;  xx[1]=x1;  xx[2]=x2                      // initial values in waves
f[0] = myfunc(x0);  f[1] = myfunc(x1);  f[2] = myfunc(x2)       // and their function values
Do
xx[k+1] = Muller_step(xx[k],xx[k-1],xx[k-2],f[k],f[k-1],f[k-2])
f[k+1] = myfunc(xx[k+1])                       // Muller requires only one functional eval /iteration
k += 1
While ((k<nmax+3) && (cabs(xx[k]-xx[k-1]) > tol))
// deletepoints k+1,nmax-k+2, xwave,fwave          // delete "unused" points and keep waves for debug
Variable /c root = xx[k]
Killwaves xx,f                                 // delete waves (normal operation)
Return root
End

Function/C Muller_step(xk,xk1,xk2,fk,fk1,fk2)                  // from Numerical Recipes, 2nd ed.
Variable/C xk,xk1,xk2,fk,fk1,fk2                           // initial values
Variable/C q,A,B,C                             // use complex quantities (/c flag)
q = (xk-xk1)/(xk1-xk2)
A = q*fk-q*(1+q)*fk1+q^2*fk2
B = (2*q+1)*fk - (1+q)^2*fk1 + q^2*fk2
C = (1+q)*fk
Variable/c  d = Sqrt(Cmplx(B^2-4*A*C,0))                   // it's this step that can make a real guess go complex
if (Magsqr(B+d) > Magsqr(B-d))
Return xk - (xk-xk1)*2*C/(B+d)
else
Return xk - (xk-xk1)*2*C/(B-d)
endif
End

Function/c myfunc(z)                                //  your analytic function
Variable/c z
//  Return z^3+1                                // FindRoots could actually do this one
Return Sqrt(z+Cmplx(1,1))-2                 // but it can't do this one
End

Muller's method solver:
Bech has posted a useful snippet to find complex roots. I have takent the liberty of adding to its flexibility and making it more consistent with other Igor operation syntax:
(1) added a parameter wave to argument list
(2) passed a user specified function to Muller() using FUNCREF definition (with added prototype function)
The latter usage means that you may not have to alter either the Muller() code or the function code, if you wish to use a different function.

Bech's earlier Muller_step() code remains unchanged.

I have also shown an example function (TM mode of a dielectric waveguide) that Muller is capable of solving. This example is used within other functions, so a History snippet example is not informative.

Function/C Muller(pwcplx, x0, x1, x2, fin)  //  modified from Igor Exchange code snippet by Bech
wave/C pwcplx       // ADDED complex parameter wave for input function
Variable/C x0,x1,x2     // initial values for root iteration
FUNCREF myprotofunc fin // ADDED flexibility in changing input function without changing code
Variable tol=1e-9,nmax=100,k=2                // tol = tolerance; nmax = max iterations (adjust to taste)
Make /C/D/O/N=(nmax+3) xx=0,wf=0    // temporary waves [keep when debugging]; renamed 'wf'
xx[0]=x0;  xx[1]=x1;  xx[2]=x2      // initial values in waves
wf[0] = fin(pwcplx,x0);  wf[1] = fin(pwcplx,x1);  wf[2] = fin(pwcplx,x2)    // and their function values
Do
xx[k+1] = Muller_step(xx[k],xx[k-1],xx[k-2],wf[k],wf[k-1],wf[k-2])
wf[k+1] = fin(pwcplx,xx[k+1])   // Muller requires only one functional eval /iteration
k += 1
While ((k<nmax+3) && (cabs(xx[k]-xx[k-1]) > tol))
// deletepoints k+1,nmax-k+2, xwave,fwave    // delete "unused" points and keep waves for debug
Variable /C root = xx[k]
Killwaves xx,wf             // delete waves (normal operation)
Return root
End
//------------------------------------------------------------------------------------------------------------------------------------
function/C myprotofunc(parmw, cx)
wave/C      parmw
variable/C  cx
end
//------------------------------------------------------------------------------------------------------------------------------------
function/C fTM(pwC, xi)               //    new parameter wave argument
wave/C pwC                //    complex parameter wave
variable/C xi                 //    xi = sigma*d
variable/C  Kb  =   pwC[0]
variable/C  Ka  =   pwC[1]
variable/C  Kd  =   pwC[2]
variable    rho =   real(pwC[3])
variable    d   =   real(pwC[4])
return  atan( (Ka/Kb)*sqrt( (Ka-Kb)*(rho/xi)^2 - 1 ) )  + atan( (Ka/Kd)*sqrt( (Ka-Kd)*(rho/xi)^2 - 1 ) ) - xi
end

Forum

Support

Gallery