Main Content

使用初始条件阶跃函数求解 PDE 方程组

此示例说明如何求解初始条件中使用步函数的偏微分方程组。

以如下 PDE 为例

nt=x[dnx-ancx]+Srn(N-n),ct=2cx2+S(nn+1-c).

这些方程涉及常量参数 daSrN,并为 0x1t0 定义。这些方程源于描述肿瘤相关血管生成的初始步骤的数学模型 [1]。n(x,t) 表示内皮细胞的细胞密度,c(x,t) 表示因为肿瘤而释放的蛋白质的浓度。

当出现以下情况时,此问题实现常数稳态

[n0c0]=[10.5].

然而,稳定性分析预测方程组会演化出非齐次解 [1]。因此,需要使用阶跃函数作为初始条件,以扰动稳态和促进方程组演化。

边界条件要求两个解分量在 x=0x=1 处通量为零。

xn(0,t)=xn(1,t)=0,xc(0,t)=xc(1,t)=0.

要在 MATLAB® 中求解此方程组,您需要对方程、初始条件和边界条件编写代码,然后在调用求解器 pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

在编写方程代码之前,您需要确保它的形式符合 pdepe 求解器的要求:

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

由于该 PDE 方程组中有两个方程,PDE 方程组可以重写为

[1001]t[nc]=x[dnx-ancxcx]+[Srn(N-n)S(nn+1-c)].

则方程中系数的值为

m=0

c(x,t,u,ux)=[11](仅对角线值)

f(x,t,u,ux)=[dnx-ancxcx]

s(x,t,u,ux)=[Srn(N-n)S(nn+1-c)]

现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = angiopde(x,t,u,dudx)

  • x 是独立的空间变量。

  • t 是独立的时间变量。

  • u 是关于 xt 微分的因变量。它是二元素向量,其中 u(1)n(x,t)u(2)c(x,t)

  • dudx 是偏空间导数 u/x。它是二元素向量,其中 dudx(1)n/xdudx(2)c/x

  • 输出 cfs 对应于 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 的任何值提供 n(x,t0)c(x,t0) 的值。使用函数签名 u0 = angioic(x) 编写函数。

当出现以下情况时,此问题实现常数稳态

[n0c0]=[10.5].

然而,稳定性分析预测方程组会演化出非齐次解[1]。因此,需要使用阶跃函数作为初始条件,以扰动稳态和促进方程组演化。

u(x,0)=[n0c0],u(x,0)={1.05u10.3x0.61.0005u20.3x0.6

对应的函数是

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

编写边界条件代码

现在,编写计算以下边界条件的函数

xn(0,t)=xn(1,t)=0,xc(0,t)=xc(1,t)=0.

对于在区间 axb 上提出的问题,边界条件应用于所有 t 以及 x=ax=b。求解器所需的边界条件的标准形式是

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

对于 x=0,边界条件方程为

[00]+[11][dnx-ancxcx]=0.

因此系数为:

  • pL(x,t,u)=[00],

  • qL(x,t)=[11].

对于 x=1,边界条件是相同的,因此 pR(x,t,u)=[00]qR(x,t)=[11]

边界函数应使用函数签名 [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t),其中:

  • 对于左边界,输入 xlul 对应于 ux

  • 对于右边界,输入 xrur 对应于 ux

  • t 是独立的时间变量。

  • 对于左边界,输出 plql 对应于 pL(x,t,u)qL(x,t)(对于此问题,x=0)。

  • 对于右边界,输出 prqr 对应于 pR(x,t,u)qR(x,t)(对于此问题,x=1)。

此示例中的边界条件由以下函数表示:

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

选择解网格

要了解方程的限制行为,需要很长的时间区间,因此使用 10 个位于区间 0t200 中的点。此外,在 0x1 区间内,c(x,t) 的限值分布仅变化约 0.1%,因此具有 50 个点的相对精细的空间网格是合适的。

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

求解方程

最后,使用对称性值 m、PDE 方程、初始条件、边界条件以及 xt 的网格来求解方程。

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

pdepe 以三维数组 sol 形式返回解,其中 sol(i,j,k) 是在 t(i)x(j) 处计算的解 uk 的第 k 个分量的逼近值。将解分量提取到单独的变量中。

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

对解进行绘图

创建基于所选的 xt 网格点绘制的解分量 nc 的曲面图。

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

Figure contains an axes object. The axes object with title c(x,t): Concentration of Fibronectin, xlabel Distance x, ylabel Time t contains an object of type surface.

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

Figure contains an axes object. The axes object with title n(x,t): Density of Endothelial Cells, xlabel Distance x, ylabel Time t contains an object of type surface.

现在,仅绘制在 tf=200 处的解的最终分布。这些图对应于 [1] 中的图 3 和 4。

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

Figure contains an axes object. The axes object with title Final distribution of n(x,t indexOf f baseline ) contains an object of type line.

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

Figure contains an axes object. The axes object with title Final distribution of c(x,t indexOf f baseline ) contains an object of type line.

参考

[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
% ---------------------------------------------

另请参阅

相关主题