В этом примере показано, как вычислить режимы вибрации круговой мембраны.
Вычисление режимов вибрации требует решения дифференциального уравнения с частными производными собственного значения. Этот пример сравнивает решение, полученное при помощи 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.11, New conv eig= 0 Basis= 13, Time= 0.14, New conv eig= 2 Basis= 16, Time= 0.15, New conv eig= 2 Basis= 19, Time= 0.16, New conv eig= 2 Basis= 22, Time= 0.20, New conv eig= 3 Basis= 25, Time= 0.21, New conv eig= 3 Basis= 28, Time= 0.24, New conv eig= 5 End of sweep: Basis= 28, Time= 0.25, New conv eig= 5 Basis= 15, Time= 0.28, New conv eig= 0 End of sweep: Basis= 15, Time= 0.29, 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: 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')