抛物形扩散方程的有限差分法及数值实例
扩散方程的差分解法
（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随时间变化图
抛物型方程的差分方法
抛物型方程的差分方法抛物型方程是描述物理现象中的薄膜振动、热传导、扩散等过程的方程，具有非常重要的应用价值。
差分方法是一种常用的数值计算方法，用于求解微分方程，对于抛物型方程的数值求解也是非常有效的方法之一、本文将介绍抛物型方程的差分方法，并具体讨论用差分方法求解抛物型方程的一些具体问题。
首先，我们来介绍一下抛物型方程的一般形式。
抛物型方程一般可以表示为：∂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）边界条件等。
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）是一种常用的数值求解抛物型方程的方法。
它基于离散化的思想，将偏微分方程转化为差分方程，然后利用中心差分近似代替偏导数，得到差分方程的离散形式。
(十二章)抛物型方程有限差分法
（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
,我们需要求解这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) 逆风差分格式
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
有限差分求解扩散的数值解
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
扩散模型的原理范文
扩散模型的原理范文扩散模型是一种模拟物质传播和扩散过程的数学模型，广泛应用于社会科学、物理学、生物学、工程学等领域。
其基本原理是描述和预测物质在空间和时间上的分布变化规律。
以下将详细介绍扩散模型的原理。
1.扩散过程扩散是指物质由高浓度区域向低浓度区域的自发传播过程。
这种传播可以是由分子、粒子、信息等传递激发的，也可以是由温度、浓度梯度等驱动的。
扩散过程遵循费克定律，即物质扩散的速率与浓度梯度成正比。
2.扩散方程∂C/∂t=D∂²C/∂x²其中，C是物质浓度，t是时间，x是空间位置，D是扩散系数。
该方程描述了物质浓度随时间变化的速率和空间扩散的强度。
对于多维情况，扩散方程的形式会有所不同。
3.初始条件与边界条件扩散模型还需要指定初始条件和边界条件。
初始条件是指在初始时刻，物质浓度在空间上的分布情况。
边界条件则是指在边界位置上，物质浓度的变化规律。
这些条件通常通过实验测量或数值模拟获得。
4.数值求解由于扩散方程是偏微分方程，对于大多数情况，无法求得解析解。
因此，需要使用数值方法来近似求解扩散方程。
常见的数值方法包括有限差分法、有限元法等。
这些方法将空间和时间分割成小的网格，通过在网格上进行计算来模拟物质扩散过程。
5.扩散系数和扩散机制扩散系数是一个重要的参数，它反映了物质在扩散过程中的扩散速率。
扩散系数的大小取决于物质本身的性质，以及扩散过程中的环境条件。
不同的扩散机制也会导致不同的扩散系数。
例如，在气体中，扩散主要由分子碰撞引起；在液体中，扩散主要由颗粒之间的扩散引起。
6.扩散模型的应用扩散模型在实际中有广泛的应用。
在社会科学中，扩散模型可以用于描述和预测信息、流行病等的传播过程。
在物理学中，扩散模型可以用于研究热传导、原子扩散等现象。
在生物学中，扩散模型可以用于模拟细胞内物质的运输过程。
在工程学中，扩散模型可以用于设计和优化化学反应器、催化剂等。
总结：扩散模型是一种用于模拟物质传播和扩散过程的数学模型。
第四章 抛物型方程的有限差分方法
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
抛物型方程的有限差分解法及其在复杂电磁环境中的应用
II
with the results of Fourier method and AREPS. The research achievements are the use of FDM to solve parabolic equation with accurate and available results, which is demonstrated through the comparison with other algorithms. The new ideas are proposed as follows: a method to improve the existing absorbing layer, establishment of some GPS initial field models which are suitable in different conditions with circular polarization wave included. In addition, conclusions of low gazing GPS signals propagation on sea surface are summarized. Key Words: Parabolic wave equation, Finite difference method, Inverse algorithm, Atmospheric duct, GPS initial field
研究生姓名 指导教师 姓名 单位名称 申请学位级别 论文提交日期 学位授予单位 答辩委员会主席
430070 无线电物理 2010 年 5 月 2010 年 5 月 王嘉赋 徐晓英
学科专业名称 论文答辩日期 学位授予日期 评阅人
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 格式。
有限差分法求解抛物型方程说明
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例.§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 层上的值求出。
抛物型方程的有限差分方法
抛物型方程的有限差分方法一，求解问题考虑一维非齐次热传导方程的定解问题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 +≤≤。
偏微分方程数值解法抛物型方程差分法PPT课件
C-N 格式矩阵形式
[(1 ra2 )I ra 2 C]uk1
2 [(1 ra2 )I
ra 2
C ]uk
(fk
f k1 )
2
2
H [(1 ra2 )I ra2 C]1[(1 ra2 )I ra2 C]
2
2
特征值
j
(1 (1
ra2 ) ra2 )
ra 2 ra 2
cos(j cos(j
j
/
2(n
1))
|
(H) 1
C-N 格式是无条件稳定的.
第13页/共17页
数值实验题 用三种差分格式求 初边值问题数值解
ut uxx , 0 x 1, t 0
u( x,0) si nx, 0 x 1
u(0, t ) u(1, t ) 0, t 0
并与准确解比较
u(x, t) exp( 2t)sin x
第16页/共17页
感谢您的观看！
第17页/共17页
显格式 隐格式
1
[ukj
1
ukj ]
a2 h2
2 x
ukj
1
[ukj
1
ukj ]
a2 h2
u2 k1 xj
C-N格式
1
[ukj
1
ukj ]
a2 2h2
2 x
(ukj
1
ukj )
第14页/共17页
数值计算实验
显格式: input T:=1 error = 7.9443e-006 k = 200
其中 M 与 无关. k > k0
第4页/共17页
简单显式差分格式
uk1 j
(1
- 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 (常数)是扩散系数。
取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 网格边界点的集合。
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， k u 0=k N u =0 (1.2.3)计算后得：111(12)k k k k j j j j j u 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 ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨⎪⎪⎪=+-++⎩ (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 φ+⎧=+=-⎨=⎩u Au f u (1.2.6)其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=r r r r r r r rr r21002100210021A而对于向前差分格式，当网比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 N hτ++---+==-=边值条件：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 t数值解精确解相对误差0.1 0.1 0.11510.115241.170010-⨯0.2 0.2 0.08150.08164⨯1.658010-0.3 0.3 0.04180.04194⨯1.275210-0.4 0.4 0.01830.01845⨯7.445610-0.5 0.45 0.01170.01185⨯5.375510-0.6 0.35 0.03000.030141.067410-⨯0.7 0.25 0.06840.06864⨯1.741010-0.8 0.15 0.50840.50855⨯7.589710-0.9 0.05 0.36540.36545⨯1.740710-表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在410-以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。
下面，我们比较在不同的分割下对有限差分算法精度的影响。
在扩散系数1a=不变的情况下，讲时间和空间进行更加细密的分割，取==，其中，M表示空间上的分割，N表示时间上的分割。
观察数M N50,10000值解与精确解在节点,x t处的绝对误差值，如下图所示：j k图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在410-和510-以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。
最后，我们比较网比1/2r=以及1r=时扩散方程的收敛情况。
当网比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_V547zVS6UwenT8rGKsbwNf9CuDmh2qmsy5K8eQ 32PzhYazZ_sWiHfz_Pj3LA7ufHH4O3t8NIlUCnXNUiyYhssqmNBeArsrLwQG[3]未知.偏微分方程的Matlab解法[OL]./link?url=WBlR0q9n02YsNg6UOLkCQ4Z5QOCefDKNdEglrNR2TJ8Vlxa J9IkCrLlaEDlJ_OHCLehf19UJU_B7PRe1skQTYaaLcAa1UWcwlv7GYsWNtG7[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');。