초기 조건으로 계단 함수를 사용하는 연립 PDE 풀기

이 예제에서는 초기 조건에 계단 함수를 사용하는 연립편미분방정식을 푸는 방법을 보여줍니다.

다음과 같은 PDE가 있다고 가정해 보겠습니다.

$$\begin{array}{l}\frac{\partial \mathit{n}}{\partial \mathit{t}}=\frac{\partial }{\partial \mathit{x}}\left[\mathit{d}\frac{\partial \mathit{n}}{\partial \mathit{x}}-\mathit{a}\mathit{n}\frac{\partial \mathit{c}}{\partial \mathit{x}}\right]+\mathit{S}\mathit{r}\mathit{n}\left(\mathit{N}-\mathit{n}\right),\\ \frac{\partial \mathit{c}}{\partial \mathit{t}}=\frac{{\partial }^2\mathit{c}}{\partial {\mathit{x}}^2}+\mathit{S}\left(\frac{\mathit{n}}{\mathit{n}+1}-\mathit{c}\right).\end{array}$$

이 방정식은 상수 파라미터 $\mathit{d}$, $\mathit{a}$, $\mathit{S}$, $\mathit{r}$ 및 $\mathit{N}$을 포함하고 $0\le \mathit{x}\le 1$ 및 $\mathit{t}\ge 0$에 대해 정의됩니다. 이 방정식은 종양 관련 혈관 형성의 첫 번째 단계에 대한 수학 모델 [1]에서 기인합니다. $\mathit{n}\left(\mathit{x},\mathit{t}\right)$는 내피 세포의 세포 밀도를 나타내고, $\mathit{c}\left(\mathit{x},\mathit{t}\right)$는 내피 세포가 종양에 반응하여 방출하는 단백질 농도를 나타냅니다.

이 문제는 다음과 같은 경우에 일정한 정상 상태입니다.

$\left[\begin{array}{c}{\mathit{n}}_0\\ {\mathit{c}}_0\end{array}\right]=\left[\begin{array}{c}1\\ 0.5\end{array}\right]$.

하지만 안정성 분석에서는 시스템이 비동차해로 변화할 것으로 예측합니다 [1]. 따라서 계단 함수는 초기 조건으로 사용되어 정상 상태를 섭동하고 시스템의 변화를 촉진합니다.

경계 조건의 경우 두 해 성분 모두 $\mathit{x}=0$ 및 $\mathit{x}=1$에서 플럭스가 0이어야 합니다.

$$\begin{array}{l}\frac{\partial }{\partial \mathit{x}}\mathit{n}\left(0,\mathit{t}\right)=\frac{\partial }{\partial \mathit{x}}\mathit{n}\left(1,\mathit{t}\right)=0,\\ \frac{\partial }{\partial \mathit{x}}\mathit{c}\left(0,\mathit{t}\right)=\frac{\partial }{\partial \mathit{x}}\mathit{c}\left(1,\mathit{t}\right)=0.\end{array}$$

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

방정식 코딩하기

방정식을 코딩하려면 먼저 방정식이 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에 방정식이 2개 있기 때문에 연립 PDE를 다음과 같이 다시 작성할 수 있습니다.

$$\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\frac{\partial }{\partial \mathit{t}}\left[\begin{array}{c}\mathit{n}\\ \mathit{c}\end{array}\right]=\frac{\partial }{\partial \mathit{x}}\left[\begin{array}{c}\mathit{d}\frac{\partial \mathit{n}}{\partial \mathit{x}}-\mathit{a}\mathit{n}\frac{\partial \mathit{c}}{\partial \mathit{x}}\\ \frac{\partial \mathit{c}}{\partial \mathit{x}}\end{array}\right]+\left[\begin{array}{c}\mathit{S}\mathit{r}\mathit{n}\left(\mathit{N}-\mathit{n}\right)\\ \mathit{S}\left(\frac{\mathit{n}}{\mathit{n}+1}-\mathit{c}\right)\end{array}\right].$$

방정식의 계수 값은 다음과 같습니다.

$$\mathit{m}=0$$

$$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\left[\begin{array}{c}1\\ 1\end{array}\right]$$(대각 값만 해당)

