상태 종속성이 강한 질량 행렬에 대한 ODE 풀기

이 예제에서는 이동 메시 기법 [1]을 사용하여 버거스 방정식을 푸는 방법을 보여줍니다. 문제에 질량 행렬이 포함되고, 질량 행렬의 강한 상태 종속성과 희소성을 고려하도록 옵션이 지정되어 풀이 과정이 더 효율적입니다.

문제 설정

버거스 방정식은 PDE로 지정되는 대류-확산 방정식입니다.

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}=ϵ\frac{{\partial }^2\mathit{u}}{{\partial \mathit{x}}^2}-\frac{\partial }{\partial \mathit{x}}\left(\frac{{\mathit{u}}^2}{2}\right),0<\mathit{x}<1,\mathit{t}>0,ϵ=1{\times \mathrm{10}}^{-4}.$$

좌표 변환([1]의 방정식 18)을 적용하면 다음과 같이 좌변에 항이 추가됩니다.

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}-\frac{\partial \mathit{u}}{\partial \mathit{x}}\frac{\partial \mathit{x}}{\partial \mathit{t}}=ϵ\frac{{\partial }^2\mathit{u}}{{\partial \mathit{x}}^2}-\frac{\partial }{\partial \mathit{x}}\left(\frac{{\mathit{u}}^2}{2}\right).$$

$\mathit{x}$에 대한 편도함수를 근사하기 위해 유한 차분을 사용하여 PDE를 일변수 ODE로 변환할 수 있습니다. 유한 차분이 $\Delta$로 작성되었다면 PDE는 $\mathit{t}$에 대한 도함수만 포함하는 ODE로 다시 작성할 수 있습니다.

$$\frac{\mathit{du}}{\mathit{dt}}-\Delta \mathit{u}\frac{\mathit{dx}}{\mathit{dt}}=ϵ{\Delta }^2\mathit{u}-\Delta \left(\frac{{\mathit{u}}^2}{2}\right).$$

이 형식에서 ODE 솔버(예: ode15s)를 사용하여 시간 경과에 따른 $\mathit{u}$와 $\mathit{x}$의 해를 구할 수 있습니다.

이 예에서 문제는 $\mathit{N}$개의 점으로 구성된 이동하는 메시 상에서 정식화되며 [1]에 설명된 이동 메시 기법은 각 시간 스텝에서 메시 점이 변경 영역에 집중되도록 배치합니다. 경계 조건과 초기 조건은 다음과 같습니다.

$$\begin{array}{l}\mathit{u}\left(0,\mathit{t}\right)=\mathit{u}\left(1,\mathit{t}\right)=0,\\ \mathit{u}\left(\mathit{x},0\right)=\sin \left(2\pi \mathit{x}\right)+\frac{1}{2}\sin \left(\pi \mathit{x}\right).\end{array}$$

주어진 $\mathit{N}$개 점으로 구성된 초기 메시의 경우, 풀어야 할 $2\mathit{N}$개의 방정식이 있습니다. 버거스 방정식에 해당하는 $\mathit{N}$개 방정식과 각 메시 점의 움직임을 결정하는 $\mathit{N}$개 방정식입니다. 따라서 최종 연립방정식은 다음과 같습니다.

$$\begin{array}{l}\frac{\mathit{d}{\mathit{u}}_1}{\mathit{dt}}-\Delta {\mathit{u}}_1\frac{\mathit{d}{\mathit{x}}_1}{\mathit{dt}}=ϵ{\Delta }^2{\mathit{u}}_1-\Delta \left(\frac{{{\mathit{u}}_1}^2}{2}\right),\\ ⋮\\ \frac{\mathit{d}{\mathit{u}}_{\mathit{N}}}{\mathit{dt}}-\Delta {\mathit{u}}_{\mathit{N}}\frac{\mathit{d}{\mathit{x}}_{\mathit{N}}}{\mathit{dt}}=ϵ{\Delta }^2{\mathit{u}}_{\mathit{N}}-\Delta \left(\frac{{{\mathit{u}}_{\mathit{N}}}^2}{2}\right),\\ \frac{{\mathit{d}}^2\dot{{\mathit{x}}_1}}{\mathit{d}{\mathit{t}}^2}=\frac{1}{\tau }\frac{\mathit{d}}{\mathit{dt}}\left(\mathit{B}\left({\mathit{x}}_1,\mathit{t}\right)\frac{\mathit{d}{\mathit{x}}_1}{\mathit{dt}}\right),\\ ⋮\\ \frac{{\mathit{d}}^2\dot{{\mathit{x}}_{\mathit{N}}}}{\mathit{d}{\mathit{t}}^2}=\frac{1}{\tau }\frac{\mathit{d}}{\mathit{dt}}\left(\mathit{B}\left({\mathit{x}}_{\mathit{N}},\mathit{t}\right)\frac{\mathit{d}{\mathit{x}}_{\mathit{N}}}{\mathit{dt}}\right).\end{array}$$

