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

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

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

2ut2-u=0.

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

m2ut2-(cu)+au=f.

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

Figure contains an axes. The axes with title Geometry With Edge Labels Displayed contains 5 objects of type line, text.

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

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

Figure contains an axes. The axes contains 2 objects of type line.

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

  • u(x,0)=arctan(cos(πx2)).

  • ut|t=0=3sin(πx)exp(sin(πy2)).

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

Figure contains an axes. The axes contains an object of type patch.

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

Для просмотра документации необходимо авторизоваться на сайте