pdepe

1차원 포물선 및 타원 PDE 풀기

구문

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

설명

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)은 하나의 공간 변수 x와 시간 t에서 포물선 및 타원 PDE로 구성된 연립방정식을 풉니다. 최소 하나 이상의 방정식이 포물선 방정식이어야 합니다. 스칼라 m은 문제의 대칭성을 나타냅니다(슬래브, 원통 또는 구면). 풀려는 방정식은 pdefun에 코딩되고 초기값은 icfun에 코딩되며 경계 조건은 bcfun에 코딩됩니다. 공간을 이산화할 때 생성되는 상미분 방정식(ODE)은 적분되어, tspan에 지정한 시간에서 근사해를 얻게 됩니다. pdepe 함수는 xmesh에서 제공한 메시에 대한 해의 값을 반환합니다.

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)는 odeset 함수를 사용하여 생성된 options로 정의된 적분 설정도 사용합니다. 예를 들어, 절대 허용오차와 상대 허용오차를 지정하려면 AbsTol 옵션과 RelTol 옵션을 사용하십시오.

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)는 이벤트 함수라고 하는 (t,u(x,t))의 함수가 0인 위치도 찾습니다. 출력값에서 te는 이벤트 발생 시간이고, sole는 이벤트 발생 시 계산된 해이며, ie는 트리거된 이벤트의 인덱스입니다. tsol은 tspan에 지정한 시간 중에서 첫 번째 종료 이벤트 이전에 해당하는 시간의 열 벡터입니다.

각 이벤트 함수에 대해, 0에서 적분을 종료할지 여부와 영점교차의 방향을 고려할지 여부를 지정합니다. 이를 수행하려면 odeset의 'Events' 옵션을 함수(예: @myEventFcn)로 설정하고 대응 함수 [value,isterminal,direction] = myEventFcn(m,t,xmesh,umesh)를 생성합니다. xmesh 입력값은 공간 메시를 포함하고 umesh는 메시 점에서 구한 해입니다.

예제

원통 좌표에서 열 방정식 풀기

pdepe를 사용하여 원통 좌표에서 열 방정식을 풀고 그 해를 플로팅합니다.

