1求解 PDE 方程组
此示例说明由两个偏微分方程构成的方程组的解的构成，以及如何对解进行计算和绘图。
以如下 PDE 方程组为例
要在 MATLAB® 中求解该方程，您需要对方程、初始条件和边界条件编写代码，然后在调用求解器pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
1.1编写方程代码
在编写方程代码之前，您需要确保它的形式符合 pdepe
求解器的要求：
因此，此示例中的方程可由以下函数表示：
function [c,f,s] = pdefun(x,t,u,dudx)
c = [1; 1];
f = [0.024; 0.17] .* dudx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
end1.2编写初始条件代码
接下来，编写一个返回初始条件的函数。初始条件应用在第一个时间值处，并为 x 的任何值提供 u (x
,
t
0 )的 值。初始条件的数量必须等于方程的数量，因此对于此问题，有两个初始条件。使用函数签名
u0 =
pdeic(x)
编写函数。
对应的函数是
function u0 = pdeic(x)
u0 = [1; 0];
end1.3编写边界条件代码
现在，编写计算以下边界条件的函数
此示例中的边界条件由以下函数表示：
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end1.4选择解网格
当 t
较小时，此问题的解会快速变化。虽然
pdepe 选择了适合解析急剧变化的时间步，但要在输出绘图中 显示该行为，您需要选择适当的输出时间。对于空间网格，在
0 ≤
x
≤ 1 两端的解中都存在边界层，因此 您需要在那里指定网格点来解析急剧变化。
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];1.5求解方程
最后，使用对称性值
m
、PDE 方程、初始条件、边界条件以及
x
和
t
的网格来求解方程。
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
pdepe 以三维数组
sol
形式返回解，其中
sol(i,j,k)
是在
t(i)
和
x(j)
处计算的解
u
k
的第
k 个分量的逼近 值。将每个解分量提取到一个单独变量中。
u1 = sol(:,:,1);
u2 = sol(:,:,2);
1.6 对解进行绘图
创建在 x
和
t
的所选网格点上绘制的
u
1
和
u
2
的解的曲面图。
surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')
surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')
1.6 局部函数
此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve
c = [1; 1];
f = [0.024; 0.17] .* dudx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
end
% ---------------------------------------------
function u0 = pdeic(x) % Initial Conditions
u0 = [1; 0];
end
% ---------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end2使用初始条件阶跃函数求解 PDE 方程组
此示例说明如何求解初始条件中使用步函数的偏微分方程组。以如下 PDE 为例
要在 MATLAB® 中求解此方程组，您需要对方程、初始条件和边界条件编写代码，然后在调用求解器 pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
2.1编写方程代码
在编写方程代码之前，您需要确保它的形式符合 pdepe
求解器的要求：
因此，此示例中的方程可由以下函数表示
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))];
end2.2编写初始条件代码
然而，稳定性分析预测方程组会演化出非齐次解。因此，需要使用阶跃函数作为初始条件，以扰动稳态和促进方程组演化。
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
end2.3 编写边界条件代码
function [pl,ql,pr,qr] = angiobc(xl,ul,xr,ur,t)
pl = [0; 0];
ql = [1; 1];
pr = pl;
qr = ql;
end2.4选择解网格
要了解方程的限制行为，需要很长的时间区间，因此使用 10 个位于区间 0 ≤
t
≤ 200 中的点。此外，在0 ≤
x
≤ 1
区间内，
c(x
,
t
)的限值分布仅变化约 0.1%，因此具有 50 个点的相对精细的空间网格是合适的。
x = linspace(0,1,50);
t = linspace(0,200,10);2.5求解方程
最后，使用对称性值
m
、PDE 方程、初始条件、边界条件以及
x
和
t
的网格来求解方程。
m = 0;
sol = pdepe(m,@angiopde,@angioic,@angiobc,x,t);
pdepe 以三维数组
sol
形式返回解，其中
sol(i,j,k)
是在
t(i)
和
x(j)
处计算的解
u
k
的第
k 个分量的逼近值。将解分量提取到单独的变量中。
n = sol(:,:,1);
c = sol(:,:,2);2.6 对解进行绘图
创建基于所选的 x
和
t
网格点绘制的解分量
n
和
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')
