有限差分法求解抛物型方程说明
有限差分法
两端都要给定边界条件(双程坐标) 。
9
(C) 双曲型方程:适当的边界条件和初始条件,与波动传 播的性质有关 如:一维对流方程
∂u ∂u +c =0 ∂t ∂x u (x ,0) = f (x )
解为 u (x , t ) = f (x − ct ) ,代表一个向右(c > 0 时)或向左 ( c < 0 时)传播的波形。必须在波形传来的一侧提供边界条 件(单程坐标) 。
10
不适定的例子:
utt + u xx = 0 u (x ,0) = u t (x ,0) = 0
拉普拉斯方程+非闭域边界条件,解为 u (x , t ) ≡ 0 。 然而,若定解条件为 u (x ,0) = 0, ut (x ,0) =
u (x , t ) = 1 sin nx ,解为 n
1 sinh nt sin nx n
(
)
n n um+1 = um −
cτ n n um +1 − um −1 2h
(
)
设计算到第 n 步时的累积误差
n ~n εn = 计算值um − 差分法精确解um m
反之
n ~n um = εn + um m
15
则第 n+1 步的计算值
~n ~ n cτ u n − u n ~ ~ um+1 = um − m +1 m −1 2h cτ n cτ n n n = um − um +1 − um −1 + εn − εm +1 − εn −1 m m 2h 2h n = um+1 + εn +1 m
uin +1 − uin −1 uin+1 − uin +1 − uin −1 − uin−1 −α =0 Lh u = τ h2 ατ 2 ⎛ ∂ 2u ⎞ τ 2 ⎛ ∂ 3u ⎞ Ti = Lh u − Lu (x i , t n ) = 2 ⎜ 2 ⎟ + ⎜ 3 ⎟ − L 截断误差 6 ⎜ ∂t ⎟i h ⎜ ∂t ⎟i ⎝ ⎠ ⎝ ⎠
抛物型方程
前言抛物型方程解的估计及其应用1前言数学物理方程主要指从物理学及其它各门自然科学、技术科学中所产生的偏微分方程(有时也包括积分方程、微分方程等),它们反映了有关的未知变量关于时间的导数和关于空间变量的导数之间的制约关系.连续介质力学、电磁学、量子力学等等方面的基本方程属于数学物理方程的范围.它以具有物理背景的偏微分方程(组)作为研究的主要对象.它与其他数学分支及物理、化学等自然科学和工程技术的很多领域都有着广泛的联系,因此,无论在历史上还是在今天的现实生活中,它对于推动数学理论的发展,加强理论与实际的联系,帮助人们认识世界和改造世界都起着重要大的作用.微积分产生以后,人们就开始把力学中的一些问题,归结为偏微分方程进行研究.早在18世纪初,人们已经将弦线振动问题归结为弦振动方程,并探讨了它的解法.随后,人们又陆续了解了流体的运动、弹性体的平衡和振动、热传导、电磁相互作用、原子核和电子的相互作用、化学反应过程等等自然现象的基本规律,把它们写成偏微分方程的形式,并且求出了典型问题的解答,从而能通过实践,验证这些基本规律的正确性,显示了数学物理方程对于认识自然界基本规律的重要性.有了基本规律,人们还要利用这些基本规律来研究复杂的自然现象和解决复杂的工程技术问题,这就需要求出数学物理方程中许多特定问题的解答,随着计算机的出现及计算技术的发展,即使是相当复杂的问题,也可以计算出足够精确的数值来,这对于预测自然现象的变化(如气象预报)和进行各种工程设计(如机械强度的计算)都有着很重要的作用.在研究数学物理方程的同时,人们对偏微分方程的性质也了解得越来越多、越来越深入,形成数学中的一门重要分支——偏微分方程理论.它既有悠久的历史,又不断地更新着它的对象、内容和方法.它直接联系着众多自然现象和实际问题,不断地提出或产生需要解决的新课题和新方法.它所面临的数学问题多样而复杂,不断地促进着许多相关数学分支(如泛函分析、复变函数、微分几何、计算数学等)的发展,并从它们之间引进许多有力的解决问题的工具.因此,数学物理方程又是纯粹数学的抛物型方程解的估计及其应用许多分支和自然科学各部门及工程技术等领域之间的一个重要的桥梁.2 选题背景2.1 题目类型及来源题目类型:研究论文题目来源:专题研究2.2研究目的和意义数学物理方程将数学与物理紧密地结合在一起,用精微的数学思想和方法应用于实际的物理研究中,通过物理过程建立数学模型(偏微分方程),通过求解和分析模型,对实际物理过程进一步深入理解,提出解决实际问题的途径和方法.而抛物型方程是偏微分方程中的一类,且在研究热传导、扩散等物理现象时都会遇到,具有巨大的理论价值,同时,抛物型方程的概念和性质也决定了它在工程数学,物理等方向的巨大实用价值.研究抛物型方程解的估计及其应用,有助于我们更好的理解和掌握偏微分方程理论,并在认识和了解抛物型方程广泛的使用价值的基础上,能够探索抛物型方程更广泛的应用.随着物理科学所研究的现象在广度和深度两方面的扩展,偏微分方程的应用范围更广泛.从数学自身的角度看,抛物型方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展.2.3国内外现状和发展趋势与研究的主攻方向自18世纪以来,偏微分方程理论在得到广泛应用的同时,也得到不断发展和完善,内容也越来越丰富,它直接联系着众多自然现象和实际问题,不断地提出或产生需要解决的新课题和新方法.它所面临的数学问题多样而复杂,不断地促进着许多相关数学分支的发展,如泛函分析、复变函数、微分几何、计算数学等,并从它们之中引进许多有力的解决问题的方法.关于抛物型方程的解,有经典解、广义解、数值解等方面的研究.经典解在求解区域中具有方程中所出现的连续偏导数,并按通常意义满足方程与定解条件,也就是热传导方程的一些知识说,将解代入方程及定解条件后即可使其化为恒等式.求经典解的方法有分离变量法、Fourier变换法.经典解容易理解,且应用方便,但实际求解抛物型方程的定解问题时,往往不一定能得到经典解.于是就提出了广义解[1]的理论,即先寻求一个正则性较低的函数(广义解),它按较弱的意义满足方程与定解条件,然后再进一步证明这个函数实际上就是原来问题的经典解.广义解可按逼近过程来定义,也可按分部积分的方法来定义.由于抛物型方程的广义解是在“较差”的函数类中寻求相应定解问题按拓广意义下的解,因而又提出了广义函数的概念、性质、应用等方面的研究.在实际求解抛物型方程的定解条件时,除了一些特殊的情况下可以方便地求得其精确解外,在一般的情况下,当方程或定解条件具有比较复杂的形式,或求解区域具有比较复杂的形状时,往往求不到或不易求到其精确解,我们就只得去寻求抛物型方程定解问题的近似解,特别是数值近似解,于是就有了抛物型方程数值解的理论研究.求抛物型数值解的方法主要有:有限差分法、有限元素法.随着社会文明的发展,我国与其它国家的文化交流沟通很全面,偏微分方程的研究方向基本上也是一致的.在微分方程定性理论中有着重要应用的时滞积分不等式以及在差分方程定性理论中有重要应用的时滞离散不等式和关于时间尺度的动力方程理论,在研究了抛物型方程及抛物型方程系统、双曲型方程及双曲型方程系统的强迫振动性理论以及某些时滞脉冲偏微分方程在特定边值条件下解的振动的若干充要条件和多时滞脉冲抛物型微分方程系统解的强迫振动性,推广了已有的结果,建立了若干新定理,促进了偏泛函微分方程理论的发展,并主要创新提出并建立了偏泛函微分方程系统的强迫振动性理论以及具有连续分布变元的双曲型偏泛函微分方程系统和抛物型偏泛函微分方程系统的振动性理论;获得了变系数抛物型偏泛函微分方程解的振动的若干充分必要条件和时滞脉冲抛物型微分方程解的振动的若干充分必要条件和多时滞脉冲抛物型微分系统解的强迫振动性;研究了时滞积分不等式理论和关于时间尺度的动力不等式理论.这些理论的建立,为偏泛函微分方程理论和积分不等式理论的进一步发展起到了非常大的促进作用.3热传导方程的一些知识3.1 热传导方程的导出若物体内各点的温度不相同,则其热量就会从温度较高的地方向温度较低的地方抛物型方程解的估计及其应用传递,这就是常说的热传导现象.由于热量的传递,所以物体内温度随时间和点的位置不同而不同,因此热传导问题可归结为研究物体内部的温度分布情况.下面我们考察一个均匀的、各向同性的物体G 在Ω内部的温度变化规律. 设以(),,u x y z 表示物体G 在Ω内任一点(),,M x y z 处在时刻t 的温度.在Ω内任取一小块区域V ,使V -⊂Ω,并且其边界Γ是光滑的闭曲面,Γ上面积元素的单位外法向量记作n .根据传热学中的傅里叶实验定律[2],物体在无穷小时段dt 内,从V 内经过dS 流出的热量dQ 与时间dt ,流经面积dS 以及温度沿dS 的外法向量的方向导数un∂∂成正比,即udQ k dSdt k u ndSdt n∂=-=-∇⋅∂其中0k >是物体的热传导系数,,,x y z ⎛⎫∂∂∂∇= ⎪∂∂∂⎝⎭.上式中的负号表示热流的方向与温度梯度的方向相反(因为热量总是由温度高处流向温度低处),因此从时刻1t 到时刻2t 经过Γ流入V 内的全部热量 211t t Q dt k u ndsdt Γ=∇⋅⎰⎰⎰若物体Ω内有热源,且热源强度为(),,,F x y z t (即在时刻t 点(),,x y z 处的单位面积在单位时间内发出的热量),则在[]12,t t 内,V 从热源上吸收的热量为 ()212,,,t t VQ F x y z t dxdydzdt =⎰⎰⎰⎰另一方面,在[]12,t t 内,V 内温度从()1,,,u x y z t 升高到()1,,,u x y z t 所需吸收的热量为 ()()321,,,,,,VQ c u x y z t u x y z t dxdydz ρ=-⎡⎤⎣⎦⎰⎰⎰其中为c 物体的比热,ρ为物体的密度. 根据能量守恒,有热传导方程的一些知识123Q Q Q +=若(),,,u x y z t 关于,,x y z 具有二阶连续偏导数,则由高斯公式得 22111t t t t VQ dt k u ndsdt k udxdydzdt Γ=∇⋅=∆⎰⎰⎰⎰⎰⎰⎰这里 ∆ 是laplace 算子,222222x y z∂∂∂∆=++∂∂∂若(),,,u x y z t 关于t 具有一阶连续偏导数,则由Newton-Leibniz 公式有 213t t t VQ dt c u dxdydz ρ=⎰⎰⎰⎰因此有()2211t t t t t VVdt c u dxdydz dt k u F dxdydz ρ=∇+⎰⎰⎰⎰⎰⎰⎰⎰由于时间段[]12,t t 及区域V 是任意取定的,并且被积函数是连续的,则 2t u a u f -∆= 其中2k a c ρ=,Ff c ρ=,并且当0f ≥时,表示Ω内有热源;当0f ≤时,表示Ω内有冷源(即热汇).在适当情况下,方程中描述空间坐标的独立变量的数目还可以减少.例如当物体是均匀细杆时,假如它的侧面是绝热的,也就是说不产生热交换,又假定温度的分布在同一截面是相同的,则温度函数u 仅与坐标x 及时间t 有关,我们就得到一维热传导方程222u u a t x∂∂=∂∂ 同样,如考虑薄片的热传导,薄片的侧面绝热,可得二维热传导方程22222u u u a t x y ⎛⎫∂∂∂=+ ⎪∂∂∂⎝⎭抛物型方程解的估计及其应用3.2 定解问题的提法方程描述的是同类物理现象的共性,但是每一具体的物理现象都是处在各自特定条件之下的,这就需要我们把它所处的特定条件也用数学形式表达出来,我们称这些特定条件为定解条件.定解条件分为初始条件和边界条件.初始条件是说明初始状态的条件,边界条件是描述边界状态的条件,边界条件可分为三类,第一类边界条件(又称Dirichlet 边界条件)是直接给出未知函数在研究区域Ω的边界∂Ω上的值;第二类边界条件(又称Neumann 边界条件)是在∂Ω上给出未知函数u 沿∂Ω沿外法方向n 的方向导数;第三类边界条件(又称为Robin 条件)是在边界∂Ω上给出未知函数u 及其沿∂Ω的外法方向导数的某种线性组合的值.从物理学角度来看,如果知道了物体在边界上的温度状况(或热交换状况)和物体在初始时刻的温度,就可以完全确定物体在以后时刻的温度.因此热传导方程最自然的一个定解问题就是在已给的初始条件与边界条件下求解问题的解.初始条件的提法显然为()(),,,0,,u x y z x y z ϕ=其中(),,x y z ϕ为已知函数,表示物体在0t =时的温度分布第一边界条件:在3R 中的有界区域Ω的导热问题中,若Ω的边界∂Ω处于恒温0u 的环境下,则边界条件为0u u ∂Ω|=若边界温度按已知规律(),,,g x y z t 变化,则(),,,u g x y z t ∂Ω|=第二边界条件:若热量在边界曲面∂Ω各点的流速为(),,,G x y z t ,则由Fourier 定律,边界条件可写成(),,,ug x y z t n ∂=∂ 其中Gg k =-,若0G =,则0u n ∂Ω∂=∂,此时称之为绝热边界条件.定解问题的求解第三边界条件:如果物体内部与周围的介质通过边界∂Ω有热量交换,物体外介质的温度为2u ,物体表面的温度为1u ,内外两种介质间的热交换系数为()110k k >,根据Newton 定律,从物体内部流到外部的热量与两介质间的温度差成正比,即有()112dQ k u u dsdt =-另一方面,由Fourier 定律[3],在时间间隔内从边界曲面上面积元流出的热量为udQ k dsdt n∂=-∂从而有()112uk u u dsdt k dsdt n∂-=-∂即(),,,u u g x y z t n σ∂Ω∂⎛⎫+=⎪∂⎝⎭ 其中1k kσ=, ()1,,,u g x y z t σ= 4 定解问题的求解4.1 初值问题的求解我们可以利用傅里叶变换来求解热传导问题的初值问题.其思想是把原函数变换到另一类函数中去,经过变换,使热传导方程变为常微分方程,从而可以找出一个解,再经过Fourier 的逆变换,得到原热传导方程的解.()()()()2,,,,0,t xx yy u a u u f x y t u x y x y ϕ⎧-+=⎪⎨=⎪⎩ (1) 视t 为参数,先求解齐次热传导方程的初值问题()()()2,,0,txx yy u a u u u x y x y ϕ⎧=+⎪⎨=⎪⎩ (2)对,x y 进行Fourier 变换,记()()12,,,,F u x y t U t λλ=⎡⎤⎣⎦,抛物型方程解的估计及其应用()()12,,F x y ϕλλ=Φ⎡⎤⎣⎦在(1)式两边关于,x y 进行Fourier 变换,原问题变为()()()()()()()222121122121212,,,,,,,,0,d U t a i U t i U t dtU λλλλλλλλλλλλ⎧⎡⎤=+⎪⎣⎦⎨⎪=Φ⎩(3) (2)式是带参数12,λλ的常微分方程的柯西问题,它的解为()()()2121212,,,a tU t e λλλλλλ-+=Φ (4)函数()212a teλλ-+的Fourier 逆变换[4]为()()()()()()()2222221212122222112211221221F 21=2a a i x y a t i xa t i xe t ete d d ed ed λλλλλλλλλλλλπλλπ+∞-+-++-∞----+∞+∞-∞-∞⎡⎤=⎢⎥⎣⎦⎰⎰⎰-()222222111122111111+11cos sin =2cos a t i xa ta ta ted exd i exd exd λλλλλλλλλλλλ----+∞+∞+∞-∞-∞-∞∞-=+⎰⎰⎰⎰令()221+110cos a tI x exd λλλ∞-=⎰()()221222211+/111111202sin 1 =sin cos 2 =2a t a t a t I x e xd e x x xe d a t xI x a tλλλλλλλλλ∞-+∞--+∞0=-⎡⎤∣-⎢⎥⎣⎦-⎰⎰ 解得()224x a tI x ce-=又()2212+1+00 a t y I e d e dy λλ∞-∞-===⎰定解问题的求解则有()222222121421F 4x y aa tet e a tλλπ+--+⎡⎤=⎢⎥⎣⎦-由(4)可得初值问题(2)的解为()()()()222421,,,4x y a tu x y t e d d a tξηϕξηξηπ-+--+∞+∞-∞-∞=⎰⎰ (5)再求解非齐次热传导方程具有齐次初始条件的柯西问题()()()2,,,.00t xx yy u a u u f x y t u x y ⎧=++⎪⎨=⎪⎩ (6) 由齐次化原理[5],此柯西问题的解可写为()()0,,,,;tu x y t x y t d ωττ=⎰而(),,;x y t ωωτ=为下述柯西问题的解:()()()2,,,,,t xx yy a t x y f x y ωωωτωττ⎧=+>⎪⎨=⎪⎩于是,利用(5)式,易知柯西问题(6)的解为()()()()()222420,,1,,4x y ta t u x y t ed d d a t ξητϕξητξητπτ-+--+∞+∞--∞-∞=-⎰⎰⎰ (7)由叠加原理[6],由(5)及(7)就得到柯西问题(1)的解为()()()()()()()()222222424201,,,4,,1 4x y a tx y ta t u x y t ed d a t ed d d a t ξηξητϕξηξηπϕξητξητπτ-+--+∞+∞-∞-∞-+--+∞+∞--∞-∞=+-⎰⎰⎰⎰⎰在上面的推导中,由于预先不知道是否满足进行傅里叶变换及有关计算的条件,所得的解还只是形式解.为证明上式确实是柯西问题(1)的解,还得进行验证.抛物型方程解的估计及其应用4.2 初边值问题的求解热传导方程的初边值问题20 t xx u a u -= (8)00x x l u u ==∣=∣= (9) ()0 t u x ϕ=∣= (10) 令()()() ,u x t X x T t = (11)并要求它满足齐次边界条件(9),这里()X x 及()T t 分别表示仅与x 有关及仅与t 有关的特定函数.将(11)代入方程(8)中,得到()()()()///0 X x T t X x T t -= (12) 将上式分离变量,有()()()()///2T t X x a T t X x λ==- (13)由于在(13)式中,左边仅是t 的函数,右边仅是x 的函数,左右两端要相等,只有等于同一个常数才可能.记次常数为λ-(其值待定),就得到()()/2T aT 0t t λ+= (14)()()//0Xx X x λ+= (15)这样方程(13)就被分离为两个常微分方程,其中一个含有自变量t ,另一个仅含有自变量x ,我们可以通过求解这两个方程来决定()T t 及()X x ,从而得到方程(8)的特解(11)为了使此解是满足齐次边界条件(9)的非平凡解,就必须找到方程(8)满足边界条件定解问题的求解()()00,0X X l == (16) 的非平凡解.方程(15)的通解随0λ>,0λ=以及0λ<而不同,下面分三种情况讨论:情形1 当0λ<时,方程(15)的通解可写成 ()12X x C C e =+要使它满足边界条件(16),就必须1200C C e +=⎧⎪⎨+=⎪⎩由于110e≠只能120C C ==.故在0λ<的情况得不到非平凡解.情形2 当0λ=时,方程(15)的通解可以写成 ()12X x C C X =+ 要满足边界条件(16),()X x 也只能恒等于零.情形3 当0λ>时,方程(15)的通解具有如下形式: ()12X x C C =+ 由边界条件()00X =知10C =,再由()20X l C ==可知,为了使20C ≠,就必须sin 0=.于是222(1,2,)k k k lπλλ===⋯这样就找到了一族非零解()sin(1,2,)k k k X x A x k lπ==⋯ 将固有值代入方程(14)中,可得到其通解为()2222(1,2,)a k tl k k T t B ek π-==⋯这样就得到方程(8)的满足齐次边界(9)的下列分离变量形式的特解:()()()2222k ,sin(1,2,)a k tl k k k k u x t X x T t a ex k lππ-===⋯ 现在我们设法作这种特解的适当的线性组合,以得出初边值问题的解,也就是说,要决定常数k a 使()22221,sina k tl k k k u x t a ex lππ∞-==∑ (17) 满足初始条件(10). 故由初始条件(10)应有()1sink k k x a x lπϕ∞==∑ 由于 1,sink x l π⎧⎫⎨⎬⎩⎭在[]0,l 上正交,因此,k a 是在[]0,l 区间中正弦展开的傅里叶级数的系数,即()02sin l k k a d l lπϕξξξ=⎰ (18) 故()()222201,sin sina k tll k k k u x t d ex l lπππϕξξξ∞-==⋅∑⎰ (19) 是用级数形式表示的初边值问题的形式解.为了考察由分离变量法得到的形式解是否是混合问题的经典解,还得进行验证. 当1C ϕ∈,且()()00l ϕϕ==,()x ϕ是有界函数,(18)式确定的函数(),u x t 是混合问题的解.分析:在求解过程中,级数(17)中的每一项都满足方程(8),因此只要证明级数(17)可以逐项求导两次就好了.也就是说,如果证明了级数(17)求导两次后仍是一致收敛的,那么它一定满足方程(8),此时边界条件(9)和初始条件(10)的满足也是显然的推论了.证明:由于式(19)中含有因子2222a k tl eπ-,因此对于任意0δ>,当0t >时,对任意的0p >,级数22221p a k tl k k el ππ∞-=⎛⎫ ⎪⎝⎭∑均是一致收敛的,而由ϕ是有界函数的假设(()x M ϕ<),可得()0sinlk d Ml lπϕξξξ≤⎰故(19)式中列举的所有级数是一致收敛的,因而,由式(19)表示的级数,当0t >时,关于x 及t 是无穷次可导的,并且求导与求和可以交换.由于级数的每一项都满足方程(8)及边界条件(9)、(10),从而式(19)式表示的级数在0t >时确实满足方程及边界条件.当加上条件()()00l ϕϕ==时,当0t →时,对任意[]0,x l ∈,由式(19)给出的级数趋于初值()x ϕ,即得到式(19)给出的级数确实是初边值问题(8)~(10)的经典解.5 抛物型方程解的估计及其应用先验估计是偏微分方程理论研究中的一个常用的方法.其特点是在假设定解问题解存在的前提下导出解所应当满足的估计,而常用的估计有最大模估计[7],能量估计[8]等等.一般地,我们可以根据先验估计得到定解问题解的唯一性和稳定性,并且可结合其他一些分析方法推导出解的存在性,此外,作为对解的一种估计,先验估计还可能提供关于解的某种性态(如有界性等)方面的信息.5.1 极值原理考虑热传导方程()()2,,,t xx Lu u a u f x t x t Q ≡-=∈其中(){},0,0Q x t x l t T =<<<≤,Q 的侧边和底边统称为Q 的抛物边界,记作Γ,即(){}(){}(){},0,0,,0,0,0x t x t T x t x l t T x t t x l Γ==<≤⋃=<≤⋃=≤≤在热传导过程中,如果物体内部无热源,则热量总是由温度高处向其它地方扩散,而温度最低处的温度会逐渐上升.因此物体的最高温度和最低温度总是在初始时刻或物体的边界上达到.这就是热传导方程的“极值原理”.定理 1(弱极值原理) 设函数()()()2,1,C u x t Q C Q ∈⋂满足Lu f =. (1) 若0f ≤,则u 在Q 上的最大值必在抛物边界Γ上达到,即 ()()max ,max ,Qu x t u x t Γ=(2) 若0f ≥,则()()min ,min ,Qu x t u x t Γ=(3) 若0f =,则()()max ,max ,Qu x t u x t Γ=, ()()min ,min ,Qu x t u x t Γ=同时成立,这里()2,1C Q 表示在Q 内关于x 二次连续可微,且关于t 一次连续可微的函数全体.证明:(1)不妨先考虑0f <情形. 反设存在点()00,x t Q ∈,使得()()00,max ,Qu x t u x t =则在该点处0x u =,0xx u ≤,0t u ≥(如果0t T <,则0t u =;如果0t T =,则0t u ≥).因此()()()00200,,0t xx x t f x t u a u =-≥,这与0f <的假设相矛盾.故(),u x t 不能在Q 内达到最大值,从而有 ()()max ,max ,Qu x t u x t Γ=当 (),0f x t ≤时,设法将其转化为前面的情形.为此构造辅助函数 ()(),,v x t u x t t ε=- 其中ε是任意小的正数.因为0Lv Lu f εε=-=-< 所以()()max ,max ,Qv x t v x t Γ=于是()()()()max ,max ,max ,max ,QQu x t v x t t v x t T u x t T εεεΓΓ=+≤+≤+⎡⎤⎣⎦令0ε→,得()()max ,max ,Qu x t u x t Γ=(2)若0f ≥,则对u -应用情形(1)的结论即可.(3)结合前面两种情况,若0Lu =,则u 在Q 的上的最大值与最小值都在抛物边界Γ上达到.下面我们将弱极值原理推广到稍一般的热传导方程()()()21,,,t xx x Lu u a u b x t u c x t u f x t ≡-++=定理 2 函数()()2,1u C Q C Q ∈⋂满足10L u f =≤,则u 在Q 上的正最大值必在抛物边界Γ上达到,即()()max ,max ,Qu x t u x t +Γ≤由于其证明与定理1的证明方式类似,这里不再赘述.定理3 设()0,c x t c ≥-,其中0c 为正常数.若函数()()()2,1,u x t C Q C Q ∈⋂满足10L u f =≤,且()max ,0u x t Γ≤,则必有()max ,0Qu x t ≤证明 令()()0,,c t v x t e u x t -=,则(),v x t 满足方程 ()0200c t t xx x v a v bv c c v fe --+++=≤ 由于00c c +≥,根据定理2,得()()()0max ,max ,max ,0c t Qv x t v x t e u x t -++ΓΓ≤≤≤因此结论得证.利用定理3,不难得到下列推论:推论1(比较原理) 设()()00,0c x t c c ≥-≥,又设()()2,1,u v C Q C Q ∈⋂,且11L u L v ≤,u v ΓΓ≤,则对任意的(),x t Q ∈,有()(),,u x t v x t ≤5.2 初边值问题解的最大模估计设Ω是n R 中的有界开集,0T >.记(0,]T Q T =Ω⨯,(){}()[0,)0T T Γ=∂Ω⨯⋃Ω⨯这里的T Γ称为T Q 的抛物边界.我们先在T Q 中研究抛物型方程记 []()()1,,int i x i A u u u b x t uf x t ==-∆+=∑[]()()()1,,,int ix i B u u u b x t uc x t u f x t ==-∆++=∑考察第一初边值问题[]()()()()()()()()[]1,, ,,0 ,, ,0,i nt i x Ti A u u u b x t u f x t x t Q u x x x u x t g x t x t T ϕ=⎧=-∆+=∈⎪⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎪⎩∑ (20)定理4 设()()2,1T T u C Q C Q ∈⋂是问题(20)的解,则TQ max u FT B ≤+其中sup TQ F f =,()[]{}0,max max ,max T B x g ϕ∂Ω⨯Ω=证明 令v tF B =+,与u ±作比较.因为 [][]A u F f A u =≥±=± ,(),T x t Q ∈ ()()(),0,0v x B x u x ϕ=≥±=± , x ∈Ω v B g u ∂Ω∂Ω∂Ω≥≥±=± , 0t T ≤≤ 由比较原理知,u v FT B ±≤≤+,即 ()TQ max ,u x t FT B ≤+推论 2 第一初边值问题(20)的解在函数类()()2,1T T C Q C Q ⋂中是唯一的,且连续地依赖于f ,ϕ和g .证明 当0f g ϕ==≡时,对应的解u 满足TQ max 0u =,故0u ≡,从而解是唯一的.假设i u 是对应于{},,i i i f g ϕ的解,1,2i =,则12u u -是对应于{}121212,,f f g g ϕϕ---的解.于是[]{}TT121212120,Q Q max max ax max ,max T u u T f f g g ϕϕ∂Ω⨯Ω-≤-+--所以当{}111,,f g ϕ与{}222,,f g ϕ充分接近时,1u 与2u 也充分接近,这说明问题(20)的解连续地依赖于f ,ϕ和g .现在考察第一初边值问题[]()()()()()()()[], ,,0 ,, ,0,TB u f x t x t Q u x x x u x t g x t x t T ϕ⎧=∈⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎩ (21) 定理5 设()0,c x t c ≥-,()()2,1T T u C Q C Q ∈⋂是问题(21)的解,则 ()0TQ max c T u e FT B ≤+其中sup TQ F f =,()[]{}0,max max ,max T B x g ϕ∂Ω⨯Ω=.证明 不妨认为00c ≥,令()0c t v e FT B =+,与u ±作比较.因为[]()()()()()()[]()00000000, =, ,c t c t c t c t c t c t T B u Fe c e Ft B c x t e Ft B Fe e c c x t Ft B Fe F f B u x t Q =+++++++≥≥≥±=±∈()()(),0,0v x B x u x ϕ=≥±=± , x ∈Ω v B g u ∂Ω∂Ω∂Ω≥≥±=± , 0t T ≤≤由比较原理知,()0c T u v e FT B ±≤≤+,即()()0TQ max , c T u x t e FT B ≤+5.3 初值问题解的最大模估计记[]T D 0,n R T =⨯,[](),t C u u u c x t u =-∆+ 考察初值问题[]()()()(), ,,0 TnC u f x t x tD u x x x Rϕ⎧=∈⎪⎨=∈⎪⎩ (22) 设(),c x t 连续,()()00,0c x t c c ≥->,(),f x t 和()x ϕ有界,记 sup TD F f =, sup nR ϕΦ=如果()()2,1T T u C D C D ∈⋂是初值问题(22)的解,则 ()0sup Tc T D u e FT ≤+Φ证明 令()()0,,c t v x t u x t e -=,则v 满足[]()()()(),,,0 t nD v v v c x t v f x t v x x x R ϕ⎧=-∆+=⎪⎨=∈⎪⎩ (23) 其中()()0,,0c x t c x t c =+≥,()()0,,c t f x t e f x t -=由于解得先验估计方法不能直接用于初值问题,我们希望借助于一个有界区域上的初边值问题进行讨论,任意取定较大的常数L ,记{}](,0,L T D x L T =≤⨯.因为解u 有界,所以存在正常数K 使得u K ≤在D T 上成立,在有界区域,L T D 上考虑辅助函数()()22,2K w x t Ft x nt v L =+Φ++± 直接计算知,在,L T D 上w 满足[]()()()()()()002,22220 ,,0 ,,0c t L T c tx L x L K D w F c Ft x nt e f x t D L K w x x x x LL w x t K u x t e ϕ--==⎧⎧⎫=++Φ++±≥∈⎨⎬⎪⎩⎭⎪⎪=Φ+±≤⎨⎪⎪≥Φ+±>⎪⎩利用比较原理知,(),0w x t ≥在,L T D 上成立对于D T 内的任一点()00,x t ,取L 充分大使得()00,,L T x t D ∈,于是()00,0w x t ≥ 即()()2000002,2K v x t Ft x nt L≤+Φ++ 令L →∞得()000,v x t Ft Ft ≤+Φ≤+Φ从而()()()000000,,c t c T u x t v x t e e Ft =≤+Φ由()00,T x t D ∈的任意性知,估计式(23)成立.推论3 初值问题(23)的解在函数类()()2,1T T C D C D ⋂中是唯一的,且连续地依赖于f ,ϕ.由于其证明与推论3的证明方式类似,这里不再赘述.5.4 初边值问题的能量估计设Ω是n R 中的一个光滑区域,在](0,T Q T =Ω⨯上考察第一初边值问题()()()()()[], ,,0 0 ,0,t T u u f x t x t Q u x x x u x t T ϕ-∆=∈⎧⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎩ (24) 定理6 设()()1,02,1T T u C Q C Q ∈⋂是问题(23)的解,则存在正常数()C C T =使得()222200max ,2TT t Tux t dx u dxdt C dx f dxdt ϕΩΩΩΩ≤≤⎛⎫+∇≤+ ⎪⎝⎭⎰⎰⎰⎰⎰⎰ (25) 证明 问题(24的方程两边乘以u 并在T Q 上积分,得000tttt uu dxdt u udxdt f udxdt ΩΩΩ-∆=⎰⎰⎰⎰⎰⎰(26)对(26)式左端第一项中关于t 的积分利用分部积分以及初值条件,可知()()22011,22t t uu dt u x t x ϕ=-⎰ (27)对(26)式左端第二项关于x 的积分利用散度定理以及边界条件,推出22u u udx u dS u dx u dx n Ω∂ΩΩΩ∂∆=-∇=-∇∂⎰⎰⎰⎰ (28)将(27)式和(28)式代入(26)式,得2220022ttu dx u dxdt f udxdt dx ϕΩΩΩΩ+∇=+⎰⎰⎰⎰⎰⎰ (29)利用不等式222ab a b ≤+可知 220002t ttf udxdt f dxdt u dxdt ΩΩΩ≤+⎰⎰⎰⎰⎰⎰将上式代入(29)式,得222220002tttu dx u dxdt f dxd u dxdt dx ϕΩΩΩΩΩ+∇≤++⎰⎰⎰⎰⎰⎰⎰⎰ (30)记 ()20tY t u dxdt Ω=⎰⎰,()220t F t f dxd dx ϕΩΩ=+⎰⎰⎰那么不等式蕴含()()()Y t Y t F t '≤+ 利用Gronwall 不等式[9]推出()()()()()()2022001 tt t t ttu dxdtF t Y t Y e e F t e F t e f dxd dx ϕΩΩΩ=≤+-⎛⎫≤=+ ⎪⎝⎭⎰⎰⎰⎰⎰将上式代入(30)式知()22220021tt tu dx u dxdt e f dxd dx ϕΩΩΩΩ⎛⎫+∇≤++ ⎪⎝⎭⎰⎰⎰⎰⎰⎰ 此式两边关于t 在[]0,T 上取上确界,就得到估计式(25).下面我们将讨论一般形式的二阶抛物型初边值问题.设Ω为n R 中的有界区域,且有光滑边界()0,T Q T =Ω⨯,在区域中讨论一般形式的二阶抛物型初边值问题()()()(),11,,,,i j i nnij x x i x i j i u a x t u u b x t u c x t u f x t t ==∂-++=∂∑∑ (31) ()0 t u x x ϕ==∈Ω (32)0T u ∑= (33)解的性质.式中,()0,T T ∑=Γ⨯为区域的侧边界;()12,,n x x x x =∈Ω为方便讨论,作如下假设:(1) 系数ij a 、i b 、c 及右端项f 都是T Q 上的连续函数,并且ij a 在T Q 上还具有一阶连续偏导数. (2) 对一切,1,2,i j n =;ij ji a a =且存在正常数0α>,使得对一切(),T x t Q ∈及任意给定的实向量()12,,,n ξξξ,有:()2,11,nnijiji i j i a x t ξξαξ==≥∑∑成立.对于初边值问题的解,定义能量函数:()212E t u dx Ω=⎰ (34)定理7 若(),u x t 为初边值问题(31)~(33)的解,能量函数()E t 按式(34)定义,则能量估计式:()()200 0t Ct Ct E t E e Ce f dxdt t T Ω≤+≤≤⎰⎰(35)成立.其中,C 为一个不依赖于u 的正常数.证明 用u 乘以式(31),并在Ω上关于x 积分,就得到:()(),11,,i j i n n t ij x x i x i j i u udx a x t u u dx b x t u u cu dx fudx ΩΩΩΩ==⎛⎫⎛⎫-++= ⎪ ⎪⎝⎭⎝⎭∑∑⎰⎰⎰⎰ []0,t T ∈ (36) 式左端的第一项可以写成212d u dx dt Ω⎧⎫⎨⎬⎩⎭⎰;当3n ≥时,记12,,,n ααα为侧边界T ∑法向量的方向角,dS 为广义面积微元.令(),1,2,,i ij ij x p a uu i j n ==,固定i ,让1,2,,j n =,利用高维高斯公式[10],并注意边界条件(它隐含着0Tu ∑=),边界积分项为零,可得()()()()()()12121122121212120cos cos cos = = =Ti i i n i ni i in n i i in n i x x i x x in x x i i in x x x x p p p dSp p p dx xx x a u a u a u udx a u a u a u u dxααα∑ΩΩΩ=++⎛⎫∂∂∂+++ ⎪∂∂∂⎝⎭+++++⎰⎰⎰⎰故对固定的i ,有:()()()()()12121212=i i i n i ni x x i x x in x x i i in x x x x a u a u a u udx a u a u a u u dx ΩΩ-+++++⎰⎰(37)成立,对式(37)关于i 从1到n 求和.式(36)左端的第二项可以写成:(),1,1,1i j i j i i n n n ij x x ij x x ij x x i j i j i j a u u dx a u u dx a u u dx ΩΩΩ===⎛⎫⎛⎫⎛⎫-=+ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭∑∑∑⎰⎰⎰ (38)将上式的第二项,连同式右端的第三、四项移至等式右边,并将其和记为(),t Q u u dx Ω⎰则有()()()1,1,,i i i n n ti x ij x x i i j Q u u dx fudx b x t u u cu a u u dx ΩΩΩ==⎛⎫=-++ ⎪⎝⎭∑∑⎰⎰⎰则由于系数的可微性假设(1)可得,对一切0t T ≤≤成立()21,i n t T x i Q u u dx C u u u dx ΩΩ=⎛⎫≤+ ⎪⎝⎭∑⎰⎰ (39)其中T C 为一个不依赖于T 的正常数,但与u 无关.对任意给定的0ε>,有2211122i innnx x i i i nuu dx udx u dx εεΩΩΩ===≤+∑∑∑⎰⎰⎰ (40)取TC αε=,由式(40)就得到()22111,2inntx i i Q u u dx udx C u dx αΩΩΩ==≤+∑∑⎰⎰⎰(41)其中212TT nC C C α=+,将式(41)代入式(36),容易得到2221,11111222i j i n n n ij x x x i j i i dE a u u dx u dx C u dx f dx dt αΩΩΩΩ===⎛⎫⎛⎫+≤+++ ⎪ ⎪⎝⎭⎝⎭∑∑∑⎰⎰⎰⎰ (42) 再注意到由假设(2)有2,11i j i n n ij x x x i j i a u u dx u dx αΩΩ==⎛⎫≥ ⎪⎝⎭∑∑⎰⎰ 就可得到()22dEC E t f dx dtΩ≤+⎰ (43)其中2121C C =+.在式(43)两边乘以2C t e -再对t 积分,,并放大被积函数,即可得 ()()200t Ct Ct E t E e Ce f dxdt Ω≤+⎰⎰定理证毕.5.5 能量不等式的应用5.5.1 初边值问题解的唯一性热传导方程是抛物型方程的典型代表.下面考虑二维热传导方程的初边值问题()2t xx yy u a u u f =++ (44)()0,t u x y ϕ== (45) (),,u x y t μΓ= (46) 这里,Γ表示Ω的边界,应用能量不等式可得如下定理.定理8 若热传导方程的初边值问题的解存在,则其解唯一.证明 设1u ,2u 是该定解问题的两个解,则其差12u u u =-满足相应的齐次方程及齐次初始条件和齐次边界条件.此时的齐次方程满足假设(1)、(2),有(34)式定义的能量函数知,在初始时刻有()00E =,故由能量不等式(35)得:()()22220x y E t u a u u dxdy Ω⎡⎤=++=⎣⎦⎰⎰ 即0x y u u u ===,从而可推出(),,u x y t const =.又由于在初始时刻0u =,故得(),,0u x y t ≡.即12u u =.这样就证明了初边值问题(44)~(46)解的唯一性. 5.5.2 初边值问题解的稳定性为了记号简单起见,对于定义在区域Ω上的函数ϕ和定义在区域上()0,T ⨯Ω的函数f ,常以()2L ϕΩ和()()20,L T f ⨯Ω分别表示()122dxdy ϕΩ⎰⎰和()1220Tf dxdydt Ω⎰⎰⎰.定理9 热传导方程的初边值问题:()2t xx yy u a u u f =++()0,t u x y ϕ== 0u Γ=的解(),,u x y t ,在下述意义下关于初始值ϕ与方程右端项f 是稳定的:对任何给定的0ε>,一定可以找到仅依赖于ε和T 的0η>,只要 ()212L ϕϕηΩ-≤ ()212x xL ϕϕηΩ-≤()212y yL ϕϕηΩ-≤ ()()2120,L T f f η⨯Ω-≤ (47)那么以1ϕ为初值、1f 为右端项的解1u 与以2ϕ为初值、2f 为右端项的解2u 之差在上满足()212L u u εΩ-≤ ()212x xL u u εΩ-≤ ()212y yL u u εΩ-≤ (48)证明 记12u u u =-,12ϕϕϕ=-,1f f f =-,则u 满足()2t xx yy u a u u f =++ (49)()0,t u x y ϕ== (50) 0u Γ= (51) 方程(49)满足假设(1)、(2),从而利用能量不等式(35),可得:()()()()222000tTCt Ct E t E e Ce f dxdydt C E f dxdydt ΩΩ≤+≤+⎰⎰⎰⎰⎰⎰[]0,t T ∈ (52)式中,2C 为一个仅依赖于T 的正常数.记。
10_抛物型方程的有限差分方法
10_抛物型方程的有限差分方法抛物型方程是一类常见的偏微分方程,广泛应用于自然科学和工程学的领域中。
有限差分方法是一种常用的数值求解抛物型方程的方法之一、本文将介绍抛物型方程的有限差分方法(II)。
有限差分方法主要基于离散化的思想,将偏微分方程转化为差分方程,进而求解差分方程的数值解。
对于抛物型方程,其一般形式可以表示为:∂u/∂t=Δu+f(x,t)其中,u(x, t)是未知函数,表示空间位置x和时间t上的解,Δu表示Laplace算子作用于u的结果,f(x, t)是已知函数。
有限差分方法的基本思想是将空间和时间域进行离散化,将连续的空间和时间划分为有限个网格点,然后使用差分近似代替偏导数,得到差分方程。
假设空间域被划分为Nx个网格点,时间域被划分为Nt个网格点,对于每个网格点(i,j),可以表示为(x_i,t_j),其中i=0,1,...,Nx,j=0,1,...,Nt。
在有限差分方法中,我们使用中心差分近似来代替偏导数。
对于时间导数,可以使用向前差分或向后差分,这里我们使用向前差分,即:∂u/∂t≈(u_i,j+1-u_i,j)/Δt对于空间导数,可以使用中心差分,即:∂^2u/∂x^2≈(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2将上述差分近似代入抛物型方程中,可以得到差分方程的离散形式:(u_i,j+1-u_i,j)/Δt=(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2+f_i,j其中,f_i,j=f(x_i,t_j)。
重排上式,可以得到递推关系式:u_i,j+1=αu_i-1,j+(1-2α)u_i,j+αu_i+1,j+Δt*f_i,j其中,α=Δt/Δx^2通过设置初始条件和边界条件,可以利用以上递推关系式得到抛物型方程的数值解。
总结来说,抛物型方程的有限差分方法(II)是一种常用的数值求解抛物型方程的方法。
它基于离散化的思想,将偏微分方程转化为差分方程,然后利用中心差分近似代替偏导数,得到差分方程的离散形式。
抛物方程的有限差分法
抛物方程的有限差分法作者:李娜来源:《科技视界》2014年第32期【摘要】抛物方程是描述物理现象的一类重要方程,其中差分方法和有限元方法是求其数值解的两类主要方法。
本文主要介绍有限元方法中的向前差分法,首先简单介绍向前差分法,给出稳定性和收敛性的概念,然后以一维热传导方程为例进行求解,同时给出收敛性和稳定性分析,并利用Matlab软件做出了误差分析图。
【关键词】抛物方程;有限元方法;向前差分法;误差分析0 引言由于抛物型方程与时间t有关,称为非驻定问题。
非驻定问题可用差分法,也可用有限元法求解。
热传导方程式(或称热方程)是一个重要的偏微分方程,它描述一个区域内的温度如何随时间变化。
热传导在三维等方向均匀介质里的传播可用方程式u■=kΔu表示,其中u■=u (t,x,y,z)表示温度,它是时间变量t与空间变量(x,y,z)的函数,■是空间中一点的温度对时间的变化率,uxx、uyy和uzz是温度对三个空间坐标轴的二次导数。
k决定于材料的热传导率、密度与热容。
求解方程时,如果考虑的介质不是整个空间,则为了得到方程唯一解,必须指定的边界条件。
如果介质是整个空间,为了得到唯一解,必须假定解的增长速度有一个指数型的上界,并且此假定与实验结果相吻合。
1 本文研究的方程本文主要研究一维热传导方程的有限差分解法,下面给出了各向同性介质中无热源的一维热传导方程及初始条件:■=a(x,t)■a>0 0<x<1,0<t<Tux,0=?覫x=sin(πx) 0<x<1u0,t=u(1,t)=0 0≤t≤T (1)在此,本文利用有限元方法中的向前差分法求解偏微分方程式(1),首先需要建立差分格式,而在建立差分格式时通常取空间步长和时间步长为常量。
下面介绍向前差分的概念以及如何利用该方法对其进行收敛性、精确性和稳定性分析。
1.1 向前差分格式有限差分法和有限元方法是求解偏微分方程的两种主要的数值方法。
二维抛物方程的有限差分法
二维抛物方程的有限差分法摘要二维抛物方程是一类有广泛应用的偏微分方程,由于大部分抛物方程都难以求得解析解,故考虑采用数值方法求解。
有限差分法是最简单又极为重要的解微分方程的数值方法。
本文介绍了二维抛物方程的有限差分法。
首先,简单介绍了抛物方程的应用背景,解抛物方程的常见数值方法,有限差分法的产生背景和发展应用。
讨论了抛物方程的有限差分法建立的基础,并介绍了有限差分方法的收敛性和稳定性。
其次,介绍了几种常用的差分格式,有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等,重点介绍了古典显式格式和交替方向隐式格式。
进行了格式的推导,分析了格式的收敛性、稳定性。
并以热传导方程为数值算例,运用差分方法求解。
通过数值算例,得出古典显式格式计算起来较简单,但稳定性条件较苛刻;而交替方向隐式格式无条件稳定。
关键词:二维抛物方程;有限差分法;古典显式格式;交替方向隐式格式FINITE DIFFERENCE METHOD FORTWO-DIMENSIONAL PARABOLICEQUATIONAbstractTwo-dimensional parabolic equation is a widely used class of partial differential equations. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensional parabolic equation.Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direction implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differential method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable.Keywords:Two-dimensional Parabolic Equation; Finite-Difference Method; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme目录摘要 (I)Abstract (II)1绪论 (1)1.1课题背景 (1)1.2发展概况 (1)1.2.1抛物型方程的常见数值解法 (1)1.2.2有限差分方法的发展 (2)1.3差分格式建立的基础 (3)1.3.1区域剖分 (3)1.3.2差商代替微商 (3)1.3.3差商代替微商格式的误差分析 (4)1.4本文主要研究容 (5)2显式差分格式 (7)2.1常系数热传导方程的古典显式格式 (7)2.1.1古典显式格式格式的推导 (7)2.1.3古典显式格式的算法步骤 (8)3隐式差分格式 (10)3.1古典隐式格式 (10)3.2 Crank-Nicolson隐式格式 (12)3.3 Douglas差分格式 (13)3.4加权六点隐式格式 (14)3.5交替方向隐式格式 (15)3.5.1 Peaceman-Rachford格式 (15)3.5.2 Rachford-Mitchell格式 (15)3.5.3 Mitchell-Fairweather格式 (15)3.5.4交替方向隐式格式的算法步骤 (16)4实例分析与结果分析 (17)4.1算例 (17)4.1.1已知有精确解的热传导问题 (17)4.1.2未知精确解的热传导问题 (19)4.2结果分析 (20)5稳定性探究与分析 (21)5.1稳定性问题的提出 (21)5.2 几种分析稳定性的方法 (21)5.3 r变化对稳定性的探究 (23)5.3.1 古典显式格式的稳定性 (23)5.3.2 P-R格式格式的稳定性 (24)结语 (26)参考文献 (27)附录P-R格式的C++实现代码 (28)致谢 (30)1绪论1.1课题背景抛物方程是一类特殊的偏微分方程,二维抛物方程的一般形式为u Lu t∂=∂ (1-1) 其中1212((,,))((,,))(,,)(,,)(,,)u u u u u u L a x y t a x y t b x y t b x y t C x y t x x y y x y∂∂∂∂∂∂=++++∂∂∂∂∂∂ 120,0,0a a C >>≥。
有限差分法求解抛物型方程说明
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下,才能求得定解问题解的解析式,对比较复杂的问题要找到解的解析表达式是困难的,因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法,只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数,得到相应的差分方程,通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想,并给出一些具体的数值实例.§1 差分方法的基本思想有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格,主要采用Taylor 级数展开等方法,在每个网格节点上用有限差分近似公式代替方程中的导数,从而建立以网格节点上的函数值为未知数的代数方程组.有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑,可分为中心格式和逆风格式.考虑时间因子的影响,差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式,主要是上述几种格式的组合,不同的组合构成不同的差分格式.泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式.首先考虑单变量函数()u x ,如图1把区域x 离散为一批结点,记0()(), =0,1,2,i i u x u x ih u i =+=图1 单变量函数离散化函数()u x 在点i x 处的泰勒展开式为23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''+=++++ (1)或23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''-=-+-+ (2)式(1)和(2)重新整理可得2()()()()()2!3!i i i i i u x h u x u x u x u x h h h '''''+-'=---(3)和2()()()()()2!3!i i i i i u x u x h u x u x u x h h h '''''--'=+++(4)于是给出在点i x 处函数u 的一阶导数的两个近似公式1()()()i i i ii u x h u x u u u x h h ++--'≈= (5)1()()()i i i i i u x u x h u u u x h h----'≈= (6)因为级数被截断,这两个近似公式肯定要产生误差,此误差与h 同阶,形式分别为()(), ,2()(), .2i i i i i i hE u O h x x h hE u O h x h x ξξξξ''=-=≤≤+''==-≤≤ 若把式(3)和(4)相加并求()i u x ',可得11()()()22i i i i i u x h u x h u u u x h h+-+---'≈= (7)其截断误差与2h 同阶,形式为22()(), ,6i i i h E u O h x h x h ξξ''=-=-≤≤+若把式(3)和(4)相减并求()i u x '',可得1122()2()()2()i i i i i i i u x h u x u x h u u u u x h h +-+-+--+''≈= (8)其截断误差与2h 同阶,其形式为22()(), ,12i i i h E u O h x h x h ξξ''=-=-≤≤+我们可继续用这种方式来推导更复杂的公式,类似的公式还有很多,这里不再一一列举.公式(5)、(6)分别称为一阶向前、向后差分格式,这两种格式具有一阶计算精度,公式(7)、(8)分别称为一阶、二阶中心差分格式,这两种格式具有二阶计算精度.图2 二维区域网格剖分上面的结果可直接推广使用于导出二元函数(,)u x y 的许多有限差分近似公式.如图7.2,把求解区域进行网格剖分,使12(,)(,), ,=0,1,2,i j ij u x y u ih jh u i j ==其中x 方向的网格间距为1,h y 方向的网格间距为2,h 整数i 和j 分别表示函数(,)u x y 沿x 坐标和y 坐标的位置.二元函数(,)u x y 对x 求偏导时y 保持不变,对y 求偏导时x 保持不变,根据向前差分公式(7.5)可以给出在点(,)i j x y 处函数(,)u x y 的一阶偏导数的两个近似公式1,,1(,)i j i j i ju x y u u xh +∂-≈∂ (9),1,2(,)i j i j i ju x y u u yh +∂-≈∂ (10)相类似地,根据二阶中心差分格式(8)可以得到函数(,)u x y 的二阶偏导数的近似公式21,,1,221(,)2i j i j i j i ju x y u u u x h +-∂-+≈∂ (11)2,1,,1222(,)2i j i j i j i j u x y u u u yh+-∂-+≈∂ (12)下面我们推导函数(,)u x y 的二阶混合偏导数2ux y∂∂∂在(,)i j x y 的有限差分表达式.根据一阶中心差分格式(7),112111,11,11,11,122121221,11,1(,)(,)(,)1()21 ()()222 i j i j i j i j i j i j i j i j i j i u x y u x y u x y O h x y h y y u u u u O h O h h h h u u u +-+++--+--+++-∂∂∂⎡⎤⎡⎤∂=-+⎢⎥⎢⎥∂∂∂∂⎣⎦⎣⎦--⎡⎤=-++⎢⎥⎣⎦--≈1,11,1124j i j u h h -+--+二维有限差分近似可以直接推广到三维空间或三维空间加一维时间的情形.定义1 当步长趋于零时,差分方程的截断误差趋于零,则称差分格式与微分方程是相容的.定义2 当步长趋于零时,差分方程的解收敛于微分方程的解,则称差分格式是收敛的. 定义3 当差分方程的解由于舍入误差的影响,所产生的偏差可以得到控制时,则称差分格式是稳定的.§2 抛物型方程的有限的差分法为了说明如何使用有限差分法来求解偏微分方程,本节我们给出以下几个数值实例.算例1 考虑一维非齐次热传导方程的初边值问题:2212(,), 01,01,(,0)(), 01,(0,)(), (1,)(), 0 1.u ua f x t x t t x u x q x x u t g t u t g t t ⎧∂∂=+<<<≤⎪∂∂⎪⎪=≤≤⎨⎪==<≤⎪⎪⎩(7.13),其中2,a =函数11(,)[cos()2sin()],22xf x t e t t =--+-初始条件1()sin,2xq x e =左、右边界条件分别为11()sin(),2g t t =-21()sin()2g t e t =-.该定解问题的解析解为1(,)sin(),(,)[0,1][0,1].2xu x t e t x t =-∈⨯将求解区域{(,)|,0}x t a x b t T Ω=≤≤≤≤进行网格剖分,[,]a b 作m 等分,[0,]T 作n 等分,记,,b a Th m nτ-==则 ,0,,0i k x a ih i M t k k n τ=+≤≤=≤≤对该问题建立如下向前差分格式:11122, 11, 11,k kk k k k i i i i i i u u u u u a f i m k n hτ+-+--+=+≤≤-≤≤-(14) (,0)(),1,i i u x q x i m =≤≤ (15) 12(,)(), (,)(),1.k k k k u a t g t u b t g t k n ==≤≤ (16)令2r ah τ=,差分格式(7.14)整理得111(12), 11, 1 1.k k k k k i i i i i u ru r u ru f i m k n τ+-+=+-++≤≤-≤≤- (17)显然时间在1k t +上的每个逼近值可独立地由k t 层上的值求出。
3-抛物型方程的有限差分法
中都是精确的,则初始 误差的传递情况如表 1:
表1 r=1/2时Richardson格式的误差传播
-4
0 1 2 3 4 5 6
-3
-2
-1
0
1
2
3
4
0 0 0 0
0 0 0
0 0
0
0
0 0
0 0 0
0 0 0 0
2
4 7 4
24 17 6
6 17
于( x j , t
1 k 2
)(t
k
1 2
1 ( k ) )展开,则得 2
2 2 Rk ( u ) 0 ( h ). (1.9) j
(四)
Richardson 格式,即 a
k k uk 2 u u j 1 j j 1 2
1 k 1 uk u j j
2 h 1 k k k k 1 或u k 2 r ( u 2 u u ) u 2f j .(1.10) j j 1 j j 1 j
例1 写出向前差分格式的矩 阵形式。
1 k uk u j j
解
a
k k uk 2 u u j 1 j j 1
h
2
fj
1 k k k uk ru ( 1 2 r ) u u j j 1 j j 1 f j
显然A I (( N 1)阶单位矩阵, B (1 2r ) I rS , 其中 0 1 0 1 0 0 S 0 0 0 1 0 0 1 0 故C (1 2r ) I rS .
( 2.4)
同理:对于向后差分格 式,即
抛物方程的有限差分法
图1
,我们需要求解这1/h +1()×T/τ+1()个点对应的函数值实上由已知的初边值条件蓝色标记附近的点可直接得到,所以只要确定微分方程的解在其它点上的取值即可,可记为u []
k j
=u (x j ,t k )。
建立差分格式
j =1, (1)
-1;k =0,1,…,T τ-1,用向前差分代替关于时间的
一阶偏导数,用二阶中心差分代替关于空间的二阶偏导数,则可定义最简显格式:
-u k j =u k j+1-2u k j +u k
j-1
h
2
变形有:
(上接第50页)极大值理论,检测初始行波、故障点反射波和对端母线反射波到达测量端的时间,测量故障点距离,从测试结果看,该方案有效弥补传统行波测距的不足之处,提高了故障测距的精确度。
【参考文献】
[1]陈靖.行波法故障测距的理论研究及其实现方案[D].武汉:武汉大学,2004.数值解的剖分图如图2:
图2
真解与数值解的误差剖分图如图3:
图3
3数值实验及结果分析
我们对所求解的初边值问题(1)进行算法精度的数值实验,当
u 0
(x )sin πx 时,边界值仍然为u (0,t )=u (1,t )=0,其精确解为:u (x ,t )
从表中我们可以看出。
. All Rights Reserved.。
经典偏微分方程课后习题答案
第四章 抛物型微分方程有限差分法1设已知初边值问题22, 01, 0<(,0)sin , 01(0,)(1,)0, 0 u ux t t x u x x x u t u t t T π⎧∂∂=<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩T ≤, 试用最简显格式求上述问题的数值解。
取h=0.1,r=0.1.0 1/10 2/10 … 1 T 2τ τt解: 1.矩形网格剖分区域. 取空间步长1, 时间2510h =0.00τ=以及0.01τ=的矩形网格剖分区域, 用节点)表示坐标点(,j k (,)(,)j k x t jh k τ=, 0,1,...1/; 0,1,...,/j h k T τ==, 如图所示.显然, 我们需要求解这(1/1)(/1)h T τ+×+个点对应的函数值. 事实上由已知初边界条件蓝标附近的点可直接得到, 所以只要确定微分方程的解在其它点上的取值即可. 沿用记号[]k(,)j j k u x t =。
u 2. 建立差分格式, 对于11,...1; 0,1,...,1Tj k hτ=−=−, 用向前差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式:1122k k k k k1jj j j u u u u u h ++−+=. 变形j τ−−有:1112(12) (k k k kj j j j u ru r u ru r h τ+−+=+−+=(4.1)用向后差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式最简隐格式:111122k k k k k j jj j j u u u u u h τ++++−−+=11+−1kj +,变形有:1111(12) k k k j j j ru r u ru u ++−−−++−= (4.2)(4.1)*0.5+(4.2)*0.5得CN 格式为:111112222k k k k k k k k j jj j j j j j u u u u u u u u h τ+++−+−−++−+=111++−1kj +x x变形有:111111(22)(22) k k k k k j j j j j ru r u ru ru r u ru ++−−+−−++−=+−+ (4.3)3 初边界点差分格式处理.对于初始条件u x (,0)sin , 01=π≤≤h 离散为(4.4)0sin 0,1,...1/j u jh j π==对于边界条件离散为(0,)(1,)0, 0 u t u t t T ==≤≤00 0,1,.../k k N u u k T τ===(4.5)总结: 联立方程(4.1)(4.4)(4.5)得到已知问题的最简显格式差分方程组:11100(12)1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N u ru r u ru T j k h u jh j h u u k T τπτ+−+⎧=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.2)( 4.4)( 4.5)得到已知问题的最简隐格式差分方程组:1111100(12) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N ru r u ru u T j k h u jh j h u u k T τπτ++−−+⎧−++−=⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.3)( 4.4)( 4.5)得到已知问题的CN 格式差分方程组:11111100(22)(22) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k k j j j j j jk k N ru r u ru ru r u ru T j k h u jh j h u u k T τπτ++−−+−⎧−++−=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩1k j + 4 求解并显示结果利用软件计算(Matlab)如上最简显格式差分方程组.h=1/10;tau=0.0025;T=0.5; r=tau/h^2;M=1/h+1;N=T/tau+1; u=zeros(M,N);for m=1:Mu(m,1)=sin((m-1)*h*pi); endu(1,1:N)=0;u(M,1:N)=0;for n=1:N-1for m=2:M-1u(m,n+1)=r*(u(m+1,n)+u(m-1,n))+(1-2*r)*u(m,n); end end u=u’ 这样我们就计算出不同时刻不同位置k t j x 对应的函数值(,)j k u x t 取tau=0.0025, 即r=0.25绘图, 取tau=0.01, r=1再绘图,如图()图4.2 习题1数值解图示(左r=0.25, 右r=1)2.试构造初边值问题 ()()()()(), 0.51, 0,,0, 0.51,0.5,0, 1,0.51,, 0u u x x x T t x x u x x x u ⎪∂u t t u t t T x ϕ⎧∂∂∂⎛⎞=<<<≤⎜⎟⎪∂∂∂⎝⎠⎪⎪=≤≤⎨⎪==−≤≤⎪∂⎩的显格式,并给出其按最大范数稳定的充分条件。
抛物型方程的计算方法
分类号:O241.82本科生毕业论文(设计)题目:一类抛物型方程的计算方法作者单位数学与信息科学学院作者姓名专业班级2011级数学与应用数学创新2班指导教师论文完成时间二〇一五年四月一类抛物型方程的数值计算方法(数学与信息科学学院数学与应用数学专业2011级创新2班)指导教师摘要: 抛物型方程数值求解常用方法有差分方法、有限元方法等。
差分方法是一种对方程直接进行离散化后得到的差分计算格式,有限元方法是基于抛物型方程的变分形式给出的数值计算格式。
本文首先给出抛物型方程的差分计算方法,并分析了相应差分格式的收敛性、稳定性等基本理论问题.然后,给出抛物型方程的有限元计算方法及理论分析。
关键词:差分方法,有限元方法,收敛性,稳定性Numerical computation methods for a parabolic equationYan qian(Class 2, Grade 2011,College of Mathematics and Information Science)Advisor: Nie huaAbstract:The common methods to solve parabolic equations include differential method,finite element method etc。
The main idea of differential method is to construct differential schemes by discretizing differential equations directly. Finite element scheme is based on the variational method of parabolic equations。
In this article, we give some differential schemes for a parabolic equation and analyze their convergence and stability. Moreover,the finite element method and the corresponding theoretical analysis for parabolic equation are established.Key words:differential method,finite element method, convergence,stability1 绪 论1。
偏微分方程的离散化方法4重要知识
u(x x1, y, t) u(x, y, t)
,
u
u(x, y, t) u(x x2 , y, t)
x x x1
x1
x x x1
x2
2
2
k u
k x x1 2
u(x
x1,
y, t) x1
u( x,
y,
t)
k x x2 2
u(x, y, t) u(x x2 , y, t) x2
(1 2 )Pin
(
Pn i1
Pn i1
)
,
t x
2
,截断误差:O(t
x2 )
从方程可以看出:如果已知第 n(本步时间)的值 Pin ,就可以求得第 n+1
时刻(下步时间)的值
P n1 i
。因此如初始条件,即
n=0
时各网格的
P
值已给定,
就可以依次求得以后各时间的 P 值。这种差分格式是显式差分格式。在显式差分
re
re 0.14 x 2 y 2
2)、对于各向异性地层,步长: x, y ,其等值供给半径:
re
0.28
Ky Kx
1/ 2
x 2
Kx Ky
1/ 2
1/ 2
y 2
K K
y x
1/ 4
Kx Ky
1/ 4
重点辅导
28
ai,
j
P n1 i 1, j
ei,
j
P n1 i, j
bi,
j
P n1 i 1
di,
j
P n1 i. j 1
fi, j (形成五对
角矩阵)
重点辅导
23
2.2 抛物型方程的差分解法
u ( j 1, n) 2u ( j , n) u ( j 1, n) u 2h 4 ( j , n) u ( j , n) 2 2 4 h 2 t 4! x
n
(8)
0
Lu j
n
Lh, u j R j n
式中:
2 4 2 2 h 2 Rn u ( j , n ) u ( j , n ) O ( h ) j 4 2 4! x 2 t
(backward space difference) (backward time difference)
u n j
(3)一阶中心差分(central difference)
hu
n j
un 1 un
j 2
j
1 2
h
u
n j
uj
n
1 2
uj
n
1 2
1 n 1 un u j j
n
(22)
n+1 n
j-1
j
j+1
注意:
① 泰勒展开点在格边上,不是在结点上,但在格式中未出现格边量。 ② ③
O( 2 h2 ) ——全二阶精度。 1 在 ( j, n ) 点展开时,用到了周围6个结点上的量,该格式又称为六点格式。 2 Rj
2u idea:是将微分方程中的 2 项以 u ( x, t ) x
u j n1 u j n 1 2
u j 1n 2u j n u j -1n h2
0 (23)
(24)
u j n1 2r(u j 1n - 2u j n u j 1n ) u j n1
第四章 抛物型方程的有限差分方法
2 h 称为Du Fort -Frankel格式,仍为三层显式格式.
2
a
n 1 n 1 n un ( u u ) u j 1 j j j 1
0
截断误差: T x j , tn a u x j , tn u x j , tn 2 u x j h, tn u x j , tn u x j , tn u x j h, tn h2
1 2a G , k 0
0 4a cos kh 1 2a 1 1 0 4a cos kh 1 2a 1 2a 1 2a 0 1
2
1
4a cos kh 2a 1 G , k 的特征方程: 0 1 2a 1 2a
修正 Richardson:无条件不稳定显格式
Du Fort Frankel:无条件稳定的三层显格式. 但后者的相容性是有条件的.事实上, 显格式中,无条件相容和无条件稳定是无法同时成立的.
4 三层隐式格式
先考虑
n 1 n u u 3 j j n n 1 u u 1 j j 1 n1 n1 un 2 u u j 1 j j 1
引理1.1实系数二次方程 2 b c 0的根: c 1. 模 1 b 1 c, " "设1 , 2是方程两根,且 i 1 i 1, 2 证: c b 则12 c1 2 b a a 12 c c 1 2 1 1 2 b 1 c b 1 12 1 2 1 12 1 2 1 1 1 2 0, 若 1 2 0 1 12 1 2 1 1 1 2 0, 若1 2 0 b 1 c
抛物型方程有限差分法
抛物型方程有限差分法1. 简单差分法考虑一维模型热传导方程(1.1) )(22x f xua t u +∂∂=∂∂,T t ≤<0 其中a 为常数。
)(x f 是给定的连续函数。
(1.1)的定解问题分两类:第一,初值问题(Cauchy 问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件:(1.2) ()()x x u ϕ=0,, ∞<<∞-x第二,初边值问题(也称混合问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件:()13.1 ()()x x u ϕ=0,,l x l <<-及边值条件()23.1 ()()0,,0==t l u t u ,T t ≤≤0假定()x f 和()x ϕ在相应的区域光滑,并且于()0,0,()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
现在考虑边值问题(1.1),(1.3)的差分逼近 取 N l h =为空间步长,MT=τ为时间步长,其中N ,M 是自然数, jh x x j ==, ()N j ,,1,0Λ=; τk y y k ==, ()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。
其中 ()j i y x ,表示网格节点;h G 表示网格内点(位于开矩形G 中的网格节点)的集合; h G 表示位于闭矩形G 中的网格节点的集合;h Γ表示h G -h G 网格边界点的集合。
k j u 表示定义在网点()k i t x ,处的待求近似解,N j ≤≤0,M k ≤≤0。
注意到在节点()k i t x ,处的微商和差商之间的下列关系((,)kj k ju u x t t t ∂∂⎛⎫≡ ⎪∂∂⎝⎭):可得到以下几种最简差分格式 (一) 向前差分格式()24.1 ()j j j x u ϕϕ==0, k u 0=kN u =0其中1,,1,0-=N j Λ,1,,1,0-=M k Λ。
04有限差分法.ppt
n Rj
O t x
2
无条件稳定
2.一维混合问题
u 2u 2 0 t x u x ,0 F x u a, t t u b, t t
0 x b, t 0, 0
对于[a,b]区间的内点,可以构造以上各种格式。 如四点显式
例:驱动腔内的流体流动。
3.网格划分
x h y l xi ih
-----称为步长。
u x, y u i , j
xi , y j i, j
y j jl
4.差分格式 将u在(i,j)附近展成Taylor级数
ui 1, j ui , j ui 1, j ui , j 1 2u 1 3u u h 2 h 2 3 h 3 ... 2 x 3! x x i , j i, j i, j 1 2u 1 3u u h 2 h 2 3 h 3 ... 2 x 3! x x i , j i, j i, j
-----中心差分式
O h 表示具有二阶精度。
2
两Taylor展式相加
2u 1 ui 1, j 2ui , j ui 1, j O h 2 x 2 h2 i, j
抛物型方程差分法
2m
从而要求 4rsin2i2, 1im1
2m
a 1
易见,只要 r h 2 2 就可以保证数值格式稳定。 称为稳定性条件
对于非齐次方程、非零边界条件的情形,其稳定性 分析仿上,只是差分格式现在变成
u r k 1 A u r k b r k r 其中向量 b k 依赖于方程的右端项和边界条件。
u ( 0 ,tk )( tk ) ,u ( 1 ,tk )( tk ) , 0kn.
3.处理方程 u
t
2u ax2
f(xi, tk)中的偏导数
(xi,tk)
(xi,tk)
关于时间的一阶偏导数用向前差商近似,
u u(xi,tk1)u(xi,tk)
t (xi,tk)
r 12r
0
r O
O r
0
12r r
1r2ruuuum m kkM 12kk12Auuuum m kkM 12kk12
也可以简写成 u rk1A u rk ,从而有
u r k 1 A ( A u r k 1 ) L A k 1 u r 0
( xi , tk )
( xi , tk ) — 网格节点
用
u
k i
表示温度分布函
数 u( x , t ) 在点 ( xi , tk )
处的网格函数 , 相当于
x x i 1 x i x i 1
u( x , t ) 在该点的近似 .
2. 原方程弱化为节点处的离散方程
连续方程 离散方程
u 2u tax2f(x,t), 0x1 , 0tT
将数值解 u
k i
代替精确解 u( xi , tk )
一类非线性偏微分方程的数值求解
课例研究引言:偏微分方程可以用来描述真实世界的实际问题。
简单的拋物型偏微分方程即热传导方程有效地表征了物体内温度随着时间的演化过程与温度分布。
对于具有简单边界条件的偏微分方程,解析解可以通过分离变量法或拉普拉斯变化得到[1]。
由于问题本身的复杂性,非线性偏微分方程目前主要采用数值方法求解且没有统一的求解方法。
因此,针对非线性微分方程的特点选取合理的求解方法是十分重要的。
有限差分法是求解偏微分方程普遍采用的数值方法之一。
基于有限差分法,目前已有很多学者针对偏微分方程的数值求解展开了相关研究[2-4],如:二维波动方程的差分方法[5],以及有限差分法在求解一类非线性微分方程时的稳定性问题[6]。
1 一类非线性偏微分方程的化简设函数(),f x t ,其自变量为x 和t ,考虑下面的非线性偏微分方程()22,mf f f h x t t x x ∂∂∂ =− ∂∂∂ ,a x b <<,0t <<∞ (1)()(),0=f x x ϕ(2)()()1,=f a t g t ,()()2,=f b t g t(3)其中m 为幂指数。
上述的非线性偏微分方程若直接按照有限差分法进行离散求解会出现不收敛的情形。
由于该方程为非线性二阶偏微分方程,利用复合函数求导的关系,可以将上式右边进行化简。
化简后可以写为()11,1m f f h x t t m x x + ∂∂∂ =− ∂+∂∂(4)非线性微分方程由于收敛较为困难,目前普遍采用隐式方法求解。
在下面的计算中,将针对非线性微分方程采用显式求解。
为求解上述非线性偏微分方程,分别在时间和空间进行离散.假定变量x 的区间为a x b ≤≤,将x 划分为n 个网格,则i b ax ia n−=−,0,1,i n = (5)其中每一网格的宽度b ax n−∆=。
同理,可以对时间进行离散k t k t =∆,t ∆为时间步长,k 为时间步数。
利用时间的向前差分和空间的中心差分法可以将方程(4)进行离散()1,k kiik i f f f tt+−∂=∂∆ (6)()22,i x i x k k k i f f f x x+∆−∆−∂=∂∆ (7)其中i 表示节点编号。
抛物型方程的有限差分方法
抛物型方程的有限差分方法一,求解问题考虑一维非齐次热传导方程的定解问题22(,),0,0(,0)(),0(0,)(),(1,)(),0u ua f x t x l t T t xu x t x l u t t u t t t T ϕαβ∂∂-=<<<≤∂∂=≤≤==<≤......(1)..................(2) (3)其中α为正长数,(,)f x t ,()t ϕ,()t α,()t β为已知函数,(0)(0),(1)(0)ϕαϕβ==,式(2)为初值条件,(3)为边值条件。
二,网格剖分取空间步长/h l M =和时间步长/T N τ=,其中M 、N 都是整数。
用两族平行直线,(0,1,,)i x x ih i M ===和(0,1,,)k t t k i N τ===将矩形域{0;0}Gx l t T =≤≤≤≤分割成矩形网格,网格结点为(,)i k x t 。
以h G 表示网格内点集合,即位于开矩形G 的网点集合;h G 表示所有位于闭矩形G 的网点集合;h h G G Γ=-是网格界点集合。
其次,用ki u 表示定义在网点(,)i k x t 的函数,11,01i Mk N ≤≤-≤≤-。
用适当的差商代替方程(1)中相应的偏微商。
三, 差分格式 1, 向前差分 向前差分格式111202()(),11,01k kk k kiii i i ii i kki i i M u u u u u af hf f x u x u u i M k N ττϕϕ++---+=+====≤≤-≤≤-以2/ra h τ=为网比。
将上式改写为便于计算的形式,则得以下向量形式111(12)()(,)11,01k k k kii i i i k u r u r u u f x t i M k N τ+-+=-+++≤≤-≤≤-上式表示第k 层的值显示表示出来。
已知第k 层的值{|1}k i u i M ≤≤,则可以直接得到第k+1的值1{|1}k i u i M +≤≤。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下,才能求得定解问题解的解析式,对比较复杂的问题要找到解的解析表达式是困难的,因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法,只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数,得到相应的差分方程,通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想,并给出一些具体的数值实例.§1 差分方法的基本思想有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格,主要采用Taylor 级数展开等方法,在每个网格节点上用有限差分近似公式代替方程中的导数,从而建立以网格节点上的函数值为未知数的代数方程组.有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑,可分为中心格式和逆风格式.考虑时间因子的影响,差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式,主要是上述几种格式的组合,不同的组合构成不同的差分格式.泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式.首先考虑单变量函数()u x ,如图1把区域x 离散为一批结点,记0()(), =0,1,2,i i u x u x ih u i =+=图1 单变量函数离散化函数()u x 在点i x 处的泰勒展开式为23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''+=++++ (1)或23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''-=-+-+ (2)式(1)和(2)重新整理可得2()()()()()2!3!i i i i i u x h u x u x u x u x h h h '''''+-'=---(3)和2()()()()()2!3!i i i i i u x u x h u x u x u x h h h '''''--'=+++(4)于是给出在点i x 处函数u 的一阶导数的两个近似公式1()()()i i i ii u x h u x u u u x h h ++--'≈= (5)1()()()i i i i i u x u x h u u u x h h----'≈= (6)因为级数被截断,这两个近似公式肯定要产生误差,此误差与h 同阶,形式分别为()(), ,2()(), .2i i i i i i hE u O h x x h hE u O h x h x ξξξξ''=-=≤≤+''==-≤≤ 若把式(3)和(4)相加并求()i u x ',可得11()()()22i i i i i u x h u x h u u u x h h+-+---'≈= (7)其截断误差与2h 同阶,形式为22()(), ,6i i i h E u O h x h x h ξξ''=-=-≤≤+若把式(3)和(4)相减并求()i u x '',可得1122()2()()2()i i i i i i i u x h u x u x h u u u u x h h +-+-+--+''≈= (8)其截断误差与2h 同阶,其形式为22()(), ,12i i i h E u O h x h x h ξξ''=-=-≤≤+我们可继续用这种方式来推导更复杂的公式,类似的公式还有很多,这里不再一一列举.公式(5)、(6)分别称为一阶向前、向后差分格式,这两种格式具有一阶计算精度,公式(7)、(8)分别称为一阶、二阶中心差分格式,这两种格式具有二阶计算精度.图2 二维区域网格剖分上面的结果可直接推广使用于导出二元函数(,)u x y 的许多有限差分近似公式.如图7.2,把求解区域进行网格剖分,使12(,)(,), ,=0,1,2,i j ij u x y u ih jh u i j ==其中x 方向的网格间距为1,h y 方向的网格间距为2,h 整数i 和j 分别表示函数(,)u x y 沿x 坐标和y 坐标的位置.二元函数(,)u x y 对x 求偏导时y 保持不变,对y 求偏导时x 保持不变,根据向前差分公式(7.5)可以给出在点(,)i j x y 处函数(,)u x y 的一阶偏导数的两个近似公式1,,1(,)i j i j i ju x y u u xh +∂-≈∂ (9),1,2(,)i j i j i ju x y u u yh +∂-≈∂ (10)相类似地,根据二阶中心差分格式(8)可以得到函数(,)u x y 的二阶偏导数的近似公式21,,1,221(,)2i j i j i j i ju x y u u u x h +-∂-+≈∂ (11)2,1,,1222(,)2i j i j i j i j u x y u u u yh+-∂-+≈∂ (12)下面我们推导函数(,)u x y 的二阶混合偏导数2ux y∂∂∂在(,)i j x y 的有限差分表达式.根据一阶中心差分格式(7),112111,11,11,11,122121221,11,1(,)(,)(,)1()21 ()()222 i j i j i j i j i j i j i j i j i j i u x y u x y u x y O h x y h y y u u u u O h O h h h h u u u +-+++--+--+++-∂∂∂⎡⎤⎡⎤∂=-+⎢⎥⎢⎥∂∂∂∂⎣⎦⎣⎦--⎡⎤=-++⎢⎥⎣⎦--≈1,11,1124j i j u h h -+--+二维有限差分近似可以直接推广到三维空间或三维空间加一维时间的情形.定义1 当步长趋于零时,差分方程的截断误差趋于零,则称差分格式与微分方程是相容的.定义2 当步长趋于零时,差分方程的解收敛于微分方程的解,则称差分格式是收敛的. 定义3 当差分方程的解由于舍入误差的影响,所产生的偏差可以得到控制时,则称差分格式是稳定的.§2 抛物型方程的有限的差分法为了说明如何使用有限差分法来求解偏微分方程,本节我们给出以下几个数值实例.算例1 考虑一维非齐次热传导方程的初边值问题:2212(,), 01,01,(,0)(), 01,(0,)(), (1,)(), 0 1.u ua f x t x t t x u x q x x u t g t u t g t t ⎧∂∂=+<<<≤⎪∂∂⎪⎪=≤≤⎨⎪==<≤⎪⎪⎩(7.13),其中2,a =函数11(,)[cos()2sin()],22xf x t e t t =--+-初始条件1()sin,2xq x e =左、右边界条件分别为11()sin(),2g t t =-21()sin()2g t e t =-.该定解问题的解析解为1(,)sin(),(,)[0,1][0,1].2xu x t e t x t =-∈⨯将求解区域{(,)|,0}x t a x b t T Ω=≤≤≤≤进行网格剖分,[,]a b 作m 等分,[0,]T 作n 等分,记,,b a Th m nτ-==则 ,0,,0i k x a ih i M t k k n τ=+≤≤=≤≤对该问题建立如下向前差分格式:11122, 11, 11,k kk k k k i i i i i i u u u u u a f i m k n hτ+-+--+=+≤≤-≤≤-(14) (,0)(),1,i i u x q x i m =≤≤ (15) 12(,)(), (,)(),1.k k k k u a t g t u b t g t k n ==≤≤ (16)令2r ah τ=,差分格式(7.14)整理得111(12), 11, 1 1.k k k k k i i i i i u ru r u ru f i m k n τ+-+=+-++≤≤-≤≤- (17)显然时间在1k t +上的每个逼近值可独立地由k t 层上的值求出。
该格式为显格式,采用显格式时,应注意时间步长和空间步长的选取,当12r ≤时向前差分格式是稳定的.我们采用步长0.1h =和0.0025τ=,选取0.25t =时的数据进行比较,得到精确解与近似解的最大误差是9.00670e-005.表1 算例1 时节点处数值解、精确解和误差的绝对值(显式格式)数值解与精确解的比较见下图:图3 t =0.25精确解与近似解的比较 图4 t =0.25精确解与近似解的绝对误差对该问题还可建立如下向后差分格式:11122, 11, 11,k k k k k k i i i i i i u u u u u a f i m k n hτ--+--+=+≤≤-≤≤-(18) 差分格式(14)整理得111(12), 11, 1 1.k k k k k i i i i i ru r u ru u f i m k n τ--+-++-=+≤≤-≤≤- (19)显然时间在k t 层上的逼近值需要通过求解一个三对角线性方程组得到。
该格式为隐格式,该隐格式对于任意网格比r 都是稳定的.我们采用步长0.1h =和0.0025τ=,选取0.25t =时的数据进行比较,得到精确解与近似解的最大误差是1.274828e-005.数值解与精确解的比较结果如下图:图5 t =0.25精确解与近似解的比较 图6 t =0.25精确解与近似解的绝对误差算例2 考虑一维非线性Chaffee-Infante 方程的初边值问题:11(())(), (,)[0,1][0,1],(,0)(), [0,1],(,)(), (,)(), [0,1].u u k x f u x t t x x u x q x x u a t g t u b t g t t ∂∂∂⎧=+∈⨯⎪∂∂∂⎪=∈⎨⎪==∈⎪⎩(20) 其中扩散系数()1,k x =3(),f u u u =-+初始条件11()1),22q x =-+左右边界条件分别为:1113()1),224g t t =--+2113()1)224g t t =--+.该问题的解析解为113(,)1),(,)[,][0,].224u x t t x t a b T =--+∈⨯. 将求解区域进行剖分,方法同上,对该问题建立如下向前差分格式:1111122(),11,11,k kk k k k k k i i i i i i i i i i i u u k k u u u u u k f u i m k n h h hτ+++-+----+=++≤≤-≤≤- (,0)(),1,i i u x q x i m =≤≤12(,)(), (,)(),1.k k k k u a t g t u b t g t k n ==≤≤采用步长0.1h =和0.005,τ=我们选取0.5t =时的数据进行比较,得到精确解与近似解的最大误差是6.866286e-005.表2 算例2 时节点处数值解、精确解和误差的绝对值运行结果如下图,图7 t =0.5精确解与近似解的比较 图8 t =0.5精确解与近似解的绝对误差。