抛物形扩散方程的有限差分法及数值实例(完整资料).doc
扩散方程的差分解法
（24）
则
， （25）
则误差方程为
（26）
误差放大因子为
（27）
要满足稳定性条件，则要求对所有的k值均有 。从（28）式中可以看出，当 （即 ）时， 恒成立。因此，全隐格式是无条件稳定的。
4.4收敛性
如果差分方程的解为 ，微分方程的解为 ，若当 ， 时，差分方程的解与微分方程的解之差
扩散方程的差分解法
在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
write(2,*) 'x=',x,'m'
do n=1,nt,200
write(2,*) (n-1)*dt,u(j,n)
enddo
if
enddo
!-----------!
end
5.3.2全隐格式
!----------------------------------------全隐格式求解扩散方程-----------------------------------------------!
由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。
在本文计算中，取 ， 。
初始条件（ 时）
（29）
边界条件为
（30）
其初始时刻（ ）时的u分布如图1所示，x=0m处u随时间变化情况如图2所示，x=10m处u随时间变化情况如图3所示。
图1初始时刻u分布图
图2 x=0处u随时间变化图
抛物型方程的有限差分法
令
L(h3
)
u
k j
uk1 j
u
k j
a
[
uk1 j1
2
2u
k j
1
h2
uk 1 j 1
uk j1
2u
k j
h2
uk j 1
]
将截断误差
R
k j
(u)
L(h3)u( x j , tk
)
[Lu]kj
于(
x
j
,
t
k
1
2
)(tk 1 2
(k
1 )
2
)展 开 ， 则 得
R
k j
(u)
0(
2
h2
).
(1.9)
j
k1
k1
k1
k
k
k
j1
j
j1
j1
j
j1
2
h2
h2
j
u0 ( x ), uk uk 0 ,
j
j
j
0
N
( 1.8 ) 2
( 1.8 ) 1
将(1.8)1改 写 为
r 2
uk1 j1
(1
r
)u
k j
1
r 2
uk 1 j 1
r 2
uk j1
(1
r
)u
k j
1
r 2
uk j 1
fj
(1.8)1
u(
t x,0)
a x2
(x),
f (x), 0 x
0 l
t
T (1.3)1
u(0,
t)
u(l, t)
0,
抛物型方程的差分方法
抛物型方程的差分方法抛物型方程是描述物理现象中的薄膜振动、热传导、扩散等过程的方程，具有非常重要的应用价值。
差分方法是一种常用的数值计算方法，用于求解微分方程，对于抛物型方程的数值求解也是非常有效的方法之一、本文将介绍抛物型方程的差分方法，并具体讨论用差分方法求解抛物型方程的一些具体问题。
首先，我们来介绍一下抛物型方程的一般形式。
抛物型方程一般可以表示为：∂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）边界条件等。
扩散方程的差分解法
扩散方程的差分解法在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。
热传导方程中的自变量中包括时间t ，它是描述一种随时间变化的物理过程，即所谓不定常现象。
这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。
1.扩散方程一维扩散方程为：22u u t xα∂∂=∂∂ （1）式中，u 为因知量，α为扩散系数，x 为坐标，t 为时间。
其定解条件如下： 初始条件： (,0)() 0x u x f x L =≤≤（2）边界条件： 12(0,)() , (,)()u t f t u L t f t ==（3） 一般假定函数()f x ，1()f t ，2()f t 满足连接条件，即1(0)(0) f f =，2()(0) f L f =。
2.有限差分法有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。
其控制方程中的导数用离散点上函数值的差商代替。
差分格式可以分为显格式和隐格式。
所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。
由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。
隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。
因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。
为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。
因此，隐格式总是伴随着求解巨大的代数方程组。
10_抛物型方程的有限差分方法
10_抛物型方程的有限差分方法抛物型方程是一类常见的偏微分方程，广泛应用于自然科学和工程学的领域中。
有限差分方法是一种常用的数值求解抛物型方程的方法之一、本文将介绍抛物型方程的有限差分方法（II）。
有限差分方法主要基于离散化的思想，将偏微分方程转化为差分方程，进而求解差分方程的数值解。
对于抛物型方程，其一般形式可以表示为：∂u/∂t=Δu+f(x,t)其中，u(x, t)是未知函数，表示空间位置x和时间t上的解，Δu表示Laplace算子作用于u的结果，f(x, t)是已知函数。
有限差分方法的基本思想是将空间和时间域进行离散化，将连续的空间和时间划分为有限个网格点，然后使用差分近似代替偏导数，得到差分方程。
假设空间域被划分为Nx个网格点，时间域被划分为Nt个网格点，对于每个网格点(i,j)，可以表示为(x_i,t_j)，其中i=0,1,...,Nx，j=0,1,...,Nt。
在有限差分方法中，我们使用中心差分近似来代替偏导数。
对于时间导数，可以使用向前差分或向后差分，这里我们使用向前差分，即：∂u/∂t≈(u_i,j+1-u_i,j)/Δt对于空间导数，可以使用中心差分，即：∂^2u/∂x^2≈(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2将上述差分近似代入抛物型方程中，可以得到差分方程的离散形式：(u_i,j+1-u_i,j)/Δt=(u_i-1,j-2u_i,j+u_i+1,j)/Δx^2+f_i,j其中，f_i,j=f(x_i,t_j)。
重排上式，可以得到递推关系式：u_i,j+1=αu_i-1,j+(1-2α)u_i,j+αu_i+1,j+Δt*f_i,j其中，α=Δt/Δx^2通过设置初始条件和边界条件，可以利用以上递推关系式得到抛物型方程的数值解。
总结来说，抛物型方程的有限差分方法（II）是一种常用的数值求解抛物型方程的方法。
它基于离散化的思想，将偏微分方程转化为差分方程，然后利用中心差分近似代替偏导数，得到差分方程的离散形式。
抛物型对流扩散方程
抛物型对流扩散方程
抛物型对流扩散方程是水力学中一个重要的基本方程，它描述了液体中湍流运
动的数学表达形式。
抛物型对流扩散方程公式可由下式得到：
∂u∂t+u⋅∇u=−g⋅∇h+(∇⋅Δ)u-k∇2η,其中u是几何位移，t是时间，g是重力
加速度，h是重力场，Δ是拉普拉斯算子，k是拉格朗日运动等弦水动力系数，η
是密度。
抛物型对流扩散方程的应用很广泛，它可以用来分析流体的动态特性，并有助
于求解海洋涡场、各种湍流模式、源汇问题等。
举例来说，该方程可用来研究气候变化中河流流动物理过程，也可用来研究表面温带对于对流层等层结构、平流变化等关键过程中的影响。
此外，它还能够提供关于机械装置的流动特性的精确模拟。
抛物型对流扩散方程的求解不是一件容易的事情，它要求求解方法具有较高的
计算效率和求解准确度，尤其是人工网格的定义。
现阶段，多流变技术和网格技术均在快速发展，为使抛物型对流扩散方程能够尽可能反映实际环境中湍流流动特性，给求解方法提供更多可能。
总之，抛物型对流扩散方程是一个非常重要的基础性方程，它可以帮助我们深
入探究水力过程的机制，为水力学的研究和设计提供更为丰富的软件工具，从而满足现代水力学研究题目的需要。
有限差分法求解抛物型方程说明
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例.§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.。
有限差分求解扩散的数值解
Jeff
C2
无孔结构
J0 =-D ∂C ∂z
多孔介质
ε ∂C Jeff =- D τ ∂z
其中ε为孔隙率
曲折因子
τ= J0 ε Jeff
3.扩散模拟表征多孔介质物质输运性能（相转换流延制备） SRCT构建孔结构数组
灰度图
二值图
黑色：孔 白色：实体 孔隙率：0.35
三维结构重构图A
3.扩散模拟表征多孔介质物质输运性能
1.有限差分的基本原理（ 一维扩散例子） 计算结果
稳态下
∂C ∂2 C =D 2 = 0 ∂t ∂x
2.边界条件处理（ 二维柱坐标扩散例子）
扩散进口
扩散出口
绝缘 环
圆柱上底面和侧面为高浓度面（红色） 圆柱下底面同心小圆为低浓度面（白色） 圆柱下底面圆环为绝缘面 （蓝色）
2.边界条件处理（ 二维柱坐标扩散例子） 0 Z
r
扩散微分方程： ∂C ∂2 C 1 ∂C ∂2 C =D( 2 + r + 2) ∂t ∂r ∂r ∂Z
圆心对称
柱坐标下取体积微元由质量衡算推导扩散偏微分方程
r+ △r
△Z
r
体积微元质量变化量
∆m=-2π r+∆r ∆Z∆t*J -π r+∆r 2 -r2 ∆t*J
移项化简得
r+∆r+2πr∆Z∆t*J r Z
2 2 + π r+ ∆r -r ∆t*J Z+∆Z
∆m 1 -J = * 2πr∆Z∆r∆z r
r+∆r *
即
r+∆r -J r *r J Z+∆Z -J Z ∆r ∆Z ∂C ，将 J=-D 代入可得 ∂x
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
第四章 抛物型方程的有限差分方法
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
第三章_抛物型方程的有限差分法(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 Λ。
抛物型方程差分法
2m
从而要求 4rsin2i2, 1im1
2m
a 1
易见，只要 r h 2 2 就可以保证数值格式稳定。 称为稳定性条件
对于非齐次方程、非零边界条件的情形，其稳定性 分析仿上，只是差分格式现在变成
u r k 1 A u r k b r k r 其中向量 b k 依赖于方程的右端项和边界条件。
u ( 0 ,tk )( tk ) ,u ( 1 ,tk )( tk ) , 0kn.
3.处理方程 u
t
2u ax2
f(xi, tk)中的偏导数
(xi,tk)
(xi,tk)
关于时间的一阶偏导数用向前差商近似，
u u(xi,tk1)u(xi,tk)
t (xi,tk)
r 12r
0
r O
O r
0
12r r
1r2ruuuum m kkM 12kk12Auuuum m kkM 12kk12
也可以简写成 u rk1A u rk ，从而有
u r k 1 A ( A u r k 1 ) L A k 1 u r 0
( xi , tk )
( xi , tk ) — 网格节点
用
u
k i
表示温度分布函
数 u( x , t ) 在点 ( xi , tk )
处的网格函数 , 相当于
x x i 1 x i x i 1
u( x , t ) 在该点的近似 .
2. 原方程弱化为节点处的离散方程
连续方程 离散方程
u 2u tax2f(x,t), 0x1 , 0tT
将数值解 u
k i
代替精确解 u( xi , tk )
3-抛物型方程的差分方法
，则退化为古典隐式格式；
(3)取 1/ 2 ，则退化为Crank-Nicholson六点格式
为了提高对时间的截断误差，可用中心差分
u
n 1 j
u
n 1 j
2
a
u
n j 1
2u u h
n j 2
n j 1
0
Richardson格式，它是二阶精度的三层显式格式。 通过将其化为等价的二层差分格式，可给出其增 长矩阵为
n u1n u1n 1 au0 n n 1 u2 0 u2 n n 1 u u 0 3 3 n n 1 a u N 2 u N 2 u n u n 1 au n 1 2a N 1 N 1 N
u 2u 0 x 1, t 0 t a x 2 , u ( x, 0) ( x), 0 x 1 u / x u t0 x0 g1 (t ), t0 u / x u x 1 g 2 (t ),
古典显式格式
u
截断误差是 增长因子是
n 1 j
u
n j
a
u
n j 1
2u u h
n j 2
n j 1
0
T O( h2 )
kh G( , k ) 1 4a sin 2 其中网格比 / h2
2
稳定性条件是
1 a 2
古典隐式格式
n 1 un u j j
0
a 0 0
a 0
如用Crank-Nicholson六点格式 n n n n 1 n 1 n 1 n 1 1 1 1 a u (1 a ) u a u u a ( u 2 u u j 1 j j 1 j j 1 j j 1 ) 2 2 2 可得如下代数方程组
抛物形扩散方程的有限差分法与数值实例
偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生：向聘抛物形扩散方程的有限差分法及数值实例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 网格边界点的集合。
教案4(抛物型方程的有限差分方法)
第四章 抛物型方程的有限差分方法1 常系数扩散方程22u ua t x∂∂=∂∂，x R ∈，0t > （1.1） 其中a 为正常数。
如果给点初始条件(,0)()u x g x =，x R ∈ （1.2）（1.1）式和（1.2）式构成了一个初值问题。
1.1 向前差分格式，先后差分格式以上初值问题的向前差分格式111220n n n n nj jj j j u u u u u ahτ++---+-= （1.3）0()j j u g x = （1.4）其截断误差为2()O h τ+，易其增长因子为2(,)14sin 2khG k a τλ=- 其中2h τλ=，如果有12a λ≤，则有(,)1G k τ≤，即von Neumann 条件满足。
因此（1.3）稳定的条件是12a λ≤。
（1.1） 的先后差分格式 111220n n n n nj jj j j u u u u u ah τ-+---+-= （1.5）是无条件稳定的格式。
1.2 加权隐式格式把（1.3）改写为111111220n n n n n j jj j j u u u u u ah τ----+---+-= (1.3)' 用(1.5)(1)(1.3)θθ'⋅+-⋅得，111111112222[(1)]0n n n n nn n n j jj j j j j j u u u u u u u u a hhθθτ----+-+---+-+-+-= （1.6）其中01θ≤≤，称差分格式（1.6）为加权隐式格式。
其节点分布如图（4.1）把（1.6）改写为便于计算的形式11(12)n n nj j j a u a u a u λθλθλθ+--++-=11111(1)[12(1)](1)n n n j j j a u a u a u λθλθλθ---+--+--+-（1.7）其中2h τλ=。
(1.6)式中，取12θ=，用1n n +→，得 111111112[(2)(2)02n n j jn n n n n nj j j j j j u u a u u u u u u hτ+++++-+----++-+= （1.8） 此格式称为Crank-Nicolson 格式。
抛物型方程的有限差分方法
抛物型方程的有限差分方法一，求解问题考虑一维非齐次热传导方程的定解问题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 +≤≤。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
此文档下载后即可编辑偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：向聘抛物形扩散方程的有限差分法及数值实例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 (常数)是扩散系数。
取Nlh =为空间步长，MT =τ为时间步长，其中N ，M 是自然数，用两族平行直线jhx x j ==, ()N j ,,1,0Λ=和k t t k τ==,()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。
其中(),jkx t 表示网格节点；hG表示网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。
k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。
现在对方程进行差分近似：(一) 向前差分格式=-+τk jk j u u 11122(())k k kj j j j j j u u u af f f x h +--++=(1.2.2)()j j j x u ϕϕ==0，ku 0=kNu =0(1.2.3) 计算后得：111(12)k k k k j j j j ju ru r u ru f τ++-=+-++(1.2.4) 其中，2a r h τ=，1,,1,0-=N j Λ，1,,1,0-=M k Λ。
显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。
方程组如下：10001210110002321210003432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨⎪⎪⎪=+-++⎩M(1.2.5)若记()TkN k k k u u u 121,,,-=Λu ，()()()()T N x x x 121,,,-=ϕϕϕϕΛ，()()()()T N x f x f x f 121,,,-=τττΛf则显格式(1.2.4)可写成向量形式10,0,1,,1k k k M φ+⎧=+=-⎨=⎩L u Au f u(1.2.6) 其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=r r r r r r r rr r21002100210021ΛO MO O O M O ΛA 而对于向前差分格式，当网比12r ≤时稳定，当12r >时不稳定。
这就意味着给定空间步长h 以后，时间步长τ必须足够小，才能保证稳定。
1.3抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取1a =已知方程的精确解为2sin t u e x ππ-=：22(,0)sin()(0,)(1,)001,00.5u ut x u x x u t u t x t π⎧∂∂=⎪∂∂⎪⎪=⎨⎪==⎪≤≤≤≤⎪⎩(1.3.1)设空间步长1/h M =，时间步长为0.5/N τ=，网格比2/r g h =。
向前格式：11122,1,...,1,1,...,k k k k kj jj j j u u u u u j M k Nh τ++---+==-=边值条件：110010111220,,k kk k kk k u u u u u j u u h τ++----+===， 11111122,,k k k k k k k M MM M M M M u u u u u j M u u hτ+++-+---+===. 初值条件：(,0)sin ,0,1,...,j j u x x j Mπ==对时间和空间进行分割，令M=40，N=1600，通过Matlab 计算得到该方程的解析解，数值解以及相对误差如下：图(1)解析解的图像图(2)数值解的图像图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行比较，结果如表(1)数值解精确解相对误差x t0.1 0.1 0.1151 0.1152 4⨯1.170010-0.2 0.2 0.0815 0.0816 4⨯1.658010-0.3 0.3 0.0418 0.0419 4⨯1.275210-0.4 0.4 0.0183 0.0184 57.445610-⨯0.5 0.45 0.0117 0.0118 5⨯5.375510-0.6 0.35 0.0300 0.0301 4⨯1.067410-0.7 0.25 0.0684 0.0686 4⨯1.741010-0.8 0.15 0.5084 0.5085 5⨯7.589710-0.9 0.05 0.3654 0.3654 5⨯1.740710-表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在410-以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。
下面，我们比较在不同的分割下对有限差分算法精度的影响。
在扩散系数1a=不变的情况下，讲时间和空间进行更加细密的分割，取50,10000==，其中，M表示空间上的分割，N表M N示时间上的分割。
观察数值解与精确解在节点,x t处的绝对误差j k值，如下图所示：图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在410-10-和5以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。
最后，我们比较网比1/2r=时扩散方程的收敛情况。
r=以及1当网比1r=时，此时我们取M=10，N=50，这时，方程的数值解与解析解还有相对误差图如下：图(5)M=10,N=50的解析解的图像图(6)M=10,N=50的数值解的图像图(7)M=10,N=50的绝对误差的图像此时，我们观察绝对误差发现，扩散方程(1.3.1)时不收敛不稳定的。
而前面我们已经知道，到网格比为12r =时，方程是收敛稳定的。
所以，我们可以验证，当网比12r ≤时稳定，当12r >时不稳定。
[参考文献][1]李荣华，刘播.微分方程数值解法[M].北京.高等教育出版社.2009.1.[2]王曰朋.偏微分方程数值解[OL]. /link?url=-UAAXRI8_V547zVS6UwenT8rG KsbwNf9CuDmh2qmsy5K8eQ32PzhYazZ_sWiHfz_Pj3LA7ufHH4O 3t8NIlUCnXNUiyYhssqmNBeArsrLwQG[3]未知.偏微分方程的Matlab 解法[OL]./link?url=WBlR0q9n02YsNg6UOLkCQ4Z5 QOCefDKNdEglrNR2TJ8VlxaJ9IkCrLlaEDlJ_OHCLehf19UJU_B7P Re1skQTYaaLcAa1UWcwlv7GYsWNtG7[4]周品，何正风.MATLAB数值分析.[M].北京.机械工业出版社.2009.1.附录：L=1;M=40;N=1600;alfa=1;lambda=0.5;%网格比%**********************************************%h=L/M;%空间步长x=0:h:L;x=x';tao=lambda*h^2/alfa;%时间步长tm=N*tao;%热传导的总时间tm%tm=0.1;t=0:tao:tm;t=t';%计算初值和边值U=zeros(M+1,N+1);U(:,1)=sin(pi*x);U(1,:)=0;U(M+1,:)=0;%*************用差分法求出温度U，与杆长L，时间T的关系****************%for k=1:Nj=2;while j<=MU(j,k+1)=lambda*U(j+1,k)+(1-2*lambda)*U(j,k)+lambda*U(j-1,k);j=j+1;endendlength(U);%***************设置立体网格****************%for i=1:N+1X(:,i)=x;endfor j=1:M+1Y(j,:)=t;endmesh(X,Y,U);legend('数值解');xlabel('X');ylabel('T');zlabel('U');z=zeros(M+1,N+1);for j=1:M+1for k=1:N+1z(j,k)=exp(-pi*pi*t(k))*sin(pi*x(j));endend%mesh(x,t,z')legend('解析解');xlabel('X');ylabel('T');zlabel('Z');for j=1:M+1for k=1:N+1y(j,k)=abs(z(j,k)-U(j,k));endendmesh(x,t,y');legend('绝对误差');xlabel('X');ylabel('Y');zlabel('error');。