이동 메시에 대한 항은 [1]의 MMPDE6에 해당합니다. 파라미터 $\tau$는 메시가 등분포가 되게 강제하는 시간 스케일을 나타냅니다. $\mathit{B}\left(\mathit{x},\mathit{t}\right)$ 항은 [1]의 방정식 21으로 주어진 모니터 함수입니다. [1]에 있는 21:

$$\mathit{B}\left(\mathit{x},\mathit{t}\right)=\sqrt{1+{\left(\frac{\mathit{d}{\mathit{u}}_{\mathit{i}}}{\mathit{d}{\mathit{x}}_{\mathit{i}}}\right)}^2}.$$

이동하는 메시 점을 사용하여 버거스 방정식을 풀기 위해 이 예제에서 사용된 접근 방식은 몇 가지 기법을 보여줍니다.

• 연립방정식은 질량 행렬 정식화 $\mathit{M}{\mathit{y}}^{\prime }=\mathit{f}\left(\mathit{t},\mathit{y}\right)$를 사용하여 표현됩니다. 질량 행렬은 ode15s 솔버에 함수로 제공됩니다.

• 도함수는 버거스 방정식에 대한 수식뿐만 아니라 이동하는 메시의 선택을 제어하는 일련의 수식도 포함합니다.

• 야코비 행렬 $\mathrm{dF}/\mathrm{dy}$의 희소성 패턴과 벡터 $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$를 곱한 질량 행렬의 도함수는 솔버에 함수로 제공됩니다. 이러한 희소성 패턴을 제공하면 솔버가 보다 효율적으로 연산할 수 있습니다.

• 유한 차분은 여러 편도함수를 근사하는 데 사용됩니다.

MATLAB®에서 이 방정식을 풀려면 도함수, 질량 행렬 함수, 야코비 행렬 $\mathrm{dF}/\mathrm{dy}$의 희소성 패턴 함수 및 $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$의 희소성 패턴 함수를 작성하십시오. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.

질량 행렬 코딩하기

연립방정식의 좌변은 1계 도함수의 선형 결합을 포함하므로 모든 항을 표현하려면 질량 행렬이 필요합니다. 연립방정식의 좌변을 $\mathit{M}{\mathit{y}}^{\prime }$과 같도록 설정하여 질량 행렬의 형식을 추출합니다. 질량 행렬은 4개의 블록으로 구성되며 각각은 차수가 $\mathit{N}$인 정사각 행렬입니다.

$\left[\begin{array}{c}\frac{\partial {\mathit{u}}_1}{\partial \mathit{t}}-\frac{\partial {\mathit{u}}_1}{\partial {\mathit{x}}_1}\frac{\partial {\mathit{x}}_1}{\partial \mathit{t}}\\ ⋮\\ \frac{\partial {\mathit{u}}_{\mathit{N}}}{\partial \mathit{t}}-\frac{\partial {\mathit{u}}_{\mathit{N}}}{\partial {\mathit{x}}_{\mathit{N}}}\frac{\partial {\mathit{x}}_{\mathit{N}}}{\partial \mathit{t}}\\ \frac{{\partial }^2\dot{{\mathit{x}}_1}}{\partial {\mathit{t}}^2}\\ ⋮\\ \frac{{\partial }^2\dot{{\mathit{x}}_{\mathit{N}}}}{\partial {\mathit{t}}^2}\end{array}\right]=\mathit{M}{\mathit{y}}^{\prime }=\left[\begin{array}{cc}{\mathit{M}}_1& {\mathit{M}}_2\\ {\mathit{M}}_3& {\mathit{M}}_4\end{array}\right]\left[\begin{array}{c}\dot{{\mathit{u}}_1}\\ ⋮\\ \dot{{\mathit{u}}_{\mathit{N}}}\\ \dot{{\mathit{x}}_1}\\ ⋮\\ \dot{{\mathit{x}}_{\mathit{N}}}\end{array}\right]$.

