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

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

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

Создайте модель УЧП.

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.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

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