Azzera filtri
Azzera filtri

Why do I receive Interval might contain additional eigenvalues / with pdetool ?

6 visualizzazioni (ultimi 30 giorni)
Hello !
I am trying to solve partial differential equation with neumann and dirichlet condition, on a rectangle pierced by another rectangle, but I have this error after launching my algorithm.
Error using pde.PDEModel/solvepdeeig (line 69)
Interval might contain additional eigenvalues.
Error in rectangle_in_rectangle (line 53)
result2=solvepdeeig(model2,evr2);
evr2 is in the interval where I am searching the eigenvalues. In order to try to solve the problem, I compute a new interval n-times, which is growing from [0 0.1] to [0 1] for instance. But it did not solve the problem.
Does someone had the same problem, or any idea which can help me ?
Thanks in advance !
I post my code below if you want to check :
function rectangle_in_rectangle(n)
% cercle troué d'un carré
mus=zeros(1,n);
lambdas=zeros(1,n);
for i = 1:n
R1 = [3,4,-i,i,i,-i,i,i,-i,-i]';
c_r2=sqrt((2*i)^2-1); %côté du second carré pour que l'aire soit égale à 1
R2 = [3,4,-c_r2/2,c_r2/2,c_r2/2,-c_r2/2,c_r2/2,c_r2/2,-c_r2/2,-c_r2/2]';
gm = [R1,R2];
sf = 'R1-R2';
ns = char('R1','R2');
ns = ns';
g = decsg(gm,sf,ns);
sz=size(g); %afin de savoir le nombre de edges pour appliquer les conditions sur tous les bords
% c'est le nombre de colonnes qui nous intéresse, car chaque colonne
% correspond à une arrête
%model pour vp
model1=createpde(1);
geometryFromEdges(model1,g);
applyBoundaryCondition(model1,'neumann','edge',1:sz(2),'g',0,'q',0);
generateMesh(model1);
pdeplot(model1);
specifyCoefficients(model1,'m',0,'f',0,'c',1,'a',0,'d',1);
evr=[0.01,3];
result=solvepdeeig(model1,evr);
mu=result.Eigenvalues;
while isempty(mu)
disp("pas de vp trouvée dans l'intervalle");
evr(2)=evr(2)*2;
result=solvepdeeig(model1,evr);
mu=result.Eigenvalues;
end
mus(i)=mu(1);
%dirichlet
model2=createpde(1);
geometryFromEdges(model2,g);
applyBoundaryCondition(model2,'dirichlet','edge',1:sz(2),'u',0);
generateMesh(model2);pdeplot(model2);
specifyCoefficients(model2,'m',0,'f',0,'c',1,'a',0,'d',1);
evr2=[0,0.05];
result2=solvepdeeig(model2,evr2);
lambda=result2.Eigenvalues;
while isempty(lambda)
disp("pas de vp trouvée dans l'intervalle");
evr2(2)=evr2(2)*1.2;
result2=solvepdeeig(model2,evr2);
lambda=result2.Eigenvalues;
end
lambdas(i)=lambda(1);
end
%tracé du graphe
plot(lambdas,mus,'.')
hold on
%calcul du premier 0 de la fonction de bessel
f=@(x) besselj(0,x);
x=fzero(f,1);
%Inequalities from theorem 3.3 (The diagram (λ1, µ1))
%Inequality 1
plot(lambdas,(lambdas).^(-1)*(pi^4)/4,'.')
%Inequality 2
plot(lambdas,(lambdas).^(-1)*((pi*x)^2)*9,'.')
%Rayleigh-Faber-Krahn inequality
plot([pi*x^2 pi*x^2],[1 30])
legend('Diag cercle troué','first inequality','second inequality','Rayleigh-Faber-Krahn')
hold off
xlabel('1*λ(Ω)')
ylabel('1*µ(Ω)')
title("Diagramme (λ,µ), pour des un rectangle troué d'un carré");
end

Risposte (1)

Jaynik
Jaynik il 25 Gen 2024
Hi Thibault,
I tried the code you have provided and was successfully able to reproduce the error in R2023a. Upon debugging and stepping into "solvepdeeig", I found that the function internally calls "sptarn" on line 66. The error message "Interval might contain additional eigenvalues." is triggered when the variable "ires" is less than zero. On opening this function, the documentation suggests to "Try again with a smaller interval".
I see that you have tried increasing the range till results are found but this skips some upper bound for which the eigen values are actually found. Using binary search, I was able to find the upper bound to be around 2480 to 2483 where the solver successfully identifies eigenvalues, but increasing the upper bound beyond this range leads to the error. You can use this code to obtain the result:
evr2 = [0, 2483]; % Can try with 2480, 2481, 2482
result2 = solvepdeeig(model2, evr2);
lambda = result2.Eigenvalues;
I tested the some code in R2023b and found that the issue does not occur. It's possible that this behavior was addressed in the latest release. You can try upgrading MATLAB to resolve the issue.
Hope this helps!

Prodotti


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by