Finding roots of a complex function

Hi,

I'm trying to find the roots of a complex function z=x+iyz in MATLAB. However, it seems that the root-finding routine is missing some roots, which leads to inaccurate results. Could you advise on how I can improve the code to ensure more reliable root detection? I have the two matlab files.

Thanks,

clearvars
close all

% set parameter values
pars.gamma1=0.1093;
pars.alpha3=-0.1104e-2;
pars.K1=6e-12;
pars.d=0.2e-3;
pars.eta1=0.240e-1;
pars.chia=1.219e-6;
pars.alpha=1-pars.alpha3^2/(pars.gamma1*pars.eta1);
pars.Ha=pi*sqrt(pars.K1/pars.chia)/pars.d;

% set lists of u (field) and xi (activity) values
uvals=0:0.05:3;
xivals=-0.3:0.01:0.3;

nu=length(uvals);
nxi=length(xivals);

% initiate arrays for output
taumin=ones(nu,nxi);
wavenummin=ones(nu,nxi);

% start timer
tic
disp('Starting u and xi loops');

% start loop around u values
for i=1:nu
    
    pars.H=uvals(i)*pars.Ha;
    
    % start loop around xi value
    for j=1:nxi
        
        xi=xivals(j);
        
        % set initial tau values for root finding, tau is a complex
        % variable
        tauRvals=-50:0.1:50;
        tauIvals=0.1*ones(size(tauRvals));
        ntau=length(tauRvals);
        
        % plot equation (projected onto real line) to solve for these values of u and xi
        fig1 = figure(1);
        y=zeros(size(tauRvals));
        for ii=1:ntau
            tau=tauRvals(ii);
            y(ii) = (pars.H ^ 2 * sin(pars.K1 ^ (-0.1e1 / 0.2e1) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) * pars.d) * pars.chia * pars.d * tau ^ (0.3e1 / 0.2e1) * sqrt(pars.eta1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) + sin(pars.K1 ^ (-0.1e1 / 0.2e1) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) * pars.d) * pars.d * pars.gamma1 * sqrt(pars.eta1) * sqrt(tau) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) - 0.2e1 * sqrt(pars.K1) * cos(pars.K1 ^ (-0.1e1 / 0.2e1) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) * pars.d) * pars.alpha3 * tau ^ 2 * xi + 0.2e1 * sqrt(pars.K1) * cos(pars.K1 ^ (-0.1e1 / 0.2e1) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) * pars.d) * pars.alpha3 ^ 2 * tau + 0.2e1 * sqrt(pars.K1) * pars.alpha3 * tau ^ 2 * xi - 0.2e1 * sqrt(pars.K1) * pars.alpha3 ^ 2 * tau) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.3e1 / 0.2e1) * (pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) ^ (-0.1e1 / 0.2e1) / pars.alpha3 / sin(pars.K1 ^ (-0.1e1 / 0.2e1) * pars.eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) * pars.d);
        end
        plot(tauRvals,y);
        hold on
        plot(tauRvals,zeros(size(tauRvals)),'-k');
        xline(0);
        yline(0);
        xlabel('tau');
        ylabel('tau equation');
        axis([min(tauRvals) max(tauRvals) -1e-2 1e-2]);
        drawnow
        hold off
        
       
        % loop around initial tau values for root finding
        tausol=zeros(1,ntau);
        flag=zeros(1,ntau);
        wavenum=zeros(1,ntau);
        for k=1:ntau
            
            % set function, options and inital tau value
            fun = @(x)rootsolver_complex(x,xi,pars);
            options = optimset('TolFun',1e-15,'MaxFunEvals',1e5,'Maxiter',1e5,'Display','none');
            tauinit=[tauRvals(k),tauIvals(k)];
            
            % find root of tau equation
            [x,fval,exitflag,output] = fsolve(fun,tauinit,options);
            % save complex tau solution
            tausol(k)=complex(x(1),x(2));
            % set solve flag (if exitflag>0 the root finder has solved)
            flag(k)=(exitflag>0);
            % calulate wavenumber (real part of) using equation from Maple
            % file
            wavenum(k)=real(sqrt((pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau * xi - pars.alpha3 ^ 2 + pars.eta1 * pars.gamma1) / pars.K1 / pars.eta1 / tau) * pars.d / pi / 0.2e1);
        end
        
        tauflag=[tausol',flag'];
        tausol_found=tauflag(flag==1);
        tausolR=real(tausol_found);
        tauRvals=tauRvals(flag==1);
        wavenum=wavenum(flag==1);
        
        % plot real part of tau solution versus initial tau
        fig2 = figure(2);
        plot(tauRvals,tausolR,'g-')
        hold on
        plot(tauRvals,tauRvals,'k-')
        xlabel('initial $\tau$', 'Interpreter','latex');
        ylabel('$tau$ solution', 'Interpreter','latex');
        
        hold off
        
        % calculate min value of real part of tau and the wavenumber at
        % that min value of tau
        taumin(i,j)=min(tausolR);
        wavenummins=wavenum(tausolR==min(tausolR));
        wavenummin(i,j)=wavenummins(1);
        
        % filled contour plot of minimum tau value (negative tau means instability)
        fig3 = figure(3);
        [Xi,U] = meshgrid(xivals,uvals);
        N=[0:0.1:1];
        map = [0.95*(1-N') 0.95*(1-N') N'];
        contourf(U,Xi,taumin,[-100:10:100])
        colormap(map)
        colorbar
        xlabel('Orienting field, $u$', 'Interpreter','latex');
        ylabel('Activity, $xi$ [Pa]', 'Interpreter','latex');
        %title('minimum tau');
        drawnow
        
        % filled contour plot of minimum tau value (negative tau means instability)
        fig4 =figure(4);
        contourf(U,Xi,wavenummin,10)
        colormap(map)
        colorbar
        xlabel('Orienting field, $u$', 'Interpreter','latex');
        ylabel('Activity, $xi$ [Pa]', 'Interpreter','latex');
        %title('wavenumber at minimum tau');
        drawnow
        
        % stability domain in (u,xi) plane
        fig5 = figure(5);
        S = 25; % size of symbols in pixels
        % normalize colouring vector to go from zero to 1
        normtau = (taumin>0);
        normtau=reshape(normtau,nu*nxi,1);
        C = [0.95*(1-normtau) 0.95*(1-normtau) normtau];
        scatter(reshape(U,nu*nxi,1),reshape(Xi,nu*nxi,1),S,C,'filled','Marker','o')
        xlabel('Orienting field, $u$', 'Interpreter','latex');
        ylabel('Activity, $xi$ [Pa]', 'Interpreter','latex');
        title('Blue = stable, Yellow = unstable', 'Interpreter','latex');
        drawnow
        
    end
    
    % display time taken and percentage complete
    toc
    disp(['Progress: ' num2str(round(100*(i*(j-1)+j)/(nu*nxi))) ' % completed']);  
end

function F = rootsolver_complex(x,xi,pars)
% function to provide right-hand-side of the equation for tau

gamma1=pars.gamma1;
alpha3=pars.alpha3;
K1=pars.K1;
d=pars.d;
eta1=pars.eta1; 
chia=pars.chia;
alpha=pars.alpha;
H=pars.H;

tau=complex(x(1),x(2));

% equation taken directly from Maple file eq.mw
y = (H ^ 2 * sin(K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) * d) * chia * d * tau ^ (0.3e1 / 0.2e1) * sqrt(eta1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) + sin(K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) * d) * d * gamma1 * sqrt(eta1) * sqrt(tau) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) - 0.2e1 * sqrt(K1) * cos(K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) * d) * alpha3 * tau ^ 2 * xi + 0.2e1 * sqrt(K1) * cos(K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) * d) * alpha3 ^ 2 * tau + 0.2e1 * sqrt(K1) * alpha3 * tau ^ 2 * xi - 0.2e1 * sqrt(K1) * alpha3 ^ 2 * tau) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.3e1 / 0.2e1) * (H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) ^ (-0.1e1 / 0.2e1) / alpha3 / sin(K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1) * sqrt(H ^ 2 * chia * eta1 * tau + alpha3 * tau * xi - alpha3 ^ 2 + eta1 * gamma1) * d);
F=[real(y),imag(y)];