이 정식화는 ${\mathit{M}}_1$ 및 ${\mathit{M}}_2$가 버거스 방정식의 좌변(연립방정식의 첫 번째 $\mathit{N}$개 방정식)을 구성하고 ${\mathit{M}}_3$ 및 ${\mathit{M}}_4$는 메시 방정식의 좌변(연립방정식의 마지막 $\mathit{N}$개 방정식)을 구성하는 것을 보여줍니다. 블록 행렬은 다음과 같습니다.

$$\begin{array}{l}{\mathit{M}}_1={\mathit{I}}_{\mathit{N}},\\ {\mathit{M}}_2=-\frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial {\mathit{x}}_{\mathit{i}}}{\mathit{I}}_{\mathit{N}},\\ {\mathit{M}}_3=0_{\mathit{N}},\\ {\mathit{M}}_4=\frac{{\partial }^2}{\partial {\mathit{t}}^2}{\mathit{I}}_{\mathit{N}}.\end{array}$$

${\mathit{I}}_{\mathit{N}}$은 $\mathit{N}\times \mathit{N}$ 단위 행렬입니다. ${\mathit{M}}_2$의 편도함수는 유한 차분을 사용하여 계산되고 ${\mathit{M}}_4$의 편도함수는 라플라시안 행렬을 사용합니다. 메시 움직임에 대한 방정식은 $\dot{\mathit{u}}$의 영향을 받지 않기 때문에 ${\mathit{M}}_3$은 0만 포함한다는 점을 주목하십시오.

이제 질량 행렬을 계산하는 함수를 작성할 수 있습니다. 함수는 시간 $\mathit{t}$ 및 해 벡터 $\mathit{y}$에 대한 두 개의 입력을 받아야 합니다. 해 벡터 $\mathit{y}$는 절반의 $\dot{\mathit{u}}$ 성분과 절반의 $\dot{\mathit{x}}$ 성분을 포함하고 있기 때문에 함수는 이 값들을 먼저 추출합니다. 그런 다음 함수는 문제의 경계값을 고려하여 모든 블록 행렬을 생성하고 4개의 블록을 사용하여 질량 행렬을 조합합니다.

function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.

도함수 코딩하기

이 문제의 도함수는 $2\mathit{N}$개의 요소를 가진 벡터를 반환합니다. 첫 번째 $\mathit{N}$개 요소는 버거스 방정식에 해당하고 마지막 $\mathit{N}$개 요소는 이동 메시 방정식에 대한 것입니다. 함수 movingMeshODE는 다음 단계를 거쳐 연립방정식에서 모든 방정식의 우변을 계산합니다.

1. 유한 차분을 사용하여 버거스 방정식을 계산합니다(첫 번째 $\mathit{N}$개 요소).

2. 모니터 함수를 계산합니다(마지막 $\mathit{N}$개 요소).

3. 모니터 함수에 공간 평활화를 적용하고 이동 메시 방정식을 계산합니다.

도함수의 첫 번째 $\mathit{N}$개 방정식은 버거스 방정식의 우변을 인코딩합니다. 버거스 방정식은 다음 형식의 공간 도함수를 포함하는 미분 연산자로 간주될 수 있습니다.

$\mathit{f}\left(\mathit{u}\right)=ϵ\frac{{\partial }^2\mathit{u}}{{\partial \mathit{x}}^2}-\frac{\partial }{\partial \mathit{x}}\left(\frac{{\mathit{u}}^2}{2}\right)$.

참고 문헌 [1]에서는 다음의 중심 유한 차분을 사용하여 미분 연산자 $\mathit{f}$를 근사하는 과정을 설명합니다.

$${\mathit{f}}_{\mathit{i}}=ϵ\left[\frac{\left(\frac{{\mathit{u}}_{\mathit{i}+1}-{\mathit{u}}_{\mathit{i}}}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}}}\right)-\left(\frac{{\mathit{u}}_{\mathit{i}}-{\mathit{u}}_{\mathit{i}-1}}{{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{i}-1}}\right)}{\frac{1}{2}\left({\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}\right)}\right]-\frac{1}{2}\left(\frac{{\mathit{u}}_{\mathit{i}+1}^2-{\mathit{u}}_{\mathit{i}-1}^2}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}}\right).$$