$$\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\left[\begin{array}{c}\mathit{d}\frac{\partial \mathit{n}}{\partial \mathit{x}}-\mathit{a}\mathit{n}\frac{\partial \mathit{c}}{\partial \mathit{x}}\\ \frac{\partial \mathit{c}}{\partial \mathit{x}}\end{array}\right]$$

$$\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\left[\begin{array}{c}\mathit{S}\mathit{r}\mathit{n}\left(\mathit{N}-\mathit{n}\right)\\ \mathit{S}\left(\frac{\mathit{n}}{\mathit{n}+1}-\mathit{c}\right)\end{array}\right]$$

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

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

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

• u는 x와 t에 따라 미분되는 종속 변수입니다. 요소를 2개 가진 벡터로, 여기서 u(1)은 $\mathit{n}\left(\mathit{x},\mathit{t}\right)$이고 u(2)는 $\mathit{c}\left(\mathit{x},\mathit{t}\right)$입니다.

• dudx는 편미분 공간 도함수 $\partial \mathit{u}/\partial \mathit{x}$입니다. 요소를 2개 가진 벡터로, 여기서 dudx(1)은 $\partial \mathit{n}/\partial \mathit{x}$이고 dudx(2)는 $\partial \mathit{c}/\partial \mathit{x}$입니다.

• 출력값 c, f 및 s는 pdepe에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다.

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

function [c,f,s] = angiopde(x,t,u,dudx)
d = 1e-3;
a = 3.8;
S = 3;
r = 0.88;
N = 1;

c = [1; 1];
f = [d*dudx(1) - a*u(1)*dudx(2)
     dudx(2)];
s = [S*r*u(1)*(N - u(1));
     S*(u(1)/(u(1) + 1) - u(2))];
end

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

초기 조건 코딩하기

이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값에서 적용되고 모든 x 값에 대한 $\mathit{n}\left(\mathit{x},{\mathit{t}}_0\right)$ 및 $\mathit{c}\left(\mathit{x},{\mathit{t}}_0\right)$ 값을 제공합니다. 함수 시그니처 u0 = angioic(x)를 사용하여 함수를 작성하십시오.

이 문제는 다음과 같은 경우에 일정한 정상 상태입니다.

$\left[\begin{array}{c}{\mathit{n}}_0\\ {\mathit{c}}_0\end{array}\right]=\left[\begin{array}{c}1\\ 0.5\end{array}\right]$.

하지만 안정성 분석에서는 시스템이 비동차해로 변화할 것으로 예측합니다 [1]. 따라서 계단 함수는 초기 조건으로 사용되어 정상 상태를 섭동하고 시스템의 변화를 촉진합니다.

