I'm trying to solve a nonlinear eigenvalue problem in MATLAB, still without success. It's a problem about graphene plasmonics.
The nonlinear eigenvalue problem is given below:
\begin{equation} \frac{\epsilon_1}{\kappa_{1||n}}E_{x||n}+\frac{\epsilon_2}{\kappa_{2||n}}E_{x||n}+ \frac{i}{\omega \epsilon_0} \sum_{p=-M}^{M}\hat{\sigma}_{n-p} E_{x||p}=0 \end{equation}
\begin{equation} \kappa_{m||n}=\sqrt{(q+nG)^2-\epsilon_m (\omega/c)^2} \end{equation}
The parameters $\epsilon_1$, $\epsilon_2$, $\epsilon_0$, $G$, $c$, $\hat{\sigma}_{j}$ are known. The problem stands for determining $q$ (wavevector), $\omega$ (frequency), so as for the system of the first equation to has nontrivial solution. This system of course is a $(2M+1) \times (2M+1)$ square system of the form $Ax=0$, where $A=A(\omega, q)$. As far as I know, I should give a specific value to $\omega$, and then solve the equation $det(A)=0$ to find $q$. Here is my MATLAB 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)
The main function (heart of the code) is 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=0.5*10^-12; %[sec]
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 large cutoff integer $M$, the determinant becomes extremely big (or small, if I scale properly the equations with a constant number). Anyway, the problem seems unstable. Is there a way to fix my code? Or is there a way to treat the problem in a different mathematical way? My knowledge in computational linear algebra is poor.
Thank you in advance.
See corresponding post on scicomp stackexchange: https://scicomp.stackexchange.com/questions/29710/nonlinear-eigenvalue-problem-matlab-code/30668#30668