Уравнение волны в квадратной области

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

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

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

Поэтому стандартное волновое уравнение имеет коэффициенты $$m\mathrm{=}$$1$$,\mathrm{c}$$= $$1\mathrm{,}$$a = 0 $$\mathrm{}$$и f = 0.

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

Решите проблему на квадратном домене. 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

Укажите коэффициенты PDE.

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);

Создайте и просмотрите сетку конечных элементов для проблемы.

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) команда.