end

[Torsten]

Maybe you could explain a little more what you try to do in your code. It's hard to interprete all those figures coming up.

In the line

wavenum(k)=real(sqrt((pars.H ^ 2 * pars.chia * pars.eta1 * tau + pars.alpha3 * tau *...    

shouldn't "tau" be replaced by "tausol(k)" ?

[University]

The main figure of concern is:

% filled contour plot of minimum tau value (negative tau means instability)
fig3 = figure(3);
[Xi,U] = meshgrid(xivals,uvals);
N=[0:0.1:1];
map = [0.95*(1-N') 0.95*(1-N') N'];
contourf(U,Xi,taumin,[-100:10:100])
colormap(map)
colorbar
xlabel('Orienting field, $u$', 'Interpreter','latex');
ylabel('Activity, $xi$ [Pa]', 'Interpreter','latex');
%title('minimum tau');
drawnow

I am trying to solve two equations equations simultaneously:

Find the roots q (which is y in my MATLAB codes), use the results to solve for \tau. Then plot the minimum of real part of tau for each xi and u values.

[Matt J]

No, there is no reliable way to find all roots with a numerical solver. Numerical solvers will find one root, if all goes well.

[John D'Errico]

Is there a good reason why you are doing things like this in multiple places?

pars.K1 ^ (-0.1e1 / 0.2e1)

It seems to me you might as well have written

1/sqrt(pars.K1)

Of course, the former is much more impressive. And you certainly use sqrt in other places.

And then we see:

K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1)

which might as well have been written as:

1/sqrt(K1*eta1*tau)

Again, far less impressive though.

As @Matt J has already said though, there is no general way to insure you find all the roots of a nonlinear model using numerical methods. In fact, it is easy to show (by way of example) that any such numerical method will provably fail on at least some problems.

[University]

K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1)

This is taken directorly from MAPLE.

[John D'Errico]

The crappiest code is often that generated by computer algebra.