Вибрация круговой мембраны

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

Вычисление режимов вибрации требует решения дифференциального уравнения с частными производными собственного значения. Этот пример сравнивает решение, полученное при помощи 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);

Сгенерируйте mesh.

generateMesh(model,'Hmax',0.2);

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

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.06,  New conv eig=  0
              Basis= 13,  Time=   0.08,  New conv eig=  2
              Basis= 16,  Time=   0.09,  New conv eig=  2
              Basis= 19,  Time=   0.10,  New conv eig=  2
              Basis= 22,  Time=   0.13,  New conv eig=  3
              Basis= 25,  Time=   0.14,  New conv eig=  3
              Basis= 28,  Time=   0.16,  New conv eig=  5
End of sweep: Basis= 28,  Time=   0.17,  New conv eig=  5
              Basis= 15,  Time=   0.20,  New conv eig=  0
End of sweep: Basis= 15,  Time=   0.20,  New conv eig=  0

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

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

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

Max difference in eigenvalues from solvepdeeig and eigs: 1.080025e-12

Обе функции вычисляют те же собственные значения. Для любого собственного значения можно умножить собственный вектор на произвольный скаляр. 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')