Волновое уравнение в квадратной области

Этот пример показывает, как решить волновое уравнение с помощью solvepde функция.

Стандартное волновое уравнение второго порядка

Чтобы выразить это в форме тулбокса, обратите внимание, что solvepde функция решает задачи вида

Таким образом, стандартное волновое уравнение имеет коэффициенты $$m=1$$, $$c=1$$, $$a=0$$, и $$f=0$$.

c = 1;
a = 0;
f = 0;
m = 1;

Решите задачу в квадратной области. The squareg функция описывает эту геометрию. Создайте model объект и включать геометрию. Постройте график геометрии и просмотрите метки ребер.

numberOfPDE = 1;
model = createpde(numberOfPDE);
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on'); 
ylim([-1.1 1.1]);
axis equal
title 'Geometry With Edge Labels Displayed';
xlabel x
ylabel y

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

specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);

Установите нулевые граничные условия Дирихле на левой (ребро 4) и правой (ребро 2) и нулевые граничные условия Неймана на верхней (ребро 1) и нижней (ребро 3).

applyBoundaryCondition(model,'dirichlet','Edge',[2,4],'u',0);
applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);

Создайте и просмотрите конечный элемент, mesh для задачи.

generateMesh(model);
figure
pdemesh(model);
ylim([-1.1 1.1]);
axis equal
xlabel x
ylabel y

Установите следующие начальные условия:

u0 = @(location) atan(cos(pi/2*location.x));
ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));
setInitialConditions(model,u0,ut0);

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

Задайте время решения как 31 с равными интервалами точек во времени от 0 до 5.

n = 31;
tlist = linspace(0,5,n);

Установите SolverOptions.ReportStatistics от model на 'on'.

model.SolverOptions.ReportStatistics ='on';
result = solvepde(model,tlist);

457 successful steps
38 failed attempts
992 function evaluations
1 partial derivatives
113 LU decompositions
991 solutions of linear systems

u = result.NodalSolution;

Создайте анимацию, чтобы визуализировать решение для всех временных шагов. Сохраните фиксированную вертикальную шкалу путем вычисления сначала максимального и минимального значений u в любое время, и масштабировать все графики, чтобы использовать эти $$z$$- пределы по осям.

figure
umax = max(max(u));
umin = min(min(u));
for i = 1:n
    pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous','Mesh','off');
    axis([-1 1 -1 1 umin umax]); 
    caxis([umin umax]);
    xlabel x
    ylabel y
    zlabel u
    M(i) = getframe;
end

Чтобы воспроизвести анимацию, используйте movie(M) команда.