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

1、解偏微分方程的有限差分法摘要：本文主要利用有限差分法求解计算力学的系统数学模型，推导有限差分法的理论基础，并在此基础上，通过部分有限差分法求解偏微分方程的计算实例验证了推导的正确性和操作可行性。关键词：计算力学、偏微分方程、有限差分法abstract:this disseration mainly focuss on solving the mathematic model of computation mechanics with finite-difference method。the theoretical basis Of finite-difference is derived in

2、 the second part Of the dissertation，And then I use MATLAB to program the algorithms to second partKey words: computation mechanics，partial differential equations，finite-difference method引言1机械系统设计经常需要从机械角度进行结构设计和结构分析，这种分析的前提是建立工程问题的数学模型。将自然基本定律和原理应用于机械系统，得到构成系统数学模型的相关边界条件和初始条件的微分积分方程。解决这种数学模型的方法主要分为

3、分析方法和数值方法两种，据悉，分析方法的局限性在系统的边界条件和载荷条件复杂的情况下，往往找不到问题的解析或近似解决方案。另一方面，计算机技术的发展使计算更准确、更快。因此，对大多数工程问题进行数值解法研究具有更实用的价值。对于微分方程，主要分为差分法和积分法两种，本文主要讨论差分法。2有限差分法的理论基础2.1有限差分法的基本思想建立系统的数学模型时，我们面临的主要问题是微分积分方程的解法。基本思想是用有限未知量的离散差分方程近似替换连续变量的微分方程和有限解条件，并将差分方程的解作为微分方程解问题的近似解。通过将原始方程和边界条件的微分近似为“差”，将方程的积分近似为求和或机器求积公式，将

4、原始微分积分方程和边界条件转换为差分方程。求解偏微分方程的有限差分方法的步骤主要为：n区域是不连续的。也就是说，将给定偏微分方程的解析区域细分为称为网格中节点的有限格。用有限差分公式替换每个晶格的导数的n近似替换；可以把n近似解，也就是偏微分方程的解看作插值多项式及其微分法的替换过程。原则上，该方法仍然可以达到任意满意的计算精度。方程式的连续数值分析是因为您可以减少独立变数的离散值，或透过插补不连续点的函数值来近似。理论上，当网格阶段接近0时，差分方程的解必须收敛为正确的解，但由于机器字节的限制，网格阶段也不必减小到无穷大，差分方法的收敛或算法的稳定性非常重要。因此，使用有限差分法时，不仅要保

5、证精度，还要保证其收敛性。2.2系统微分方程的一般形式(1)大多数工程问题是二维问题，所以结果微分方程一般是偏微分方程，对于一维问题，与常微分方程类似的解是偏微分方程，所以这里只讨论非一般的偏微分方程。本文提出了二阶偏微分方程的一般形式，并对高阶偏微分进行了类似的推广。二阶偏微分方程的一般形式如下：Axx bxy cyy=f (x，y，x，y)其中是弹性体特定特性的物理量(连续函数)。如果a、b和c都是常数，则(1)表达式称为次线性，它有三种次线性表达式格式：N =B2-4ac 0时，称为椭圆方程。N =B2-4ac=0时，称为抛物线方程式。N =B2-4ac 0时，称为双曲方程。椭圆方程主要

6、用于处理热传导等稳态或静态问题。抛物线方程主要用于处理渗透、扩散等瞬态问题。双曲型方程主要用于处理弦振动、薄膜振动等振动问题。除了上述微分方程外，一般还需要提供以下三类解条件：n第一类边界条件(Dirichlet条件):| = (x，y)；n第二类边界条件(Neumann条件):n | = 1 (x，y)；n第三个边界条件(Robin条件):n (x，y) | = (x，y)；其中是解决域的边界，n是的单位外部法向矢量， (x，y) | 0。第二种和第三种类型的边界条件统称为衍生边界条件。2.3有限差分方程的数学基础2.3.1一元函数微分的差分公式x点的函数，大致可以表示为两个相邻点的函数值之

7、间的差值。在X=x0处，函数FX的“泰勒”卷展栏如下所示：Fx=n=01n！Fn(x0)x-x0n(2)=fx0 fxx-x0 12！Fxx-x02 13！F (3) xx-x03(3)对于单变量函数FX，用步骤x=h将a，b部分除以相等的间隔。x0=a，x1=x0 x=a h，x2=x1 x=a 2h，Xi=Xi-1 x=a ih，xn=b (n=b-ah)，这些与节点xi相邻的节点具有xi-h和Xi h，因此可以在点Xi中配置以下形式的扩展：Fxi-h=fxi-fxih 12！Fxih2 R2(x)(4)Fxi h=fxi fxih 12！Fxih2 R2(x)在格式(3)和样式(4)中，

