有限差分法求解偏微分方程M A T L A B
南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 115104000545 成绩： 有限差分法求解偏微分方程
一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程： 22(,)()u u f x t t x αα??-=??其中为常数 具体求解的偏微分方程如下： 22001 (,0)sin()(0,)(1,)00 u u x t x u x x u t u t t π???-=≤≤?????? =??? ==≥??? 2.推导五种差分格式、截断误差并分析其稳定性； 3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析； 4.结论及完成本次实验报告的感想。 二、推导几种差分格式的过程： 有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2 100000000()()()()()()()......()(()) 1!2!! n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1) 求解区域的网格划分步长参数如下：
一、差分的概念 二、差分方程的概念 三、三、常系数线性差分方程解的结构 常系数线性差分方程解的结构第六节第六节 差分与差分方程的概念 常系数线性差分方程解的结构 常系数线性差分方程解的结构第九章
一、差分的概念 1.差分的定义 . Δ,)1()()1()0(: ).(111210x x x x x x x y y y y y y y y y y y x f x f f f x x f y ?=?+=+++也称为一阶差分，记为的差分，为函数称函数的改变量，，，，，将之简记为 ，，，，，列函数值可以排成一个数取非负整数时，当设函数… …??
x x x x x x x x x x x y y y y y y y y y y y y x f y +?=???=?===++++++1211212 2)()()(Δ)Δ(ΔΔ,)(即 差分的一阶差分的的二阶差分为函数函数. 以上的差分高阶差分：二阶及二阶) (),(3 4 2 3 x x x x y y y y ??=???=?差分：同样可定义三阶、四阶??
例1 1 求 求(),(),(2 3222x x x ???. 解 ，则 设2 x y =1 2)1()(2 22+=?+=?=?x x x x y x []2 )12(1)1(2) 12()(2 2 2 =+?++=+?=?=?x x x x y x 0 22)(2 3 3 =?=?=?x y x
解例2 2 求下列函数的差分求下列函数的差分求下列函数的差分 y x y a sin )2(; log )1(== ); 1 1(log log )1(log )1(1x x x y y y a a a x x x +=?+=?=?+. 2 sin )21 (cos 2sin )1(sin Δ)2(a x a ax x a y x ?+=?+=
课程论文
有限差分法在微分方程中的应用 本学期学习了《微分方程数值解》，本书中有限差分法给我留下的印象比较深刻，下边说说自己在方面的一点理解，请老师指正。 1.有限差分法的基本思想： 当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。 2.有限差分法求解偏微分方程的步骤： 区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格，这些离散点称作网格的节点； 近似替代，即采用有限差分公式替代每一个格点的导数。 逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。 从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小，那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此，在运用有限差分法时，除了要保证精度外，还必须要保证其收敛性。 3.构造差分法的几种形式： 主要草用的是泰勒级数展开的方法。其基本差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等。其中前两种形式为一阶计算精度，后一种为二阶计算精度。
1、常系数线性差分方程的解 方程（ 8）其中为常数，称方程（8）为常系数线性方程。 又称方程（9） 为方程（8）对应的齐次方程。 如果（9）有形如的解，带入方程中可得： （10） 称方程（10）为方程（8）、（9）的特征方程。 显然，如果能求出（10）的根，则可以得到（9）的解。 基本结果如下： （1）若（10）有k个不同的实根，则（9）有通解： ， （2）若（10）有m重根，则通解中有构成项： （3）若（10）有一对单复根，令：，，则（9）的通解中有构成项： （4）若有m 重复根：，，则（9）的通项中有成项：
综上所述，由于方程（10）恰有k 个根，从而构成方程 （9）的通解中必有k个独立的任意常数。通解可记为： 如果能得到方程（8）的一个特解：，则（8）必有通解： + （11） （1）的特解可通过待定系数法来确定。 例如：如果为n 的多项式，则当b不是特征根时，可设成形如形式的特解，其中为m次多项式；如果b是r重根时，可设特解：，将其代入（8）中确定出系数即可。 2、差分方程的z变换解法 对差分方程两边关于取Z变换，利用的Z 变换F（z）来表示出的Z变换，然后通过解代数方程求出F（z），并把F(z)在z=0的解析圆环域中展开成洛朗级数，其系数就是所要求的 例1设差分方程，求 解：解法1：特征方程为，有根： 故：为方程的解。 由条件得： 解法2：设F（z）=Z(),方程两边取变换可得：
由条件得 由F（z）在中解析，有 所以， 3、二阶线性差分方程组 设，，形成向量方程组 （12）则 （13）（13）即为（12）的解。 为了具体求出解（13），需要求出，这可以用高等代数的方法计算。常用的方法有： （1）如果A为正规矩阵，则A必可相似于对角矩阵，对角线上的元素就是A的特征值，相似变换矩阵由A的特征向量构成：。 （2）将A 分解成为列向量，则有 从而，
z 二阶线性齐次差分方程012=++++n n n cx bx ax 的特征根法求解： 令形式解 ,代入方程得特征方程: , 根: n n x λ=02=++c b a λλ(1) βα,为实根, 对应有解： 和 ; n n x α=)1(n n x β=)2((2) αα,为重根, 对应有解： 和n n x α=) 1(1) 2(lim ?→=??=n n n n n x αα βαβαβ ，或者 n n n x α=)2((3) , ?βαλi e r i ±?=±=()()??λ?λn i n e e e x r n i r n n n n sin cos ln ln ln ±====±?, 对应有解： 和. ?n e x r n n cos ln )1(=?n e x r n n sin ln )2(=(4) 关于解的结构理论与线性微分方程类似，由此得一般解: )2(2)1(1n n n x c x c x +=1. (98) 求差分方程的一般解。 (n y y n n 51021=++()72 51255?+?=n C y n n ) 解：齐次方程的通解为，设非齐次方程的特解为：()n n C y 5?=b an y n +=~，代入求。 b a ,2. 斐波拉契数( ???==+=++11012x x x x x n n n ??? ???????????????????????+=++1125125151n n n x ) 3． 银行实行贷款购房业务，A 贷元，月利r ，n 个月本利还清，在这个月内按复利计息，每月连本带息还n x 元。 (1) 求的关系； (2) 记个月的平均利息(r n A f x ,,=)n n A x n v ?=，求r v n ∞→lim . 设第i 个月欠元，则 i A (),101???=?+=?A A x r A A i i 齐次方程的通解为 ();1n n r C A +=非齐次方程的特解为r x A n =~； 非齐次方程的通解为：();1r x r C A n n ++= 代入初始条件得非齐次方程的特解为()();111r r x r A A n n n ?+?+= 0=n A 得x 值。。。。。
差分方程的解法分析及MATLAB 实现（程序） 摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言 线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域 法[1]. 1 迭代法 例1 已知离散系统的差分方程为)1(3 1)()2(81)1(43)(-+=-+--n x n x n y n y n y ,激励信号为)()4 3()(n u n x n =,初始状态为21)2(4)1(=-=-y y ，.求系统响应. 根据激励信号和初始状态,手工依次迭代可算出24 59)1(,25)0(==y y . 利用MATLAB 中的filter 函数实现迭代过程的m 程序如下: clc;clear;format compact; a=[1,-3/4,1/8],b=[1,1/3,0], %输入差分方程系数向量,不足补0对齐 n=0:10;xn=(3/4).^n, %输入激励信号 zx=[0,0],zy=[4,12], %输入初始状态 zi=filtic(b,a,zy,zx),%计算等效初始条件 [yn,zf]=filter(b,a,xn,zi),%迭代计算输出和后段等效初始条件 2 时域经典法 用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项,根据自由项形 式求特解；然后根据边界条件求完全解[3].用时域经典法求解例1的基本步骤如下. （1）求齐次解.特征方程为081432=+-αα,可算出4 1 , 2121==αα.高阶特征根可用MATLAB 的roots 函数计算.齐次解为. 0 , )4 1()21()(21≥+=n C C n y n n h （2）求方程的特解.将)()4 3()(n u n x n =代入差分方程右端得自由项为 ?????≥?==-?+-1,)4 3(9130 ,1)1()43(31)()43(1n n n u n u n n n 当1≥n 时,特解可设为n p D n y )4 3()(=,代入差分方程求得213=D . （3）利用边界条件求完全解.当n =0时迭代求出25)0(=y ,当n ≥1时,完全解的形式为 ,)4 3(213 )41()21()(21n n n C C n y ?++=选择求完全解系数的边界条件可参考文[4]选)1(),0(-y y .根据边界条件求得35,31721=-=C C .注意完全解的表达式只适于特解成立的n 取值范围,其他点要用 )(n δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为
南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 成绩： 有限差分法求解偏微分方程
一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程： 22(,)()u u f x t t x αα??-=??其中为常数 具体求解的偏微分方程如下： 22001 (,0)sin()(0,)(1,)00 u u x t x u x x u t u t t π???-=≤≤?????? =??? ==≥??? 2.推导五种差分格式、截断误差并分析其稳定性； 3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析； 4.结论及完成本次实验报告的感想。 二、推导几种差分格式的过程： 有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2100000000()()()()()()()......()(()) 1!2!! n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1) 求解区域的网格划分步长参数如下：
使用有限差分方法解边值问题： 由两点边值问题的一般形式： 根据差分方程： 当网格划分均匀，即有，化简差分方程： 代入再次化简： 用方程组展开写成矩阵形式： MATLAB编程：
运行后算出的结果：0 0.00376645934479969 0.00752341210586145 0.0112613555020809 0.0149707943560995 0.0186422448923756 0.0222662385306948 0.0258333256736017 0.0293340794862392 0.0327590996670822 0.0360990162080584 0.0393444931425513 0.0424862322797872 0.0455149769241112 0.0484215155776656 0.0511966856249889 0.0538313769980622 0.0563165358203363 0.0586431680282822 0.0608023429690169
0.0627851969725639 0.0645829368973219 0.0661868436473210 0.0675882756598612 0.0687786723621374 0.0697495575954688 0.0704925430057619 0.0709993313988528 0.0712617200593841 0.0712716040318917 0.0710209793627865 0.0705019463019362 0.0697067124625652 0.0686275959382091 0.0672570283754778 0.0655875580013963 0.0636118526041142 0.0613227024657904 0.0587130232464804 0.0557758588178718 0.0525043840457360 0.0488919075199819 0.0449318742312199 0.0406178681927653 0.0359436150070336 0.0309029843752992 0.0254899925498146 0.0196988047273101 0.0135237373829146 0.00695926054356603 0 与精确解比较：
--以有限差分法为例偏微分方程数值求解 1. 偏微分方程求解问题的描述 教材P653[12.1.1]椭圆型 教材P653[12.1.2] 教材P664[12.2.1]双曲型 教材P665[12.2.4]拉普拉斯泊松 对流 波动 教材P684[12.3.1]抛物型 教材P685[12.3.6]扩散 对流扩散 教材P686[12.3.8]二维扩散 教材P678[12.2.23]二维对流
??????????????????????≥≤≤==≥≤≤==≤≤=>≥≤≤≤≤???? ????+??=??0,0, ),(),,(),(),0,(0,0,),(),,(),(),,0(,0,),()0,,(0,0 , 0 , 0 21212222t L x t x v t L x u t x v t x u t L y t y t y L u t y t y u L y x y x y x u b t L y L x y u x u b t u μμ?Ω 求解域初值条件边值条件) ,,(t y x u 未知函数
????? ? ????????????????????≥<<-==≥<<==≥≤≤-==≥≤≤==≤≤==≤≤≤≤≤≤???? ????+??=??0 , 50 , sin 255sin ),(),5,(0 , 50 , 0),(),0,(0 , 50 , 5sin sin 25),(),,5(0 , 50 , 0),(),,0(5,0,0),()0,,( 10000 , 50 , 50 001.022********t x x x t x v t x u t x t x v t x u t y y y t y t y u t y t y t y u y x y x y x u t y x y u x u t u μμ?Ω 求解域初值条件边值条件以具体问题为例演示具体的求解过程) ,,(t y x u 未知函数
例题：已知差分方程51 (2)(1)()(+1)+0.5()66 x k x k x k r k r k +-++=，其中r (k )=1，k ≥0，x (0)=1， x (1)=2。 (1) 试由迭代法求其全解的前5项； (2) 分别由古典法求其零输入解、零状态解，以及全解； (3) 用Z 变换法求解差分方程。 解：注：解题过程中出现的下标“zi ”和“zs ”分别表示零输入条件和零状态条件。 1. 迭代法 题目中给出的条件仅仅是零输入初始条件，进行迭代求解时的初始条件应该是全解初始条件。 (1) 零输入初始条件 本题已给出零输入时的两个初始条件x zi (0)=1，x zi (1)=2。 (2) 零状态初始条件 取k =-2时，则51 (0)(1)(2)(1)0.5(2)66x x x r r --+-=-+-，得x zs (0)=0； 取k =-1 时，则51 (1)(0)(1)(0)0.5(1)66 x x x r r -+-=+-，求得x zs (1)=1。 (3) 全解初始条件 x (0)= x zi (0)+ x zs (0)=1； x (1)= x zi (1)+ x zs (1)=3。 (4) 根据求出的全解x (0)和x (1)，利用迭代法求解 取k =0时，则51(2)(1)(0)(1)0.5(0)66x x x r r -+=+，求得23(2)6x =； 取k =1时，则51(3)(2)(1)(2)0.5(1)66x x x r r -+=+，求得151 (3)36x =； 取k =2时，则51(4)(3)(2)(3)0.5(2)66x x x r r -+=+，求得941 (4)216 x =。 2. 古典法 (1) 零输入解 令输入为零，则得齐次方程 51 (2)(1)()066 x k x k x k +-++= (a) 根据差分方程定义的算子()()n d x k x k n =+，可得它的特征方程251 066 d d -+= 求得特征根为： 112d = ，21 3 d =
有限差分法 有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散 点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函 数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差 分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便 可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原 微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和 计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分 格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格 式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过 程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致 差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以 控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能 任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是 数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的 微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用 待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法 将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从 而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数 问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分 的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目 前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分 方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
第三节 差分方程常用解法与性质分析 1、常系数线性差分方程的解 方程)(...110n b x a x a x a n k k n k n =+++-++ （ 8） 其中k a a a ,...,,10为常数，称方程（8）为常系数线性方程。 又称方程0...110=+++-++n k k n k n x a x a x a （9） 为方程（8）对应的齐次方程。 如果（9）有形如 n n x λ=的解，带入方程中可得： 0 ...1110=++++--k k k k a a a a λλλ （10） 称方程（10）为方程（8）、（9）的特征方程。 显然，如果能求出（10）的根，则可以得到（9）的解。 基本结果如下： （1） 若（10）有k 个不同的实根，则（9）有通解： n k k n n n c c c x λλλ+++=...2211， （2） 若（10）有m 重根λ，则通解中有构成项： n m m n c n c c λ )...(121----+++
（3）若（10）有一对单复根 βαλi ±=，令：?ρλi e ±=， αβ?βαρarctan ,22=+=，则（9）的通解中有构成项： n c n c n n ?ρ?ρsin cos 21--+ （4） 若有m 重复根：βαλi ±=，φρλi e ±=，则（9）的通项中有成 项： n n c n c c n n c n c c n m m m m n m m ?ρ?ρsin )...(cos )...(1221121---++---+++++++ 综上所述，由于方程（10）恰有k 个根，从而构成方程 （9）的通解中必有k 个独立的任意常数。通解可记为：-n x 如果能得到方程（8）的一个特解：*n x ，则（8）必有通解： =n x -n x +* n x （11） （1） 的特解可通过待定系数法来确定。 例如：如果)(),()(n p n p b n b m m n =为n 的多项式，则当b 不是特征 根时，可设成形如)(n q b m n 形式的特解，其中)(n q m 为m 次多项式；如 果b 是r 重根时，可设特解：r n n b )(n q m ，将其代入（8）中确定出系 数即可。
有限差分法解偏微分方程综述 绪论 有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 从差分的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。 目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式， 目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。 有限差分法求解偏微分方程 在采用数值计算方法求解偏微分方程时，若将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。有限差分法求解偏微分方程的步骤如下： 1、区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格； 2、近似替代，即采用有限差分公式替代每一个格点的导数； 3、逼近求解。换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程 有限差分法的应用 抛物型方程的差分方法 1. 简单差分法
有限差分法求解偏微分方程 摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。 关键词：计算力学，偏微分方程，有限差分法 Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method. Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Method
1 引言 机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。 求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，主要分为差分法和积分法两种，本论文主要讨论差分法。 2 有限差分法理论基础 2.1 有限差分法的基本思想 当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的步骤主要有以下几步： ?区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格， 这些离散点称作网格的节点； ?近似替代，即采用有限差分公式替代每一个格点的导数； ?逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代 替偏微分方程的解的过程。 从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小，
线性差分方程 在数学大家庭中，线性差分方程是重要的一员。如离散状态下的数学建模一般会产生差分方程，将微分方程离散化仍然会产生差分方程。线性差分方程与线性方程组及线性常微分方程有许多相似的性质，下面让我们讨论其解性质。 考虑如下一般的线性差分方程 （1） ()()()110k m k k m k m m a m u a m u a m u b +-+-+++= 0,1,2,m = 称为k 阶线性差分方程，其中(),j m a m b 为给定的关于m 的函数，并且()()00k a m a m ≠。当0m b =时方程 （2） ()()()1100k m k k m k m a m u a m u a m u +-+-+++= 0,1,2,m = 称为齐次差分方程。如果整变量j 的函数j u 使（1）和（2）成为等式，则称j u 为相应方程的解。容易看出，如果021,,k u u u - （称为初始值）已给定，则由（1）或（2）可以逐次地定出(),1,j u j k k =+ 。 与线性方程组及线性常微分方程类似，对上述差分方程有 命题1 如果() () () 12,,,k m m m u u u 是齐次差分方程（2）的解，则它们的任意组合 （3） () () () 1212k m m m k m v c u c u c u =+++ 也是（2）的解，其中()1,2,,j c j k = 为常数。 命题2 设() () () 12,,,k m m m u u u 是k 阶齐次差分方程（2）的解，且行列式 （4） () ()() ()()()()()()12111122 22 120k k k k k k u u u u u u u u u ≠ 则齐次差分方程（2）的任何解均可表成（3）的形式，此时称（3）为（2）的通解。 如果() () () 12,,,k m m m u u u 满足条件（4），则称() () () 12,,,k m m m u u u 线性无关，故命题2可叙述
第八节 二阶常系数线性差分方程 二阶常系数齐次线性差分方程的一般形式 12= ++++x x x by ay y （1） 二阶常系数非齐次线性差分方程的一般形式 ) x f by ay y x x x =++++12（2） 一、二阶常系数齐次线性差分方程的求解 令x x y λ= 代入方程（1）得特征方程 02=++b a λλ （3） 根据特征方程（3）的根的三种情形写出通解 (1)第一种情形：特征方程（3）有两个不同的 实根 21λλ≠ ，通解为 () 为任意常数212 211,C C C C y x x x λλ+= (2)第二种情形：特征方程（3）有两个相同的 实根 λλλ==21 ，通解为 ()() 为任意常数2121,C C x C C y x x λ+= (2)第三种情形：特征方程（3）有一对共轭 复根 i βαλ±=2,1 ，通解为 () () 为任意常数2121,sin cos C C x C x C r y x x θθ+= 其中 () πβα β θβα<<>=+=0, 0arctan ,22r r 例1 求差分 方程 0612=--++x x x y y y 的通解。 解 特征方程 062 =--λλ 的根为 2,321-==λλ 原方程的通解为
()()为任意常数2121,, 23C C C C y x x x -+= 例2 求差分 方程 04412=+-++x x x y y y 的通解。 解 特征方程 0442 =+-λλ 的根为 221==λλ 原方程的通解为 ()()为任意常数2121,, 2C C x C C y x x += 例3 求差分方程 041 2=++x x y y 的通解。 解 特征方程 0412=+λ 的根为 i 21 2,1±=λ 原方程的通解 () 为任意常数2121,, 2sin 2cos 21C C x C x C y x x ??? ? ? +??? ??=ππ 例4 求差分 方程 016412=+-++x x x y y y 的通解。 解 特征方程 01642 =+-λλ 的根为 i 3222,1±=λ 原方程的通解 () 为任意常数2121,, 3sin 3cos 4C C x C x C y x x ??? ? ? +=ππ 二、二阶常系数非齐次线性差分方程的求解 ) x f by ay y x x x =++++12（2） 由上节定理3知道，差分方程（2）的通解应由对应齐次差分方程的通解 （前面已学过）和非齐次差分方程的特解两部分组成。我们只学习后部分。 二阶常系数非齐次线性差分方程的特解求法——待定系数法 1. 非齐 次项 ()()x P x f n = 型 （1）1不是特征方程的根，即1+a +b ≠0, 设
有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 从差分的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差 分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
计算机与信息工程学院验证性实验报告 一、实验目的 1．学习并掌握系统的差分方程表示方法以及差分方程的相关概念。 2．熟练使用filter 函数对差分方程进行数值求解。 3．掌握差分方程的求解及MATLAB 实现方法。 二、实验原理及方法 1．一LTI 系统可以用一个线性常系数差分方程表示： ()()N M k m k m a y n k b x n m ==-=-∑∑，任意n 如果N a ≠0，那么这个差分方程就是N 阶的，已知系统的输入序列，用这个方程可以根据当前输入x(n)和以前M 点的输入x(n-m ),…,x(n-1),以及以前N 点的输出y(n-N),…,y(n -1)来计算当前输出y(n)。在实际中这个方程在时间上是从n =-∞到n =+∞朝前计算的，因此该方程的另一种形式是： ()()()M N m k m k y n b x n m a y n k ===---∑∑ 方程的解能以下面形式求得：()()()H p y n y n y n =+分别为方程的齐次解跟特解部分。已知输入和差分方程的稀疏，可用filter 对差分方程进行数值求解。最简单形式为：y=filter(b,a,x) 其中b=[b0,b1,…,bM];a=[a0,a1,…,aN]; 2．上面差分方程解的形式为齐次解和特解，另外还可以求零输入解和零状态解理论计算中要用到z 变换，请好好掌握z 变换的内容。用MATLAB 实现时，若已知初始条件，则应用y=filter(b,a, x, xic)来求完全响应。这里xic 是初始状态输入数组。MATLAB 还提供一种filtic 函数来得到xic 。
收稿日期:2009-09-01 第一作者简介:贾新民(1956-),男,四川邻水人,新疆昌吉学院计算机工程系,副教授,研究方向:计算机程序设计及其 语言教学和理论物理研究。 有限差分法求解拉普拉斯方程 贾新民1 严 文2(1.昌吉学院计算机工程系新疆昌吉831100;2.昌吉学院物理系 新疆 昌吉831100)摘 要:以极板上具有半圆截面沟槽的电容器内的电势分布为例,介绍了综合应用计算机软件利用有限 差分法求解复杂边界的拉普拉斯方程数值解的方法。并利用数值解的结果讨论了沟槽表面的电场分布和电荷分布。 关键词:拉普拉斯方程;有限差分法;五点差分格式 中图分类号:O411.2 文献标识码:A 文章编号:1671-6469(2009)05-0105-05 1 引言 无源空间的引力场、静电场、稳定的温度分布等问题都满足拉普拉斯(Laplace )方程 2u (x ,y ,z )=0(1)但由于方程(1)是偏微分方程,只有在问题具有高度对称的情况下,才能求出解析解,而这种情形是极少 的。有些情形看上去很简单,但却求不出解析解。对于这些情况,只能寻求数值解。 2 计算机数值解法方案文献[1][2][3]给出了拉普拉斯方程数值解的方法———有限差分法。 有限差分法的思想是用差分Δu (x +Δx ,y )Δx ,Δu (x ,y +Δy )Δy 代替导数5u 5x ,5u 5y ,用网格将求解区域覆盖,对于平面拉普拉斯方程,第i 行第j 列小格的电势由Laplace 方程的五点差分格式给出。 u ij =14(u ij -1+u ij +1+u i -1j +u i +1j )(2) 考虑图1所示的具有半圆形截面的槽的电容器内部的电势和电场分布。 为了能够对坑(槽)内部的电场进行比较细致的观察, 应该将半径R 取的大些,为了满足无限远的条件,应该使 求解区域尽量大些。 我们选Excel 为计算工具,因为Excel 具有不用编写 程序和直观的优点。Excel 的一个单元格代表求解区域的 一个网格,单元格的值表示该网格处的电势。根据Excel 的能力,取R =28,极板宽度为16×28=448,两极板距离 d 取8×28=224。由于Excel 的最大列数为256(A - IG ),最大行数为65536,为了能够用Excel 计算,需要将问题的行列置换,即将224+28=252行换成252列,而将448列换成448行。 5 01
