Nonlinear eigenvalue problem - Problem with fsolve

4 次查看（过去 30 天）
Nikolaos Xenidis
Nikolaos Xenidis 2018-3-28
编辑： Nikolaos Xenidis 2018-3-28
Hello everybody. I'm trying to solve a problem about graphene plasmonics, and I'm dealing with a non-linear eigenvalue problem. More specifically, the eigenvalue problem is decribed by these equations:
It is convenient for numerical puproses to scale these equations into this form (that is, dividing the whole set with h_bar*c):
These equations form an infinite system (as n runs from -infinite to infinite). By limiting n to -M<=n<=M and cutting the other terms, we have a finite system (N+1) X (N+1) of equations. The physical interpetation of the problem is that for a given frequency ω, there is a matced wavenumber q, which is given by the above system, that is setting to zero the determinant of the system. All other parameters are known, but there is also some kind of complexity in the σ terms, as they are function of the frequency (the σ terms are Fourier series expansion of the optical conductivity of graphene, as in this problem, the conductivity is periodically modulated).
The whole system has to be solved for a spectrum of frequencies (that is, multiple times for different frequencies ω). The idea is to form the matrix of the system, and set the determinant to zero, searching for the wavenumber q by using fsolve. This process is repeated for every frequency.
Here is my code
global w % angular frequency - it is global because I pass it in the function fun
c0=299792458; %speed of light
f=linspace(25*10^12,45*10^12,1000); %frequency spectrum
guess=32*2*pi*25*10^12/c0;; % initial guess for the wavenumber
for x=1:1000 % scanning for every frequency in the spectrum
w=2*pi*f(x);
q(x)=fsolve(@fun,guess); % find the wavenumber q via fsolve
guess=q(x); % next guess is the previous solution (wavenumber)
end
plot(q,f)
and the function fun (that is the main code) is given below
function F=fun(q)
global w
hbc=3.161*10^-26; % hbar*c0, SI units, [J*m]
c0=299792458; %speed of light
e0=8.854*10^-12; %vacuum permittivity [F/m]
e1=1; %dielectric permittivity, air
e2=2.1; %dielectric permittivity, spacer
elcharge=1.602*10^-19; % elementary charge, [C]
kbol=1.38*10^-23; %Boltzmann constant, SI units [J/K]
hbar=1.054*10^-34; %reduced Planck constant, SI units, [J*s]
h=hbar*2*pi; % Planck constant
T=300; % Temperature [K]
mc1=0.2*elcharge; %chemical potential of 1st region, [J]
mc2=0.5*elcharge; %chemical potential of 2nd region, [J]
w1=10*10^-9; %width of 1st region, [m]
w2=60*10^-9; %width of 2nd region, [m]
R=w1+w2; %spatial period
G=2*pi/R; %reciprocal lattice wavevector
M=1; % M is the cutoff of the infinite terms of the system, -M<=n<=M
hbw=hbar*w; %hbar * w
sg1=i*((kbol*T*elcharge^2)/(pi*(w+i/t)*hbar^2))*((mc1/(kbol*T))+2*log(exp(-mc1/(kbol*T))+1))+i*(elcharge^2/(4*pi*hbar))*log((2*mc1-hbar*(w+i/t))/(2*mc1+hbar*(w+i/t))); %1st region's conductivity (Kubo formula for graphene conductivity)
sg2=i*((kbol*T*elcharge^2)/(pi*(w+i/t)*hbar^2))*((mc2/(kbol*T))+2*log(exp(-mc2/(kbol*T))+1))+i*(elcharge^2/(4*pi*hbar))*log((2*mc2-hbar*(w+i/t))/(2*mc2+hbar*(w+i/t))); %2nd region's conductivity (Kubo formula)
for n=-M:M
k1(n+M+1)=sqrt((hbc*(q+n*G))^2-e1*hbw^2);
k2(n+M+1)=sqrt((hbc*(q+n*G))^2-e2*hbw^2);
end
sf0=(sg1*w1/R)+(sg2*w2/R);; %zeroth Fourier component
%fill the matrix of Fourier components, sf[]
for l=-2*M:2*M
if l~=0
sf(l+2*M+1)=sg1*sin(l*G*w1)/(R*l*G)+j*sg1*(cos(l*G*w1)-1)/(R*l*G)-sg2*sin(l*G*w1)/(R*l*G)-j*sg2*(cos(l*G*w1)-1)/(R*l*G);
else
sf(l+2*M+1)=sf0;
end
end
%forming the main matrix of the problem, of which the determinant we want to set to zero
for n=-M:M
v0(n+M+1)=(e1/k1(n+M+1))+(e2/k2(n+M+1))+j*sf0/(hbw*e0*c0); %diagonal elements
end
A0=diag(v0);
A=A0;
for n=1:2*M
vu=ones(1,2*M+1-n)*j*sf((4*M+2)/2-n)/(hbw*e0*c0); %upper diagonal elements
vd=ones(1,2*M+1-n)*j*sf((4*M+2)/2+n)/(hbw*e0*c0); %down diagonal elements
Au=diag(vu,n);
Ad=diag(vd,-n);
A=A+Au+Ad;
end
F=det(A);
end
For the specific problem, I need many terms (let's say up to 10) to have good convergence to the right result. Fsolve "somewhat" works with maximum M=2 and seems unable to work for further terms (up to M=3). Usually, it says "fsolve stopped because the problem appears to be locally singular". Any ideas on how to solve this issue?

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

回答（0 个）

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by