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 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 10

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More