정사각형 영역에서의 파동 방정식

이 예제에서는 solvepde 함수를 사용하여 파동 방정식을 푸는 방법을 보여줍니다.

표준 2차 파동 방정식은 다음과 같습니다.

$$\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot \nabla u=0.$$

이를 solvepde 함수가 풀 수 있는 문제 형식인, 툴박스의 다음 형식으로 표현합니다.

$$m\frac{{\partial }^{2}u}{\partial t^2}-\nabla \cdot (c\nabla u)+au=f.$$

따라서 표준 파동 방정식의 계수는 $$m=1$$, $$c=1$$, $$a=0$$, $$f=0$$입니다.

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

정사각형 영역에서 문제를 풉니다. squareg 함수는 이 지오메트리를 설명합니다. model 객체를 만들고 지오메트리를 포함시킵니다. 지오메트리를 플로팅하고 모서리 레이블을 표시합니다.

model = createpde;
geometryFromEdges(model,@squareg);
pdegplot(model,EdgeLabels="on"); 
xlim([-1.1 1.1]);
ylim([-1.1 1.1]);
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)에는 디리클레 경계 조건을 0으로 설정하고 위쪽(모서리 1)과 아래쪽(모서리 3)에는 노이만 경계 조건을 0으로 설정합니다.

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

다음 초기 조건을 설정합니다.

• $$u(x,0)=\arctan \left(\cos \left(\frac{\pi x}{2}\right)\right)$$.

• $${\frac{\partial u}{\partial t}|}_{t=0}=3\sin (\pi x)\exp \left(\sin \left(\frac{\pi y}{2}\right)\right)$$.

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

이렇게 설정하면 높은 진동 모드로 에너지가 들어가는 것을 피하고 적절한 시간 스텝 크기를 사용할 수 있습니다.

풀이 시간을 0 ~ 5의 시간 범위에서 균일한 간격으로 배치된 점 31개로 지정합니다.

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

문제를 풉니다.

result = solvepde(model,tlist);
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]); 
    xlabel x
    ylabel y
    zlabel u
    M(i) = getframe;
end

애니메이션을 재생하려면 movie(M) 명령을 사용합니다.