$$\begin{array}{l}\mathit{u}\left(\mathit{x},0\right)=\left[\begin{array}{c}{\mathit{n}}_0\\ {\mathit{c}}_0\end{array}\right],\\ \mathit{u}\left(\mathit{x},0\right)=\{\begin{array}{cc}1.05{\mathit{u}}_1& 0.3\le \mathit{x}\le 0.6\\ 1.0005{\mathit{u}}_2& 0.3\le \mathit{x}\le 0.6\end{array}\end{array}$$

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

function u0 = angioic(x)
u0 = [1; 0.5];
if x >= 0.3 && x <= 0.6
  u0(1) = 1.05 * u0(1);
  u0(2) = 1.0005 * u0(2);
end
end

경계 조건 코딩하기

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

$$\begin{array}{l}\frac{\partial }{\partial \mathit{x}}\mathit{n}\left(0,\mathit{t}\right)=\frac{\partial }{\partial \mathit{x}}\mathit{n}\left(1,\mathit{t}\right)=0,\\ \frac{\partial }{\partial \mathit{x}}\mathit{c}\left(0,\mathit{t}\right)=\frac{\partial }{\partial \mathit{x}}\mathit{c}\left(1,\mathit{t}\right)=0.\end{array}$$

구간 $\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{x}=0$인 경우, 경계 조건 방정식은 다음과 같습니다.

$$\left[\begin{array}{c}0\\ 0\end{array}\right]+\left[\begin{array}{c}1\\ 1\end{array}\right]\cdot \left[\begin{array}{c}\mathit{d}\frac{\partial \mathit{n}}{\partial \mathit{x}}-\mathit{a}\mathit{n}\frac{\partial \mathit{c}}{\partial \mathit{x}}\\ \frac{\partial \mathit{c}}{\partial \mathit{x}}\end{array}\right]=0.$$

계수는 다음과 같습니다.

• $${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)=\left[\begin{array}{c}0\\ 0\end{array}\right],$$

• ${\mathit{q}}_{\mathit{L}}\left(\mathit{x},\mathit{t}\right)=\left[\begin{array}{c}1\\ 1\end{array}\right]$.

$\mathit{x}=1$인 경우, 경계 조건이 동일하므로 ${\mathit{p}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)=\left[\begin{array}{c}0\\ 0\end{array}\right]$ 및 ${\mathit{q}}_{\mathit{R}}\left(\mathit{x},\mathit{t}\right)=\left[\begin{array}{c}1\\ 1\end{array}\right]$이 됩니다.

경계 함수는 함수 시그니처 [pl,ql,pr,qr] = angiobc(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] = angiobc(xl,ul,xr,ur,t)
pl = [0; 0];
ql = [1; 1];
pr = pl;
qr = ql;
end

해 메시 선택하기

방정식의 제한 동작을 확인하려면 긴 시간 구간이 필요하므로, 구간 $0\le \mathit{t}\le 200$에서 10개의 점을 사용하십시오. 또한 $\mathit{c}\left(\mathit{x},\mathit{t}\right)$의 극한 분포는 구간 $0\le \mathit{x}\le 1$에서 약 0.1%만 변화하므로, 50개의 점으로 구성된 상대적으로 세밀한 공간 메시가 적합합니다.

x = linspace(0,1,50);
t = linspace(0,200,10);

방정식 풀기

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

m = 0;
sol = pdepe(m,@angiopde,@angioic,@angiobc,x,t);

pdepe는 3차원 배열 sol로 해를 반환합니다. 여기서 sol(i,j,k)는 t(i)와 x(j)에서 구한 해의 k번째 성분 ${\mathit{u}}_{\mathit{k}}$에 대한 근삿값을 계산합니다. 해 성분을 별도의 변수로 추출하십시오.

n = sol(:,:,1);
c = sol(:,:,2);

해 플로팅하기

$\mathit{x}$와 $\mathit{t}$에 대해 선택한 메시 점에서 플로팅된 해 성분 $\mathit{n}$과 $\mathit{c}$의 곡면 플롯을 만듭니다.

surf(x,t,c)
title('c(x,t): Concentration of Fibronectin')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,n)
title('n(x,t): Density of Endothelial Cells')
xlabel('Distance x')
ylabel('Time t')

이제 ${\mathit{t}}_{\mathit{f}}=200$일 때 해의 최종 분포를 플로팅합니다. 이 플롯은 [1]의 Figure 3 및 4에 대응합니다.

plot(x,n(end,:))
title('Final distribution of n(x,t_f)')

plot(x,c(end,:))
title('Final distribution of c(x,t_f)')

참고 문헌

[1] Humphreys, M.E.와 M.A.J. Chaplain. "A mathematical model of the first steps of tumour-related angiogenesis: Capillary sprout formation and secondary branching." IMA Journal of Mathematics Applied in Medicine & Biology, 13 (1996), pp. 73-98.

로컬 함수(Local Function)

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

function [c,f,s] = angiopde(x,t,u,dudx) % Equation to solve
d = 1e-3;
a = 3.8;
S = 3;
r = 0.88;
N = 1;

c = [1; 1];
f = [d*dudx(1) - a*u(1)*dudx(2)
     dudx(2)];
s = [S*r*u(1)*(N - u(1));
     S*(u(1)/(u(1) + 1) - u(2))];
end
% ---------------------------------------------
function u0 = angioic(x) % Initial Conditions
u0 = [1; 0.5];
if x >= 0.3 && x <= 0.6
  u0(1) = 1.05 * u0(1);
  u0(2) = 1.0005 * u0(2);
end
end
% ---------------------------------------------
function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [0; 0];
ql = [1; 1];
pr = pl;
qr = ql;
end
% ---------------------------------------------

참고 항목

pdepe

관련 항목

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