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

April 14, 2010 at 06:19 am - Permalink

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More