偏微分方程组的MATLAB解法.pdf
 
第第3737 卷第卷第23 期期··学术学术 湖南农机湖南农机 年年35 月月
 
第37 卷第3 期 张向利：户县拖拉机及驾驶员年度检审现状及对策
 
Vol.37Vol.37 No.2No.3 HUNANHUNAN AGRICULTURALAGRICULTURAL MACHINERYMACHINERY May.2010Mar.2010
 
偏微分方程的MATLAB解法
 
李 明
 
(长江大学信息与数学学院，湖北荆州 434023)
 
摘 要：一般来说，在实际中偏微分方程的通解是不容易求出的，用定解条件确定函数更是比较困难。本文在运
 
用MATLAB 解偏微分方程时，列出两种方法：pdepe 函数和PDE 工具箱，并应用实例展现出两种方法的实现过程。结
 
果表明：MATLAB 对解偏微分方程带来极大的方便，并且在此基础上可以解决更多更复杂的问题。
 
关键词：偏微分方程；MATLAB ；PED Toolbox
 
： ： ：
 
中图分类号 TE121 文献标识码 A 文章编号 1007- 8320(2010)03- 0089- 02
 
MATLAB Solution partial differential equations
 
LI Ming
 
( ， ， ， ， )
 
School of Mathematics Yangtze University Jingzhou Hubei 434023 China
 
： ， ， ，
 
Abstract Generally speaking in practice the general solution of partial differential equations is not easy to find and
 
use the boundary conditions to determine the function is more difficult. This articl using MATLAB for solving partial
 
， ：
 
differential equations two methods are listed pdepe function and the PDE Toolbox.On the basis of describing the two
 
，
 
methods in detail this article uses examples

 
文章目录
 
(1)偏微分方程的类型(二阶)
(2)抛物线型
1.显式法
2.Crank-Nicholson隐式算法 
  
(3)双曲线型
(4)椭圆型

 
(1)偏微分方程的类型(二阶)
 
$$a\frac{{\partial }^{2}u}{\partial x^2}+b\frac{{\partial }^{2}u}{\partial y\partial x}+c\frac{{\partial }^{2}u}{\partial x^2}+d\frac{\partial u}{\partial x}+e\frac{\partial u}{\partial y}+fu+g=0$$
 
$b^2-4ac<0$ 椭圆
$b^2-4ac=0$ 抛物线
$b^2-4ac>0$ 双曲线
 
(2)抛物线型
 
 
  
 
 
  
  
 
1.显式法
 
求解思想：通过差分的方法一排一排向上推。
做划分并代入方程$$\frac{u_{i,j+1}-u_{i,j}}{k}=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}(\Delta x=h,\Delta t=k)$$
通过化简得到$$u_{i,j+1}=r{u}_{i-1,j}+(1-2r)u_{i,j}+r{u}_{i+1,j}(r=\frac{k}{h^2})$$
具体推的步骤大概如下： 
  
由于已知 $u(x,0)=f(x)$，因此相当于知道 $u_{0,0},u_{1,0},u_{2,0}\dots$。
通过上面的公式就可以推出来 $u_{1,1},u_{2,1},u_{3,1}\dots$，注意由于已知左边界和右边界，因此$u_{0,1}$也知道，所以第二排就可以全部推出来。
通过上面的方式可以求出区域内全部的数值解。
 
 
2.Crank-Nicholson隐式算法 
 
求解思想：也是一排一排向上推，但是这次是使用线性方程组一次性求出一排。
这里采用相同的划分方式，但是代入不同的差分方程
 
 
  
 
通过化简得到
 
 
  
 
具体推的步骤大概如下： 
  
由于已知 $u(x,0)=f(x)$，因此相当于知道 $u_{0,0},u_{1,0},u_{2,0}\dots$。
通过上面的公式就可以推出来方程$$-r{u}_{i-1,1}+(2+2r)u_{i,1}-r{u}_{i+1,1}=c(c是一个常数)$$，注意由于已知左边界和右边界，所以这个其实就转化成在一维上的差分问题，最后列出全部的方程构成方程组求解即可。
通过上面的方式可以求出区域内全部的数值解。
 
一个例子：
 
 
  
 
做划分并且代入差分方程
 取$k=0.01,h=0.1$
 $$-u_{i-1,j+1}+4u_{i,j+1}-u_{i+1,j+1}=u_{i-1,j}+u_{i+1,j}$$
进行求解(这里利用了对称性，在 $x=0.5$ 两边是对称的，将$j=0$隐去，并根据对称性将$u_6$替换成$u_4$)
 
 
  
 
(3)双曲线型
 
 
  
 
得到的差分方程为 ($r=\frac{k}{h}$注意和之前的定义不同)：
 
 
  
 
划分需要满足一定的条件$$\frac{k}{h}\le \frac{1}{c}$$
具体求解按照之前类似的方法即可。
 
