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

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

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

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

Таким образом, стандартное уравнение волны имеет коэффициенты $$m=1$$, $$c=1$$, $$a=0$$, и $$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

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

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