8、可以获得：(5)n第一次正向差异：Fxi=fxi h-fxih(6)n第一次反向差异：Fxi=fxi-fxi-hh(7)n主要中心差：Fxi=fxi h-fxi-h2h(fi=f(xi)、样式(5)、(6)、(7)可以分别缩写为：(8)n第一次正向差异：Fi=fi 1-fih(9)n第一次反向差异：Fi=fi-fi-1h(10)n主要中心差：Fi=fi 1-fi-12h次要差异可用于型式(8)、型式(9)和型式(10)。(11)n次正向差异：Fi=fi 1-fih=fi 2-2fi 1 fih2(12)n次反向差异：Fi=fi-fi-1h=fi-2fi-1 fi-2h2(13)n次中心差：Fi

9、=fi 1-fi-12h=fi 2-2fi-24h 2差分公式(13)将中间节点上的一阶导数表示为2h的两个节点上的函数值，可以说是中间微分公式。型式(11)和型式(12)将一个端点的主要导出值表示为相邻的三节点函数值，也可以称为结束导出公式。必须指出，中点微分公式比端点微分公式精度高。因为前者反映节点两侧函数的变化，而后者仅反映节点一侧函数的变化。因此，我们必须总是尽可能地应用前者，只有在不能应用前者的情况下才应用后者。(14)然而，表达式(11)中的每个阶导数使用前向差将节点I中使用的主顺序转换为逆顺序差，同时，确保使用的节点不相邻(11)修改如下：fi=fi1-fih=fi1-fih-f

10、i-1hh=fi1-2fi-1 H2同样，根据上述推导过程，可以得到随机阶差公式。(15)N n次正向差异：Fi(n)=fi 1(n-1)-fi(n-1)h(16)n阶反向差异：Fi(n)=fi(n-1)-fi-1(n-1)h(17)n次中心差：Fi(n)=fi 1(n-1)-fi-1(n-1)2h在上述公式中，每个节点上先前主要微分的替代可能不一致，这表明最终公式在系数上可能存在差异，因为它可能是正向差异、反向差异或中心差异。当然，对于每个相邻节点，也可以通过扩展需要阶的泰勒，创建直接寻找每个阶导数的方程。将2.3.2微分方程转换为线性方程Ym(18)以椭圆微分方程为例，三类微分方程被转换为

11、代数方程，双曲线和抛物线方程按该顺序求解。2u=uxx uyy (2称为拉普拉斯运算符)，fx，y和g(x，y)是解决方案域的连续函数。假设解决方案面积为R=x，y 33600xa，0yb，ba=m/n，则解决方案面积为(n-1)(m-1)其中a=nh，b=mh如图1所示。记住Fi，j=f(xi，yj)，可以根据格式(14)获得以下内容：2u=uxx uyy=ui1、j-2ui、Jui-1、jh2ui、J1-2ui、Jui、j-1 H2=ui 1、j ui-1、j ui、j 1 ui、j-1-4ui、jh2o (H2)Yj 1YjYj-1Y1X1 xi-1 Xi Xi 1 xn图1 5点差公式

12、如图(2)所示，也称为五点差分公式的表达式(18)和表达式(13)可以分别得到正向差分公式(19)和反向差分公式(20)。(19)n向前差分2u=uxx uyy=ui 2、j-2ui 1、j ui、jh2ui、J2-2ui、J1 ui和jh2=ui2、Jui、j2ui、j-2ui1、j-2ui、J1 H2O (H2)YmYm(20)n后向差异2u=uxx uyy=ui、j-2ui-1、j ui-2、jh2ui、j-2ui、j-1 ui、j-2 H2=ui-2、j ui、j-2 2ui、j-2ui-1、j-2ui、j-1 h2oh2X1Xi-2Xi-1诗诗1诗2诗XnXnX1Y1Yj-2Yj-1

13、YjYjY1Yj 1Yj 2Ui-2、j2ui、jUi，j 2Ui，j 1图2前向差异(左)和后向差异(右)-2ui-1，j-2ui，j-1Ui，j-2-2ui、j 12ui、j-2ui 1，jUi 2，j-4ui、jUi，j-1Ui 1，jUi-1，j图3中心差、前差和后差的拉普拉斯算子表示(21)使用中心差分公式(18)，点x，y=(xi，yi)具有二次精度(Oh2)，因此可以近似地替换样式(18)，如下所示：2ui、j-ui 1、Jui-1、Jui、J1 ui、j-1-4ui和jh2根据椭圆方程的具体形式，可分为以下三种形式：n拉普拉斯方程：2u=0n Poison方程式：2u=g (x

14、，y)n Helmholtz方程式：2u FX，Yu=g (x，y)根据样式(21)，您可以产生三种不同形式的椭圆方程式：n拉普拉斯方程：2u=00=2ui、j-ui 1、Jui-1、Jui、J1 ui、j-1-4ui和jh2简化后得到拉普拉斯方程的计算公式：(22)Ui 1、j ui-1、j ui、j 1 ui、j-1-4ui、j=0(23)n泊松方程：2u=GX，yUi 1、Jui-1、Jui、J1 ui、j-1-4ui、j-h2gi、j=0n海尔墨水方程式：2u FX，Yu=g(x，y)(24)Ui 1、j ui-1、j ui、j 1 ui、j-1 (h2fi、j-4) ui、j-h2g

15、i、j=02.3.3建立有限差分方程方程式系统是根据型式(22)(24)产生的，但是必须知道其边界条件才能求解方程式。2.2表明，边界条件一般分为Dirichlet边界条件和派生边界条件两种，两种边界条件的有限差分方程的生成过程如下：N Dirichlet边界条件：| = (x，y)对于降阶条件，给出了边界上每个节点的函数计算公式，并直接替换节点值(xi，yi)以计算矩形区域的边界点，如下所示：U(x1，yj)=u1，j=(x1，yj)(1jm)(左边界)(25)U(xn，yj)=un，j=(xn，yj)(1jm)(右边界)U(xj，y1)=uj，1=(xj，y1)(1jn)(下边界)U(xj

16、，yn)=uj，n=(xj，yn)(1jn)(上边界)n衍生边界条件：u (x，y) n=0(26)例如，对于右侧边界点x=xn，根据Neumann条件，您可以：U (x，y) n=u (xn，yj) x=uxxn，yj=0对于拉普拉斯方程式，根据计算公式(22)，您可以对边界上的点(xn，yj)执行以下作业：(27)Un 1、j un-1、j un、j 1 un、j-1-4un、j=0(28)显然，上面的un 1、j是解决方案域之外的未知量。根据中心差分公式(10)，您可以获得：Uxxn、yj根据表达式(28)，可以得到近似表示。un 1、j-un-1、j和二次近似(27)可以使用以下公式：(29)