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)];
endHi,
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 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 roots finding, complex finding MATLAB Answers — New Questions
