불연속을 갖는 PDE 풀기

이 예제에서는 물질의 경계면이 있는 PDE를 푸는 방법을 보여줍니다. 물질의 경계면으로 인해 $x = 0.5$ 에서 문제에 불연속이 발생하고 초기 조건은 오른쪽 경계 $x = 1$ 에서 불연속을 갖습니다.

다음과 같은 조각별 PDE를 보겠습니다.

${\begin{matrix} \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} 5 \frac{\partial u}{\partial x}) - 1000 e^{u} & (0 \leq x \leq 0.5) \\ \frac{\partial u}{\partial t} = x^{- 2} \frac{\partial}{\partial x} (x^{2} \frac{\partial u}{\partial x}) - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

초기 조건은 다음과 같습니다.

$\begin{array}{l} u (x, 0) = 0 (0 \leq x < 1), \\ u (1, 0) = 1 (x = 1) . \end{array}$

경계 조건은 다음과 같습니다.

$\begin{array}{l} \frac{\partial u}{\partial x} = 0 (x = 0), \\ u (1, t) = 1 (x = 1) . \end{array}$

MATLAB®에서 이 방정식을 풀려면 방정식, 초기 조건 및 경계 조건을 코딩하고 적합한 해 메시를 선택한 후에 솔버 pdepe를 호출해야 합니다. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.

방정식 코딩하기

방정식을 코딩하려면 먼저 방정식이 pdepe 솔버에서 요구하는 형식으로 되어 있는지 확인해야 합니다. pdepe에서 요구하는 표준 형식은 다음과 같습니다.

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x})$ .

이 경우, PDE가 적절한 형식이므로 계수 값을 파악해 낼 수 있습니다.

플럭스 항 $f (x, t, u, \frac{\partial u}{\partial x})$ 와 소스 항 $s (x . t, u, \frac{\partial u}{\partial x})$ 의 값은 $x$ 값에 따라 변합니다. 계수는 다음과 같습니다.

$m = 2$

$c (x, t, u, \frac{\partial u}{\partial x}) = 1$

${\begin{matrix} f (x, t, u, \frac{\partial u}{\partial x}) = 5 \frac{\partial u}{\partial x} & (0 \leq x \leq 0.5) \\ f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x} & (0.5 \leq x \leq 1) \end{matrix}$

${\begin{matrix} s (x, t, u, \frac{\partial u}{\partial x}) = - 1000 e^{u} & (0 \leq x \leq 0.5) \\ s (x, t, u, \frac{\partial u}{\partial x}) = - e^{u} & (0.5 \leq x \leq 1) \end{matrix}$

이제 방정식을 코딩하는 함수를 만들 수 있습니다. 이 함수에는 시그니처 [c,f,s] = pdex2pde(x,t,u,dudx)가 있어야 합니다.

x는 독립적 공간 변수입니다.
t는 독립적 시간 변수입니다.
u는 x와 t에 따라 미분되는 종속 변수입니다.
dudx는 편미분 공간 도함수 $\partial u / \partial x$ 입니다.
출력값 c, f 및 s는 pdepe에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다. 이 계수는 입력 변수 x, t, u 및 dudx를 사용해 코딩됩니다.

그 결과, 이 예제의 방정식은 다음과 같은 함수로 표현할 수 있습니다.

function [c,f,s] = pdex2pde(x,t,u,dudx)
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end

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

초기 조건 코딩하기

이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값에서 적용되고 모든 $x$ 값에 대한 $u (x, t_{0})$ 값을 제공합니다. 함수 시그니처 u0 = pdex2ic(x)를 사용하여 함수를 작성하십시오.

초기 조건은 다음과 같습니다.

$\begin{array}{l} u (x, 0) = 0 (0 \leq x < 1), \\ u (1, 0) = 1 (x = 1) . \end{array}$

대응하는 함수는 다음과 같습니다.

function u0 = pdex2ic(x)
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end

경계 조건 코딩하기

이제 경계 조건을 계산하는 함수를 작성합니다. 구간 $a \leq x \leq b$ 에 있는 문제의 경우, 경계 조건은 모든 , $x = a$ 또는 $x = b$ 에 적용됩니다. 솔버에서 요구하는 경계 조건의 표준 형식은 다음과 같습니다.

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

이 예제에는 구 대칭( $m = 2$ )이 있기 때문에, pdepe 솔버는 왼쪽 경계 조건을 자동으로 적용하여 원점에서 해의 경계를 지정하고 경계 함수의 왼쪽 경계에 지정된 모든 조건을 무시합니다. 그러므로 왼쪽 경계 조건에 $p_{L} = q_{L} = 0$ 을 지정하면 됩니다. 오른쪽 경계 조건에 대해서는 표준 형식으로 경계 조건을 다시 작성하고 $p_{R}$ 과 $q_{R}$ 에 대한 계수 값을 파악해 낼 수 있습니다.

$x = 1$ 인 경우 방정식은 $u (1, t) = 1 \to (u - 1) + 0 \cdot \frac{\partial u}{\partial x} = 0$ 입니다. 계수는 다음과 같습니다.

$p_{R} (1, t, u) = u - 1$
$q_{R} (1, t) = 0$

경계 함수는 함수 시그니처 [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)를 사용해야 합니다.

입력값 xl과 ul은 왼쪽 경계의 u와 x에 대응합니다.
입력값 xr과 ur은 오른쪽 경계의 u와 x에 대응합니다.
t는 독립적 시간 변수입니다.
출력값 pl과 ql은 왼쪽 경계(이 문제의 경우 $x = 0$ )의 $p_{L} (x, t, u)$ 와 $q_{L} (x, t)$ 에 대응합니다.
출력값 pr과 qr은 오른쪽 경계(이 문제의 경우 $x = 1$ )의 $p_{R} (x, t, u)$ 와 $q_{R} (x, t)$ 에 대응합니다.

이 예제의 경계 조건은 다음과 같은 함수로 표현됩니다.

function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = 0; 
ql = 0; 
pr = ur - 1;
qr = 0;
end

해 메시 선택하기

공간 메시는 불연속 접점을 반영하기 위해 $x = 0.5$ 근처에 있는 여러 개의 값을 포함해야 하며, 일관되지 않은 초기값( $u (1, 0) = 1$ ) 및 경계값( $u (1, t) = 0$ )으로 인해 $x = 1$ 근처에 있는 점들 역시 포함해야 합니다. 작은 $t$ 에 대해 해가 급격하게 변하므로 이러한 급격한 변화에 대응할 수 있는 시간 스텝을 사용하십시오.

x = [0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 0.95 0.975 0.99 1];
t = [0 0.001 0.005 0.01 0.05 0.1 0.5 1];

방정식 풀기

마지막으로 대칭 m, PDE 방정식, 초기 조건, 경계 조건, x와 t에 대한 메시를 사용하여 방정식을 풉니다.

m = 2;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

pdepe는 3차원 배열 sol로 해를 반환합니다. 여기서 sol(i,j,k)는 t(i)와 x(j)에서 구한 해의 k번째 성분 $u_{k}$ 에 대한 근삿값을 계산합니다. u0이 각 해 성분에 초기 조건을 지정하므로, sol의 크기는 length(t)×length(x)×length(u0)입니다. 이 문제의 경우 u에 성분이 하나뿐이므로 sol은 8×18 행렬이지만, 일반적으로 명령 u = sol(:,:,k)를 사용하여 k번째 해 성분을 추출할 수 있습니다.

sol에서 첫 번째 해 성분을 추출합니다.

u = sol(:,:,1);

해 플로팅하기

$x$ 와 $t$ 에 대해 선택한 메시 점에서 플로팅된 해 $u$ 의 곡면 플롯을 만듭니다. $m = 2$ 이기 때문에 이 문제는 구 대칭이 있는 구체 기하학적 구조에 대해 정의되며, 이에 따라 해는 방사형 $x$ 방향으로만 변합니다.

surf(x,t,u)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution u')

Figure contains an axes object. The axes object with title Numerical solution with nonuniform mesh, xlabel Distance x, ylabel Time t contains an object of type surface.

이제 곡면 플롯의 등고선을 측면에서 볼 수 있도록 $x$ 와 $u$ 만 플로팅합니다. 물질의 경계면 효과를 강조 표시하려면 $x = 0.5$ 에서 선을 추가하십시오.

plot(x,u,x,u,'*')
line([0.5 0.5], [-3 1], 'Color', 'k')
xlabel('Distance x')
ylabel('Solution u')
title('Solution profiles at several times')

Figure contains an axes object. The axes object with title Solution profiles at several times, xlabel Distance x, ylabel Solution u contains 17 objects of type line. One or more of the lines displays its values using only markers

로컬 함수(Local Function)

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

function [c,f,s] = pdex2pde(x,t,u,dudx) % Equation to solve
c = 1;
if x <= 0.5
    f = 5*dudx;
    s = -1000*exp(u);
else
    f = dudx;
    s = -exp(u);
end
end
%----------------------------------------------
function u0 = pdex2ic(x) %Initial conditions
if x < 1
    u0 = 0;
else
    u0 = 1;
end
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) % Boundary conditions
pl = 0;
ql = 0;
pr = ur - 1;
qr = 0;
end
%----------------------------------------------

참고 항목

pdepe