有限差分法求解偏微分方程

有限差分法求解偏微分方程摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。关键词：计算力学，偏微分方程，有限差分法Abstract：Thisdissertationmainlyfocusesonsolvingthemathematicmodelofcomputationmechanicswithfinite-differencemethod.Thetheoreticalbasisoffinite-differenceisderivedinthesecondpartofthedissertation,andthenIuseMATLABtoprogramthealgorithmstosolvesomepartialdifferentialequationstoconfirmthecorrectnessofthederivationandthefeasibilityofthemethod.Keywords：ComputationMechanics,PartialDifferentialEquations,Finite-DifferenceMethod1引言机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，主要分为差分法和积分法两种，本论文主要讨论差分法。2有限差分法理论基础2.1有限差分法的基本思想当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的步骤主要有以下几步：区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格，这些离散点称作网格的节点；近似替代，即采用有限差分公式替代每一个格点的导数；逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小，那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此，在运用有限差分法时，除了要保证精度外，还必须要保证其收敛性。2.2系统微分方程的一般形式由于大多数工程问题都是二维问题，所以得到的微分方程一般都是偏微分方程，对于一维问题得到的是常微分方程，解法与偏微分方程类似，故为了不是一般性，这里只讨论偏微分方程。由于工程中高阶偏微分较少出现，所以本文仅仅给出二阶偏微分方程的一般形式，对于高阶的偏微分，可进行类似地推广。二阶偏微分方程的一般形式如下：Aϕ𝑥𝑥+𝐵ϕ𝑥𝑦+𝐶ϕ𝑦𝑦=𝑓(𝑥,𝑦,ϕ,ϕ𝑥,ϕ𝑦)其中，ϕ为弹性体上的某一特征物理量（连续函数）。当A、B、C都是常数时，（1）式称为准线性，有三种准线性方程形式：如果Δ=𝐵2−4𝐴𝐶0，则称为椭圆型方程；如果Δ=𝐵2−4𝐴𝐶=0，则称为抛物型方程；如果Δ=𝐵2−4𝐴𝐶0，则称为双曲型方程。椭圆型方程主要用来处理稳态或静态问题，如热传导等问题；抛物线方程主要用来处理瞬态问题，如渗透、扩散等问题；双曲型方程主要用来处理振动问题，如玄震动、薄膜震动等问题。除了上述微分方程外，必须给出定解条件，通常有如下三类：第一类边界条件（Dirichlet条件）：ϕ|Γ=𝜑(𝑥,𝑦)；第二类边界条件（Neumann条件）：∂ϕ∂n|Γ=𝜑1(𝑥,𝑦)；第三类边界条件（Robin条件）：[∂ϕ∂n+λ(x,y)ϕ]|Γ=𝜓(𝑥,𝑦)；其中，Γ为求解域Ω的边界，𝒏⃗⃗为Γ的单位外法矢，λ(x,y)|Γ≢0。第二类和第三类边界条件统称为导数边界条件。2.3有限差分方程的数学基础2.3.1一元函数导数的差分公式一个函数在x点上的导数，可以近似地用它所临近的两点上的函数值的差分（1）来表示。函数𝑓(𝑥)在𝑥=𝑥0处的泰勒展式如下：𝑓(𝑥)=∑1𝑛!𝑓(𝑛)(𝑥0)∞𝑛=0(𝑥−𝑥0)𝑛=𝑓(𝑥0)+𝑓′(𝑥)(𝑥−𝑥0)+12!𝑓′′(𝑥)(𝑥−𝑥0)2+13!𝑓(3)(𝑥)(𝑥−𝑥0)3+⋯对一个单变量函数𝑓(𝑥)，以步长∆𝑥=ℎ将[𝑎,𝑏]区间等距划分，我们得到一系列节点：𝑥0=𝑎，𝑥1=𝑥0+∆𝑥=𝑎+ℎ，𝑥2=𝑥1+∆𝑥=𝑎+2ℎ，⋯，𝑥𝑖=𝑥𝑖−1+∆𝑥=𝑎+𝑖ℎ，⋯，𝑥𝑛=𝑏（𝑛=𝑏−𝑎ℎ），然后求出𝑓(𝑥)在这些节点上的近似值。与节点𝑥𝑖相邻的节点有𝑥𝑖−ℎ和𝑥𝑖+ℎ，因此在点𝑥𝑖处可以构造如下形式的展开式：𝑓(𝑥𝑖−ℎ)=𝑓(𝑥𝑖)−𝑓′(𝑥𝑖)ℎ+12!𝑓′′(𝑥𝑖)ℎ2+𝑅2(𝑥)𝑓(𝑥𝑖+ℎ)=𝑓(𝑥𝑖)+𝑓′(𝑥𝑖)ℎ+12!𝑓′′(𝑥𝑖)ℎ2+𝑅2(𝑥)由式（3）和式（4）可得到：一阶向前差分：𝑓′(𝑥𝑖)=𝑓(𝑥𝑖+ℎ)−𝑓(𝑥𝑖)ℎ一阶向后差分：𝑓′(𝑥𝑖)=𝑓(𝑥𝑖)−𝑓(𝑥𝑖−ℎ)ℎ一阶中心差分：𝑓′(𝑥𝑖)=𝑓(𝑥𝑖+ℎ)−𝑓(𝑥𝑖−ℎ)2ℎ不妨，记𝑓𝑖=𝑓(𝑥𝑖)，则式（5）、（6）、（7）分别简写为：一阶向前差分：𝑓𝑖′=𝑓𝑖+1−𝑓𝑖ℎ一阶向后差分：𝑓𝑖′=𝑓𝑖−𝑓𝑖−1ℎ一阶中心差分：𝑓𝑖′=𝑓𝑖+1−𝑓𝑖−12ℎ（2）（3）（4）（5）（6）（7）（8）（9）（10）根据式（8）、式（9）和式（10），可得二阶差分：二阶向前差分：𝑓𝑖′′=𝑓𝑖+1′−𝑓𝑖′ℎ=𝑓𝑖+2−2𝑓𝑖+1+𝑓𝑖ℎ2二阶向后差分：𝑓𝑖′′=𝑓𝑖′−𝑓𝑖−1′ℎ=𝑓𝑖−2𝑓𝑖−1+𝑓𝑖−2ℎ2二阶中心差分：𝑓𝑖′′=𝑓𝑖+1′−𝑓𝑖−1′2ℎ=𝑓𝑖+2−2𝑓𝑖+𝑓𝑖−24ℎ2差分公式（13）是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。式（11）和式（12）是以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。但是，由于式（11）中的各阶导数均使用的是向前差分，导致用到的节点不相邻，同时为了均衡误差，将节点𝑖处用到的一阶差分换成向后差分，则式（11）修正为：𝑓𝑖′′=𝑓𝑖+1′−𝑓𝑖′ℎ=𝑓𝑖+1−𝑓𝑖ℎ−𝑓𝑖−𝑓𝑖−1ℎℎ=𝑓𝑖+1−2𝑓𝑖+𝑓𝑖−1ℎ2同理，根据上述推导过程，可得到任意阶的差分公式：n阶向前差分：𝑓𝑖(𝑛)=𝑓𝑖+1(𝑛−1)−𝑓𝑖(𝑛−1)ℎn阶向后差分：𝑓𝑖(𝑛)=𝑓𝑖(𝑛−1)−𝑓𝑖−1(𝑛−1)ℎn阶中心差分：𝑓𝑖(𝑛)=𝑓𝑖+1(𝑛−1)−𝑓𝑖−1(𝑛−1)2ℎ说明，上述公式中各节点处前一阶导数的代入可能存在不一致，可能是向前差分、向后差分或者中心差分，从而使最终的公式在系数上存在差别。当然，也（11）（12）（13）（15）（16）（17）（14）可以对各相邻节点进行需要阶数的泰勒展开，从而建立方程组直接求各阶导数。2.3.2微分方程转化为线性方程由于三种类型的微分方程解法类似，故这里仅以椭圆型微分方程为例，将微分方程转化为代数方程，对于双曲型和抛物型方程依次类推即可。不妨记：∇2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦（∇2称为拉普拉斯算子），𝑓(𝑥,𝑦)和𝑔(𝑥,𝑦)是求解域上的连续函数。假设求解区域为：𝑅={(𝑥,𝑦):0≤𝑥≤𝑎,0≤𝑦≤𝑏,𝑏𝑎=𝑚/𝑛}，将求解区域划分成(𝑛−1)×(𝑚−1)个网格，其中：𝑎=𝑛ℎ,𝑏=𝑚ℎ，如图1所示。记𝑓𝑖,𝑗=𝑓(𝑥𝑖,𝑦𝑗)，则根据式（14）可得到：∇2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦=𝑢𝑖+1,j−2𝑢𝑖,𝑗+𝑢𝑖−1,𝑗ℎ2+𝑢𝑖,𝑗+1−2𝑢𝑖,𝑗+𝑢𝑖,𝑗−1ℎ2=𝑢𝑖+1,j+𝑢𝑖−1,𝑗+𝑢𝑖,𝑗+1+𝑢𝑖,𝑗−1−4𝑢𝑖,𝑗ℎ2+𝑂(ℎ2)𝑥1𝑥𝑖−1𝑥𝑖𝑥𝑖+1𝑥𝑛图1五点差分公式式（18）也称为五点差分公式，同理根据式（12）和式（13）可分别得到向前差分公式（19）和向后差分公式（20），如图（2所示）。向前差分∇2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦=𝑢𝑖+2,j−2𝑢𝑖+1,𝑗+𝑢𝑖,𝑗ℎ2+𝑢𝑖,𝑗+2−2𝑢𝑖,𝑗+1+𝑢𝑖,𝑗ℎ2=𝑢𝑖+2,j+𝑢𝑖,𝑗+2+2𝑢𝑖,𝑗−2𝑢𝑖+1,𝑗−2𝑢𝑖,𝑗+1ℎ2+𝑂(ℎ2)（18）（19）𝑦𝑚𝑦𝑗+1𝑦𝑗𝑦𝑗−1𝑦1向后差分∇2𝑢=𝑢𝑥𝑥+𝑢𝑦𝑦=𝑢𝑖,j−2𝑢𝑖−1,𝑗+𝑢𝑖−2,𝑗ℎ2+𝑢𝑖,𝑗−2𝑢𝑖,𝑗−1+𝑢𝑖,𝑗−2ℎ2=𝑢𝑖−2,j+𝑢𝑖,𝑗−2+2𝑢𝑖,𝑗−2𝑢𝑖−1,𝑗−2𝑢𝑖,𝑗−1ℎ2+𝑂(ℎ2)图2向前差分（左）和向后差分（右）图3中心差分、向前差分和向后差分的拉普拉斯算子表示利用中心差分公式（18），由于式（18）在点(𝑥,𝑦)=(𝑥𝑖,𝑦𝑖)处具有二阶精度（𝑂(ℎ2)），所以式（18）可近似改写成下式：∇2𝑢𝑖,𝑗≈𝑢𝑖+1,j+𝑢𝑖−1,𝑗+𝑢𝑖,𝑗+1+𝑢𝑖,𝑗−1−4𝑢𝑖,𝑗ℎ2根据椭圆方程的具体形式可以将其分为以下三种形式：拉普拉斯（Laplace）方程：∇2𝑢=0泊松（Poison）方程：∇2𝑢=𝑔(𝑥,𝑦)赫耳墨次（Helmholtz）方程：∇2𝑢+𝑓(𝑥,𝑦)𝑢=𝑔(𝑥,𝑦)根据式（21），可建立三种不同形式椭圆方程的代数方程如下：拉普拉斯方程：∇2𝑢=00=∇2𝑢𝑖,𝑗≈𝑢𝑖+1,j+𝑢𝑖−1,𝑗+𝑢𝑖,𝑗+1+𝑢𝑖,𝑗−1−4𝑢𝑖,𝑗ℎ2化简后得到拉普拉斯方程的计算公式：（20）𝑦𝑚𝑦1𝑦𝑗𝑦𝑗+1𝑦𝑗+2𝑦𝑚𝑦𝑗𝑦𝑗−2𝑦𝑗−1𝑦1𝑥1𝑥𝑛𝑥1𝑥𝑖𝑥𝑛𝑥𝑖+1𝑥𝑖+2𝑥𝑖−1𝑥𝑖𝑥𝑖−2𝑢𝑖−1,𝑗−4𝑢𝑖,𝑗𝑢𝑖+1,𝑗𝑢𝑖,𝑗+1𝑢𝑖,𝑗−1−2𝑢𝑖,𝑗+1𝑢𝑖,𝑗+2−2𝑢𝑖+1,𝑗𝑢𝑖+2,𝑗2𝑢𝑖,𝑗2𝑢𝑖,