Finding Imaginary roots of a function

15 次查看（过去 30 天）
University
University 2023-9-27
评论： University 2023-9-28
Ran in:
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

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

采纳的回答

Torsten
Torsten 2023-9-27
编辑：Torsten 2023-9-27
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.
  1 个评论
University
University 2023-9-28
Hi Torsen,
Thank you for pointing this out. I have changed the initial guess and is working perfectly now.

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by