각도 대칭성이 있는 원통 좌표에서 열 방정식은 다음과 같습니다.

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{1}{\mathit{x}}\frac{\partial }{\partial \mathit{x}}\left(\mathit{x}\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$$

이 방정식은 시간 $\mathit{t}\ge 0$에서 $0\le \mathit{x}\le 1$에 대해 정의됩니다. 초기 조건은 베셀 함수 ${\mathit{J}}_0\left(\mathit{x}\right)$와 그 첫 번째 영점 $\mathit{n}=2.404825557695773$으로 다음과 같이 정의됩니다.

$$\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_0\left(\mathrm{nx}\right).$$

이 문제는 원통 좌표(m = 1)에 있으므로, pdepe는 $\mathit{x}=0$에서의 대칭 조건을 자동으로 적용합니다. 오른쪽 경계 조건은 다음과 같습니다.

$$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$$

초기 조건과 경계 조건은 다음 문제에 대한 해석적 해와 일치하도록 선택됩니다.

$$\mathit{u}\left(\mathit{x},\mathit{t}\right)={\mathit{J}}_0\left(\mathrm{nx}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$$

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

방정식 코딩하기

방정식을 코딩하려면 먼저 pdepe 솔버에서 요구하는 형식으로 방정식을 다시 작성해야 합니다. pdepe에서 요구하는 표준 형식은 다음과 같습니다.

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

이 표준 형식으로 작성된 PDE는 다음과 같습니다.

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

방정식을 적절한 형식으로 작성하면 관련 항을 다음과 같이 파악할 수 있습니다.

• $\mathit{m}=1$

• $\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=1$

• $\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{\partial \mathit{u}}{\partial \mathit{x}}$

• $\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0$

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

• x는 독립적 공간 변수입니다.

• t는 독립적 시간 변수입니다.

• u는 x와 t에 따라 미분되는 종속 변수입니다.

• dudx는 편미분 공간 도함수 $\partial \mathit{u}/\partial \mathit{x}$입니다.

• 출력값 c, f 및 s는 pdepe에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다. 이 계수는 입력 변수 x, t, u 및 dudx를 사용해 코딩됩니다.

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

function [c,f,s] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

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

초기 조건 코딩하기

이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값 tspan(1)에서 적용됩니다. 이 함수에는 시그니처 u0 = heatic(x)가 있어야 합니다.

$\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_0\left(\mathrm{nx}\right)$에 대응하는 함수는 다음과 같습니다.

function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end

경계 조건 코딩하기

이제 경계 조건을 계산하는 함수를 작성합니다.

$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$

이 문제는 원통 좌표(m = 1)에 있으므로, pdepe는 $\mathit{x}=0$에서의 대칭 조건을 자동으로 적용하며, 따라서 왼쪽 경계 조건을 지정할 필요가 없습니다.

구간 $\mathit{a}\le \mathit{x}\le \mathit{b}$에 있는 문제의 경우, 경계 조건은 모든 $\mathit{t}$, 그리고 $\mathit{x}=\mathit{a}$ 또는 $\mathit{x}=\mathit{b}$에 적용됩니다. 솔버에서 요구하는 경계 조건의 표준 형식은 다음과 같습니다.

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

이 형식으로 작성된 경우 $\mathit{u}$의 편도함수에 대한 경계 조건은 플럭스 $\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$에 대해 표현해야 합니다. 따라서 이 문제에 대한 오른쪽 경계 조건은 다음과 같습니다.

$$\mathit{u}\left(1,\mathit{t}\right)-{\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

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

• 입력값 xl과 ul은 왼쪽 경계의 $\mathit{u}$와 $\mathit{x}$에 대응합니다.

• 입력값 xr과 ur은 오른쪽 경계의 $\mathit{u}$와 $\mathit{x}$에 대응합니다.

• t는 독립적 시간 변수입니다.

• 출력값 pl과 ql은 왼쪽 경계(이 문제의 경우 $\mathit{x}=0$)의 ${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$와 ${\mathit{q}}_{\mathit{L}}\left(\mathit{x},\mathit{t}\right)$에 대응합니다.

• 출력값 pr과 qr은 오른쪽 경계(이 문제의 경우 $\mathit{x}=1$)의 ${\mathit{p}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$와 ${\mathit{q}}_{\mathit{R}}\left(\mathit{x},\mathit{t}\right)$에 대응합니다.

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

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end

해 메시 선택하기

방정식을 풀기 전에 pdepe가 해를 계산할 메시 점 $\left(\mathit{t},\mathit{x}\right)$를 지정해야 합니다. 벡터 t와 x로 점을 지정합니다. 벡터 t와 x는 솔버에서 각각 다른 역할을 수행합니다. 특히, 해에 대한 비용과 정확도는 벡터 x의 길이에 크게 종속됩니다. 하지만, 이 계산은 벡터 t의 값에 훨씬 덜 민감합니다.

이 문제에서는 x와 t 모두에 대해 공간 구간 [0,1]에서 균일한 간격의 점 25개를 포함한 메시를 사용합니다.

x = linspace(0,1,25);
t = linspace(0,1,25);

방정식 풀기

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

m = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

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

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

u = sol(:,:,1);

해 플로팅하기

해의 곡면 플롯을 만듭니다. 디스크 상의 원통 좌표로 제기되는 문제이므로, x 값은 중심에서 일정한 거리에 있는 디스크의 온도를 표시하고 t는 특정 위치에서 온도가 시간 경과에 따라 어떻게 변하는지 표시합니다.

surf(x,t,u)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])

디스크 중심($\mathit{x}=0$)의 온도 변화를 플로팅합니다.

plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of disc')

로컬 함수(Local Function)

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

function [c,f,s] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end
%----------------------------------------------

이벤트 기록으로 진동 PDE 풀기

편미분 방정식을 풀고 이벤트 함수를 사용하여 진동 해에 영점교차를 기록합니다.

다음과 같은 방정식이 있다고 가정해 보겠습니다.

$$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial }{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$$

이 방정식은 $0\le \mathit{x}\le 1$ 및 $\mathit{t}\ge \mathrm{0}.1$에 대해 정의됩니다. 초기 조건은 다음과 같습니다.

$$\mathit{u}\left(\mathit{x},0.1\right)=1.$$

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

$$\mathit{u}\left(0,\mathit{t}\right)=1,$$

$$\mathit{u}\left(1,\mathit{t}\right)=\cos \left(\pi \mathit{t}\right).$$

또한 해의 영점교차가 관심 지점입니다.

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

방정식 코딩하기

방정식을 코딩하려면 먼저 pdepe 솔버에서 요구하는 형식으로 방정식을 다시 작성해야 합니다. pdepe에서 요구하는 표준 형식은 다음과 같습니다.

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

PDE 방정식은 이미 다음과 같은 형식으로 되어 있습니다.

$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial }{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$

따라서 관련 항을 다음과 같이 파악할 수 있습니다.

• $$\mathit{m}=0$$

• $$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{1}{\mathit{x}}$$

• $$\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{1}{\mathit{t}}\mathit{u}$$

• $$\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0$$

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

• x는 독립적 공간 변수입니다.

• t는 독립적 시간 변수입니다.

• u는 x와 t에 따라 미분되는 종속 변수입니다.

• dudx는 편미분 공간 도함수 $\partial \mathit{u}/\partial \mathit{x}$입니다.

• 출력값 c, f 및 s는 pdepe에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다. 이 계수는 입력 변수 x, t, u 및 dudx를 사용해 코딩됩니다.

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

function [c,f,s] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end

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

초기 조건 코딩하기

이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값 tspan(1)에서 적용됩니다. 이 함수에는 시그니처 u0 = oscic(x)가 있어야 합니다.

$\mathit{u}\left(\mathit{x},0.1\right)=1$에 대응하는 함수는 다음과 같습니다.

function u0 = oscic(x)
u0 = 1;
end

경계 조건 코딩하기

이제 경계 조건을 계산하는 함수를 작성합니다.

$$\mathit{u}\left(0,\mathit{t}\right)=1,$$

$$\mathit{u}\left(1,\mathit{t}\right)=\cos \left(\pi \mathit{t}\right).$$

구간 $\mathit{a}\le \mathit{x}\le \mathit{b}$에 있는 문제의 경우, 경계 조건은 모든 $\mathit{t}$, 그리고 $\mathit{x}=\mathit{a}$ 또는 $\mathit{x}=\mathit{b}$에 적용됩니다. 솔버에서 요구하는 경계 조건의 표준 형식은 다음과 같습니다.

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

이 문제의 경계 조건을 표준 형식으로 작성하면 다음과 같이 됩니다.

$$\mathit{u}\left(0,\mathit{t}\right)-1+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0,$$

$$\mathit{u}\left(1,\mathit{t}\right)-\cos \left(\pi \mathit{t}\right)+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

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

• 입력값 xl과 ul은 왼쪽 경계의 $\mathit{x}$와 $\mathit{u}$에 대응합니다.

• 입력값 xr과 ur은 오른쪽 경계의 $\mathit{x}$와 $\mathit{u}$에 대응합니다.

• t는 독립적 시간 변수입니다.

• 출력값 pl과 ql은 왼쪽 경계(이 문제의 경우 $\mathit{x}=0$)의 ${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$와 ${\mathit{q}}_{\mathit{L}}\left(\mathit{x},\mathit{t}\right)$에 대응합니다.

• 출력값 pr과 qr은 오른쪽 경계(이 문제의 경우 $\mathit{x}=1$)의 ${\mathit{p}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$와 ${\mathit{q}}_{\mathit{R}}\left(\mathit{x},\mathit{t}\right)$에 대응합니다.

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

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

이벤트 함수 코딩하기

이벤트 함수를 사용하여 적분에서 해의 영점교차를 기록합니다. 이벤트 함수에는 [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)라는 함수 시그니처가 있습니다.

• m은 pdepe에 대한 첫 번째 입력값으로 지정된 좌표 대칭입니다.

• t는 현재 시간(스칼라)입니다.

• xmesh는 공간 메시입니다.

• umesh는 메시 점에서 구한 해를 포함합니다.

• value는 흔히 해 umesh로 표현되는 관심 방정식입니다. value가 0인 경우 이벤트가 발생합니다.

• isterminal은 이벤트가 적분을 중지하는지 여부를 지정합니다. isterminal이 0인 경우 이벤트가 기록되지만 적분은 중지되지 않습니다. isterminal이 1인 경우 이벤트가 발생하면 적분이 중지됩니다.

• direction은 영점교차의 방향을 지정합니다. 1이면 양의 기울기를 가진 영점교차만 이벤트를 트리거합니다. -1이면 영점교차가 음의 기울기를 가져야 합니다. 0이면 어떤 영점교차든 이벤트를 트리거합니다.

적분의 각 시점에서 솔버는 이벤트 함수를 호출하여 영점교차 여부를 확인합니다. 모든 영점교차를 기록하려면 value가 해 벡터 umesh에서 부호의 변화를 찾아야 합니다. 이 예제의 이벤트는 종료 이벤트가 아니고 어떤 기울기에서든 영점교차가 발생할 수 있으므로, isterminal과 direction을 0으로 구성된 동일한 크기의 벡터로 지정합니다.

이 문제의 이벤트 함수는 다음과 같습니다.

function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end

해 메시 선택하기

방정식을 풀기 전에 pdepe가 해를 계산할 메시 점 $\left(\mathit{t},\mathit{x}\right)$를 지정해야 합니다. 이 문제의 경우 구간 $0\le \mathit{x}\le 1$ 및 $0.1\le \mathit{t}\le 1$에서 50개의 점으로 구성된 세밀한 메시를 사용합니다. 세밀한 메시를 통해 우수한 분해능의 진동 해를 구할 수 있습니다.

x = linspace(0,1,50);
t = linspace(0.1,pi,50);

방정식 풀기

마지막으로 대칭 m, PDE 방정식, 초기 조건, 경계 조건, 이벤트 함수 및 x와 t에 대한 메시를 사용하여 방정식을 풉니다. odeset을 사용하여 이벤트 함수를 참조하는 options 구조체를 만들고 구조체의 형태로 pdepe에 대한 마지막 입력 인수로 전달합니다. 이벤트 함수와 솔버 모두에서 정보를 반환할 5개 출력 인수를 지정합니다.

• sol은 pdepe로 계산한 해입니다.

• tsol은 종료 이벤트 이전의 시간으로 구성된 벡터입니다. 종료 이벤트가 없는 경우 tsol은 t와 같습니다.

• sole는 각 이벤트 시간에서의 해입니다.

• te는 각 이벤트의 시간입니다.

• ie는 각 이벤트의 인덱스입니다. 이벤트 함수에서 values = umesh이므로, ie는 각 시간 스텝에서 이벤트를 트리거한 umesh의 인덱스를 제공합니다.

m = 0;
options = odeset('Events',@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);

해를 행렬 u로 추출합니다.

u = sol(:,:,1);

해 플로팅하기

해의 곡면 플롯을 만들고 이 플롯을 위에서 봅니다.

surf(x,t,u)
view(2)

곡면 $\mathit{f}\left(\mathit{x},\mathit{t}\right)=0$을 기준으로 이벤트가 발생한 점을 플로팅합니다. 출력 인덱스 벡터 ie는 이벤트 위치를 선택하는 데 유용합니다. 표현식 x(ie)'는 이벤트가 발생한 x값을 제공하고, 표현식 sole(x==x(ie)')는 대응되는 해 값을 제공합니다.

view([39 30])
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
hold on
plot3(x(ie)',te,sole(x==x(ie)'),'r*')
surf(x,t,zeros(size(u)),'EdgeColor','flat')
hold off

로컬 함수(Local Function)

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

function [c,f,s] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end
%----------------------------------------------
function u0 = oscic(x)
u0 = 1;
end
%----------------------------------------------
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end
%----------------------------------------------
function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
%----------------------------------------------

입력 인수

m — 대칭 상수

0 | 1 | 2

대칭 상수로, 아래 표에 있는 값 중 하나로 지정됩니다. m은 pdefun으로 지정되는 방정식에서 문제의 유형을 지정합니다. pdepe에 필요한 형식으로 방정식을 다시 작성하면 방정식에서 m의 값을 확인할 수 있습니다.

값	설명	플럭스의 라플라시안
0	대칭성이 없는 1차원 카테시안 좌표	$$\Delta f=\frac{{\partial }^{2}f}{\partial x^2}$$
1	방위각 대칭성이 있는 2차원 원통 좌표	$$\Delta f=\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{\partial f}{\partial \rho }\right)+\frac{1}{{\rho }^2}\frac{{\partial }^{2}f}{\partial {\phi }^2}$$, 여기서 $$\frac{\partial f}{\partial \phi }=0$$(방위각 대칭성)입니다.
2	방위각 대칭성과 천정각 대칭성이 있는 3차원 구면 좌표	$$\Delta f=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial f}{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^{2}f}{\partial {\phi }^2}$$, 여기서 $$\frac{\partial f}{\partial \theta }=\frac{\partial f}{\partial \phi }=0$$(방위각 대칭성과 천정각 대칭성)입니다.

m > 0이면 pdepe는 왼쪽 공간 경계 a ≥ 0이어야 하고, 왼쪽 경계 조건을 자동으로 적용하여 원점에서의 좌표 불연속을 고려합니다. 이 경우, pdepe는 bcfun에서 지정된 왼쪽 경계 조건을 모두 무시합니다.

pdefun — 방정식에 대한 PDE 함수

함수 핸들

방정식에 대한 PDE 함수로, PDE의 계수를 정의하는 함수 핸들로 지정됩니다.

pdefun은 다음 형식의 PDE에 대한 계수를 인코딩합니다.

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

이 방정식의 항은 다음과 같습니다.

• $$f\left(x,t,u,\frac{\partial u}{\partial x}\right)$$는 플럭스 항입니다.

• $$s\left(x,t,u,\frac{\partial u}{\partial x}\right)$$는 소스 항입니다.

• 시간에 대한 편도함수의 결합은 대각 행렬 $$c\left(x,t,u,\frac{\partial u}{\partial x}\right)$$의 곱셈으로 제한됩니다. 이 행렬의 대각선 요소는 0이거나 양수입니다. 0인 요소는 타원 방정식에 대응되고, 다른 모든 요소는 포물선 방정식에 대응됩니다. 하나 이상의 포물선 방정식이 있어야 합니다. 포물선 방정식에 대응되는 c의 요소는 x의 고립된 값이 메시 점인 경우 그러한 x의 값에서 0이 될 수 있습니다.

• 물질의 경계면으로 인한 c 및 s에서의 불연속은 메시 점이 각 경계면에 있는 경우 허용됩니다.

PDE는 t₀ ≤ t ≤ t_f 및 a ≤ x ≤ b에서 성립합니다. 값 tspan(1)과 tspan(end)는 t₀과 t_f에 대응하고, xmesh(1)과 xmesh(end)는 a와 b에 대응합니다. 구간 [a,b]는 유한해야 합니다. m은 0, 1, 2 중 하나여야 하며, 각각 슬래브(Slab) 대칭, 원통 대칭, 구 대칭에 해당합니다. m > 0인 경우 a ≥ 0이어야 합니다.

방정식 코딩하기

계수 c, f, s의 값을 계산하는 함수 pdefun을 작성합니다. 다음과 같이 함수 시그니처를 사용합니다.

[c,f,s] = pdefun(x,t,u,dudx)

pdefun에 대한 입력 인수는 스칼라 x 및 t와 벡터 u 및 dudx입니다. 벡터 u는 해 u의 근삿값을 계산하고 dudx는 x에 대한 편도함수의 근삿값을 계산합니다. 출력 인수 c, f, s는 방정식 개수와 같은 요소 개수를 가진 열 벡터입니다. c는 행렬 c의 대각선 요소를 저장합니다.

예

열 방정식 $$\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}$$는 다음과 같이 솔버에 필요한 형식으로 다시 작성할 수 있습니다.

