初期条件のステップ関数をもつ 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$ において両方の解要素の流束がゼロであることが要求されます。

$$\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

(メモ: 関数はすべて例の最後にローカル関数として含まれます。)

初期条件のコード化

次に、初期条件を返す関数を記述します。初期条件は最初の時間値において適用され、$\mathit{n}\left(\mathit{x},{\mathit{t}}_0\right)$ と $\mathit{c}\left(\mathit{x},{\mathit{t}}_0\right)$ の値を任意の x の値に対して提供します。関数シグネチャ 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] の図 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. and 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.

ローカル関数

ここでは、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 系の求解