(4)椭圆型
 
 
  
  
 
得到的差分方程为 (这里取$k$和$h$相等)：
 
 
  
 
求解 
  
可以采用类似之前的隐式或者显式方法求解。
可以采用迭代法求解，比如雅克比迭代，转换成下面的迭代式
 
  
 
    
  
 

 
MATLAB偏微分方程求解
 
题目：用MATLAB求解偏微分方程 主讲人： 班级　： 时间　： 基础知识预习 微分方程的ＭＡＴＬＡＢ求解包含 １：常微分方程的ＭＡＴＬＡＢ求解(上　　　　　　　　　节课已经讲过)这里不再赘述。 ２：偏微分方程的ＭＡＴＬＡＢ求解(本　　　次教学内容) 偏微分方程概念 偏微分方程(Partial Differential Equation，简称PDE)指含有未知函数及其偏导数的方程。描述自变量、未知函数及其偏导数之间的关系。 偏微分方程分为①线性偏微分方程式与②非线性偏微分方程式，常常有几个解而且涉及额外的边界条件。 求解偏微分方程的方法 求解偏微分方程的数值方法： 1. 有限元法(Finite Element Method, FEM)--- hp-FEM 2. 有限体积法(Finite Volume Method, FVM) 3. 有限差分法(Finite Difference Method, FDM)。 其它：广义有限元法(Generalized Finite Element Method, FFEM)、扩展有限元法(eXtended Finite Element Method, XFEM)、无网格有限元法(Meshfree Finite Element Method)、离散迦辽金有限元法(Discontinuous Galerkin Finite Element Method, DGFEM)等。 MATLAB解偏微分方程 MATLAB提供了两种方法解决PDE 问题：①pdepe()函数，它可以求解一般的PDEs，具有较大的通用性，但只支持命令行形式调用。 ②PDE 工具箱，可以求解特殊PDE 问题，PDEtool 有较大的局限性，比如只能求解二阶PDE 问题，并且不能解决偏微分方程组，但是它提供了GUI界面，从繁杂的编程中解脱出来了，同时还可以通过File->Save As直接生成M代码 使用pdeval()直接计算某个点的函数值？？？ 一般偏微分方程组(PDEs)的MATLAB求解 直接求解一般偏微分方程(组)，它的调用格式为sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t) 【输入参数】(1) @pdefun：是PDE 的问题描述函数，它必须换成下面的标准形式 PDE 就可以编写下面的入口函数 [c,f,s]=pdefun(x,t,u,du) m,x,t就是对应于(式1)中相关参数和自变量，du是u的一阶导数，由给定的输入变量即可表示出出c,f,s这三个函数 【输入参数】(2) @pdeic：是PDE 的初值条件，必须化为下面的形式 我们使用下面的简单的函数来描述为u0=pdeic(x) 【输入参数】(3) @pdebc：是PDE的边界条件描述函数，必须先化为下面的形式 于是边值条件可以编写下面函数描述为[pa,qa,pb,qb]=pdebc(x,t,u,du)其中a 表示下边界，b 表示下边界 【输入参数】(4) m：就是对应于(式1)中相关参数 x,t：就是对应于(式1)中自变量 【输出参数】 sol：是一个三维数组，sol(:,:,i)表示ui的解，换句话说uk对应x(i)和t(j)时的解为sol(i,j,k) 实例讲解(题目) 例： 实例讲解(解法) 【解】第一步根据(1)对照给出的偏微分方程，则原方程可以改写为 输入参数(1’)目标PDE函数 %% 目标PDE函数 function [c,f,s]=pdefun (x,t,u,du) c=[1;1]; f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2); s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)); 输入参数(2’)初值条件 初值条件改写为 %% 初值条件函数 function u0=pdeic(x) u0=[1;0]; 输入参数(3’)边界条件 边界条件改写为 %% 边界条件函数 function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) %a表示左边界，b表示右边界 pa=[0;ua(2)];qa=[1;0]; pb=[ub(1)-1;0];qb=[0;1]; (4')主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0; sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); figure('numbertitle','off','name','PDE Demo——by Matlabsky')%创建个窗口，窗口名字是name后边的名字'NumberTitle','off'是关掉默认显示名字。 subplot(21

 
matlab使用杂谈4-偏微分方程求解之pdede函数使用
 
偏微分方程
求解偏微分方程的数值方法
Matlab解偏微分方程
pdepe()函数
pdepe函数使用示例
PDE方程求解格式
PDE方程初始条件格式
PDE边界条件格式
Matlab代码
   
求解偏微分方程
总结
 

 
偏微分方程
 
偏微分方程（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) 对应 t0 和 tf 。
 
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函数的用法，说实话有些复杂，我也是看了挺久才理解，但由于没有真正应用在建模实例上，现在也不知道对他的掌握情况是咋样，之后建模过程中上传一下我的应用实例给大家参考一下，嘿嘿。