偏微分方程边值问题的数值解法
- 格式:doc
- 大小:743.00 KB
- 文档页数:14
偏微分积分方程的周期边值问题偏微分方程周期边值问题可分为两大方面:解析解法和数值解法。
其中只有很少一部分偏微分方程能求得解析解,所以实际应用中,多求数值解。
数值解法又可以分为最常见的有三种:差分法、有限体积法、有限元法。
其中,差分法是最普遍最通用的方法。
(1)直接积分的方法当场源与场域的形状比较简单,位函数仅是一个坐标的函数,所求解的泊松方程和拉普拉斯方程为二阶的常微分方程,可采用直接积分的方法求解。
(2)分离变量法当位函数是两个或三个坐标的函数,但场域的边界与所选择的坐标系中坐标面相吻合时,常采用分离变量法。
先将待求的位函数如分离成两个或三个各自仅含一个坐标的函数的乘积,组成把它代入场方程,借助“分离常数”可得每一变量的常微分方程,并分别求得其通解,然后组合成偏微分方程的通解,再由边界条件决定分离常数与积分常数,得到位函数的解。
(3)复位函数法能用来处理场域边界的几何形状比较复杂的问题,如椭圆、多角形截面的电极、偏芯电缆、电机气隙及波导等电磁场问题。
它是利用复变函数中解析函数的实部与虚部在复平面的某一区域内都满足拉普拉斯方程的特性,当所求解的二维拉普拉斯场域边界与某一解析函数的图形一致时,则此解析函数的实部或虚部就是所求位函数的解。
(4)保角变换法是利用解析函数的保角变换特性,将平面上的边界形状较复杂的场域,以对应的几何方式变换到边界形状较为简单的平面,求解后再反变换到平面,获得原问题的解。
(5)镜像法是边值问题中一种间接求解法,其理论依据是场的惟一性定理。
镜像法的基本原理是在求解的场域之外用虚设的镜像电荷或镜像电流等效替代边界上复杂分布的感应电荷、极化电荷或磁化电流等,只要求解区在等效前后满足同一边值问题,则其解答是惟一的。
应用镜像法的关键是找到镜像电荷或电流的位置与大小。
二阶线性与非线性偏微分方程始终是重要的研究对象。
这类方程通常划分成椭圆型、双曲型与抛物型三类,围绕这三类方程所建立和讨论的基本问题是各种边值问题、初值问题与混合问题之解的存在性、唯一性、稳定性及渐近性等性质以及求解方法。
第十章 偏微分方程数值解法偏微分方程问题,其求解十分困难。
除少数特殊情况外,绝大多数情况均难以求出精确解。
因此,近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
§1 差分方法的基本概念1.1 几类偏微分方程的定解问题椭圆型方程:其最典型、最简单的形式是泊松(Poisson )方程),(2222y x f yu x u u =∂∂+∂∂=∆ 特别地,当0),(≡y x f 时,即为拉普拉斯(Laplace )方程,又称为调和方程2222=∂∂+∂∂=∆yux u u Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(),(),(),(),(2222y x y x u y x y x f y ux u y x ϕ 其中Ω为以Γ为边界的有界区域,Γ为分段光滑曲线,ΓΩ称为定解区域,),(y x f ,),(y x ϕ分别为Ω,Γ上的已知连续函数。
第二类和第三类边界条件可统一表示为),(),(y x u u y x ϕα=⎪⎪⎭⎫ ⎝⎛+∂∂Γ∈n 其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件, 0≠α时为第三类边界条件。
抛物型方程:其最简单的形式为一维热传导方程220(0)u ua a t x∂∂-=>∂∂ 方程可以有两种不同类型的定解问题:初值问题⎪⎩⎪⎨⎧+∞<<∞-=+∞<<-∞>=∂∂-∂∂x x x u x t x u a tu )()0,(,0022ϕ初边值问题221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x lu t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩其中)(x ϕ,)(1t g ,)(2t g 为已知函数,且满足连接条件)0()(),0()0(21g l g ==ϕϕ边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类边界条件。
第6章:偏微分方程数值解法6.1对流方程【6.1.1】考虑边值问题, 01,0(0,)0,(1,)1(,0)t x x u au x t u t u t u x x=<<>ìï==íï=î如果取:2/7x D =,(0.5),1,2,3j x j x j =-D =,8/49t D =,k t k t=D 求出111123,,u u u 【解】采用Crank-Nicolson 方法()11111111211222k k k k k k k k j j j j j j j j u u u u u u u u t x ++++-+-+éù-=-++-+ëûD D 11111113k k k k k kj j j j j j u u u u u u +++-+-+-+-=-+由边界条件:(0,)0x u t =,取100k ku u x-=D ,10,0,1,k ku u k ==L (1,)1u t =,41ku =-1 1 0 0 - (1+2s) -s 0 0 -s (1+2s) -s 0 -s (1+2s) -s 0 s L L L L 101210 0 0 0 (1-2s) s 0 0 s (1-2s) s 0 s ( 1 k n n u u s u u u +-éùéùêúêúêúêúêúêú=êúêúêúêúêúêúêúêúêúëûëûL L L L L 01211-2s) s 0 1 1kn u u u u -éùéùêúêúêúêúêúêúêúêúêúêúêúêúêúêúêúëûëûL 由初始条件:021(72j j u x j ==-,1,2,3j =,212()t s x D ==D -1 1 0 0 0-1 3 -1 0 0 0 -1 3 -1 0 -1 3 -1 0 1012340 0 0 0 01 -1 1 0 00 1 -1 1 0 1 -1 1 1 u u u u u éùéùêúêúêúêúêúêú=êúêúêúêúêúêúëûëû00123 0 1 1u u u u éùéùêúêúêúêúêúêúêúêúêúêúêúêúëûëû000117u u ==,0237u =,0357u =1112327u u -=,111000123123337u u u u u u -+-=-+=,11100234235317u u u u u -+-=-+=114591u =125191u =,136991u =6.2抛物形方程【6.2.1】分别用下面方法求定解问题22(,0)4(1)(0,)(1,)0u u t x u x x x u t u t ì¶¶=ï¶¶ïï=-íï==ïïî01,0x t <<>(1)取0.2x D =,1/6l =用显式格式计算1i u ;(2)取0.2,0.01x t D =D =用隐式格式计算两个时间步。
偏微分方程中的边值问题偏微分方程(Partial Differential Equations,简称PDE)是数学中重要的研究对象,它描述了物理、工程、生物等学科中许多实际问题的数学模型。
在解决偏微分方程的过程中,边值问题(Boundary Value Problem,简称BVP)扮演着重要的角色。
本文将探讨在偏微分方程中的边值问题及其解决方法。
一、边值问题的定义在求解偏微分方程时,我们通常需要给定一些额外的条件,这些条件被称为边界条件或边值条件。
边值问题是指在解偏微分方程时,除了给出方程本身外,还给出了在某些边界上的条件限制。
通常边界包括定解区域的整个边界以及初始时刻的条件。
二、常见类型的边值问题1. 狄利克雷边值问题狄利克雷边值问题是指在求解偏微分方程时,给定了方程在边界上的函数值。
具体而言,对于一个定义在定解区域Ω上的偏微分方程,狄利克雷边值问题给定了方程在Ω的边界∂Ω上的值,即f(x)=g(x),其中f(x)是方程的解,g(x)是边界条件给定的函数。
通过求解方程和验证边界条件,可以得到满足狄利克雷边值问题的解。
2. 诺依曼边值问题诺依曼边值问题是指在求解偏微分方程时,给定了方程在边界上的法向导数。
具体而言,对于一个定义在定解区域Ω上的偏微分方程,诺依曼边值问题给定了方程在Ω的边界∂Ω上法向导数的值,即∂f/∂n = h(x),其中f(x)是方程的解,h(x)是边界条件给定的函数。
通过求解方程和验证边界条件,可以得到满足诺依曼边值问题的解。
3. 罗宾边值问题罗宾边值问题是指在求解偏微分方程时,给定了方程在边界上的线性组合形式,即同时给定了边界上的函数值和法向导数的线性组合。
具体而言,对于一个定义在定解区域Ω上的偏微分方程,罗宾边值问题给定了方程在Ω的边界∂Ω上函数值和法向导数的线性组合,即f(x) + ∂f/∂n = k(x),其中f(x)是方程的解,k(x)是边界条件给定的函数。
通过求解方程和验证边界条件,可以得到满足罗宾边值问题的解。
高师理科学刊Journal of Science of Teachers' College and University 第41卷第1期2021年 1月Vol. 41 No.1Jan. 2021文章编号:1007-9831 (2021 ) 01-0013-04偏微分方程边值问题的数值解法研究王福章*'2,郭玄玄2,周童2(1.徐州工程学院数学与统计学院,江苏徐州221018; 2.淮北师范大学数学科学学院,安徽淮北235000)摘要:利用椭圆型偏微分方程边值问题的数学模型,对Kansa 法数值模拟进行了探讨.考虑切比雪 夫节点作为数值模拟过程中的配点,与传统Kansa 法和解析解的结果进行了对比研究.数值结果表 明,对于椭圆型偏微分方程边值问题,基于切比雪夫节点的Kansa 法数值模拟具有令人更为满意的 结果.关键词:径向基函数;偏微分方程;无网格法;Kansa 法中图分类号:O242.1 文献标识码:A doi : 10.3969/j.issn.1007-9831.2021.01.004Study on the numerical method for the boundary value problems ofpartial differential equationsWANG Fuzhang 1,2, GUO Xuanxuan 2, ZHOU Tong 2(1. School of Mathematical and Statistics, Xuzhou University of Technology, Xuzhou 221018, China;2. School of Mathematical Sciences, Huaibei Normal University, Huaibei 235000, China )Abstract : Using the mathematical model of the boundary value problems of elliptic partial differential equation , the numerical simulation of Kansa's method was discussed. Considering the Chebyshev node as the collocation point in the process of numerical simulation , it was compared with the results of traditional Kansa's method and analytical solution. The numerical results showed that the numerical simulation based on the Chebyshev node of Kansa's method has more satisfactory results for the boundary value problems of elliptic partial differential equations.Key words : radial basis function ; partial differential equation ; meshless method ; Kansa's methodKansa 法是对应于有限元法发展起来的一种区域配点型无网格法,利用Multiquadric ( MQ )函数作为径 向基函数,具有无需网格划分、原理简单、易编程和谱收敛性等优点[1-4].目前对于Kansa 法的研究大多集 中于理想条件下相关问题的数值模拟和最优参数选取.文献⑸利用MQ 径向基函数数值模拟一维和二维地下 水污染运移模型;周德亮[6-7]等用MQ 径向基函数法数值模拟二维地下水稳定流和非稳定流问题;龙玉桥冏 等深入探讨了MQ 径向基函数法对地下水稳定流问题的研究;白玉峰[9]等对MQ 函数拟插值在数值积分与微分 中的应用进行了实例研究;张继红[10]等提出了一种新的MQ 径向基函数插值算法.尽管Kansa 法在实际问题 的数值模拟中有较多的应用,然而针对Kansa 法精度的改进工作尚未见相关报道.本文考虑椭圆型偏微分方程边值问题的数学模型,用切比雪夫节点作为数值模拟过程中的配点,与传 统Kansa 法和解析解的结果进行了对比研究,分析了基于切比雪夫节点的Kansa 法对椭圆型偏微分方程边 值问题数值结果的影响.收稿日期:2020-08-27作者简介:王福章( 1984-),男,山东日照人,副教授,博士,从事微分方程数值解法研究.E-mail: ***********************14高 师 理 科 学 刊第 41 卷1 Kansa 法的基本思想为了说明Kansa 法的基本思想,考虑平面区域W u R 2上的椭圆型偏微分方程边值问题Lu(x , y ) = f (x , y ), (x , y ) wW u(x, y) = u (x , y), (x, y) wG 0"(x , y ) = q (x , y ), (x , y ) w G n(1)(2)d n3)其中:L = D =d 2 a 2为拉普拉斯型微分算子;f (x , y )为已知源函数;r D U r N = dw , r D I G = f ;msa 法中所用的MQ 径向基函数,u (x , y )为r D 上的第一类边界条件;q (x , y )为G N 上的第二类边界条件.在Kansa 法中,用U (X ) = £ a f(^X - X 』)来近似边值问题(1)~(3)的解u (X ),其中:{X ,是区 域W 中的n 个不同的节点;{勺}为待定系数;f (|| X - X j ||2 r j =| |X - X j || ^( x - x j ^ +(y - y j )为二维区域 W u R 2 中点 X = (x , y )和 X j =( x }., y j 之间的欧氏距离, c 为MQ 径向基函数中的参数.将式U (X ) = f a f(^X - X j|| )代入控制方程(1)和边界条件(2) - (3 )中,考虑N 个配点(M 个 j =1 '2 /内点,N d -M 个第一类边界点,N -N d 个第二类边界点)可得N t a L f X ,, X j ) = f (X t ), X j =( x 」,y Jw W , j = 1,…, i =1M(4)N t af X ,, Xj ) = U (X^), Xj =( X j , y JwG ,j = M +1,-i -1•, N d (5)N df (x ,, X 」/ 、 / 、t a —- - q (X,), Xj =( x , y JwG N , j - N d +1,i =1 a n …, N (6)方程(4)~(6)可改写成矩阵的形式A a =b (7)L f (X,, X 7) j -1,…,Ma 、其中:A = ( a tJ )为N x N 系数矩阵,% =■f ( X ,, X 7) j - M +1,…,Nd a -为待求系数向量;af ( x ,, x 7)、a N )—4;_- j - N d+1,…,N I a n D 丿=N D +1'b = (f 1,…,f M , U M +1,…,U N D , q N D +1 ,…,q N )-求出待定系数后,对于待求区域内任意点处的数值解,可通过将对应点处的坐标代入式U (X )= t a f (|X -X 7|2)求得,对应点处法向的数值解可由该式求法向导数得到.2非均匀配点格式在数值模拟过程中,传统配点型无网格法的配点大多采用(类似)均匀布点的格式,为了提高传统配 点型无网格法的数值模拟精度,本文提出了一种非均匀格式的切比雪夫配点法.该方法的基本思路是采用 定义在开区间(-1, 1)上的非均匀分布的切比雪夫节点,在新的配点内增加2个区间端点-1和1,从而构成 整个闭区间[-1, 1]的计算节点,其数学表达式为X o =1, X ” =-1, X j = cos f + ] 1 < j < n -18)第1期王福章,等:偏微分方程边值问题的数值解法研究15利用数学表达式(8),通过MATLAB或者其它程序语言可以编程得到切比雪夫节点图,与传统的均匀布点图相比,切比雪夫节点在所考虑物理区域的边界,尤其是角点处分布稠密,而越靠近物理区域中间部分节点越稀疏.3算例分析为了说明基于切比雪夫节点的配点型无网格法的精度,考虑正方形区域^=[-1,1]x[-1,1]上的椭圆型偏微分方程边值问题x,y)=0,(x,y)e W(9)u(x,y)=x+y,(x,y)e dW(10)为了实施方便,仅考虑第一类边界条件,通过理论推导,可以得到该问题的理论解为u(x,y)=x+y.为了与基于均匀配点的传统无网格Kansa方法进行比较,在正方形区域每个边上取21个节点,整个区域上的总配点数N=441,对于非均匀布点的切比雪夫节点,在正方形区域选取的总配点数N=437,2种计算格式选取的均匀分布计算节点数均为40 401个.关于径向基函数中的参数c有较多的研究,这不是本文研究的重点,因此本文仅考虑固定参数的情况.当参数c=0.3时,给出传统Kansa方法所得到的误差图和基于切比雪夫配点的方法所得到的误差图(见图1).a传统Kansa方法b切比雪夫配点法图1传统Kansa方法和基于切比雪夫配点法所得到的误差由图1a可以看出,传统Kansa方法在点(-1,-1)和(1,1)附近误差较大,而其它点处的误差非常小,所有计算点处的平均相对误差为3.20X10-4.由图1b可以看出,基于切比雪夫配点的方法在点(-1,-1)和(1,1)附近误差仍然较大,但与传统Kansa方法相比提高了一个精度,所有计算点处的平均相对误差为1.65x10-5,该结果与传统Kansa方法相比也提高了一个精度.算例分析表明,传统的MQ径向基函数法在有些边界节点附近数值模拟精度较低,然而基于切比雪夫节点的径向基函数法在边界处的误差有明显减小,并且可以改进整体区域上的计算结果精度.4结语将切比雪夫节点和径向基函数结合用于数值模拟偏微分方程边值问题可以得到令人满意的结果.与传统的MQ径向基函数法相比,基于切比雪夫节点的径向基函数法优越性在于计算结果精度高,在边界处的误差有明显减小.尽管本文仅考虑将切比雪夫节点与MQ径向基函数相结合,但该方法可以非常容易地推广到其它同类方法中[,,].(下转第20页)20高师理科学刊第41卷证明假设F=f m(f(k))"-j只有有限多个零点,则N I r,F j<N I r,F\=S(r,f).由式(18)可推得T(r,f(k))=S(r,f),则严是有理函数,从而f(z)也是有理函数,这与条件中f(z)是超越亚纯函数矛盾.故F有无穷多个零点.或者说F=f m(严)"-j取每一个非零有穷复数无穷多次.证毕. 3结语亚纯函数结合其导数的Picard例外值是Nevanlinna值分布理论中的一个活跃且富有成果的课题.常见的推广是用f(z)的微分单项式M[f]=f n0(/')"•••(f⑷)"k替换f m,或用更一般的微分多项式y=(/,)"•••(f(k))"代替(f(k))".近年来,随着亚纯函数Nevanlinna值分布理论差分模拟的建立,Hayman 问题中的一些经典结果也都被相应的做了差分模拟.如何把本文定理中复域微分的研究成果模拟到差分上,是一个后续值得研究的课题.参考文献:[1]Hayman W K.Meromorphic Functions[M].Oxford:Oxford University Press,1964[2]Hayman W K.Research problems function and their derivatives[J].Ann of Nath,1959,70(2):9-42[3]Chen H H,Fang M L.On the value distribution of f f[J].Science in China,1994,25(2):789-798[4]Bergweiler W,Eremenko A.On the singularities of the inverse to a meromorphic function of finite order[J].Revista MathematicalIberoamericana,1995,11(2):355-374[5]TSE C K,Yang C C.On the value distribution of f(f<k))n[J].Kodai Math,1994,17:163-169[6]Qiu G D,Lin X Q.On the Picard value of f(f<k))n[J].Journal of Mathematical Study,1997,30(3):249-252[7]刘希普,张占亮.Tse和Yang的一个结果的改进[J].山东师范大学学报:自然科学版,2000,15(4):395-397[8]金瑾.亚纯函数的一个不等式的证明及其应用[J].曲靖师范学院学报,2010,29(6):13-17[9]金瑾.关于亚纯函数j(z)厂⑵严(z)的值分布[J].纯粹数学与应用数学,2012,28(6):711-718[10]徐俊峰.一类与小函数有关的微分多项式不等式估计[J].纯粹数学与应用数学,2015,31(5):441-448[11]戴瑞芳,徐俊峰.一类微分多项式的零点的不等式估计[J].五邑大学学报,2017,31(3):1-7[12]徐俊峰,叶水超.关于微分多项式f'(f3)"-1的零点[J].五邑大学学报,2018,32(2):1-7[13]Hayman W K,Miles J.On the growth of a meromorphic function and its derivatives[J].Complex Variables,1989,12:245-248 (上接第15页)参考文献:[1]Ling L,Trummer M R.Multiquadric collocation method with integral formulation for boundary layer problems[J].Comput MathAppl,2004,48:927-941[2]乔远阳,吴技莲,冯新龙.MQ径向基函数的理论、方法及应用[J].新疆大学学报:自然科学版,2015(4):5-13[3]齐静.径向基函数插值若干问题研究[D].重庆:重庆师范大学,2016[4]杜珊.基于MQ径向基函数的逼近性能研究[D].银川:宁夏大学,2018[5]Li J,Chen Y,Pepper D.Radial basis function method for1-D and2-D groundwater contaminant transport modeling[J].ComputMech,2003,32:10-15[6]周德亮,张静.稳定渗流计算的Multiquadric方法[J].辽宁师范大学学报:自然科学版,2004,27(4):399-402[7]周德亮,程晓生.用径向基函数法计算不稳定渗流达西速度[J].辽宁师范大学学报:自然科学版,2005,28(4):396-398[8]龙玉桥,李伟,李砚阁,等.MQ点插值法在地下水稳定流计算中的应用[J].水利学报,2011(42):572-578[9]白玉峰,李艳.关于MQ函数拟插值的研究[J].内蒙古民族大学学报:自然科学版,2010 ,25(6):17-19,22[10]张继红,王瑞林.一种新的MQ径向基函数拟插值格式[J].大连交通大学学报,2017(5):122-124[11]郭玄玄,王福章.圆锥型径向基函数的配点法及其应用[J].海南热带海洋学院学报,2020,27(2):21-24。
一维抛物型偏微分方程初边值问题求解摘要:一、引言二、一维抛物型偏微分方程初边值问题概述三、求解方法四、数值模拟与分析五、结论正文:一、引言一维抛物型偏微分方程在数学和物理等领域有着广泛的应用,比如热传导方程、波动方程等。
对于这种方程的初边值问题,人们进行了大量的研究,提出了多种求解方法。
本文将对这些方法进行综述和分析。
二、一维抛物型偏微分方程初边值问题概述一维抛物型偏微分方程形式为:$$frac{partial^2 u}{partial t^2} = c^2 frac{partial^2 u}{partial x^2}$$其中,$u(x,t)$ 是未知函数,$c$ 是常数。
初边值问题要求解该方程,并满足以下条件:1.$u(x,0) = f(x)$,即$t=0$ 时的函数值已知。
2.$frac{partial u}{partial t}(x,0) = g(x)$,即$t=0$ 时的导数值已知。
三、求解方法针对一维抛物型偏微分方程的初边值问题,目前主要有以下几种求解方法:1.分离变量法:适用于$c=1$ 的情况。
该方法将方程分解为两个独立的一阶线性微分方程,可以求得解析解。
2.矩方法:适用于$ceq 1$ 的情况。
该方法将方程转化为关于矩的递推关系式,可以求得数值解。
3.有限差分法:将方程离散化,通过差分方程求解。
该方法可以得到数值解,但可能会出现数值稳定性问题。
4.有限元法:将方程转化为有限个单元的积分方程,通过插值函数求解。
该方法可以得到较高质量的数值解,但计算复杂度较高。
四、数值模拟与分析为了比较不同方法的求解效果,我们取一维抛物型偏微分方程的一个具体例子,采用以上方法进行数值模拟。
通过对比分析,我们可以得出以下结论:1.分离变量法适用于$c=1$ 的情况,可以得到解析解,但求解范围有限。
2.矩方法对于$ceq 1$ 的情况有较好的适用性,可以得到数值解,但计算复杂度较高。
3.有限差分法易出现数值稳定性问题,求解精度较低。
偏微分方程数值解法偏微分方程(Partial Differential Equations,简称PDE)是数学中重要的研究对象,其在物理学、工程学、经济学等领域有广泛的应用。
然而,对于大多数偏微分方程而言,很难通过解析方法得到精确解,因此需要借助数值解法来求解。
本文将介绍几种常见的偏微分方程数值解法。
一、有限差分法(Finite Difference Method)有限差分法是一种常见且直观的偏微分方程数值解法。
其基本思想是将偏微分方程中的导数通过差分近似来表示,然后通过离散化的方式转化为代数方程组进行求解。
对于一维偏微分方程,可以通过将空间坐标离散化成一系列有限的格点,并使用中心差分格式来近似原方程中的导数项。
然后,将时间坐标离散化,利用差分格式逐步计算每个时间步的解。
最后,通过迭代计算所有时间步,可以得到整个时间域上的解。
对于二维或高维的偏微分方程,可以将空间坐标进行多重离散化,利用多维的中心差分格式进行近似,然后通过迭代计算得到整个空间域上的解。
二、有限元法(Finite Element Method)有限元法是另一种重要的偏微分方程数值解法。
其基本思想是将求解区域分割成有限数量的子区域(单元),然后通过求解子区域上的局部问题来逼近整个求解区域上的解。
在有限元法中,首先选择适当的形状函数,在每个单元上构建近似函数空间。
然后,通过构建变分问题,将原偏微分方程转化为一系列代数方程。
最后,通过求解这些代数方程,可以得到整个求解区域上的解。
有限元法适用于各种复杂的边界条件和几何构型,因此在实际工程问题中被广泛应用。
三、谱方法(Spectral Methods)谱方法是一种基于特定基函数(如切比雪夫多项式、勒让德多项式等)展开解的偏微分方程数值解法。
与有限差分法和有限元法不同,谱方法在整个求解区域上都具有高精度和快速收敛的特性。
在谱方法中,通过选择适当的基函数,并利用其正交性质,可以将解在整个求解区域上展开为基函数系数的线性组合。
偏微分方程的数值解法偏微分方程(Partial Differential Equation,PDE)是描述物理、化学、工程学等许多科学领域中变化的方程。
由于PDE的求解通常是困难的,因此需要使用数值方法。
本文将介绍偏微分方程的数值解法。
一般来说,求解PDE需要求得其解析解。
然而,对于复杂的PDE,往往不存在解析解,因此需要使用数值解法求解。
数值解法可以分为两类:有限差分法和有限元法。
有限差分法是将计算区域分成网格,利用差分公式将PDE转化为离散方程组,然后使用解线性方程组的方法求解。
有限元法则是将计算区域分成有限数量的单元,每个单元内使用多项式函数逼近PDE的解,在单元之间匹配边界条件,得到整个区域上的逼近解。
首先讨论有限差分法。
常见的差分公式包括前向差分、后向差分、中心差分等。
以一维热传导方程为例,其偏微分方程形式为:$$ \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2} $$其中,$u(x,t)$表示物理量在时刻$t$和位置$x$处的值。
将其离散化,可得到:$$ \frac{u(x_i,t_{j+1})-u(x_i,t_j)}{\Delta t}=\frac{u(x_{i+1},t_j)-2u(x_i,t_j)+u(x_{i-1},t_j)}{\Delta x^2} $$其中,$x_i=i\Delta x$,$t_j=j\Delta t$,$\Delta x$和$\Delta t$分别表示$x$和$t$上的网格大小。
该差分方程可以通过简单的代数操作化为:$$ u_{i,j+1}=u_{i,j}+\frac{\Delta t}{\Delta x^2}(u_{i+1,j}-2u_{i,j}+u_{i-1,j}) $$其中,$u_{i,j}$表示在网格点$(x_i,t_j)$处的数值解。
由于差分方程中一阶导数的差分公式只具有一阶精度,因此需要使用两个网格点来逼近一阶导数。
偏微分方程中的边值问题解析与数值求解偏微分方程是数学中的一个重要分支,它描述了自然界中的许多现象和过程。
在实际问题中,我们通常需要求解偏微分方程的边值问题,即在给定边界条件下找到满足方程的解。
本文将探讨偏微分方程中的边值问题的解析与数值求解方法。
1. 解析方法解析方法是指通过数学分析的手段,直接求解偏微分方程的边值问题。
这种方法通常需要利用数学工具和技巧,如分离变量法、特征线法、格林函数等。
以一维热传导方程为例,假设有一根长为L的金属棒,两端分别与温度为T1和T2的热源接触。
我们需要求解该金属棒上的温度分布。
通过分离变量法,可以将该问题转化为一系列常微分方程,进而得到温度分布的解析解。
解析方法的优点是能够给出问题的精确解,从而提供了对问题本质的深入理解。
然而,解析方法通常只适用于简单的边值问题,对于复杂的问题往往难以求解。
此外,解析解往往只存在于理想化的情况下,现实问题中的边界条件往往是复杂和不确定的,这使得解析方法的应用受到限制。
2. 数值方法数值方法是指通过数值计算的手段,近似求解偏微分方程的边值问题。
这种方法通常需要将偏微分方程离散化,将连续的问题转化为离散的问题,然后利用数值计算方法求解离散问题。
常见的数值方法包括有限差分法、有限元法和谱方法等。
有限差分法是最常用的数值方法之一,它将偏微分方程中的导数用差分近似表示,从而将偏微分方程转化为一个线性方程组,进而求解出近似解。
有限元法则是将求解区域划分为若干个小区域,然后在每个小区域内构造一个适当的试验函数,通过求解试验函数的系数来得到近似解。
谱方法则是利用傅里叶级数展开,将偏微分方程转化为一个无穷维的代数方程,通过截断级数求解出近似解。
数值方法的优点是适用范围广,可以求解各种复杂的边值问题。
同时,数值方法还可以通过增加计算精度和网格分辨率来提高计算结果的精确度。
然而,数值方法也存在一些问题,如舍入误差、稳定性问题和收敛性问题等,需要仔细处理。
§4 偏微分方程的数值解法一、 差分法差分法是常用的一种数值解法.它是在微分方程中用差商代替偏导数,得到相应的差分方程,通过解差分方程得到微分方程解的近似值. 1. 网格与差商在平面 (x ,y )上的一以S 为边界的有界区域D 上考虑定解问题.为了用差分法求解,分别作平行于x 轴和y 轴的直线族.⎩⎨⎧====jhy y ihx x i i (i ,j =0,±1,±2,…,±n ) 作成一个正方形网格,这里h 为事先指定的正数,称为步长;网格的交点称为节点,简记为(i ,j ).取一些与边界S 接近的网格节点,用它们连成折线S h ,S h 所围成的区域记作D h .称D h 内的节点为内节点,位于S h 上的节点称为边界节点(图14.7).下面都在网格D h + S h 上考虑问题:寻求各个节点上解的近似值.在边界节点上取与它最接近的边界点上的边值作为解的近似值,而在内节点上,用以下的差商代替偏导数:()()[]()()[]()()()[]()()()[]()()()[]y x u h y x u y h x u h y x u hy x u h y x u y x u h y x u hy u y h x u y x u y h x u h x u y x u h y x u hyu y x u y h x u h x u ,),(,,1,,2,1,,2,1,,1,,122222222++-+-+≈∂∂∂-+-+≈∂∂-+-+≈∂∂-+≈∂∂-+≈∂∂注意, 1︒ 式中的差商()()[]y x u y h x u h ,,1-+称为向后差商,而()()[]y h x u y x u h,,1--称为向前差商,()()[]y h x u y h x u h,,21--+称为中心差商.也可用向前差商或中心差商代替一阶偏导数.2︒ x 轴与y 轴也可分别采用不同的步长h ,l ,即用直线族⎩⎨⎧====jhy y ihx x j i (i,j =0, ±1, ±2 , ) 作一个矩形网格.图14.72. 椭圆型方程的差分方法[五点格式] 考虑拉普拉斯方程的第一边值问题()()⎪⎪⎩⎪⎪⎨⎧=∈=∂∂+∂∂y x u D y x y ux u S ,,02222μ式中μ(x ,y )为定义在D 的边界S 上的已知函数.采用正方形网格,记u (x i ,y j )=u ij ,在节点(i ,j )上分别用差商 u u u h u u u h i j ij i j i j ij i j -+-+-+-+11211222,,,,,代替2222,yux u ∂∂∂∂,对应的差分方程为u u u h u u u hi j ij i j i j ij i j -+-+-++-+=112112220,,,, (1) 或u u u u u ij i j i j i j i j =+++-+-+141111,,,,即任一节点(i ,j )上u ij 的值等于周围相邻节点上解的值的算术平均,这种形式的差分方程称为五点格式,在边界节点上取()()()h j i ij S j i y x u ∈=,,**μ(2) 式中(x i *,y j *)是与节点(i ,j )最接近的S 上的点.于是得到了以所有内节点上的u ij 值为未知量的若干个线性代数方程,由于每一个节点都可列出一个方程,所以未知量的个数与方程的个数都等于节点的总数,于是,可用通常的方法(如高斯消去法)解此线性代数方程组,但当步长不很大时,用高斯消去法将会遇到很大困难,可用下面介绍的其他方法求解. 若h →0时,差分方程的解收敛于微分方程的解,则称差分方程为收敛的.在计算过程中,由于进行四则运算引起舍入误差,每一步计算的舍入误差都会影响以后的计算结果,如果这种影响所产生的计算偏差可以控制,而不至于随着计算次数的增加而无限增大,则称差分方程是稳定的.[迭代法解差分方程] 在五点格式的差分方程中,任意取一组初值{u ij },只要求它们在边界节点(i ,j )上取以已知值μ(x i *,y j *),然后用逐次逼近法(也称迭代法)解五点格式:()()()()()[]() ,2,1,0411,1,,1,11=+++=+-+-+n u u u u u n j i n j i n j i n j i n ij 逐次求出{u ij (n )}.当(i+1,j ),(i -1,j ),(i ,j -1),(i ,j+1)中有一点是边界节点时,每次迭代时,都要在这一点上取最接近的边界点的值.当n →∞时,u ij (n )收敛于差分方程的解,因此n 充分大时,{u ij (n )}可作差分方程的近似解,迭代次数越多,近似解越接近差分方程的解.[用调节余数法求节点上解的近似值] 以差商代替Δu 时,用节点(i+1,j ),(i -1,j ),(i ,j+1),(i ,j -1)上u 的近似值来表示u 在节点(i ,j )的值将产生的误差,称此误差为余数R ij ,即()()()()()ij j i j i j i j i j i R y x u h y x u h y x u y h x u y h x u =--+++-++,4,,,,设在(i ,j )上给u ij 以改变量δu ij ,从上式可见R ij 将减少4δu ij ,而其余含有u (x i ,y j )的差分方程中的余数将增加δu ij ,多次调整δu ij 的值就可将余数调整到许可的有效数字的范围内,这样可获得各节点上u (x ,y )的近似值.这种方法比较简单,特别在对称区域中计算更简捷.例 求Δu =0在内节点A ,B ,C ,D 上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4(图14.8) 解 记u (A )=u A ,点A ,B ,C ,D 的余数分别为-4u A + u B + u c +5=R A u A -4 u B + u D +7=R Bu A-4 u c + u D +3=R Cu B + u c -4u D +5=R D以边界节点的边值的算术平均值作为初次近似值,即u A (0)=u B (0)=u C (0)=u D (0)=2.5则相应的余数为:R A =0, R B =2, R C = -2, R D =0最大余数为±2.先用δu C =-0.5把R C 缩减为零,u C 相应地变为2,这时R A , R D 也同时缩减(-0.5),新余数是R A =-0.5,R B =2,0=C R , R D =-0.5.类似地再变更δu B =0.5,从而 u B 变为3,则得新余数为0====D C B A R R R R .这样便可消去各节点的余数,于是u 在各节点的近似值为:u A =2.5, u B =3, u C =2, u D =2.5现将各次近似值及余数列表如下:次数调 整 值第n 次近似值及余数u A R A u B R B u C R C u D R D 0 1 2δu C = -0.5 δu B = 0.5 2.5 2.5 2.5 0 -0.5 0 2.5 2.5 3 2 2 0 2.5 2 2 -2 0 0 2.5 2.5 2.5 0 -0.5 0 结果近似值2.5322.5[解重调和方程的差分方法] 在矩形D (x 0≤x ≤x 0+a ,y 0≤y ≤y 0+a )中考虑重调和方程024*******=∂∂+∂∂∂+∂∂=yuy x u x u u ∆取步长h an=,引直线族图14.8⎩⎨⎧+=+=jh y y ihx x 00 (i , j = 0, 1, 2,, n ) 作成一个正方形网格.用差商代替偏导数()()()()()[]{()()()()[]()()()()[]}h y x u h y x u y h x u y h x u h y h x u h y h x u h y h x u h y h x u h y x u h y x u y h x u y h x u y x u 2,2,,2,2,,,,2,,,,8201,-+++-++---++-+-++++--+++-++=上式表明了以(x ,y )为中心时,u (x ,y )的函数值与周围各点函数值的关系,但对于邻近边界节点的点(x ,y ),如图14.9中的A ,就不能直接使用上式,此时将划分网格的直线族延伸,在延伸线上定出与边界距离为h 的点,称这些点为外邻边界节点,如图14.9以A 为中心时,点E ,C 为边界节点,点J ,K 为E ,C 的外邻边界节点,用下法补充定义外邻边界节点J 处函数的近似值u J ,便可应用上面的公式.1︒ 边界条件为()()()S P P x uP u SS ∈==21,μ∂∂μ 时,定义u J =u A -2μ2(E )h .2︒ 边界条件为()()()S P P xuP u SS ∈=∂∂=2221,μμ时,定义u J =2μ1(E )-u A -h 2μ2(E ). [其他与Δu 有关的网格] 1︒ 三角网格(图14.10(a ))取P 0(x ,y )为中心,它的周围6个邻近节点分别为:()()⎪⎪⎭⎫⎝⎛-+⎪⎪⎭⎫ ⎝⎛---⎪⎪⎭⎫⎝⎛+-⎪⎪⎭⎫⎝⎛+++h y h x P h y h x P y h x P h y h x P h y h x P y h x P 23,2,23,2,,23,223,2,,654321则 R u h u u u h i i +∆+∆=⎪⎭⎫⎝⎛-∑=226102161632式中u i =u (P i ), u 0=u (P 0),R 表示余项. 2︒ 六角网格(图14.10(b ))取P 0(x ,y )为中心,它的三个邻近节点分别为图14.9()⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛++h y h x P y h x P h y h x P 23,2,,23,2321则 R u u u h i i +∆=⎪⎭⎫⎝⎛-∑=0312334.图14.103︒ 极坐标系中的网格(图14.10(c ))取P 0(r ,θ)为中心,它的四个邻近节点分别为()()()()l r P h r P l r P h r P ++--θθθθ,,,,,,4321而拉普拉斯方程01122222=∂∂+∂∂+∂∂=θ∆u r r u r ru u的相应的差分方程为()()()011221110222134222312=⎪⎭⎫ ⎝⎛+--++++u l r h u u rh u u l r u u h 3. 抛物型方程的差分方法 考虑热传导方程的边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=><<=∂∂-∂∂0,,,,00,0,0,0,021222t t t b u t t u bx x x u t b x x u a t u μμϕ 将[0,b ]分为n 等份,每段长为∆x bn=.引两族平行线(图14.11)图14.11x =x i =i ∆x (i =0,1,2,, n )y =y j =j ∆t (j =0,1,2,, ∆t 取值见后)作成一个长方形的网格,记u (x i ,t j )为u ij ,节点(x i ,t j )为(i ,j ),在节点(i ,j )上分别用(),2,1,1,,2,1Δ2,Δ2,1,11,=-=+---++j n i x u u u t u u ji ij j i ij j i 代替22,xut u ∂∂∂∂,于是边值问题化为差分方程()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-===-==+----++ ,2,1,0,Δ,Δ1,,2,1,Δ,2,1,0,1,,2,10Δ2Δ21002,1,121,j t j u t j u n i x i u j n i x u u u a tu u nj j i j i ij j i ijj i μμϕ 记()22x ta ∆∆=λ,差分方程可写成 () ,2,1,1,,2,121,1,11,=-=+-+=-++j n i u u u u ji ij j i j i λλλ (1)由此可按t 增加的方向逐排求解.在第0排上u i 0的值由初值ϕ(i ∆x )确定,j +1排u i ,j +1的值可由第j 排的三点(i +1,j ),(i ,j ),(i -1,j )上的值u i +1,j , u ij ,u i -1,j 确定,而u 0,j +1,u n ,j +1已由边界条件μ1((j +1)∆t )及μ2((j +1)∆t )给定,于是可逐排计算一切节点上的u ij 值.当ϕ(x ), μ1(x )和μ2(x )充分光滑,且λ≤12时,差分方程收敛而且稳定.所以利用差分方程(1)计算时,必须使λ≤12,即()22Δ21Δx at ≤.热传导方程还可用差分方程()0Δ2Δ21,11,1,121,=+---+-++++x u u u a t u u j i j i j i ij j i 代替,此时如已知前j 排u ij 的值,为求第j +1排的u i ,j +1 必须解包含n -1个未知量u u j n j 1111,,,,+-+ 的线性代数方程组,这种差分方程称为隐式格式的差分方程,前面所提的差分方程称为显式格式差分方程.隐式格式差分方程对任意的λ都是稳定的.4. 双曲型方程的差分方法 考虑弦振动方程的第一边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=∂∂=><<=∂∂-∂∂0,,,,00),()0,(,0,0,0,02122222t t t b u t t u b x x t x u x x u t b x x u a tu μμψϕ 用矩形网格,列出对应的差分方程:()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-=∆=∆-==-==+--+--+-+ ,2,1,0,Δ,Δ1,,2,1),(,Δ,2,1,1,,2,1,0Δ2)(Δ22100102,1,1221,1,j t j u t j u n i x i t u u x i u j n i x u u u a t u u u nj j i i i j i ij j i j i ij j i μμψϕ 记ω=a tx∆∆与上段一样,利用u u n 022,和在第0排及第1排的已知数值(初始条件)u i 0 , u i 1可计算u i 2,然后用已知的u i 1 , u i 2及u u n 033,可计算u i 3,类似地可确定一切节点上的u ij 值.当ϕ(x ),ψ(x ),μ1(x )和μ2(x )充分光滑,且ω≤1时,差分方程收敛且稳定,所以要取∆∆t ax ≤1.二、 变分方法1. 自共轭边值问题将§3定义的共轭微分算子的概念推广到一般方程.设D 是n E 中的有界区域,S 为其边界,在D 上考虑2k 阶线性微分方程()x f x x uaLu km mi i i ni m m i i n n n=∂∂≡∑∑==++201111 ∂ 的齐次边值问题()r j u l Sj ,,2,10==式中f (x )是D 内的已知函数,l j u 是线性微分算子. 将 ⎰DvLud Ω分部积分k 次得()⎰∑⎰⎪⎪⎭⎫ ⎝⎛+=Ω=S j j j D S v R u R v u vLu d ~,Λd k 1 式中Λ(u ,v )是一个D 上的积分,其被积函数包含u ,v 的k 阶导数;R j 和 R j 是定义在边界S 上的两个线性微分算子.再将Λ(u ,v )分部积分k 次得()()⎰∑⎰⎪⎪⎭⎫⎝⎛-Ω=Λ=S k j j j D S u R v R v uL v u d ~d ,1***式中L*是一个2k 阶的微分算子,称为L 的共轭微分算子.若L=L*,则称L 为自共轭微分算子.从上面可推出格林公式()()⎰∑⎰=-=Ω-Skj jjjjDS u R v R v R u R v uL vLu 1***d ~~d 如从l j u |S =l j v |S =0可推出在边界S 上()∑==-kj jjjju R v R v R u R 1**0~~ 则称l j u |S =0为自共轭边界条件.如果微分算子及边界条件都是自共轭的,则称相应的边值问题为自共轭边值问题,此时有()0d ][=Ω-⎰DuLv vLu每个边值问题对应于某希尔伯特空间H (例如L 2(D ),见第九章§7)中的一个算子A ,其定义域M A 是H 中一线性稠密集合,它由足够次连续可微且满足边界条件的函数组成,在M A 上,Au 的数值与Lu 的数值相同,从而求解边值问题化为解算子方程Au f = 的问题.设A 为定义在实的希尔伯特空间H 中的某线性稠密集合M A 上的线性算子.若对于M A 的任意非零元素,,v u 成立(Au ,v )=(u ,Av )则称A 为对称算子.若对任意非零元素u 成立()0,>u Au 则称A 为正算子.如成立更强的不等式(Au ,u )≥r ||u ||2 (r>0)则称A 为正定算子.此处(u ,v )表示希尔伯特空间的内积,||u ||2=(u ,u ). 2. 变分原理与广义解定理 设A 是正定算子,u 是方程Au =f 在M A 上的解的充分必要条件是: u 使泛函F (u )=(Au ,u )-2(f ,u )取极小值.上述将边值问题化为等价的求泛函极值问题的方法称为能量法.在算子的定义域不够大时,泛函F (u )的极值问题可能无解.不过对于正定算子,可以开拓集合M A ,使在开拓了的集合上,泛函的极值问题有解.为开拓M A ,在M A 上引进新的内积[u ,v ]=(Au ,v ),定义模||u ||2=[u ,u ]=(Au ,u ),在模||u ||的意义下,补充极限元素,得到一个新的完备希尔伯特空间H 0,在H 0上,泛函F (u )仍然有意义,而泛函的极值问题有解.但必须注意,此时使泛函F (u )取极小的元素u 0不一定属于M A ,因此它不一定在原来的意义下满足方程Au=f 及边界条件.称u 0为广义解. 3. 极小化序列与里兹方法在处理变分问题中,极小化序列起着重要的作用.考虑泛函F (u )=(Au ,u )-2(f ,u )以d 表示泛函的极小值.设在希尔伯特空间中存在一列元素{u n } (n =1,2 ,),使()d u F n n =∞→lim则称{u n }为极小化序列.定理 若算子A 是正定的,则F (u )的每一个极小化序列既按H 空间的模也按H 0的模收敛于使泛函F (u )取极小的元素.这个定理不但指出利用极小化序列可求问题的解,而且提供一种近似解的求法,即把极小化序列中的每一个元素当作问题的近似解.设算子A 是正定的,构造极小化序列的里兹方法的主要步骤是:(1) 在线性集合M A 中选取H 0中完备的元素序列{ϕi } , (i =1,2 ,) 并要求对任意的n ,ϕ1,ϕ2,…,ϕn 线性无关.称这样的元素为坐标元素.(2) 令u a n k k k n==∑ϕ1 ,其中a k 为待定系数.代入泛函F (u ),得自变量a 1,a 2,…,a n 的函数()()()∑∑==-=nj jjn k j kjkj n f a A a a u F 11,,2,ϕϕϕ(3) 为使函数F (u n )取极小,必须()()n j a u F jn ,,2,10 ==∂∂,从而求出a k (k =1,2,…,n ).序列{u n }即为极小化序列,u n 可作为问题的近似解. 4. 里兹方法在特征值问题上的应用 算子方程Au -λu =0的非零解λ称为算子A 的特征值,对应的非零解u 称为λ所对应的特征函数. 对线性算子A ,若存在常数K ,使对任何M A 的元素ϕ成立(A ϕ,ϕ)≥K ||ϕ||2则称A 为下有界算子,正定算子是下有界的(此时K =0).记(A ϕ,ϕ)/||ϕ||2的下确界为d . 定理1 设A 为下有界对称算子,若存在不为零的元素ϕ0∈M A ,使()d A =200,ϕϕϕ则d 就是A 的最小特征值,ϕ0为对应的特征函数.于是求下有界对称算子的最小特征值问题化为变分问题,即在希尔伯特空间中求使泛函(A ϕ,ϕ)/||ϕ||2取极小的元素,或在||ϕ||=1的条件下求使泛函(A ϕ,ϕ)取极小的元素.定理2 设A 是下有界对称算子,λ1≤λ2≤…≤λn 是它的前n 个特征值,ϕ1,ϕ2,…,ϕn 是对应的标准正交特征函数,如果存在不为零的元素1+n ϕ,在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0, …, (ϕ,ϕn )=0下使泛函(A ϕ,ϕ)取极小,则ϕn +1是算子A 的特征函数,对应的特征值()11,++=n n A ϕϕλ就是除λ1 ,,λn 外的最小的一个特征值.于是求第n +1个特征值就化为变分问题,即在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0 ,, (ϕ,ϕn )=0 下求使泛函(A ϕ,ϕ)取极小的元素.为了利用里兹方法求特征值,在M A 中选取一列在H 0中完备的坐标元素序列{ϕi },(i =1,2 ,), 令u a n k k k n==∑ϕ1,确定a k ,使在条件 (u n ,u n )=1下,(Au n ,u n )取极小,这个问题化为求n个变元a 1,a 2,…,a n 的函数()()∑==nm k m k k m n n a a A u Au 1,,,ϕϕ在条件()()∑===nm k m k m k n n a a u u 1,1,,ϕϕ下的极值问题,一般可用拉格朗日乘数法解(见第九章§3,t ),此时()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n A A A A A A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的最小的根即为特征值的近似值,如果将上式的根按大小排列,就依次得后面的特征值的近似值,但精确度较差. 对一般算子方程Au -λBu=0如果A 为下有界对称算子,B 为正定算子,则()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n B A B A B A B A B A B A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的根就是特征值的近似值. 5. 迦辽金方法用里兹方法解数学物理问题有很多限制,最主要的限制是要求算子正定,但很多问题不一定满足这个条件,迦辽金方法弥补了这个缺陷. 迦辽金方法的主要步骤是:(1) 在M A 中选取在空间H 中完备的元素序列{ϕi } (i =1,2 ,),其中任意n 个元素线性无关,称{ϕi } (i =1,2,…)为坐标元素序列. (2) 把方程的近似解表示为u a n k k k n==∑ϕ1式中a k 是待定常数,把u n 代入方程Au=f 中的u ,两端与ϕj (j =1,2,…,n )求内积,得 a k 的n 个代数方程()()()n j f A a j n k j k k,,2,1,,1 ==∑=ϕϕϕ(3) 求出a k ,代回u n 的表达式,便得方程的近似解u n .在自共轭边值问题中,当算子是正定时,由迦辽金方法和里兹方法得到的关于a k 的代数方程组是相同的.。
求解偏微分方程的边值问题
本实验学习使用MATLAB 的图形用户命令pdetool 来求解偏微分方程的边值问题。
这个工具是用有限元方法来求解的,而且采用三角元。
我们用内个例题来说明它的用法。
一、MATLAB 支持的偏微分方程类型
考虑平面有界区域D 上的二阶椭圆型PDE 边值问题:
()c u u f α-∇∇+= (1.1)
其中 (1) , (2) a,f D c x y ⎛⎫∂∂∇=⨯ ⎪∂∂⎝⎭
是上的已知函数(3)是标量或22的函数方阵
未知函数为(,) (,)u x y x y D ∈。
它的边界条件分为三类:
(1)Direchlet 条件:
hu f = (1.2)
(2)Neumann 条件: ()n c u qu g ∇+= (1.3)
(3)混合边界条件:在边界D ∂上部分为Direchlet 条件,另外部分为Neumann 条件。
其中,,,,h r q g c 是定义在边界D ∂的已知函数,另外c 也可以是一个2*2的函数矩阵,n 是沿边界的外法线的单位向量。
在使用pdetool 时要向它提供这些已知参数。
二、例题
例题1 用pdetool 求解 22D 1 D: 10u x y u ∂⎧-∆=+≤⎪⎨=⎪⎩ (1.4)
解 :首先在MATLAB 的工作命令行中键入pdetool ,按回牟键确定,于是出现PDE Toolbox 窗口,选Genenic Scalar 模式.
( l )画区域圆
单击椭圆工具按钮,大致在(0,0)位置单击鼠标右键,拖拉鼠标到适当位置松开。
为了保证所绘制的圆是标准的单位园,在所绘园上双击,打开 Object Dialog 对话框,精确地输入圆心坐标X-center 为0 、Y-center 为0 及半径Radius 为l ,然后单击OK 按钮,这样单位画已画好.
( 2 )设置边界条件
单击工具边界模式按钮 ,图形边界变红,逐段双击边界,打开Boundary condition 对话框.输入边界条件.对于同一类型的边界,可以按Shift 键,将多个边界同时选择,统一设边界条件.本题选择Dirichlet 条件,输入h 为1 , r 为0。
,然后单击OK 按钮.也可以单击Boundary 菜单中Spocify Boundary Condition …选项,打开Boundary Condition 对话框输入边界条件.
( 3 )设置方程
单击偏微分方程按钮,打开PDE Specification 对话框,选择方程类型·本题选Ellintic (椭圆型),输入c为1 , a 为O , f 为1 ,然后单击OK 按钮.
( 4 )网格剖分
单击网格工具,或者单击Mesh 菜单中Initialize Mesh项,可进行初始网格剖分.这时在PDE Toolbox 窗口下方的状态栏内显示出初始网格的节点数和三角形单元数.本题节点数为144 个,三角形单元数为254 个(图?? )。
如果要细化网格,单击细化工具,或者单击Mesh 菜单中Refine Mesh 选项,节点数成为541 个,三角形单元数为1016 个。
( 5 )解方程
单击解方程工具,或者单击S olve菜单中Solve PDE 选项,可求得方程数值解并用彩色图形显示。
单击作图工具,或者单击Plot 菜单中Parameter…选项,出现Plot selection 对话框.从中选择于Height ( 3-D plot) ,然后单击Plot 按钮,方程的图形解如图?? 所示。
除了作定解问题解u的图形外,也可以作grad u等图形·
(6)输出网格节点的编号、单元编号以及节点坐标
单击Mesh 菜单中Show Node Labels选项,再单击网格工具,即可显示节点编号(图??) 。
若要输出节点坐标,只需单击Mesh 菜单中Export Mesh …选项,这时打开的Export对话框中的默认值为p e t,这里p、e、t 分别表示point (点)、edges(边)、triangles(三角形)数据变量,单击OK按钮,然后在MATLAB 命令行键入p,即可以显示按节点编号排列的坐标;键入e再回车则显示边界数据矩阵(7维数组);键入t按回车则显示三角形单元数据矩阵(4维数组)。
点、边、单元的部分输出为:
p =
Columns 1 through 11
-1.0000 0.0000 1.0000 0.0000 -0.7071 0.7071 0.7071 -0.7071 -0.9808 -0.9239 -0.8315
-0.0000 -1.0000 0 1.0000 -0.7071 -0.7071 0.7071 0.7071 -0.1951 -0.3827 -0.5556
e =
Columns 1 through 11
1.0000 9.0000 10.0000 11.0000 5.0000 1
2.0000 1
3.0000 1
4.0000
2.0000 15.0000 16.0000
9.0000 10.0000 11.0000 5.0000 12.0000 13.0000 14.0000 2.0000 15.0000 16.0000 17.0000
0 0.1250 0.2500 0.3750 0.5000 0.6250 0.7500 0.8750 0 0.1250 0.2500
0.1250 0.2500 0.3750 0.5000 0.6250 0.7500 0.8750 1.0000 0.1250 0.2500 0.3750
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
2.0000 2.0000 2.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0 0 0 0 0 0 0 0 0 0 0
t =
Columns 1 through 18
32 14 20 26 29 17 23 100 11 66 89 1 9 94 5 12 13 97
1 2 3 4 8 6 7 28 5 32 1 9 10 5 12 13 14 2
89 97 81 98 84 92 99 127 94 89 119 119 95 118 118 90 70 126
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
(7)输出近似解数值
单击Solve 菜单中Export Solution 。
选项,在Export 对话框中输入u 再单击OK 按钮,再在MATLAB 命令行中输入u 并回车,就会显示按节点编号排列的解u 的数值。
(8)近似解和准确解的比较
方程(1.4)的准确解为: 22
1(,)4x y u x y --= (1.5)
为了与准确解比较,单击Plot 菜单中Parameters …选项,打开Plot Selection 对话框,在Height (3-D plot )行的Property 下拉框中选User Entry,并且输入
u-(1-x.^2-y.^2)/4 ,单击Plot 按钮,就可以看到误差曲面,其数量级为410-。
练习1 用pdetool工具求解课本P128的第2、3、4这几题的解,并作出图形。
练习2 用pdetool工具求解课本P418的第2、4这两题的解(三角单单元形状不限),并作出图形。