$$1·\frac{\partial u}{\partial t}=x^0\frac{\partial }{\partial x}\left(x^0\frac{\partial u}{\partial x}\right).$$

이 형식에서 계수의 값이 m = 0, c = 1, f = ∂u/∂x, s = 0임을 알 수 있습니다. 이 방정식에 대한 함수는 다음과 같습니다.

function [c,f,s] = heatpde(x,t,u,dudx)  
  c = 1;
  f = dudx;
  s = 0;
end

데이터형: function_handle

icfun — 초기 조건 함수

함수 핸들

초기 조건 함수로, 초기 조건을 정의하는 함수 핸들로 지정됩니다.

t = t₀과 모든 x에 대해, 해 성분은 다음 형식의 초기 조건을 충족합니다.

$$u\left(x,t_0\right)=u_0\left(x\right).$$

초기 조건 코딩하기

초기 조건을 정의하는 함수 icfun을 작성합니다. 다음과 같이 함수 시그니처를 사용합니다.

function u0 = icfun(x)

인수 x를 사용하여 호출된 경우, 함수 icfun은 x에서의 해 성분의 초기값을 계산하여 열 벡터 u0에 반환합니다. u0 내 요소 개수는 방정식의 개수와 동일합니다.

예

상수 초기 조건 $$u\left(x,0\right)=1/2$$은 다음 함수에 대응됩니다.

function u0 = heatic(x)
  u0 = 0.5;
end

데이터형: function_handle

bcfun — 경계 조건 함수

함수 핸들

경계 조건 함수로, 경계 조건을 정의하는 함수 핸들로 지정됩니다.

모든 t와 x = a 또는 x = b에 대해, 해 성분은 다음 형식의 경계 조건을 충족합니다.

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

q의 요소는 0이거나 절대로 0이 아닙니다. 경계 조건은 ∂u/∂x에 대해서가 아니라 플럭스 f에 대해 표현됩니다. 또한, 두 개의 계수 중 p만 u에 종속될 수 있습니다.

경계 조건 코딩하기

경계 조건의 항 p와 q를 정의하는 함수 bcfun을 작성합니다. 다음과 같이 함수 시그니처를 사용합니다.

function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

• ul은 왼쪽 경계 xl = a에서의 근사해이고, ur은 오른쪽 경계 xr = b에서의 근사해입니다.

• pl과 ql은 xl에서 계산된 p와 q에 대응되는 열 벡터입니다. 마찬가지로, pr과 qr은 xr에 대응됩니다.

• 출력값 pl, ql, pr, qr의 요소 개수는 방정식의 개수와 동일합니다.

m > 0이고 a = 0일 경우, x = 0 근방에서의 해의 유계성(Boundedness)에 따라 플럭스 f는 a = 0에서 0이 되어야 합니다. pdepe는 이 경계 조건을 자동으로 적용하고 pl과 ql에 반환되는 값을 무시합니다.

