您好,欢迎光临本网站![请登录][注册会员]  
文件名称: 偏微分方程数值解的matlab实现.pdf
  所属分类: 讲义
  开发工具:
  文件大小: 11mb
  下载次数: 0
  上传时间: 2019-06-29
  提 供 者: chung******
 详细说明:偏微分方程数值解的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最新版进行解压.
  • 如果您发现内容无法下载,请稍后再次尝试;或者到消费记录里找到下载记录反馈给我们.
  • 下载后发现下载的内容跟说明不相乎,请到消费记录里找到下载记录反馈给我们,经确认后退回积分.
  • 如下载前有疑问,可以通过点击"提供者"的名字,查看对方的联系方式,联系对方咨询.
 输入关键字,在本站1000多万海量源码库中尽情搜索: