偏微分方程的数值解(五)： 二维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解系列博文：

偏微分方程的数值解(一):定解问题 & 差分解法

偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

偏微分方程的数值解(五)： 二维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法

---

目录

1 方程类型

（i）椭圆型偏微分方程                                 （ii）抛物型偏微分方程          

（iii）双曲型偏微分方程                                （iv）特征值问题

（v）非线性椭圆偏微分方程                         （vi）方程组

2 边界条件

（i）Dirichlet 条件：​                     （ii）Neumann 条件： ​      

（iii）对于偏微分方程组，混合边界条件为

3 求解偏微分方程

例 6 求解泊松方程                                                                 例7 考虑最小表面问题

例8 求解正方形区域上的热传导方程                                  例9 求解正方形区域上的波方程

例 10  求解泊松方程                                                           习题

---

MATLAB 中的偏微分方程（PDE）工具箱是用有限元法寻求典型偏微分方程的数 值近似解，该工具箱求解偏微分方程具体步骤与用有限元方法求解偏微分方程的过程是 一致的，包括几个步骤，即几何描述、边界条件描述、偏微分方程类型选择、有限元划 分计算网格、初始化条件输入，最后给出偏微分方程的数值解（包括画图）。 下面我们讨论的方程是定义在平面上的有界区域Ω 上，区域的边界记作∂Ω 。

1 方程类型

MATLAB 工具箱可以解决下列类型的偏微分方程：

（i）椭圆型偏微分方程

其中c, a, f 和未知的u 可以是Ω 上的复值函数。

（ii）抛物型偏微分方程

其中c,a, f ,d 可以依赖于时间t 。

（iii）双曲型偏微分方程

（iv）特征值问题

其中λ 是未知的特征值，d 是Ω 上的复值函数。

（v）非线性椭圆偏微分方程

其中c, a, f 可以是u 的函数。

（vi）方程组

2 边界条件

边界条件有如下三种：

（i）Dirichlet 条件：

（ii）Neumann 条件：

这里  为区域的单位外法线，h,r, q, g 是定义在∂Ω 上的复值函数。

对于二维方程组情形，Dirichlet 边界条件为    

Neumann 边界条件为：

（iii）对于偏微分方程组，混合边界条件为

这里 μ 的计算是使得满足 Dirichlet 边界条件。

3 求解偏微分方程

例 6 求解泊松方程 

                      求解区域为单位圆盘，边界条件为在圆盘边界上u = 0。

解 它的精确解为

下面求它的数值解，编写程序如下：

%(1)问题定义
g='circleg'; %单位圆
b='circleb1'; %边界上为零条件
c=1;a=0;f=1;
%（2）产生初始的三角形网格
[p,e,t]=initmesh(g);
%（3）迭代直至得到误差允许范围内的合格解
error=[]; err=1;
while err > 0.01,
[p,e,t]=refinemesh(g,p,e,t);
u=assempde(b,p,e,t,c,a,f); %求得数值解
exact=(1-p(1,:).^2-p(2,:).^2)/4;
err=norm(u-exact',inf);
error=[error err];
end
%结果显示
subplot(2,2,1),pdemesh(p,e,t);
subplot(2,2,2),pdesurf(p,t,u)
subplot(2,2,3),pdesurf(p,t,u-exact')

---

例7 考虑最小表面问题

g='circleg';
b='circleb2';
c='1./sqrt(1+ux.^2+uy.^2)';
rtol=1e-3;
[p,e,t]=initmesh(g);
[p,e,t]=refinemesh(g,p,e,t);
u=pdenonlin(b,p,e,t,c,0,0,'Tol',rtol);
pdesurf(p,t,u) 

---

例8 求解正方形区域上的热传导方程

求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的热传导方程

边界条件为Dirichlet条件u = 0。

解 这里是抛物型方程，其中c = 1, a = 0, f = 0, d = 1。编写程序如下：

%(1)问题定义
g='squareg'; %定义正方形区域
b='squareb1'; %边界上为零条件
c=1;a=0;f=0;d=1;
%（2）产生初始的三角形网格
[p,e,t]=initmesh(g);
%(3)定义初始条件
u0=zeros(size(p,2),1);
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);
u0(ix)=1
%(4)在时间段为0到0.1的20个点上求解
nframe=20;
tlist=linspace(0,0.1,nframe);
u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);
%(5)动画图示结果
for j=1:nframe
 pdesurf(p,t,u1(:,j));
 mv(j)=getframe;
end
movie(mv,10) 

---

例9 求解正方形区域上的波方程

              求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的波方程

解 这里是双曲型方程，其中 c = 1, a = 0, f = 0, d = 1。编写程序如下：

%(1)问题定义
g='squareg'; %定义正方形区域
b='squareb3'; %定义边界
c=1;a=0;f=0;d=1;
%（2）产生初始的三角形网格
[p,e,t]=initmesh(g);
%(3)定义初始条件
x=p(1,:)';y=p(2,:)';
u0=atan(cos(pi*x));
ut0=3*sin(pi*x).*exp(cos(pi*y));
%(4)在时间段为0到5的31个点上求解
n=31;
tlist=linspace(0,5,n);
uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);
%(5)动画图示结果
for j=1:n
 pdesurf(p,t,uu(:,j));
 mv(j)=getframe;
end
movie(mv,10) 

---

例 10  求解泊松方程

求解区域为单位圆盘，边界条件为在圆盘边界上u = 0。

解 它的精确解为

下面求它的数值解，编写程序如下：

g='circleg';
b='circleb1';
c=1;a=0;f='circlef';
[p,e,t]=initmesh(g);
[p,e,t]=refinemesh(g,p,e,t);
u=assempde(b,p,e,t,c,a,f);
exact=-1/(2*pi)*log(sqrt(p(1,:).^2+p(2,:).^2));
subplot(2,2,1),pdemesh(p,e,t);
subplot(2,2,2),pdesurf(p,t,u)
subplot(2,2,3),pdesurf(p,t,u-exact') 

---

习题

1. 求二维拉普拉斯方程

2. 求初边值问题

在0 ≤ t ≤ 3范围内的数值解。

3. 求热传导方程初边值问题

---

偏微分方程的数值解系列博文：

偏微分方程的数值解(一):定解问题 & 差分解法

偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

偏微分方程的数值解(五)： 二维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法