Finding Imaginary roots of a function

Hi,

I have been trying to find a roots roots of a function. I received an error: Error using fsolve Objective function is returning undefined values at initial point. FSOLVE cannot continue.

Please can help out? I have attached my m files.

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.5:3;
xivals=-0.3:0.1: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');

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=-50:0.1:50;
        %tauIvals=0.1*ones(size(tauIvals));
        ntau=length(tauIvals);
        
        % plot equation (projected onto imag line) to solve for these values of u and xi
        figure(1)
        y=zeros(size(tauIvals));
        for ii=1:ntau
            tau= tauIvals(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(tauIvals,y);
        hold on
        plot(tauIvals,zeros(size(tauIvals)),'-k');
        xline(0);
        yline(0);
        xlabel('tau');
        ylabel('tau equation');
        %axis([min( tauIvals) max( tauIvals) -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=501 %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=[tauIvals(k),tauIvals(k)];
            tauinit
            fun(tauinit)
            % 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 (imag part of) using equation from Maple
            % file
            wavenum(k)=imag(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=imag(tausol_found);
        tauIvals=tauIvals(flag==1);
        wavenum=wavenum(flag==1);
        
        % plot imag part of tau solution versus initial tau
        figure(2)
        plot(tauIvals,tausolR,'g-')
        hold on
        plot(tauIvals,tauIvals,'k-')
        xlabel('initial tau');
        ylabel('tau solution');
        hold off
        
        % calculate min value of imag 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)
        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('u');
        ylabel('xi');
        title('minimum tau');
        drawnow
        
        % filled contour plot of minimum tau value (negative tau means instability)
        figure(4)
        contourf(U,Xi,wavenummin,10)
        colormap(map)
        colorbar
        xlabel('u');
        ylabel('xi');
        title('wavenumber at minimum tau');
        drawnow
        
        % stability domain in (u,xi) plane
        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('u');
        ylabel('xi');
        title('Blue = stable, Yellow = unstable');
        drawnow
        
    end
    
    % display time taken and percentage complete
    toc
    disp(['Progress: ' num2str(round(100*(i*(j-1)+j)/(nu*nxi))) ' % completed']);
    
end

tauinit = 1×2

0 0

ans = 1×2

NaN 0

Error using fsolve
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

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

Answer 1

[Torsten]

As you can see from your code above, you choose tauinit = [0 0] for k = 501, and your function returns NaN for the real part of your function. So obviously, this is not a good initial guess.

[University]

Hi Torsen,

Thank you for pointing this out. I have changed the initial guess and is working perfectly now.