Документация

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

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

Вычисление режимов вибрации требует решения дифференциального уравнения с частными производными (PDE) собственного значения. В этом примере решение задачи о собственных значениях выполняется с помощью обоих Тулбокс УЧП solvepdeeig решатель и eigs базового MATLAB eigensolver.

Основная цель этого примера состоит в том, чтобы показать как eigs может использоваться с Тулбоксом УЧП. Обычно собственные значения, вычисленные solvepdeeig и eigs практически идентичны. Однако иногда просто более удобно использовать eigs чем solvepdeeig. Один пример этого - когда он желаем, чтобы вычислить конкретное количество собственных значений около заданного пользователями целевого значения. solvepdeeig требует, чтобы была задана нижняя и верхняя граница, окружающая это целевое значение. eigs требует только, чтобы целевое собственное значение и желаемое количество собственных значений были заданы.

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

numberOfPDE = 1;
model = createpde(numberOfPDE);

Создайте геометрический объект и включайте его в модель PDE. Геометрия является кругом.

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 Displayed';

Задайте коэффициенты.

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)*.99 max(eigenvaluesEigs)*1.01];
result = solvepdeeig(model,r);

              Basis= 10,  Time=   0.13,  New conv eig=  0
              Basis= 13,  Time=   0.16,  New conv eig=  2
              Basis= 16,  Time=   0.18,  New conv eig=  2
              Basis= 19,  Time=   0.20,  New conv eig=  2
              Basis= 22,  Time=   0.23,  New conv eig=  3
              Basis= 25,  Time=   0.24,  New conv eig=  3
              Basis= 28,  Time=   0.27,  New conv eig=  5
End of sweep: Basis= 28,  Time=   0.27,  New conv eig=  5
              Basis= 15,  Time=   0.32,  New conv eig=  0
End of sweep: Basis= 15,  Time=   0.32,  New conv eig=  0

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

Сравните решения, вычисленные eigs и solvepdeeig.

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: 9.094947e-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');