查看: 2327|回复: 0

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

[复制链接]

字体大小: 正常放大

浅夏110

542 主题	15 听众	1万积分

TA的每日心情

	开心 2020-11-14 17:15

签到天数: 74 天

[LV.6]常住居民II

群组: 2019美赛冲刺课程

群组: 站长地区赛培训

群组: 2019考研数学桃子老师

群组: 2018教师培训（呼伦贝

群组: 2019考研数学站长系列

电梯直达

1^#

发表于 2020-6-10 10:32 |只看该作者 |倒序浏览

|招呼Ta 关注Ta |邮箱已经成功绑定

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

————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711854

zan