예

구간 $$0\le x\le 1$$에서 경계 조건 $$u\left(0,t\right)=0,u\left(1,t\right)=1$$을 고려합니다. 솔버에 필요한 형식의 경계 조건은 다음과 같습니다.

$$\begin{array}{l}u\left(0,t\right)+0·f\left(0,t,u,\frac{\partial u}{\partial x}\right)=0,\\ u\left(1,t\right)-1+0·f\left(0,t,u,\frac{\partial u}{\partial x}\right)=0.\end{array}$$

이 형식을 보면 q 항이 모두 0임을 분명히 확인할 수 있습니다. p 항은 u의 값을 반영하기 위해 0이 아닙니다. 대응하는 함수는 다음과 같습니다.

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

데이터형: function_handle

xmesh — 공간 메시

벡터

공간 메시로, tspan의 모든 값에 대해 수치 해를 구하는 지점을 지정하는 벡터 [x0 x1 ... xn]으로 지정됩니다. 해에 대한 2계 근삿값은 xmesh에 지정된 메시를 바탕으로 계산됩니다. pdepe는 x에서 메시를 자동으로 선택하지 않습니다. 사용자가 적절한 고정 메시를 xmesh에 제공해야 합니다.

xmesh의 요소는 x0 < x1 < ... < xn을 충족해야 합니다. 일반적으로, 해가 급격히 변하는 곳에서는 좁은 간격의 메시 점을 사용하는 것이 좋습니다. xmesh의 길이는 >=3이어야 하고 해의 계산 비용은 xmesh의 길이에 크게 종속됩니다.