대신 메시의 모서리($\mathit{i}=1$ 및 $\mathit{i}=\mathit{N}$)에서는 단방향 차분만 사용됩니다. 이 예제에서는 $ϵ=1\times {10}^{-4}$를 사용합니다.

메시를 제어하는 방정식(도함수의 마지막 $\mathit{N}$개 방정식을 구성하는)은 다음과 같습니다.

$$\frac{{\partial }^2\dot{\mathit{x}}}{\partial {\mathit{t}}^2}=\frac{1}{\tau }\frac{\partial }{\partial \mathit{t}}\left(\mathit{B}\left(\mathit{x},\mathit{t}\right)\frac{\partial \mathit{x}}{\partial \mathit{t}}\right).$$

버거스 방정식과 마찬가지로 유한 차분을 사용하여 모니터 함수 $\mathit{B}\left(\mathit{x},\mathit{t}\right)$를 근사할 수 있습니다.

$$\mathit{B}\left(\mathit{x},\mathit{t}\right)=\sqrt{1+{\left(\frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial {\mathit{x}}_{\mathit{i}}}\right)}^2}=\sqrt{1+{\left(\frac{{\mathit{u}}_{\mathit{i}+1}-{\mathit{u}}_{\mathit{i}-1}}{{\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}-1}}\right)}^2}.$$

모니터 함수가 계산되면 공간 평활화가 적용됩니다([1]의 방정식 14 및 15). 이 예제에서는 공간 평활화 파라미터에 $\gamma =2$ 및 $\mathit{p}=2$를 사용합니다.

연립방정식을 인코딩하는 함수는 다음과 같습니다.

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

희소성 패턴에 대한 함수 코딩하기

도함수에 대한 야코비 행렬 $\mathrm{dF}/\mathrm{dy}$는 도함수 movingMeshODE의 모든 편도함수가 포함된 $2\mathit{N}\times 2\mathit{N}$ 행렬입니다. ode15s는 options 구조체에 행렬이 제공되지 않는 경우 유한 차분을 사용하여 야코비 행렬을 추정합니다. 야코비 행렬의 희소성 패턴을 제공하여 ode15s가 좀 더 빨리 계산하게 할 수 있습니다.

야코비 행렬의 희소성 패턴에 대한 함수는 다음과 같습니다.

function out = JPat(N)
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
       S2 S2];
end

spy를 사용하여 $\mathit{N}=80$에 대한 $\mathrm{dF}/\mathrm{dy}$의 희소성 패턴을 플로팅합니다.

spy(JPat(80))

좀 더 효율적으로 계산하기 위한 또 다른 방법은 $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$의 희소성 패턴을 제공하는 것입니다. ${\mathit{u}}_{\mathit{i}}$ 및 ${\mathit{x}}_{\mathit{i}}$ 중 어떤 항이 질량 행렬 함수에서 계산된 유한 차분에 있는지 검토하여 이 희소성 패턴을 찾을 수 있습니다.

$\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$의 희소성 패턴에 대한 함수는 다음과 같습니다.

function S = MvPat(N)
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end

spy를 사용하여 $\mathit{N}=80$에 대한 $\mathit{d}\left(\mathrm{Mv}\right)/\mathrm{dy}$의 희소성 패턴을 플로팅합니다.

spy(MvPat(80))

연립방정식 풀기

$\mathit{N}=80$ 값을 가진 연립방정식을 풉니다. 초기 조건의 경우 균일 그리드를 사용하여 $\mathit{x}$를 초기화하고 그리드에서 $\mathit{u}\left(\mathit{x},0\right)$을 계산합니다.

N = 80;
h = 1/(N+1);
xinit = h*(1:N);
uinit = sin(2*pi*xinit) + 0.5*sin(pi*xinit);
y0 = [uinit xinit];

odeset을 사용하여 다음 여러 값을 설정하는 options 구조체를 만듭니다.

• 질량 행렬에 대한 함수 핸들

