(完整版)大连理工大学高等数值分析抛物型方程有限差分法
抛物型方程的差分格式
a umn 1 umn 1 2h
a
umn 1
2umn h2
umn 1
抛物型方程的古典显格式
整理得方程(2.38)的显式格式(2.39)
U n1 m
(1
2ra)U
n m
r
(a
1 2
ha)U
n m1
r
(a
1 2
ha)U
n m1
截断误差为 O(k h2 ).
tn x
抛物型方程的古典显格式
三、算子
Dx
x
为 x 方向偏导数算子
Tx为 x 方向位移算子
Txumn umn 1, Tx1umn umn 1
μ x 为 x 方向平均算子
xu
n m
1 2
un
m
1 2
un
m
1 2
抛物型方程的古典显格式
x 方向差分算子
边界条件为 u(0,t) 1(t) 0, 0 t 0.20 u(1,t) 2(t) 0, 0 t 0.20
取步长⊿x = h = 0.2 , ⊿t = k = 0.02 。
抛物型方程的古典显格式
解 r = k / h2 = 0.02 / 0.22 = 0.5, 古典显式格式为
n m
umn
h
h 2!
2u x2
n
m
h2 3!
3u x3
n
m
(完整版)大连理工大学高等数值分析偏微分方程数值解(双曲方程书稿)
双曲型方程的有限差分法线性双曲型方程定解问题： （a ）一阶线性双曲型方程()0=∂∂+∂∂xux a t u （b ）一阶常系数线性双曲型方程组0=∂∂+∂∂xt uA u 其中A ，s 阶常数方程方阵，u 为未知向量函数。
（c ）二阶线性双曲型方程（波动方程）()022=⎪⎭⎫⎝⎛∂∂∂∂-∂∂x u x a x t u()x a 为非负函数（d ）二维，三维空间变量的波动方程0222222=⎪⎪⎭⎫⎝⎛∂∂+∂∂-∂∂y u x u t u 022222222=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-∂∂z u y u xu t u §1 波动方程的差分逼近 1.1 波动方程及其特征线性双曲型偏微方程的最简单模型是一维波动方程：（1.1） 22222xu a t u ∂∂=∂∂ 其中0>a 是常数。
（1.1）可表示为：022222=∂∂-∂∂xu a t u ，进一步有0=⎪⎭⎫ ⎝⎛∂∂+∂∂⋅⎪⎭⎫ ⎝⎛∂∂-∂∂u x a t x a t由于xat ∂∂±∂∂当a dt dx ±=时为()t x u ,的全导数（=dtdu dt dx x u t u ⋅∂∂+∂∂x ua t u ∂∂±∂∂=）,故由此定出两个方向（1.3） adx dt 1±=解常微分方程（1.3）得到两族直线 （1．4） 1C t a x =⋅+ 和 2C t a x =⋅- 称其为特征。
特征在研究波动方程的各种定解问题时，起着非常重要的作用。
比如，我们可通过特征给出（1.1）的通解。
（行波法、特征线法） 将（1.4）视为),(t x 与),(21C C 之间的变量替换。
由复合函数的微分法则212211C uC u x C C u x C C u x u ∂∂+∂∂=∂∂⋅∂∂+∂∂⋅∂∂=∂∂ x C C u C u C x C C u C u C x u ∂∂⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+∂∂⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂=∂∂2212121122 222122212212C u C C u C C u C u ∂∂+∂∂∂+∂∂∂+∂∂= 2222122122C uC C u C u ∂∂+∂∂∂+∂∂= 同理可得a t t a t C -=∂∂-=∂∂1，a tC=∂∂2 ⎪⎪⎭⎫⎝⎛∂∂-∂∂=∂∂⋅∂∂+∂∂⋅∂∂=∂∂212211C u C u a t C C u t C C u t utC C u C u a C u t C C u C u a C t u ∂∂⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛∂∂-∂∂⋅∂∂+∂∂⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛∂∂-∂∂⋅∂∂=∂∂2122112122 ⎥⎦⎤⎢⎣⎡∂∂∂-∂∂+⎥⎦⎤⎢⎣⎡∂∂-∂∂∂-=21222222221222C C u C u a C u C C u a ⎥⎦⎤⎢⎣⎡∂∂+∂∂∂-∂∂=22221221222C u C C u C u a 将22x u ∂∂和22tu∂∂代入（1.1）可得： ⎥⎦⎤⎢⎣⎡∂∂+∂∂∂-∂∂22221221222C u C C u C u a ⎥⎦⎤⎢⎣⎡∂∂+∂∂∂+∂∂=22221221222C u C C u C u a 即有0212=∂∂∂C C u求其对2C 的积分得：()11C f C u=∂∂ 其中()1C f 是1C 的任意可微函数。
微分方程数值解法课程设计---抛物型方程问题的差分格式[9页].doc
![微分方程数值解法课程设计---抛物型方程问题的差分格式[9页].doc](https://img.taocdn.com/s3/m/01f50aa57cd184254a353588.png)
目录一、问题的描述 (1)二、算法设计及流程图 (1)2.1 算法设计 (1)2.2 流程图 (2)三、算法的理论依据及其推导 (2)3.1 截断误差分析 (2)3.2 稳定性分析 (3)四、数值结果及分析 (3)五、总结 (5)六、附件（源代码） (6)抛物型方程问题的差分格式一、问题的描述有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。
此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解(即收敛性)，等等。
偏微分方程边值问题的差分法是物理上的定常问题，其定解问题为各种边值问题， 即要求解在某个区域内满足微分方程，在边界上满足给定的边界条件。
常系数扩散方程的差分解法可归结为选取合理的差分网格，建立差分格式求解。
常系数扩散问题的有限差分格式求常系数扩散问题为正常数其中a ,0,,22>∈∂∂=∂∂t R x xua t u （1.1） 的近似解，其初始条件为R x x g x u ∈=),()0,(二、算法设计及流程图2.1 算法设计运用加权隐式格式求解常系数扩散问题（1.1）02)1(22111112111=⎥⎥⎦⎤⎢⎢⎣⎡+--++-------+-+-h u u u h u u u a u u n j n j n j n j n j n j n jn j θθτ，（1.6） 10≤≤θ,h τ其中分为时间步长和空间步长。
步骤1 输入初始值，确定加权隐式格式的参数；步骤2 定义向量A ，把初边值条件离散，得到0j u ，j=0,1,…,J 的值存入向量A 步骤3 利用加权隐式差分格式由第n 层计算第n+1层，建立相应线性方程组，求解并且存入向量A;步骤4 计算到t=1，输出u2.2 流程图三、算法的理论依据及其推导3.1 截断误差分析常系数扩散问题（1.1）的加权隐式格式如下：02)1(22111112111=⎥⎥⎦⎤⎢⎢⎣⎡+--++-------+-+-h u u u h u u u a u u n j n j n j n j n j n j n jn j θθτ，（1.6） 其中10≤≤θ，,h τ其中分为时间步长和空间步长。
抛物型方程的差分方法
抛物型方程的差分方法抛物型方程是描述物理现象中的薄膜振动、热传导、扩散等过程的方程，具有非常重要的应用价值。
差分方法是一种常用的数值计算方法，用于求解微分方程，对于抛物型方程的数值求解也是非常有效的方法之一、本文将介绍抛物型方程的差分方法，并具体讨论用差分方法求解抛物型方程的一些具体问题。
首先，我们来介绍一下抛物型方程的一般形式。
抛物型方程一般可以表示为：∂u/∂t=α(∂^2u/∂x^2+∂^2u/∂y^2)其中，u(x,y,t)是待求函数，t是时间，x和y是空间变量，α是常数。
这个方程描述的是物理过程中的扩散现象，如热传导过程、溶质的扩散过程等。
差分方法的基本思想是将求解区域离散化为一个个网格点，然后在每个网格点处用近似的方式来计算待求函数的值。
差分方法的求解步骤主要包括以下几个方面：1.选择适当的网格和步长。
在求解抛物型方程时，需要确定空间变量x和y所在的网格点以及步长，同时也需要确定时间变量t所在的网格点和步长。
通常，我们会选择均匀网格，步长选择合适的值。
2.建立差分格式。
差分格式是差分方法的核心部分，它包括对方程进行近似处理和离散化。
对于抛物型方程，常用的差分格式有显式差分格式和隐式差分格式等。
其中，显式差分格式的计算速度快，但是有一定的稳定性限制，而隐式差分格式的稳定性较好，但是计算量较大。
因此，在具体问题中需要根据实际情况选择适当的差分格式。
3.编写计算程序。
在建立差分格式后，需要编写计算代码来求解离散方程。
具体编写的过程包括定义初始条件、建立迭代计算过程、以及计算结果的输出等。
4.计算结果的验证与分析。
求解方程后，需要对计算结果进行验证和分析，主要包括对数值解和解析解的比较、对误差的估计和控制等。
在具体求解抛物型方程时，还会遇到一些问题，例如边界条件的处理、稳定性和收敛性的分析等。
下面将对其中一些问题进行详细讨论。
1.边界条件的处理。
边界条件对差分格式的求解结果有着重要的影响，常见的边界条件包括固定端（Dirichlet）边界条件和自由端（Neumann）边界条件等。
大连理工大学高等数值分析主要内容总结
高等数值分析主要内容总结1. 矩阵部分(1) 矩阵变换a) Householder 变换（反射变换/镜像变换）定义设ωϵC n是一个单位向量，令H(ω)=I−2ωωH(1)则称H是一个Householder矩阵或Householder变换。
H具有以下性质：H H=H (Hermite矩阵),H H H=I (酉矩阵)H2=I (对合矩阵), H−1=H (自逆矩阵)det(H)=−1, diag(I,H)也是一个Householder矩阵定理设uϵC n是一个单位向量，则对于任意的xϵC n，存在Householder矩阵H=I−2ωωH，使得Hx=au，其中|a|=‖x‖2（a不唯一）。
当x=0时，ω可任取；当x=au≠0时，取ωH x=0；当x≠au时，应取ω=x−au。
‖x−au‖2b) Givens 变换（旋转变换）定义 设c,sϵC ，|c |2+|s |2=1，记n 阶矩阵T kl =[1⋱1c̅s̅1⋱1−sc1⋱1](k ) (l)(2)称为Gives 矩阵或初等旋转矩阵。
Givens 矩阵为酉矩阵，且det (T kl )=1。
定理 对于任意向量xϵC n ，存在Givens 变换T kl ，使得y =T kl x 的第k 个分量为非负实数，第l 个分量为0，其余分量不变。
当|x k |2+|x l |2=0时，取c =1，s =0，则T kl =I 。
当|x k |2+|x l |2≠0时，取c =k √|x k |2+|x l |2，s =l√|x k |2+|x l |2。
(2) 矩阵分解-QR 分解（正交三角分解/酉三角分解）定义 设AϵC n×n ，如果存在n 阶酉矩阵（酉矩阵是正交矩阵往复数域上的推广）Q 和n 阶上三角矩阵R ，使得A =QR ，则称A =QR 为A 的QR 分解。
当AϵR n×n 时，则称为A 的正三角分解。
定理 任意一个满秩实（复）矩阵A ，都可唯一地分解为A =QR 。
(十二章)抛物型方程有限差分法
（3.3a） （3.3b）
； 其中
， LOD格式的计算步骤可以总结如下：
1） 令，。 2） 求解三对角线性方程组（3.3a）得到差分解。 3） 若，则增加1，转步骤4）。否则转4）。 4） 令。 5） 求解三对角线性方程组（3.3b）得到差分解。 6） 若，则增加1，转步骤5）。否则转7）。 7） 若，则增加1，转步骤2）。否则结束。
,
及边值条件
，
假定和在相应的区域光滑，并且于，两点满足相容条件,
则上述问题有唯一的充分光滑的解。
现在考虑边值问题（1.1），（1.3）的差分逼近
取 为空间步长，为时间步长，其中，是自然数，
,； ,
将矩形域分割成矩形网格。其中 表示网格节点；
表示网格内点（位于开矩形中的网格节点）的集合；
表示位于闭矩形中的网格节点的集合；
时我们简单地称差分格式稳定。
冯诺依曼稳定性分析估量了误差的放大或扩大。对一
种稳定的方法，必须选取步长使误差的放大因子不大于1.
前面讨论的向前差分格式（1.4）当网比时稳定，当时不稳定。这就
意味着给定空间步长以后，时间步长必须足够小，才能保证稳定。而向
后差分格式（1.6）和Grank-Nicholson格式（1.8）则对任何网比都是稳
+=+ ，==0。求出，在由取，可利用，解出，。如此下去，即可逐 层算出所有，。
如此每层必须解一个三对角线性方程组的格式称为隐格 式。并视为的近似值。
直观地说，采用显式格式进行求解既方便又省工作量。但 是，后面我们将看到，有些情况用隐式格式更为便利。
1.2.3 Grank-Nicholson法 将向前差分格式和向后差分格式做算术平均，得到的差分 格式称之为六点对称格式，也称为Grank-Nicholson格式：
有限差分法求解抛物型方程说明
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例.§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 层上的值求出。
抛物方程的有限差分法
图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.。
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
第三章_抛物型方程的有限差分法(1)
1 1 ek ek j j
a
fj
k
k e k 用 ekj u j u j 表示相应节点处差分解的误差，则 j 满足
2
a
h2
以上是误差传播方程。设误差只在初始层的原点(j=0)发生，
0 1 即 e j j 0 ( 0, 00 1, j 0 0,当j 0 ） ， e j 0 ，而在以后的
u 2u a 2 f ( x), t x 0t T
(1)
其中 a 是正常数，f(x)是给定的连续函数。 x 称作空间变量，
t 称作时间变量。
（1)的定解问题有以下几类： 初值问题(也称 Cauchy 问题)：求函数 u(x,t)，满足方程(1) 和初始条件：
u( x,0) ( x), x
k1
a u ( x j 1,tk ) 2u ( x j ,tk ) u ( x j 1,tk ) [ 2 h2 u ( x j 1,tk 1 ) 2u ( x j ,tk 1 ) u ( x j 1,tk 1 ) k1 2 ] [ u au ] t xx j 2 h
1
k 1 j
1
占优且是三对角的，方程总是可解的，可采用追赶法 求得方程的解。 截断误差：
1 Rk j (u ) Lh u ( x j , t k 1 ) Lu ( x j , t k 1 )
u ( x j 1,tk 1 ) 2u ( x j ,tk 1 ) u ( x j 1,tk 1 ) u ( x j , tk 1 ) u ( x j , tk ) k 1 a ut au xx j 2 h ut au xx j
h2
抛物型方程有限差分法
抛物型方程有限差分法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 Λ。
抛物形扩散方程的有限差分法与数值实例
偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生：向聘抛物形扩散方程的有限差分法及数值实例1.1抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。
考虑一维热传导方程：22(),0u ua f x t T t x∂∂=+<≤∂∂ (1.1.1) 其中a 是常数，()f x 是给定的连续函数。
按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件：()()x x u ϕ=0,, ∞<<∞-x (1.1.2)第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件：()()x x u ϕ=0,, 0x l << (1.1.3) 及边值条件()()0,,0==t l u t u ， T t ≤≤0 (1.1.4)假定()x f 和()x ϕ在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
1.2抛物线扩散方程的求解下面考虑如下热传导方程22()(0.)(,)0(,0)()u ua f x t x u t u L t u x x ϕ⎧∂∂=+⎪∂∂⎪⎪==⎨⎪=⎪⎪⎩(1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。
取N l h =为空间步长，MT=τ为时间步长，其中N ，M 是自然数，用两族平行直线jh x x j ==,()N j ,,1,0Λ=和k t t k τ==, ()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。
其中 (),j k x t 表示网格节点；h G 表示网格点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。
抛物型方程的有限差分方法
抛物型方程的有限差分方法一，求解问题考虑一维非齐次热传导方程的定解问题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 +≤≤。
(完整版)大连理工大学高等数值分析有限元简述-2017
椭圆与抛物微分方程的有限元法空间m H 作为例子，我们将考虑区间()0,1I =上的微分方程。
用2()L I 表示在I 上勒贝格平方可积函数的集合，()m H I 表示本身以及直到m 阶的导数都属于2()L I 的函数的集合。
我们下面用到的主要是1()H I 。
这里所说的导数准确地说是应该是广义导数，对此我们不予详细说明，只需知道比如说，连续的分片线性函数（折线函数）就属于1()H I ，其广义导数是分片常数函数。
另外，我们还用到空间}0)0(),({)(11=∈=v I H v I H E。
（空间=函数集合。
）微分方程 考虑两点边值问题(), (0,1)pu qu f x ''-+=∈（1） (0)0u = （2）0)1(='u（3）其中, , p q f 都是区间)1,0(上的光滑函数，0≥q , 并且0p p ≥，0p 是一个正常数。
用)(1I H E中任一函数v 乘（1）式两端，并在]1,0[上积分，得10[()]0pu v quv fv dx ''-+-=⎰ （4） 利用分部积分，并注意0)1(='u 和0)0(=v ，得⎰⎰''+'-=''-10110|)(dx v u p v u p vdx u p ⎰''=10dx v u p以此代入到（4）得到⎰=-+''100)(dv fv quv v u p （5） 为了方便，定义()10,w v w vdx =⋅⎰ （7）),(),(),(v qw v w p v w a +''=（8）则相应于微分方程（1）-（3）的变分方程为：求)(1I H u E∈满足),(),(v f v u a = )(1I H v E∈∀ （9）注意在（9）中不出现二阶导数。
可以证明，满足微分方程（1）-（3）的光滑解一定满足变分方程（9）。
第四章 抛物型方程的有限差分方法
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
抛物型方程的差分方法
tn nk xm mh
n
0,1,2,
,
N;
N
T k
m 0,1,2,
在 t 0上的结点称为边界结点，属于 内的结点 称为内部结点。
对于初边值问题，设 (x,t) | 0 x 1,0 t T
，则网格是
tn nk xm mh
n
0,1,2,
,
N;
N
T k
m 0,1,2, , M ; Mh 1
u(0,t) 1(t),u(1,t) 2 (t) 0 t T (2.4)
2.1 差分格式建立的基础
为了构造微分方程(2.1)的有限差分逼近，首 先将求解区域用二组平行于t 轴和 x轴的直线构 成的网格覆盖，网格边长在方向 t为 t k，在 x 方向为x h (如图2.1所示)。h,k分别称为空间方向 和时间方向的步长，网格线的交点称为网格的结 点。对初值问题来说，网格是
4x
umn
2 x
3 x
11 12
3x
umn
2 x
1 12
4 x
1 90
6 x
umn
(2.19.1)
(2.19.2)
(2.19.3)
返回 返回 35 42
对于三阶、四阶偏导数的差分表达式为
h3
(
3u x3
)nm
3x 3x
3 2
4x
3 2
4 x
7 4 7 4
5x
5 x
umn umn
(2.20.1) (2.20.2)
)U
n m
代入
2 x
的表达式，则得差分方程
(2.28)
古典显式差分格式
U
n1 m
rU
抛物型方程的差分方法
2) 逆风差分格式
1 un un j j
a
n un j u j 1
§5 抛物型方程的差分方法
本章內容: 1． 常系数抛物型方程的初值问题 2．初边值问题 3．对流扩散方程 4．Richardson外推法 ^^数值实验4（网格比的重要性） ^^数值实验5 (Richardson 外推的精度) 5．变系数方程 6．二维抛物型方程问题的计算
1． 常系数抛物型方程初值问题
1 2
返回本节
c. 三层显式关系 Richardson格式 n n 1 1 un un un j j j 1 2u j u j 1 a 0 h2 2 Du Fort－Frankel格式
1 1 un un j j
2
a
n 1 1 n un un j 1 (u j j ) u j 1
目标点：Jh 1, x j jh, N T , ( x j , t n )
1 n n un un un j j j 1 2u j u j 1 内部点的离散： a 0 h2
边值点的离散： u0 ( n ), u J (n )
n n
初值点的离散： u j f (x j ) f j
1 4(1 )a sin 2 1 4a sin 2 kh 2
kh 2
3
G ( , k ) 1 1 G ( , k ) 1 4a (1 2 )sin 2 2a (1 2 ) 1 kh 2 2 0 1 2
大连理工大学 高等数值分析 常微分方程数值解法-2017
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。
通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。
本章讨论常微分方程最常用的近似数值解法-差分方法。
先来讨论最简单的Euler 法。
为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。
在微积分课程中我们熟知，微商（即导数）是差商的极限。
反过来，差商就是微商的近似。
在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
抛物型方程有限差分法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）的差分逼近 取Nlh =为空间步长，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 j u ux t t t∂∂⎛⎫≡ ⎪∂∂⎝⎭）： ()()()ττO t u t x u t x u kj k j k j +⎪⎭⎫⎝⎛∂∂=-+,,1 ()()()2112,,ττO t u t x u t x u k jk j k j +⎪⎭⎫⎝⎛∂∂=--+()()()h O x u h t x u t x u kj k j k j +⎪⎭⎫ ⎝⎛∂∂=-+,,1()()()h O x u ht x u t x u kj k j k j +⎪⎭⎫⎝⎛∂∂=--,,1 ()()()2112,,h O x u ht x u t x u kjk j k j +⎪⎭⎫⎝⎛∂∂=--+()()()()222211,,2,h O x u ht x u t x u t x u kjk j k j k j +⎪⎪⎭⎫ ⎝⎛∂∂=+--+可得到以下几种最简差分格式（一） 向前差分格式 ()14.1 =-+τk jk j u u 1j kj k j k j f h u u u a++--+2112()()j j x f f =()24.1()j j j x u ϕϕ==0，k u 0=kN u =0其中1,,1,0-=N j Λ，1,,1,0-=M k Λ。
取2ha r τ=为网比，则进一步有()14.1'1+k j u =k j ru 1++()r 21-k j u +kj ru 1-+j f τ此差分格式是按层计算：首先，令0=k ，得到1j u =01+j ru +()r 21-0j u +01-j ru +j f τ于是，利用初值()j j j x u ϕϕ==0和边值k u 0=kN u =0，可算出第一层的1j u ，1,,1,0-=N j Λ。
再由()14.1'取1=k ，可利用1j u 和k u 0=kN u =0算出2ju ，1,,1,0-=N j Λ。
如此下去，即可逐层算出所有k ju （1,,1,0-=N j Λ，1,,1,0-=M k Λ）。
由于第()1+k 层值可以通过第()k 层值直接得到，如此的格式称为显格式。
并视k j u 为()k j t x u ,的近似值。
若记()TkN k k k u u u 121,,,-=Λu ，()()()()T N xx x 121,,,-=ϕϕϕϕΛ，()()()()T N x f x f x f 121,,,-=τττΛf则显格式()14.1'可写成向量形式⎩⎨⎧=-=+=+ϕ11,,1,0,u f Au u M k k k Λ 其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=r r r r r r r rr r21002100210021ΛO MO O O M O ΛA 若记22xua t u Lu ∂∂-∂∂=()--=+τk jk j kjh u u u L 112112h u u u akj k j k j -++-那末截断误差（1.5）()=u R kj()()[]k jk j h Lu t x u L -,1=()ττO t x t u r k j +⎪⎪⎭⎫ ⎝⎛∂∂⎪⎭⎫ ⎝⎛--)~,~(2112122=()2h O +τ 其中(,)j k x t %%是矩形11+-<<j j x x x ，1+<<j k t t t 中某一点。
事实上，()=u R kj kj x u ⎪⎪⎭⎫⎝⎛∂∂222τ+()2τO kjx u h a ⎪⎪⎭⎫ ⎝⎛∂∂⋅-442ˆ12 =k j x u ⎪⎪⎭⎫⎝⎛∂∂222τ+()2τO ()22222ˆ112τO t u a h a kj +⎪⎪⎭⎫ ⎝⎛∂∂⋅⋅⋅- =⎥⎦⎤⎢⎣⎡-⋅-211212ττa h ()222~τO t u kj+⎪⎪⎭⎫⎝⎛∂∂ =⎥⎦⎤⎢⎣⎡--21121r τ()222~τO t u kj+⎪⎪⎭⎫ ⎝⎛∂∂=()2h O +τ。
这里⎪⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂222244x u xa x u a ⎪⎭⎫ ⎝⎛∂∂⋅∂∂=t u a x a 122⎪⎭⎫⎝⎛∂∂∂∂=t u x 22tx u ∂∂∂=2322t u ∂∂⎪⎭⎫ ⎝⎛∂∂∂∂=t u t 2⎪⎪⎭⎫ ⎝⎛∂∂∂∂=22x u a t t x u∂∂∂=23故22t u ∂∂44244x u a x u a a ∂∂=⎪⎪⎭⎫ ⎝⎛∂∂⋅=，从而=∂∂44x u 221t u a ∂∂⋅（二） 向后差分格式 ()16.1 =-+τk jk j u u 1j k j k j k j f hu u u a++-+-+++2111112()()j j x f f =()26.1()j j j x u ϕϕ==0，k u 0=kN u =0其中 1,,1,0-=N j Λ，1,,1,0-=M k Λ。
取2h a r τ=为网比，则进一步有()16.1'r -k j u 1++()r 21+1+k j u r -11+-k j u =kj u +j f τ按层计算：首先，取0=k ，则利用初值()j j j x u ϕϕ==0和边值k u 0=kN u =0，来确定出第一层的1j u ，1,,1,0-=N j Λ，即求解方程组：r -11+j u +()r 21+1j u r -11-j u =0j u +j f τ1,,1,0-=N j Λ，k u 0=kN u =0。
求出1j u ，在由()14.1'取1=k ，可利用1j u ，解出2j u ，1,,1,0-=N j Λ。
如此下去，即可逐层算出所有k j u ，1,,1,0-=M k Λ。
如此每层必须解一个三对角线性方程组的格式称为隐格式。
并视k j u 为()k j t x u ,的近似值。
直观地说，采用显式格式进行求解既方便又省工作量。
但是，后面我们将看到，有些情况用隐式格式更为便利。
1.2.3 Grank-Nicholson 法将向前差分格式和向后差分格式做算术平均，得到的差分格式称之为六点对称格式，也称为Crank-Nicholson 格式： ()18.1=-+τk jk j u u 1j k j k j k j k j k j k j f h u u u h u u u a +⎥⎥⎦⎤⎢⎢⎣⎡+-++-+-+++-+211111211222()()j jx f f=()28.1()j j j x u ϕϕ==0，k u 0=kN u =0进一步， ()18.1'2r -11++k j u +()r +11+k j u 2r -11+-k j u =2r k j u 1++()r -1kju 2r +k j u 1-+j f τ 按层计算：首先，取0=k ，则利用初值()j j j x u ϕϕ==0和边值k u 0=kN u =0，来确定出第一层的1j u ，1,,1,0-=N j Λ，即求解方程组：2r -11+j u +()r +11j u 2r -11-j u =2r 01+j u +()r -10ju 2r +01-j u +j f τ 1,,1,0-=N j Λ，k u 0=k N u =0。
求出1j u ，在由()18.1',取1=k ，可利用1j u ，解出2j u ，1,,1,0-=N j Λ。
如此下去，即可逐层算出所有k j u ，1,,1,0-=M k Λ。
若记22xua t u Lu ∂∂-∂∂=()--=+τk jk j k j h u u u L 13j k j k j k j k j k j k j f h u u u h u u u a +⎥⎥⎦⎤⎢⎢⎣⎡+-++-+-+++-+211111211222 在1(,)(,()/2)j k k x t x t t +=+处作Taylor 展开，可以算出截断误差为 （1.7） ()=u R k j()()[]k j k j h Lu t x u L -,3=()22h O +τ。
+（四）Richardson 格式（1.10） =--+τ211k jk j u u 2112hu u u akj k j k j -++-+j f进一步()110.1'1+k j u =r 2（k j u 1+k j u 2-+k j u 1-）+1+k j u +2j f τ这是三层显式差分格式。
显然截断误差的阶为()22h O +τ。
为使计算能够逐层进行，除初值0j u 外，还要用到1j u 。
它可以用其他双层格式提供。
Richardson 格式的矩阵形式为：⎪⎩⎪⎨⎧=-=++=-+另算10111,,1,u u f u u C u ϕτM k k k k Λ 其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------⋅-=21001210012100122ΛO M O O O M O Λr C2 稳定性与收敛性抛物方程的两层差分格式可以统一写成向量形式: (2.1) 1k k AU BU F τ+=+其中1111(,), (,,)k k k TN N U u u F f f --==LL ，A 和B 是1N -阶矩阵。