불연속을 포함한 문제의 경우, 문제가 각 하위 구간에서 매끄러워지도록 불연속 지점에 메시 점들을 배치해야 합니다. 하지만 m > 0일 경우, x = 0 근처에서는 좌표 특이점을 고려하기 위해 세밀한 메시를 사용할 필요가 없습니다.

예: xmesh = linspace(0,5,25)는 0과 5 사이에 있는 25개의 점으로 구성된 공간 메시를 사용합니다.

데이터형: single | double

tspan — 적분의 시간 길이

벡터

적분의 시간 길이로, xmesh의 모든 값에 대해 해를 구하는 지점을 지정하는 벡터 [t0 t1 ... tf]로 지정됩니다. pdepe 함수는 시간 스텝과 식을 모두 동적으로 선택하는 ODE 솔버를 사용하여 시간에 대한 적분을 수행합니다. tspan의 요소는 단순히 답을 구하려는 위치를 지정하므로, tspan의 길이가 해의 계산 비용에 미치는 영향은 미미합니다.

tspan의 요소는 t0 < t1 < ... < tf를 충족해야 합니다. tspan의 길이는 >=3이어야 합니다.

예: tspan = linspace(0,5,5)는 0과 5 사이에서 5개의 시간 지점을 사용합니다.

데이터형: single | double