• 질량 행렬의 상태 종속성. 질량 행렬은 $\mathit{t}$ 및 $\mathit{y}$ 모두의 함수이므로 이 문제의 경우 이 값은 'strong'입니다.

• 야코비 행렬의 희소성 패턴을 계산하는 함수 핸들

• 벡터를 곱한 질량 행렬의 도함수의 희소성 패턴을 계산하는 함수 핸들

• 절대 허용오차 및 상대 허용오차

opts = odeset('Mass',@mass,'MStateDependence','strong','JPattern',JPat(N),...
   'MvPattern',MvPat(N),'RelTol',1e-5,'AbsTol',1e-4);

마지막으로, ode15s를 호출하여 movingMeshODE 도함수, 시간 길이, 초기 조건 및 options 구조체를 사용하여 구간 $\left[0,1\right]$에서 연립방정식을 풉니다.

tspan = [0 1];
sol = ode15s(@movingMeshODE,tspan,y0,opts);

결과 플로팅하기

적분의 결과는 시간 스텝 $\mathit{t}$, 메시 점 $\mathit{x}\left(\mathit{t}\right)$, 해 $\mathit{u}\left(\mathit{x},\mathit{t}\right)$를 포함하는 구조체 sol입니다. 구조체에서 이 값들을 추출합니다.

t = sol.x;
x = sol.y(N+1:end,:);
u = sol.y(1:N,:);

시간 경과에 따른 메시 점의 움직임을 플로팅합니다. 플롯을 보면 (모니터 함수로 인해) 시간 경과에 따라 메시 점이 꽤 일정한 간격을 유지하지만 해가 이동함에 따라 해의 불연속 부분 근처에 군집화할 수 있음을 확인할 수 있습니다.

plot(x,t)
xlabel('t')
ylabel('x(t)')
title('Burgers'' equation: Trajectories of grid points')

이제 일부 $\mathit{t}$ 값에서 $\mathit{u}\left(\mathit{x},\mathit{t}\right)$를 샘플링하고 시간 경과에 따른 해의 변화를 플로팅합니다. 구간의 끝에서 메시 점이 고정되므로 x(0) = 0이고 x(N+1) = 1입니다. 경계값은 u(t,0) = 0과 u(t,1) = 0이며 Figure에 대해 계산된 알려진 값에 이 경계값을 추가해야 합니다.

tint = 0:0.2:1;
yint = deval(sol,tint);
figure
labels = {};
for j = 1:length(tint)
   solution = [0; yint(1:N,j); 0];
   location = [0; yint(N+1:end,j); 1];
   labels{j} = ['t = ' num2str(tint(j))];
   plot(location,solution,'-o')
   hold on
end
xlabel('x')
ylabel('solution u(x,t)')
legend(labels{:},'Location','SouthWest')
title('Burgers equation on moving mesh')
hold off

플롯을 보면 $\mathit{u}\left(\mathit{x},0\right)$은 $\mathit{x}=1$ 쪽으로 이동하는 동안 시간 경과에 따라 급격한 기울기가 발생하는 매끄러운 파동임을 알 수 있습니다. 메시 점은 추가 계산 점이 각 시간 스텝에서 적절한 위치에 있도록 불연속 부분의 움직임을 추적합니다.

참고 문헌

[1] Huang, Weizhang, et al. “Moving Mesh Methods Based on Moving Mesh Partial Differential Equations.” Journal of Computational Physics, vol. 113, no. 2, Aug. 1994, pp. 279–90. https://doi.org/10.1006/jcph.1994.1135.

로컬 함수(Local Function)

여기 나열된 함수는 솔버 ode15s가 해를 계산하기 위해 호출하는 로컬 헬퍼 함수입니다. 또는 이러한 함수를 MATLAB 경로에 있는 디렉터리에 고유의 파일로 저장할 수도 있습니다.

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

% -----------------------------------------------------------------------
function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

% -------------------------------------------------------------------------
function out = JPat(N)  % Jacobian sparsity pattern
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
    S2 S2];
end

% -------------------------------------------------------------------------
function S = MvPat(N)  % Sparsity pattern for the derivative of the Mass matrix times a vector
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end
% -------------------------------------------------------------------------

참고 항목

ode15s | odeset

관련 항목

• ODE 솔버 선택하기
• ODE 옵션 요약
• 공중으로 던져진 바통의 운동 방정식 풀기