matlab使用杂谈4-偏微分方程求解之pdede函数使用

偏微分方程

偏微分方程（Partial Differential Equation，PDE）指含有未知函数及其偏导数的方程，自变量的个数为两个或两个以上。描述自变量、未知函数及偏导数之间的关系
偏微分方程分为
1、线性偏微分方程式
2、非线性偏微分方程式

求解偏微分方程的数值方法

1、有限元法（Finite Element Method，FEM）
2、有限体积法（Finite Volume Method，FVM）
3、有限差分法（Finite Difference Method，FDM）
还有其他FFEM、XFEM、MFEM、DGFEM等方法

Matlab解偏微分方程

pdepe()函数

可求解一般的PDEs，具有较大的通用性，但指支持命令行形式调用
matlab中关于pdepe函数的用法
语法
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)
[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

变量

参数	含义
m	与该问题的对称性相对应的参数。m 可以是平板 = 0，柱状 = 1，或球面 = 2。
pdefun	定义函数句柄
icfun	初始条件函数句柄
bcfun	边界条件函数句柄
xmesh	向量 [x0, x1, …, xn]，用于指定需要针对 tspan 中每个值求数值解的点。xmesh 的元素必须满足 x0 < x1 < … < xn。xmesh 的长度必须 >= 3。
tspan	向量 [t0, t1, …, tf]，用于指定需要针对 xmesh 中每个值求解的点。tspan 的元素必须满足 t0 < t1 < … < tf。tspan 的长度必须 >= 3。
options	基础 ODE 求解器的部分选项可以在 pdepe 中使用：RelTol、AbsTol、NormControl、InitialStep、MaxStep 和 Events。在大多数情况下，这些选项的默认值可提供满意解

分割线----------------------------------------------------------------------------------
注 ：该部分涉及数学知识可不过分关注，但需要充分理解方程格式

pdepe函数主要用来结算以下格式的PDE方程

x与t的范围需要处于有限范围
对于 t = t0 和所有 x，解分量均满足以下格式的初始条件

对于所有 t 和 x = a 或 x = b，解分量满足以下形式的边界条件

q 的元素全部为零或都不为零（至于原因是什么，我还没有搞懂…）。请注意，边界条件以通量 f 的方式而不是 ∂u/∂x 表示。同时，在这两个系数之间，只有 p 可以依赖于 u。

分割线----------------------------------------------------------------------------------
在调用 sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) 中：
m 与 m 对应。
xmesh(1) 和 xmesh(end) 对应 a 和 b。
tspan(1) 和 tspan(end) 对应 t₀ 和 t_f 。

pdefun 计算项 c、f 和 s (公式 1)。其格式为
[c,f,s] = pdefun(x,t,u,dudx)
输入参数为标量 x 和 t，以及向量 u 和 dudx，分别接近于解 u 及其相对于 x 的偏导数
icfun 计算初始条件。其格式为

u = icfun(x)
当与参数 x 一起调用时，icfun 会计算并返回列向量 u 中 x 处的解分量的初始值。

bcfun 计算边界条件 (公式 3) 的项 p 和 q。其格式为[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
ul 是在左边界 xl = a 的近似解，ur 是在右边界 xr = b 的近似解。pl 和 ql 是对应在 xl 上计算的 p 和 q 的列向量，类似地 pr 和 qr 对应 xr。当 m > 0 且 a = 0 时，x = 0 附近解的有界性需要通量 f 在 a = 0 时消失。pdepe 会自动设置此边界条件并且会忽略 pl 和 ql 中返回的值。

pdepe 以多维数组 sol 的形式返回解。ui = ui = sol(:,:,i) 是解向量 u 的第 i 个分量的近似值。元素 ui(j,k) = sol(j,k,i) 在 (t,x) = (tspan(j),xmesh(k)) 处近似于 ui。

上面的讲解主要摘抄于matlab的官方文档，还需要更加详细的讲解可以参考matlab中关于pdepe函数的用法，该文档目前也是中文文档，便于观看，不过理解起来需要挺久，也可以直接看看我下面的实例，根据实际例子为你一一解开PDE的解法

pdepe函数使用示例

由于CSDN中不好输入公式，我以图片格式上传微分方程，对于下列微分方程、初始条件及边界条件，一一转换为matlab中pdepe函数求解格式

PDE方程求解格式

对应matlab代码

 function [c,f,s]=pdex1pde(x,t,u,DuDx)
    c=pi^2;
    f=DuDx;
    s=0;
 end

PDE方程初始条件格式

对应matlab代码

function uo=pdex1ic(x)
    uo=sin(pi*x);
end

PDE边界条件格式

对应matlab代码

function [pl,ql,pr,qr]=pdex1bc(x1,u1,xr,ur,t)
    pl=u1; 
    ql=0; 
    pr=pi*exp(-t); 
    qr=1; 
    end

Matlab代码

m=0;
x=linspace(0,1,20); % 方程区间为(0,1)
t=linspace(0,2,10); % t 的范围可以随取，只需要大于0即可
sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % 10x20的矩阵与t*x的维度一致
u = sol(:,:,1); %解向量 u 的第 1 个分量的近似值
surf(x,t,u) 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

figure
plot(x,u(end,:))% t=2时，u随x的变化曲线
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')

求解偏微分方程

求解偏微分方程的过程其实与上述过程一致，只不过所有步骤都使用向量格式来存储，此处直接Copy Matlab中的实例和代码，自己按照上面研究一下即可。


对应matlab代码

m = 0;
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];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(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]; 
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [1; 0]; 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0; ul(2)]; 
ql = [1; 0]; 
pr = [ur(1)-1; 0]; 
qr = [0; 1]; 

以上代码引用自matlab官方函数文档

总结

上述就是matlab中解偏微分方程中pdepe函数的用法，说实话有些复杂，我也是看了挺久才理解，但由于没有真正应用在建模实例上，现在也不知道对他的掌握情况是咋样，之后建模过程中上传一下我的应用实例给大家参考一下，嘿嘿。