options — options 구조체

구조체형 배열

options 구조체로, 구조체형 배열로 지정됩니다. odeset 함수를 사용하여 options 구조체를 만들거나 수정합니다. pdepe는 다음 옵션을 지원합니다.

• 'RelTol'

• 'AbsTol'

• 'NormControl'

• 'InitialStep'

• 'MaxStep'

• 'Events'

대부분의 경우, 이러한 옵션에 대해 디폴트 값을 사용하면 만족스러운 해를 얻을 수 있습니다.

예: options = odeset('RelTol',1e-5)는 1e-5의 상대 허용오차를 지정합니다.

데이터형: struct

출력 인수

sol — 해 배열

3차원 배열

해 배열로, 3차원 배열로 반환됩니다.

pdepe는 다차원 배열로 해를 반환합니다. u_i = ui = sol(:,:,i)는 해 벡터 u의 i번째 성분에 대한 근삿값입니다. 요소 ui(j,k) = sol(j,k,i)는 (t,x) = (tspan(j),xmesh(k))에서 u_i에 대한 근삿값을 구합니다.

ui = sol(j,:,i)는 시간 tspan(j)와 메시 점 xmesh(:)에서 해의 성분 i에 대한 근삿값을 구합니다. pdeval을 사용하여 xmesh에 포함되지 않은 지점에서 근삿값과 근삿값의 편도함수 ∂u_i/∂x를 계산합니다. 자세한 내용은 pdeval을 참조하십시오.

