Вибрация кольцевой мембраны

В этом примере показано, как рассчитать режимы колебаний круговой мембраны.

Расчет режимов колебаний требует решения дифференциального уравнения собственных значений в частных производных. В этом примере сравнивается раствор, полученный с использованием solvepdeeig решатель из дифференциального Toolbox™ в частных производных и eigs решатель из MATLAB ®. Собственные значения, рассчитанные поsolvepdeeig и eigs практически идентичны, но в некоторых случаях один решатель удобнее другого. Например, eigs удобнее при вычислении заданного числа собственных значений вблизи конкретного целевого значения. В то время как solvepdeeig требует, чтобы нижняя и верхняя границы окружали эту цель, eigs требуется только целевое собственное значение и требуемое количество собственных значений.

Создайте модель PDE.

model = createpde;

Создайте геометрию окружности и включите ее в модель.

radius = 2;
g = decsg([1 0 0 radius]','C1',('C1')');

geometryFromEdges(model,g);

Постройте график геометрии с метками кромок.

figure
pdegplot(model,'EdgeLabels','on')
axis equal
title 'Geometry with Edge Labels'

Укажите коэффициенты.

c = 1e2;
a = 0;
f = 0;
d = 10;
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);

Укажите, что решение равно нулю на всех четырех внешних кромках окружности.

bOuter = applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',0);

Создайте сетку.

generateMesh(model,'Hmax',0.2);

Использовать assembleFEMatrices для вычисления глобальных матриц массы и жесткости конечных элементов с граничными условиями, наложенными с использованием подхода нулевого пространства.

FEMatrices = assembleFEMatrices(model,'nullspace');
K = FEMatrices.Kc;
B = FEMatrices.B;
M = FEMatrices.M;

Решите проблему собственного значения с помощью eigs функция.

sigma = 1e2; 
numberEigenvalues = 5;
[eigenvectorsEigs,eigenvaluesEigs] = eigs(K,M,numberEigenvalues,sigma);

Изменить форму диагонали eigenvaluesEigs матрица в вектор.

eigenvaluesEigs = diag(eigenvaluesEigs);

Найдите наибольшее собственное значение и его индекс в векторе собственных значений.

[maxEigenvaluesEigs,maxIndex] = max(eigenvaluesEigs);

Добавьте значения ограничений, чтобы получить полный собственный вектор.

eigenvectorsEigs = B*eigenvectorsEigs;

Теперь решите ту же задачу собственного значения с помощью solvepdeeig. Задайте диапазон для solvepdeeig быть немного больше диапазона от eigs.

r = [min(eigenvaluesEigs)*0.99 max(eigenvaluesEigs)*1.01];
result = solvepdeeig(model,r);

              Basis= 10,  Time=   0.17,  New conv eig=  0
              Basis= 13,  Time=   0.19,  New conv eig=  2
              Basis= 16,  Time=   0.21,  New conv eig=  2
              Basis= 19,  Time=   0.22,  New conv eig=  2
              Basis= 22,  Time=   0.24,  New conv eig=  3
              Basis= 25,  Time=   0.26,  New conv eig=  3
              Basis= 28,  Time=   0.29,  New conv eig=  5
End of sweep: Basis= 28,  Time=   0.30,  New conv eig=  5
              Basis= 15,  Time=   0.35,  New conv eig=  0
End of sweep: Basis= 15,  Time=   0.35,  New conv eig=  0

eigenvectorsPde = result.Eigenvectors;
eigenvaluesPde = result.Eigenvalues;

Сравните решения.

eigenValueDiff = sort(eigenvaluesPde) - sort(eigenvaluesEigs);
fprintf('Maximum difference in eigenvalues from solvepdeeig and eigs: %e\n', ...
  norm(eigenValueDiff,inf));

Maximum difference in eigenvalues from solvepdeeig and eigs: 4.121148e-13

Обе функции вычисляют одни и те же собственные значения. Для любого собственного значения можно умножить собственный вектор на произвольный скаляр. eigs и solvepdeeig функции могут выбрать другой произвольный скаляр для нормализации собственных векторов.

h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1)
axis equal
pdeplot(model,'XYData', eigenvectorsEigs(:,maxIndex),'Contour','on')
title(sprintf('eigs eigenvector, eigenvalue: %12.4e', eigenvaluesEigs(maxIndex)))
xlabel('x')
ylabel('y')
subplot(1,2,2)
axis equal
pdeplot(model,'XYData',eigenvectorsPde(:,end),'Contour','on')
title(sprintf('solvepdeeig eigenvector, eigenvalue: %12.4e',eigenvaluesPde(end)))
xlabel('x')
ylabel('y')