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

相关推荐

  1. Julia调用Matlab, Python以及R的微分方程求解

    2024-04-03 18:26:02       29 阅读

最近更新

  1. leetcode705-Design HashSet

    2024-04-03 18:26:02       8 阅读
  2. Unity发布webgl之后打开streamingAssets中的html文件

    2024-04-03 18:26:02       8 阅读
  3. vue3、vue2中nextTick源码解析

    2024-04-03 18:26:02       9 阅读
  4. 高级IO——React服务器简单实现

    2024-04-03 18:26:02       8 阅读
  5. 将图片数据转换为张量(Go并发处理)

    2024-04-03 18:26:02       7 阅读
  6. go第三方库go.uber.org介绍

    2024-04-03 18:26:02       8 阅读
  7. 前后端AES对称加密 前端TS 后端Go

    2024-04-03 18:26:02       10 阅读

热门阅读

  1. 学习总结!

    2024-04-03 18:26:02       4 阅读
  2. Vue3中props和emits的使用总结

    2024-04-03 18:26:02       6 阅读
  3. IO和NIO的主要区别在哪里?

    2024-04-03 18:26:02       5 阅读
  4. CSS 滚动条样式修改

    2024-04-03 18:26:02       5 阅读
  5. codeforces round 936 div2(a,b,c)

    2024-04-03 18:26:02       6 阅读
  6. Python程序设计 魔法函数

    2024-04-03 18:26:02       5 阅读
  7. ACI Fabric

    2024-04-03 18:26:02       5 阅读
  8. wencoo个人的博客目录索引-更新

    2024-04-03 18:26:02       4 阅读
  9. Oracle extent、segment、datafile、tablespace及存储关系

    2024-04-03 18:26:02       4 阅读
  10. Oracle密码文件

    2024-04-03 18:26:02       4 阅读
  11. SpringMVC 中实现自定义转换

    2024-04-03 18:26:02       5 阅读
  12. 开源的分布式文件系统 Fastdfs 安装入门介绍

    2024-04-03 18:26:02       5 阅读
  13. 在stable diffusion中如何分辨lora、大模型、controlnet

    2024-04-03 18:26:02       5 阅读