Finding roots of a complex function

2 次查看（过去 30 天）
University
University 2025-6-17
评论： John D'Errico 2025-6-18
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
  6 个评论
显示 4更早的评论
University
University 2025-6-17
K1 ^ (-0.1e1 / 0.2e1) * eta1 ^ (-0.1e1 / 0.2e1) * tau ^ (-0.1e1 / 0.2e1)
This is taken directorly from MAPLE.
John D'Errico
John D'Errico 2025-6-18
The crappiest code is often that generated by computer algebra.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Just for fun 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by