matlab使用教程(32)—求解偏微分方程(3)

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];
end

1.2编写初始条件代码

        接下来，编写一个返回初始条件的函数。初始条件应用在第一个时间值处，并为 x 的任何值提供 u (x , t 0 )的 值。初始条件的数量必须等于方程的数量，因此对于此问题，有两个初始条件。使用函数签名 u0 = pdeic(x) 编写函数。

        对应的函数是

function u0 = pdeic(x)
u0 = [1; 0];
end

1.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];
end

1.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];
end

2使用初始条件阶跃函数求解 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))];
end

2.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
end

2.3 编写边界条件代码

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

2.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')