tsol — 풀이 시간

열 벡터

풀이 시간으로, tspan에 지정한 시간 중에서 첫 번째 종료 이벤트 이전에 해당하는 시간의 열 벡터로 반환됩니다.

sole — 이벤트 발생 시 계산된 해

배열

이벤트 발생 시 계산된 해로, 배열로 반환됩니다. te의 이벤트 시간은 sole로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

te — 이벤트 시간

열 벡터

이벤트 시간으로, 열 벡터로 반환됩니다. te의 이벤트 시간은 sole로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

ie — 트리거된 이벤트 함수의 인덱스

열 벡터

트리거된 이벤트 함수의 인덱스로, 열 벡터로 반환됩니다. te의 이벤트 시간은 sole로 반환되는 해에 대응되며, ie는 발생한 이벤트를 지정합니다.

팁

• uji = sol(j,:,i)가 시간 tspan(j)와 메시 점 xmesh에서 구한 해의 성분 i를 근사하면, pdeval은 점으로 구성된 배열 xout에서 이러한 근사 및 그 편도함수 ∂u_i/∂x의 값을 계산해 uout과 duoutdx에 다음과 같이 반환합니다. [uout,duoutdx] = pdeval(m,xmesh,uji,xout). pdeval 함수는 플럭스가 아니라 편도함수 ∂u_i/∂x를 계산합니다. 플럭스는 연속적이지만, 물질의 경계면에서는 편도함수에 비약이 있을 수 있습니다.

알고리즘

ode15s 솔버로 시간 적분이 수행됩니다. pdepe는 ode15s의 기능을 활용하여, PDE에 타원 방정식이 포함되어 있을 때 제기되는 미분대수 방정식(DAE)을 풀고 지정된 희소성 패턴을 갖는 야코비 행렬을 처리합니다.

이산화 후에는 타원 방정식에서 대수 방정식이 생성됩니다. 타원 방정식에 해당하는 초기 조건 벡터의 요소가 이산화를 수행한 후 일관되지 않을 경우 pdepe는 시간에 대한 적분을 시작하기 전에 이러한 요소를 조정하려고 합니다. 이러한 이유로, 초기 시간에 대해 반환되는 해에는 기타 모든 시간에 비해 이산화 오차가 있을 수 있습니다. 메시가 충분히 세밀한 경우 pdepe는 지정된 값에 가까운 모순 없는 초기 조건을 구할 수 있습니다. pdepe가 모순 없는 초기 조건을 구하기 어렵다는 메시지를 표시하면 메시를 더 세밀하게 조정해 보십시오. 포물선 방정식에 해당하는 초기 조건 벡터의 요소는 조정이 필요하지 않습니다.

참고 문헌

[1] Skeel, R. D. and M. Berzins, "A Method for the Spatial Discretization of Parabolic Equations in One Space Variable," SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp.1–32.

버전 내역

R2006a 이전에 개발됨

참고 항목

pdeval | ode15s | odeset | odeget

도움말 항목

• 편미분 방정식 풀기
• 단일 PDE 풀기
• 연립 PDE 풀기

×

예제 열기

이 예제의 수정된 버전이 있습니다. 사용자가 편집한 내용을 반영하여 이 예제를 여시겠습니까?

아니요, 수정된 버전을 덮어씁니다. 예