有限差分法

1第二章流体力学中的有限差分法2.1基本概念一、差分法：用离散点上函数值的差商代替导数如：取离散点（结点）nnxxxxx−1210LL结点函数值为nnyyyyy,,,,1210−LL用差商代替y(x)在结点xi的导数iiiiiixxyyxydxdy−−=⎟⎠⎞⎜⎝⎛ΔΔ≈⎟⎠⎞⎜⎝⎛++11向前差商2由此产生的误差称为截断误差：)()(2)()()(2)()(1)()()(2hOxyhxyxyxyhxyhxyhxyhxyhxydxdyxyiiiiiiiiiii=+′′=′−⎭⎬⎫⎩⎨⎧−⎥⎦⎤⎢⎣⎡+′′+′+=′−−+=⎟⎠⎞⎜⎝⎛−⎟⎠⎞⎜⎝⎛ΔΔLL步长iixxh−=+13若11−−−−=⎟⎠⎞⎜⎝⎛ΔΔ≈⎟⎠⎞⎜⎝⎛iiiiiixxyyxydxdy向后差商截断误差为)(hO。若1111−+−+−−=⎟⎠⎞⎜⎝⎛ΔΔ≈⎟⎠⎞⎜⎝⎛iiiiiixxyyxydxdy中心差商当结点为等步长时，即11−+−=−=iiiixxxxh，截断误差为)(2hO。4二阶差商代替二阶导数等步长时2111122222hyyyhhyyhyyxyxydxydiiiiiiiiii−+−++−=−−−=⎟⎟⎠⎞⎜⎜⎝⎛ΔΔ≈⎟⎠⎞⎜⎝⎛Δ′Δ≈⎟⎟⎠⎞⎜⎜⎝⎛截断误差为)(2hO。平面区域上，等距差分网格L,2,1,0,00±±=⎩⎨⎧Δ+=Δ+=jiyjyyxixxji。5二阶偏导数2,1,1),(222xuuuxujiijjijiΔ+−≈⎟⎟⎠⎞⎜⎜⎝⎛∂∂−+21,1,),(222yuuuyujiijjijiΔ+−≈⎟⎟⎠⎞⎜⎜⎝⎛∂∂−+例1：对拉普拉斯方程02222=∂∂+∂∂yuxu得到差分格式如下02221,1,2,1,1=Δ+−+Δ+−−+−+yuuuxuuujiijjijiijji或6()⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛ΔΔ+⎥⎥⎦⎤⎢⎢⎣⎡+⎟⎟⎠⎞⎜⎜⎝⎛ΔΔ++=−+−+21,1,2,1,112yxuuyxuuujijijijiij――五点格式，菱形格式。即：某结点的函数值是四个相邻结点的加权平均。正方形网格yxΔ=Δ()1,1,,1,141−+−++++=jijijijiijuuuuu(i,j+1)(i-1,j)(i,j)(i+1,j)(i,j-1)7二、差分法的基本问题（其它方法也存在类似问题）1.适定性：方程在一定的定解条件（初始条件，边界条件）下有唯一解，并且解连续性地依赖定解条件（解随定解条件的连续变化而连续变化）8(A)椭圆型方程：要求在闭域边界上给出边界条件。如：拉普拉斯方程02222=∂∂+∂∂yuxu第一类边界条件),(11yxfu=Γ第二类边界条件),(22yxfnu=∂∂Γ第三类边界条件),(),(33yxfuyxnu=⎥⎦⎤⎢⎣⎡α+∂∂Γ9(B)抛物型方程：要求初始条件和边界条件。如：二维扩散方程⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂α=∂∂2222yuxutu边界条件基本与拉普拉斯方程类似，为闭域边界条件；初始条件),()0,,(0yxuyxu=，可以从初始条件开始步进计算。又如：一维扩散方程22xutu∂∂α=∂∂)()0,(0xuxu=；)(),0(1tftu=，)(),(2tftLu=两端都要给定边界条件（双程坐标）。10(C)双曲型方程：适当的边界条件和初始条件，与波动传播的性质有关如：一维对流方程0=∂∂+∂∂xuctu)()0,(xfxu=解为)(),(ctxftxu−=，代表一个向右(c0时)或向左(c0时)传播的波形。必须在波形传来的一侧提供边界条件（单程坐标）。11不适定的例子：0=+xxttuu0)0,()0,(==xuxut拉普拉斯方程＋非闭域边界条件，解为0),(≡txu。然而，若定解条件为nxnxuxutsin1)0,(,0)0,(==，解为nxntntxusinsinh1),(=n→∞时，0)0,(,0)0,(→=xuxut，但∞→),(txu，不适定。122．差分格式的相容性：当步长趋于零时，差分格式应收敛于原微分方程（截断误差应趋于零）。例2：方程0=∂∂+∂∂=xuctuLu（算子xctL∂∂+∂∂=）),(ninitxuu=差分方程02111=Δ−+Δ−=−++xuuctuuuLnininininih截断误差),(2xtOLuuLTnihniΔΔ=−=→0，相容。13例3：扩散方程022=∂∂α−∂∂=xutuLu()()02111111=−−−α−τ−=−−++−+huuuuuuuLnininininininih截断误差L−⎟⎟⎠⎞⎜⎜⎝⎛∂∂τ+⎟⎟⎠⎞⎜⎜⎝⎛∂∂ατ=−=ninininihnitutuhtxLuuLT33222226),(若==τrh常数，0222≠⎟⎟⎠⎞⎜⎜⎝⎛∂∂α→niniturT，不相容。143.收敛性：当步长趋于零时，差分解应收敛于微分方程的真解。Lax等价定理：对于适定的线性微分方程的初值问题，若差分格式相容，则稳定性是收敛性的充分必要条件。4.稳定性：差分法计算中所产生的误差（舍入误差，参数误差等）随时间衰减或不增大，则称离散格式是稳定的，反之则是不稳定的。稳定性分析：谐波分析法，矩阵法等15例4：一维对流方程0=∂∂+∂∂xuctu(c=常数)(1)取中心差格式()02111=−+τ−=−++nmnmnmnmnmhuuhcuuuL即()nmnmnmnmuuhcuu1112−++−τ−=设计算到第n步时的累积误差nmnmnmuu差分法精确解计算值−=ε~反之nmnmnmuu+ε=~16则第n+1步的计算值()()()11111111122~~2~~++−+−+−++ε+=ε−ετ−ε+−τ−=−τ−=nmnmnmnmnmnmnmnmnmnmnmnmuhcuuhcuuuhcuu所以()nmnmnmnmhc1112−++ε−ετ−ε=ε(误差方程形式上与差分格式一致)17Fourier分析：)1(,)(−==ε∑∞−∞=ieAxikxkkk为波数(k次谐波)，Ak=相应的振幅。取误差为任意一个谐波分量，即ikmhkikxknmeAeAm==ε则()nmikhikhikmhknmeehceAλε=⎥⎦⎤⎢⎣⎡−τ−=ε−+211放大因子()khikhhiceehcikhikhsin1sin2121β−=τ−=−τ−=λ−(hc/τ=β)18◆用范数来表示计算误差的总体大小可取212||⎟⎠⎞⎜⎝⎛ε=ε∑mnmnm或||supnmmnmε=ε稳定性要求nmnmnmε≤ελ=ε+1，即1≤λ然而1sin1sin1||22≥β+=β−=λkhkhi所以中心差格式不满足稳定性要求（计算结果出现较大的波动）。19（2）向前差商格式()011=−+τ−++nmnmnmnmuuhcuu即()nmnmnmuuu111++β−β+=误差方程为()nmnmnm111++βε−εβ+=ε令ikmhknmeA=ε，则()nmikhnmeεβ−β+=ε+11放大因子ikheβ−β+=λ1不同k值，复数ikheβ−β+1位于复平面上圆心为)0,1(β+，半径为||β的轨迹圆上。20β0(c0)时，轨迹圆在|λ|=1之外，不满足|λ|≤1。β0(c0)，且|β|≤1时，|λ|≤1，格式稳定。(3)向后差商格式()011=−+τ−−+nmnmnmnmuuhcuu即()nmnmnmuuu111++β+β−=误差方程为()nmnmnm111++βε+εβ−=εβ0β0|λ|=121令ikmhknmeA=ε，得()nmikhnmeεβ+β−=ε+11放大因子ikheβ+β−=λ1β0(c0)时，|λ|≥1，不稳定。β0(c0)，且|β|≤1时，|λ|≤1，格式稳定。结论：c0，且|β|≤1时，后差商格式稳定c0，且|β|≤1时，前差商格式稳定22――迎风格式（上风格式，upwindscheme）⎪⎪⎩⎪⎪⎨⎧−−≈∂∂+−)0()0(11时时chuucchuucxucnmnmnmnm来流的上游决定下游。|β|≤1，即cxtΔ≤Δ――Courant条件即：≤Δt波传播Δx距离所用的时间，使umn+1的依赖区域不超出[xm-1,xm+1]。23例5：一维扩散方程022=∂∂α−∂∂xutu(1)用显格式022111=+−α−τ−=−++huuuuuuLnmnmnmnmnmnmh令2hrατ=，有())(21111nmnmnmnmuururu−++++−=误差方程为())(21111nmnmnmnmrr−++ε+ε+ε−=ε用谐波分析法得放大因子)cos1(21khr−−=λ要对所有k值满足1||≤λ，须1|41|≤−r，即2/10≤≤r所以，显格式的稳定性条件为0,2|2≥ααΔ≤Δ以及xt24(2)隐格式022111111=+−α−τ−=+−++++huuuuuuLnmnmnmnmnmnmh有())(2111111+−+++++=+nmnmnmnmuuruur误差方程()nmnmnmnmrrε=ε+ε−ε++−+++)(2111111用谐波分析法，令ikmhknmeA=ε+1得[][]nmnmnmikhikhkhreerrε=ε−+=ε+−+++−11)cos1(21)(21所以放大因子)cos1(211khr−+=λ显然，当r≥0（即α≥0）时，一定满足1||≤λ所以，隐格式是无条件稳定的。25◆抛物型方程差分格式的正系数准则：1+nmu表达式右边各项的系数应与1+nmu的系数同号，否则违背扩散过程的规律。该准则对例4中的一维对流方程也是适用的。262.2典型的偏微分方程的差分解一．椭圆型方程，以泊松方程为例Ω⊂−=∂∂+∂∂区域),(),,(2222yxyxfyuxu1.常用的差分离散格式(1)菱形五点格式ijjiijjijiijjifyuuuxuuu−=Δ+−+Δ+−−+−+21,1,2,1,122或27()()ijjijijijiijfyxyxuuyxyuuyxxu)(2)(2)(22222,1,12221,1,222Δ+ΔΔΔ=+Δ+ΔΔ−+Δ+ΔΔ−−+−+hyx=Δ=Δ时()ijijjijijijiijijhfuuuuuhuuL−=−+++=◊=−+−+411,1,,1,12截断误差O(h2)28(2)角点五点格式Lhuij=□uij()ijijjijijijifuuuuuh−=−+++=−−+++−++4211,11,11,11,12截断误差O(h2)(3)九点格式Lhuij=uij=(2◇uij+□uij)/3()[()]ijijjijijijijijijijifuuuuuuuuuh−=−+++++++=−−+++−++−+−+204611,11,11,11,11,1,,1,12截断误差O(h4)，f=0时截断误差可达O(h6)29以上各式用于区域内部结点，形成各内部结点的代数方程式。在边界结点上，其周围结点不可能全部在区域内部，需要利用边界条件来导出边界结点的代数方程式。2.边界条件的处理：(1)第一类边界条件),(),(1),(1yxfyxuyx=Γ⊂如图，由插值法得13111λ+λ+=uuuO即11113111),(1λλλ+=+−yxfuuO外P23λ2Δy1ΔxOλ1Δx内Δy430或由2、4两点的插值可得()22214221,1λλλ+=+−yxfuuO如果边界结点在实际区域之外，如图中P结点，由OPuuu)1(222λ−+λ=得22212),()11(λλyxfuuOP=−+31(2)第二类边界条件),(22yxfnu=∂∂Γ取),(2RRPyxfnu≈⎟⎠⎞⎜⎝⎛∂∂且αΔ−+αΔ−≈α⎟⎟⎠⎞⎜⎜⎝⎛∂∂+α⎟⎠⎞⎜⎝⎛∂∂=⎟⎠⎞⎜⎝⎛∂∂sincossincos43yuuxuuyuxunuPPPPP得⎟⎟⎠⎞⎜⎜⎝⎛Δ+Δ=Δ−Δ−yxyxfuyuxuRRPααααsincos),(sincos243nR3Pα4323.代数方程组的求解各内点菱形五点格式方程＋各边界点方程→代数方程组统一形式：ijnbnbijbuau=−∑其中：Ω⊂集合内部结点与边界结点的),(ji；nb代表与(i