文件名称:
偏微分方程数值解的matlab实现.pdf
开发工具:
文件大小: 11mb
下载次数: 0
上传时间: 2019-06-29
详细说明:偏微分方程数值解的MATLAB实现,提供了求解一维偏微分方程的函数和求解二维偏微分方程的工具箱14.13求解一维偏微分方程
下面结合一个简单的实例介绍一维PDE的求解。
【例14-1】求解下面的PDE问题。
式中,0≤x≤1,t≥0。1=0时,解满足初始条件:
x, 0)=sin
x=0和x=1时,解满足下面的边界条件:
a(0,)=0
re-+--(,)=0
按照下面的步骤求解此方程
1.重写PDE
按照方程(14-1)的形式重写PDE,即
a(oou
x
+0
at
参数m=0,项《,《个一只a
cl xt,u
x,t,,
ax)ax
x,f,u,
0
O
2.编写PDE问题的程序
把PDE问题写成程序,用函数 pdexlpde表示。
function (c, f, s]=pdexlpde(x, t, u, DuDx)
C-PI
2
f= DuDx
s=0;
编写初始条件的函数
编写初始条件的函数 delic,如下所示
function u0= pdexlic(x)
10-sin(pi*x);
4.编写边界条件的函数
编写边界条件的函数 pdexlbc,如下所示
function [pl, ql, pr, qr]-pdexlbc(xl, ul, XI, ur, t)
pl=ul;
ql-0;
pr=pi exp(-t)
qr
287
程序中,pl和q对应于左侧的边界条件(=0),pr和q对应于右侧的边界条件(x1)
5.定义计算网格
点(x)为希望 pdepe函数求解的点。用矢量t和x指定点。这两个失量在求解过程中起不
同的作用。特别地,计算时间和精度与矢量x的长度密切相关,对矢量t中的值则很不敏感
本例中在[01取20个等间隔的点,在[02中取5个等间隔的时间值,形成计算网格,
将计算网格节点上的解。
0,1,20);
:是际港,=1053,20,中
X= lisp
space(0, 2, 5);
1界由的所不1=x0
6.使用PDE求解器
0=(10
本例调用 pdepe函数,m=0,使用函数 pdexlpde, pdexlic, pdexlbo和x,t定义的计算网格。
pdepe函数将数值解返回到三维数组sol中,其中,so(jk近似表示解中的第k个组分l,为
t〕和x)处的解。
m=0;
dOSE
sol=pdepe(m,pdexlpde, pdexlic, pdexlbc, x, t)
本例使用了函数句柄来传递 pdexlpde, delic和 pdexlbc函数。
7.查看计算结果
下面提取和显示第一个解组分。本例中,解只有一个组分,但是为了达到演示的目的
本例从三维数组中将它“提取”出来。在命令窗口键入下面的命令行。
U= so
surf(x, t, u)
title(方程的解")
xIahe
('距离
ylabel(时间t
生成图14-1。
方程的
0.8
06小
04
0.2
02
1.5
0.8
0.5
00
2
时间t
距离x3
图14-1方程的解
288
14.2二维偏微分方程的求解本器个中0
利用 MATLAB提供的偏微分方程数值解工具箱,可以求解二维PDE问题。求解的问题
包括椭圆型问题、抛物型问题、双曲型问题、特征值问题和非线性问题等。该工具箱使用有
限单元法进行数值求解。一
,5量形干,:主
注意:学习本章后面的内容,需要首先安装 MATLAB偏微分方程数值解工具箱。
示的”,,关我
s苏
1421有限元法
有限元法最早应用于结构力学。它是利用剖分插值把区域连续求解的微分方程离散成求
解线性代数方程组,以近似解代替精确解。按照所依据的原理不同,有限元法又可分为变分
有限元法、迦辽金有限元法和均衡有限元法等。:(130)p
在进行有限元计算以前,往往要进行一系列的预处理工作,如对几何形状或形体进行离
散化,即用比较简单的形状或形体来逼近和代替实际的形状或形体。例如,对于二维的情况,
可以用很多三角形或四边形来逼近一个有着复杂边界的实际图形;对于三维的情况,可以用
四面体或六面体的集合体来逼近实际形体。这样,可以将比较复杂的曲线或曲面问题转化为
相对简单的直线或平面间题。这个过程,即为建立几何模型的过程。图142中用三角形网格
将研究域离散化。直时就常的眼AB2O2A.
B
是
oe
量
1a2M业
界6
溶文的
香
a西
B甲其,果示显圆图142研究域的三角形网格
建立几何模型以后,需要根据要求解的问题赋予初始条件、边界条件和方程参数等。组
各方面的贡献以后,组成一个大型的矩阵方程,然后采用一定的方法求得方程的解
计算完毕以后,可以得到网格上各个节点(网格顶点)上的解,可以用列表的形式列出
这些值,但使用图形更直观。于是引出有限元计算的后处理问题。通常用于表现有限元计算
结果的图形有变形网格图、云图(色谱图)、等值线图、矢量图、网格图、表面图、流线图等。
②,
14.22椭圆型问题
可以用 assempde函数求解下面的标量椭圆型问题
-V·Vn)+m=f在?上1
或系统椭圆型问题
-V·(cVn)+a=f在2上
289
assenede函数是PDE工具箱中的一个基本函数。它使用有限元公式组合PDE问题。该
命令可以有选择地生成PDE问题的解。
该函数组合PDE问题的刚度矩阵和右端项,调用格式如下。凶e
1a·u= assempde(b2p, e,t, c, a,):通过在线性方程组中剔除 Dirichlet边界条件来组合和求解
PDE问题。对于标量的情况,用解的列矢量代表解天量a,列矢量中的值对应于p的对应
节点处的解。对于有m个节点的N维系统,的前m个值描述的第一个元素,接下来
的η个值描述n的第二个元素,依此类推。这样,u的元素就作为N个节点值块放到矢
量M中
东tA
〖KF= assempde(b, P, e, t, c, a, f):通过用刚性位移( f spring)近 Dirichlet边界条件来组
合PDE问题。K和F分别为刚度矩阵和右端项。PDE问题的有限元解为“KF出
[ K, F, B, ud]= LsSempde(b, p, e, t, c, a, f,:通过从线性方程组剔除 Dirichlet边界条件来组合
PDE问题。1=F返回非 Dirichlet点上的解。完整的PDE问题可以通过 MATLAB表达式
u=B*u+ud进行求解。补面图是替地
R· K, M, F,Q,GH, R= assempde(b, p, e, t. c, a, 1:给出PDE问题的分离表示( split representation
●u= assempde(k, m, F,Q,gh, r)-将PDE问题的分离表示转化为单一矩阵或矢量的形式,
然后通过从线性方程组中剔除 Dirichlet边界条件来进行求解。
面平的
K1F= assenede( S, M,F,Q,GH R):用很大的常数修改Dcet边界条件,从而将PDE
问题的分离表示转化为单一矩阵或矢量的形式。
[KF,B,udJ= assempde( k, M, F,Q,GH,R):从线性方程组中剔除 Dirichlet边界条件,从
而将PDE问题的分离表示转化为单一矩阵或矢量形式。
b描述PDE问题的边界条件。b可以是边界条件矩阵或边界M文件的文件名。边界条件
和边界M文件的格式分别在amb函数和 rebound函数中描述。
务水PDE间题的几何模型由网格数据Pt描述。 knitmesh函数对网格数据的意义进行了具体
系数cad可以通过多种方法给定。这些系数也可以与时间有关。在 assempde函数中
可以看到所有选项的列表。
【例142】求解L形薄膜的方程
△=1,2上为 Dirichlet边界条件v=0
最后绘图显示结果。其中1 shape和
Shapes分别为表示对象几何模型和边界
条件的M文件,是工具箱自带的文件。
005
knitmesh函数和 finemesh函数分别对几
何模型进行网格初始化和细化, pdesurf
函数绘制解的表面图。变
IP, e, t]=initmesh('lshapeg', "Hmax, 0.2)
Ip,e, t]
f
-retinemeshc
"lshapcg,, e,t);i-l
图14-3L形薄膜泊松方程的解
u=assempdec'lshapeb', p, e, t, 1, 0, 1);
pdesurf(p,t,u)
结果如图143所示。
2
(Y⊙2)
290·
1423抛物型问题
可以用 parabolic函数求解标量抛物型问题
eroustlvo rotor soi
O
d-V,(cVu)+a=f在Q2上m0
或系统PDE问题
的自带的
T172 TR5n010 200011
da-v.(cva)+a=在g上圆同曲.
该函数的调用格式如下
ul=parabolic(u0, tlist, g,; p,e, t, c, a, f, d:研究域为由p,e,t描述的网格,b为边界条件,初
值为u0。
对于标量的情况,解矩阵u1中的每一行都是由p中对应列给定的坐标上的解。u1中的每
列为tist中对应项给定的时间上的解。对于具有m个节点的N维系统,ul的前m行描述
n的第一个元素,接下来的叩行描述u的第二个元素,依此类推。这样,u的元素就作为N
块节点行放到中
b描述PDE问题的边界条件,它可以是个边界条件矩阵或边界M文件的文件名。边界
条件可能与时间t有关。边界条件矩阵和边界M文件分别在aemb函数和 rebound函数中
描述。
宝中q量中1,的长
PDE问题的几何模型由网格数据p,et描述。 nitmesh函数对网格数据的意义进行了具体
描述。
n,发,式
,个面数部计0n来
系数ca,可以通过多种方法给定。这些系数也可以与时间t有关。在 assenede函数中
可以看到所有选项的列表
回d
aol和rTol为绝对和相对容限。M应关回面
u1=parabolic(uo, tlist, K,F,B,ud,M):生成下面PDE问题的解。
中
du.pmt
坚黄具的回3(q
B"MB+K·l1=F,u=Bu,+
用格
dt
要的初值为u0。兴回是路
【例14-3】求解热传导方程。
△
容状lon
日电2)py(=[:
研究域为方形-1≤,y≤1( (square)在x+y2<0.42的区域内令u(0)=1,在其他区域内令
(O)=0。使用 Dirichlet边界条件=0( quare1)要求计算时间 unspaced0.1,20)处的解。其
中 square和 squarer分别为表示对象几何模型和边界条件的M文件,是工具箱自带的文件。
knitmesh函数和 finemesh函数分别对几何模型进行网格初始化和细化。。【11
Ip, e, t]=initmeshCsquareg);
LP, e, t]-refinemeshCsquareg, P,e, t);
uo=zeros(size(p, 2), 1);
ix=ind(sqrt(p(1,;)2+p(2,)2)-04
uo(ix)ones(size(ix))
tlist=linspace(0, 0. 1, 20);
ul-parabolic(u0, tlist, ' squarebl,p, e,t, 1, 0,1,0)
运行后显示:
29
96 successful steps
0 failed attempts
194 function evaluations
I partial derivatives
20 LU decompositions 2
=1+(
)。
193 solutions of linear systems
1424双曲型问题
可以用 hyperbolic函数求解标量双曲型问题
亲d,a2uV,(cV)+m=f在上n一1
at2
或系统双曲型问题
和会,秒中可由中确,南
d2-V·(eVn)+a=f在2上雪中
at2
该函数的调用格式如下。
·ul- hyperbolic(uO,uto, tlist, b,pet,c,a2f,山):研究对象为用p,e,t描述的网格,b为边
界条件,u0为初值,初始导数为uto。
是,
对于标量的情况,解矩阵u1中的行是由p中对应列给定的坐标上的解。u1中的每一列为
dis中对应项给定的时间上的解。对于具有np个节点的N维系统,u的前叩行描述u的第
1个元素,接下来的np行描述u的第2个元素,依此类推。这样,u的元素就作为N个节点
行块放到t中。
b描述PDE问题的边界条件,它可以是一个边界条件矩阵或边界M文件的文件名。边界
条件可能与时间t有关。边界条件矩阵和边界M文件分别描述在 assemb函数和 rebound函
数中。
PDE问题的几何模型由网格数据p,e,t描述。 knitmesh函数对网格数据的意义进行了具体
描述。
-=
系数cad可以通过多种方法给定。这些系数也可以与时间t有关。在 assenede函数中
可以看到所有选项的列表。
ato和rtol为绝对和相对容限。
·ul- hyperbolic(u0, uto, tlist, K, F,B,ud,M):生成下面PDE问题的解
BMB
d i
2
+K
=F
u= Bu,+
初值为u0和ut0。
spanair
【例14-4】求解波动方程。型
长管dsu
=△
以是白
研究域为方形-1≤x,y≤1( quare),当x=±1时,有 Dirichlet边界条件u=0;当y±1
矩形3)时,有 Neumann边界条件
(((
=0
Op
选择
0,,d
pe aPr,Ur
292
和
arctan(cos(r)
du(oy
dt
= 3sin(x)exp(cos (v)
计算时间等于0,16,13,…,29/6,5时的解。其中 square和 square3分别为表示对象几
何模型和边界条件的M文件,是工具箱自带的文件。 knItmesh函数对几何模型进行网格初
始化。m1a请
29
文量自具
M
[p, e, t]-initmesh('squareg")
x=p(1,);
体
回
y=p(2,)
u0-=atan( cos(pi/2*x));
8=
ut0=3*sin(pi*x). *exp(cos(pi*y)
,3
tlist=linspace(0, 5, 31)
()1、,0,1,A,中
uu=hyperbolic(u0, uto, tlist, 'squareb3' p, e, t, 1, 0, 0, 1)
运行后显示:
462 successful steps
70 failed attempts
1066 function evaluations
I partial derivatives
后个
156 LU decompositions
1065 solutions of linear systems
1425特征值问题
21y vow
可以用 pdeeig函数求解标量特征值问题
y·(cVn)+au=dn在上
或系统特征值问题
了=VwM,里=
V·(c⑧Vn)+an=adnn在gQ上
1,,=18
该函数的调用格式如下
wJ= pdeeig(b,p,eoa,dn):研究域的几何模型由p,e和t描述边界条件由b给定
r为二元素矢量,指示实数轴上的个区间(左侧可以是-nf。算法将该区间内的所有
特征值返回到1中
ν为特征值矩阵。对于标量的情况,ν中的列为p中对应节点处解的特征值。v中的每
列为tist中对应项给定的时间上的解。对于具有个节点的N维系统,的前np行描述
的第一个元素,接下来的m行描述v的第二个元素,依此类推。这样,卩的元素就作为N块
节点行放到v中。
0=5H0
b描述PDE问题的边界条件,它可以是一个边界条件矩阵或边界M文件的文件名。边界
条件可能与时间t有关。边界条件矩阵和边界M文件分别在amb函数和 rebound函数中
描述
PDE问题的几何模型由网格数据p,e,f描述。 knitmesh函数对网格数据的意义进行了具体
描述
系数c,a,d可以通过多种方法给定。这些系数也可以与时间t有关。在 assempde函数中
可以看到所有选项的列表。
293
·[]}= pdeeig(K,B,Mr):生成下面广义稀疏矩阵特征值问题的解。
Kui=namBui, u= Bui
参数的实部位于区间r内。
【例145】计算小于100的特征值及对应的特征模式。问题为
铁不的球的-V=220,,10的
研究域为L形薄膜。然后显示第1个和第16个特征模式。其中Mg和1 Shaped分别为
表示对象几何模型和边界条件的M文件,是工具箱自带的文件。 knitmesh函数和 finemesh
函数分别对几何模型进行网格初始化和细化, resurf函数绘制解的表面图。
Ip, e, t=initmeshClshapeg"),
(.c
Ip, e, t]-refinemesh(Ishapeg', p,e,
=
Ip,e, t]=refinemesh('"lshapeg', p, e, t);
aciD
[V,I=pdeeigC'lshapeb, p,e, t, 1, 0, 1,[-Inf 1001);
l(1)
%第1个特征值
desuriip,tv(∷1)
%第1个特征模式
gu
Ire
membranc(1, 20, 9, 9)
% MATLAB函数
J i2ir cap
figure
l(16)
%第16个特征值
50
desurf(p, t, v(: 16
pe
)
%第16个特征模式
运行结果为:
snguilnoeikoine
Basis=10, Timc= 9.45, New conv eig- 2
Basis=13, Time= 11.37, New conv eig= 2
Basis= 16, Time=. 13. 24, New conv eig- 3
国己.
Basis= 19, Time= 15.27, New conv eig= 4
soeg4读
Basis= 22, Time- 17.47, New conv eig- 5
Basis= 25, Time= 19.72, New conv eig 6
Basis= 28, Time= 22.08, New conv eig=7
Basis=31, Time- 24.61, New conv eig=- 7
Basis= 34, Time= 27.25, New conv eig= 9
件,BB3
Basis=37, Time= 29.94, New conv eig= 9
40, Time- 32.79, New conv eig=11 a
Bas43,Tim=35.87, New conv eig=1l量兴
Basis-46, Time- 39.16, New conv eig=16
ihR Basis-49, Time- 42.24, New conv eig=17-120628277
Basis= 52, Time= 45.54, New conv eig= 20
End of sweep: Basis=-55, Time- 49.22, New conv eig= z 8 3 00v700se
Basis=55, Time= 49. 16, New conv ci
Baas> 38, lme- 58.44. New conv eig= 0
End of sweep: Basis-, Time- 58. 50, New conv eig- o
ns
96699
ans
933166
并生成3个特征值模式图,如图144、图145和图146所示。
294·
(系统自动生成,下载前可以参看下载内容)
下载文件列表
相关说明
- 本站资源为会员上传分享交流与学习,如有侵犯您的权益,请联系我们删除.
- 本站是交换下载平台,提供交流渠道,下载内容来自于网络,除下载问题外,其它问题请自行百度。
- 本站已设置防盗链,请勿用迅雷、QQ旋风等多线程下载软件下载资源,下载后用WinRAR最新版进行解压.
- 如果您发现内容无法下载,请稍后再次尝试;或者到消费记录里找到下载记录反馈给我们.
- 下载后发现下载的内容跟说明不相乎,请到消费记录里找到下载记录反馈给我们,经确认后退回积分.
- 如下载前有疑问,可以通过点击"提供者"的名字,查看对方的联系方式,联系对方咨询.