使用初始条件阶跃函数求解 PDE 方程组
此示例说明如何求解初始条件中使用步函数的偏微分方程组。
以如下 PDE 为例
这些方程涉及常量参数 、、、 和 ,并为 和 定义。这些方程源于描述肿瘤相关血管生成的初始步骤的数学模型 [1]。 表示内皮细胞的细胞密度, 表示因为肿瘤而释放的蛋白质的浓度。
当出现以下情况时,此问题实现常数稳态
.
然而,稳定性分析预测方程组会演化出非齐次解 [1]。因此,需要使用阶跃函数作为初始条件,以扰动稳态和促进方程组演化。
边界条件要求两个解分量在 和 处通量为零。
要在 MATLAB® 中求解此方程组,您需要对方程、初始条件和边界条件编写代码,然后在调用求解器 pdepe
之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写方程代码
在编写方程代码之前,您需要确保它的形式符合 pdepe
求解器的要求:
由于该 PDE 方程组中有两个方程,PDE 方程组可以重写为
则方程中系数的值为
(仅对角线值)
现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = angiopde(x,t,u,dudx)
:
x
是独立的空间变量。t
是独立的时间变量。u
是关于x
和t
微分的因变量。它是二元素向量,其中u(1)
是 ,u(2)
是 。dudx
是偏空间导数 。它是二元素向量,其中dudx(1)
是 ,dudx(2)
是 。输出
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 的任何值提供 和 的值。使用函数签名 u0 = angioic(x)
编写函数。
当出现以下情况时,此问题实现常数稳态
.
然而,稳定性分析预测方程组会演化出非齐次解[1]。因此,需要使用阶跃函数作为初始条件,以扰动稳态和促进方程组演化。
对应的函数是
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
编写边界条件代码
现在,编写计算以下边界条件的函数
对于在区间 上提出的问题,边界条件应用于所有 以及 或 。求解器所需的边界条件的标准形式是
对于 ,边界条件方程为
因此系数为:
.
对于 ,边界条件是相同的,因此 且 。
边界函数应使用函数签名 [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t)
,其中:
对于左边界,输入
xl
和ul
对应于 和 。对于右边界,输入
xr
和ur
对应于 和 。t
是独立的时间变量。对于左边界,输出
pl
和ql
对应于 和 (对于此问题,)。对于右边界,输出
pr
和qr
对应于 和 (对于此问题,)。
此示例中的边界条件由以下函数表示:
function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t) pl = [0; 0]; ql = [1; 1]; pr = pl; qr = ql; end
选择解网格
要了解方程的限制行为,需要很长的时间区间,因此使用 10 个位于区间 中的点。此外,在 区间内, 的限值分布仅变化约 0.1%,因此具有 50 个点的相对精细的空间网格是合适的。
x = linspace(0,1,50); t = linspace(0,200,10);
求解方程
最后,使用对称性值 、PDE 方程、初始条件、边界条件以及 和 的网格来求解方程。
m = 0; sol = pdepe(m,@angiopde,@angioic,@angiobc,x,t);
pdepe
以三维数组 sol
形式返回解,其中 sol(i,j,k)
是在 t(i)
和 x(j)
处计算的解 的第 k
个分量的逼近值。将解分量提取到单独的变量中。
n = sol(:,:,1); c = sol(:,:,2);
对解进行绘图
创建基于所选的 和 网格点绘制的解分量 和 的曲面图。
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')
现在,仅绘制在 处的解的最终分布。这些图对应于 [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 % ---------------------------------------------