偏微分方程数值解52059
第十章 偏微分方程数值解法
第十章 偏微分方程数值解法偏微分方程问题，其求解十分困难。
除少数特殊情况外，绝大多数情况均难以求出精确解。
因此，近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
§1 差分方法的基本概念1.1 几类偏微分方程的定解问题椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程),(2222y x f yu x u u =∂∂+∂∂=∆ 特别地，当0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又称为调和方程2222=∂∂+∂∂=∆yux u u Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(),(),(),(),(2222y x y x u y x y x f y ux u y x ϕ 其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩ称为定解区域，),(y x f ，),(y x ϕ分别为Ω，Γ上的已知连续函数。
第二类和第三类边界条件可统一表示为),(),(y x u u y x ϕα=⎪⎪⎭⎫ ⎝⎛+∂∂Γ∈n 其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件， 0≠α时为第三类边界条件。
抛物型方程：其最简单的形式为一维热传导方程220(0)u ua a t x∂∂-=>∂∂ 方程可以有两种不同类型的定解问题：初值问题⎪⎩⎪⎨⎧+∞<<∞-=+∞<<-∞>=∂∂-∂∂x x x u x t x u a tu )()0,(,0022ϕ初边值问题221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x lu t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩其中)(x ϕ，)(1t g ，)(2t g 为已知函数，且满足连接条件)0()(),0()0(21g l g ==ϕϕ边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类边界条件。
第6章_偏微分方程数值解法
u ( x, 0) = sin( x) ， u (0, t ) = 0, u (π , t ) = 0
利用线上法数值求解 u ( x, t ) 随时间的演化关系 解：取 Δx = π /15 ，计算程序：demo_MOL.m，和结果见右图。
对于 a > 0 ，波从 k − 1 点过来， k − 1 点状态已变化， k + 1 点状态还未变化。差分只能 uk − uk −1 。同样意
n n
义可分析 a < 0 情况。见图 6.1.1。迎风格式的精度为 O(Δt , Δx) ，稳定性条件为 Δt < Δx / | a | 。
% Upwind_Method L = 15;dx = 0.1;dt = 0.05;a = -1.; x =[-L+dx:dx:0]';n=length(x); %Initial value u1=zeros(1,n-20);u2=ones(1,10);u3 = zeros(1,10); u = [u1 u2 u3]';r = a*dt/dx; u0 = u; plot(x,u','LineWidth',2);axis([-15 0 -1 2]);pause(1); for t=dt:dt:10. u(1:n-1) = (1+r)*u(1:n-1)-r*u(2:n); % u(2:n-1)=0.5*((1.-r)*u(3:n)… % +(1.+r)*u(1:n-2)); % Lax scheme hold off;plot(x,u,'LineWidth',2); axis([-15 0 -1 2]);pause(0.05) end hold on; plot(x,u0','r','LineWidth',2);axis([-15 0 -1 2]); xlabel('position');ylabel('u(x,t)'); legend('传播的波','初始方波'); title('Upwind')
偏微分方程数值解(双曲方程书稿)word精品文档14页
双曲型方程的有限差分法线性双曲型方程定解问题： （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） 22222x u a t u ∂∂=∂∂ 其中0>a 是常数。
（1.1）可表示为：022222=∂∂-∂∂x u a t u ，进一步有 0=⎪⎭⎫ ⎝⎛∂∂+∂∂⋅⎪⎭⎫ ⎝⎛∂∂-∂∂u x a t x a t由于x a t ∂∂±∂∂当a dt dx ±=时为()t x u ,的全导数（=dt du dt dx x u t u ⋅∂∂+∂∂xua 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 ∂∂+∂∂=∂∂⋅∂∂+∂∂⋅∂∂=∂∂ xC 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 u t C 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 的任意可微函数。
偏微分方程数值解PPT课件
从（1）得到：
u(ti)u(ti1)hu(ti)O(h)
精选
14
从（2）得到：
u(ti)u(ti1)hu(ti)O(h)
从（1）-（2）得到：
u(ti)u(ti1)2 hu(ti1)O (h2)
从（1）+（2）得到：
u (ti)u (ti 1) 2 u h (2 ti) u (ti 1 ) O (h 2)
精选
15
对经典的初值问题
du
dt
f (t,u )
u ( 0 ) u 0
t (0,T)
满足Lipschitz条件
4
常微分方程的数值解
大气科学中
常微分方程和偏微分方程的关系
1. 大气行星边界层(近地面具有湍流运动特性的大 气薄层,1—1.5km), 埃克曼（V.W.Ekman）（瑞典） 螺线的导出；
2. 1963年，美国气象学家Lorenz在研究热对流的 不稳定问题时，使用高截断的谱方法，由 Boussinesq流体的闭合方程组得到了一个完全确 定的三阶常微分方程组，即著名的Lorenz系统。
2. Curtis F.Gerald and Patrick O., Applied Numerical Analysis, Person Education, Inc., 2004.
3. Eugenia Kalnay, Atmospheric Modeling, Data Assimilation and Predictability, the press Syndicate of the University of Cambridge,2003.
ìïïïïïïïïïíïïïïïïïïî
x1¢=
x
¢
第5章偏微分方程值解ppt课件
t
t nt , x ix , y jy , z kz
总目录
本章目录
5.1
5.2
5.3
5.4
5.2 基本离散化公式
以3对于二阶偏导，我们可以通过对泰勒展开式处 理技术得到下面离散化计算公式：
2u t 2 2u x 2 2u y 2 2u z 2
总目录
本章目录
5.1
5.2
5.3
5.4
5.3 几种常见偏微分方程的离散化计算
例下面介绍3种迭代格式： 1 u (u u u u (1)同步迭代： 4 1 u (u u u u (2)异步迭代： 4 1 u u u ) u (u 4 (3)超松弛迭代：
（5-4） (计算实例VB程序见课本)
总目录
本章目录
5.1
5.2
5.3
5.4
5.3 几种常见偏微分方程的离散化计算
2、一维流动传热传导方程的混合问题 一维流动传热传导方程的混合问题：
2 u u 2 u b f (u, t ) a 2 t x x u t 0 (x), u 0 x x l u x 0 μ1(t)
u
x0
1 (t ),u xt 2 (t )
为初值条件 为边值条件
当该波动方程只提初值条件时，称此方程为波动 方程的初值问题，二者均提时，称为波动方程的 混合问题。
总目录 本章目录
5.1
5.2
5.3
5.4
5.3 几种常见偏微分方程的离散化计算
t t
x
0
x
0
l
(a)初值问题
偏微分方程的数值解方法及源程序
如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程 和定解条件中的已知函数） ，则此定解问题是适定的。可以证明，上面所举各种定解问 题都是适定的。 §2 偏微分方程的差分解法 差分方法又称为有限差分方法或网格法， 是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化 区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替； 通过用网格点上函数的差商代替导数， 将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组（称为差分格式） 。如果差分格式有 解， 且当网格无限变小时其解收敛于原微分方程定解问题的解， 则差分格式的解就作为 原问题的近似解（数值解） 。因此，用差分方法求偏微分方程定解问题一般需要解决以 下问题： （i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式； （iv）讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法 以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。 考虑 Poisson 方程的第一边值问题（3）
（8）
∂u ∂u +a =0 ∂t ∂x
物理中常见的一维振动与波动问题可用二阶波动方程
2 ∂ 2u 2 ∂ u =a ∂t 2 ∂x 2
（9）
（10）
描述，它是双曲型方程的典型形式。方程（10）的初值问题为
2 ⎧ ∂ 2u 2 ∂ u ⎪ 2 =a ∂t ∂x 2 ⎪ ⎪ ⎨u ( x,0) = ϕ ( x) ⎪ ∂u ⎪ = φ ( x) ⎪ ⎩ ∂t t =0
偏微分方程数值解_图文_图文
估计误差
这种误差称为“局部截断误差”，如图。
局部截断误差是以点 的精确解 而产生的误差。
为出发值，用数值方法推进到下一个点
2.整体截断误差—收敛性
整体截断误差是以点 的初始值 为出发值，用数值方法推进i+1步到点
，所得的近似值 与精确值
的偏差：
称为整体截断误差。
特例，若不计初始误差，即 则
即 3.舍入误差—稳定性
五、线性多步（Linear Multistep Method）法
1. 预备知识：插值多项式
插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况， 估算出函数在其他点处的近似值。
从几何上理解：对一维而言，已知平面上n＋1个不同点，要寻找一条n次多项式 曲线通过这些点。插值多项式一般常见的是拉格朗日插值多项式。
把
代入 中，有
经比较得到
取 为自由参数： 从而得到不同的但都是二阶的R-K方法，对应的有中点法、Heun（亨）法 以及改进的Euler法。
基于相同的过程，通过比较五次Taylor多项式，得到更加复杂的结果，给出了包含 13个未知数的11个方程。得到多组系数，其中常用的是以下四阶R-K法：
改进的Euler法、R-K法以及解析解的比较：
是待定的系数。
Euler法就是
的R-K法。
其系数的确定如下：将 展开成 的幂级数，并与微分方程的精确解
在点 的Taylor展开式相比较，使两者的前
项相同，这样确定的R-K法，
其局部截断误差为
，根据所得关于待定系数的方程组，求出它们的值后
代入公式，就成为一个 阶R-K方法。
例题 以二阶R-K法为例说明上述过程
2. Curtis F.Gerald and Patrick O., Applied Numerical Analysis, Person Education, Inc., 2004.
偏微分方程数值解
8.Solve菜单 Solve PDE 对当前的几何结构实体 CSC、三角形网格和图形解偏微分方 程。
Parameters…打开PDE对话框，输入解 PDE的参数。 Expeort Solution… 输出PDE方程的解 矢量u。如果可行，将计算特征值1输 出到主工作区间。
Parameters…打开对话框，可以输入解
在MATLAB的命令窗口中输入pdetool命令， 然后单击回车键，就将显示PDE图形界面，如图：
PDE Toolbox 菜单
1.File 菜单 New 更新或建立一个新的几何结果实体模型 Constructive Soild Geometry(CSG), 并取名为 Untitled。 Open … 从硬盘装载M文件。 Save 将在GUI内完成的成果储存到一个M文件 中。
PDE Specifficatiom…打开一个对话框，输入偏微分方程
类型和应用参数。参数的维数决定于偏微分方程的维数。
如果选择专业应用模式，那么特殊偏微分方程和参数将 代替标准偏微分方程系数。每一个参数c、a、f和d皆可 作为有效的MATLAB表达式，以作为计算三角形单元质 量中心的参数值。下面的变量是很有用的。
二阶线性偏微分方程的分类
自变量多于一个微分方程称为偏微分方程。 一般的二阶线性方程总可以写成如下的形状：
若二阶偏导数项的系数 满足： 则称方程在点
则称方程在点
则称方程在点
在区域G中某点
为双曲型的； 为抛物型的； 为椭圆型的。
二阶偏微分方程的解法的方法主要有两大类，其中 每一类又包含若干种方法。
第一类是解析法 分为分离变量法、保角变换法、镜像法和格林函 数法等。
前处理
网格剖分
取定沿x轴和y轴方向的步长 和 ，
第十章偏微分方程的数值解.doc
第十章偏微分方程的数值解第十章偏微分方程的数值解求解偏微分方程问题非常困难。
除了少数特殊情况，在大多数情况下很难找到准确的解决方案。
因此，近似解更重要。
本章只介绍求解各种典型偏微分方程定解问题的差分方法。
(1)差分方法的基本概念 1.1几种偏微分方程固定解的椭圆方程；最典型和最简单的形式是泊松方程。
特别是，在那个时候，它是拉普拉斯方程，也称为调和方程。
泊松方程的第一个边值问题是一个有边界的有界区域，一条分段光滑曲线，称为固定解区域，以及已知的连续函数。
第二类和第三类边界条件可以统一表示为边界的外法线方向。
当时，这是第二种边界条件，当时，这是第三种边界条件。
抛物线方程:在最简单的形式中，一维热传导方程可以有两种不同类型的定解:初值问题初边值问题在初边值问题中，是一个已知的函数，满足连接条件，边界条件称为第一类边界条件。
第二个和第三个边界条件在其中。
当时，它是第二类边界条件，否则它被称为第三类边界条件。
双曲线方程:的最简单形式是一阶双曲线方程。
物理学中常见的一维振动和波动问题可以用二阶波动方程来描述，二阶波动方程是双曲方程的一种典型形式。
方程的初值问题是一个边界条件，一般有三种类型。
最简单的初边值问题是1.2差分法。
差分法的基本概念也称为有限差分法或网格法。
它是求解偏微分方程定解问题数值解最广泛使用的方法之一。
其基本思想是:首先，对求解区域进行网格划分，用一组有限离散点(网格点)代替自变量的连续变化区域。
问题中出现的连续变量的函数被定义在网格点上的离散变量的函数所代替。
通过用网格点上函数的差商代替导数，将具有连续变量的偏微分方程的固定解问题转化为只有有限个未知数的代数方程(称为差分格式)。
当网格为无穷大时，如果差分格式有解，并且其解收敛于原微分方程的解，则差分格式的解作为原问题的近似解(数值解)。
因此，在用差分法求偏微分方程的定解时，通常需要解决以下问题:(1)选择网格；(2)选择微分方程和固定解条件的差分逼近，列出差分格式；(3)求解差分方案；(4)讨论微分方程差分格式解的收敛性和误差估计。
偏微分方程数值解
ei1u(ti1)ui1h 22u()O (h2)
局部截断误差是以点 t i 的精确解 u ( t i ) 为出发值，用数值方法推进到下一个点
t i 1 而产生的误差。
20
2.整体截断误差—收敛性
整体截断误差是以点 t 0 的初始值 u 0 为出发值，用数值方法推进i+1步到点
t i 1 ，所得的近似值u i 1 与精确值u ( t i 1 ) 的偏差： i1u(ti1)ui1
12
2.差分格式求解 将积分方程通过差分方程转化为代数方程求
解，一般常用递推算法。
在常微分方程差分法中最简单的方法是 Euler方法，尽管在计算中不会使用，但从 中可领悟到建立差分格式的技术路线，下 面将对其作详细介绍:
13
差分方法的基本思想“就是以差商 代替微商”
考虑如下两个Taylor公式：
x2
9 x2 3 x1x3
x3
5 x3 7 x1x2
x4
5 x4 x1x6
x5
x5 3 x1x4
4 x4x5 Re
11
欧拉法—折线法
• 常微分方程能直接进行积分的是少数，而多数是 借助于计算机来求常微分方程的近似解；
• 有限差分法是常微分方程中数值解法中通 常有效 的方法；
• 建立差分算法的两个基本的步骤： 1. 建立差分格式，包括：a. 对解的存在域剖分； b. 采用不同的算法可得到不同的逼近误差—截断 误差（相容性）；c.数值解对真解的精度—整体 截断误差（收敛性）；d.数值解收敛于真解的速 度；e. 差分算法—舍人误差（稳定性）.
ui1ui hf(ti,ui)
预估
u i 1 u i h /2 [f( ti,u i) f( ti 1 ,u i 1 ) ] 校正
数值解偏微分方程的方法和应用
数值解偏微分方程的方法和应用数值解偏微分方程（Numerical Methods for Partial Differential Equations）是一种通过离散化空间和时间域来近似解析解的方法。
它在科学、工程和计算机领域中得到广泛应用。
本文将介绍数值解偏微分方程的基本原理和一些常见的方法，并探讨其在实际问题中的应用。
一、求解偏微分方程的基本原理偏微分方程是包含未知函数及其偏导数的方程，通常用于描述动力学、传热传质、流体力学等现象。
求解偏微分方程的解析解往往十分困难，因此需要借助数值方法来近似求解。
数值解偏微分方程的基本原理是将连续的空间和时间域划分为离散的网格，通过有限差分、有限元或谱方法等离散化技术，将偏微分方程转化为代数方程组。
通过求解这个方程组，可以得到偏微分方程的数值解。
二、常见的数值解偏微分方程方法1. 有限差分法（Finite Difference Method）：有限差分法是最常见也是最简单的数值方法之一。
它通过用中心差分逼近导数，将偏微分方程转化为代数方程组。
有限差分法易于理解和实现，广泛应用于求解各类偏微分方程。
2. 有限元法（Finite Element Method）：有限元法利用有限维空间的函数空间来逼近偏微分方程的解。
它将求解域分解为离散的有限元，将偏微分方程转化为一个求解未知函数系数的代数方程组。
有限元法适用于各种复杂的几何形状和边界条件，广泛应用于结构力学、流体力学等领域。
3. 谱方法（Spectral Method）：谱方法使用一组基函数的线性组合来逼近偏微分方程的解。
它利用高阶多项式函数的收敛性质，能够获得高精度的数值解。
谱方法在求解计算流体动力学和传热传质方程等问题中具有重要的应用价值。
三、数值解偏微分方程的应用1. 流体力学：数值解偏微分方程在流体力学领域有着广泛的应用。
通过数值模拟流体的运动和变形过程，可以预测飞机、汽车等工程结构在空气或水中的流动性能，为工程设计和优化提供指导。
偏微分方程的数值解法
偏微分⽅程的数值解法偏微分⽅程的数值解法
主要总结常见椭圆形、双曲型、抛物型偏微分⽅程的数值解法
椭圆偏微分⽅程
拉普拉斯⽅程是最简单的椭圆微分⽅程
∂2u ∂x2+∂2u
∂y2=0
确定偏微分⽅程的边界条件主要采⽤固定边界条件：u|Γ=U1(x,y) 即在边界Γ上给定u的值U1(x,y)五点差分格式
五点差分格式的形式为：
u i+1,j+u i−1,j+u i,j+1+u i,j−1=4u i,j
以u i,j为中⼼向其上下左右做差分，并⽤这些近似的代替u i,j
运⽤五点差分法可以求出下列边值问题
∂2u ∂x2+∂2u
∂x2=0
u(x1,y)=g1(x),u(x2,y)=g2(x)
u(x,y1)=f1(y),u(x,y2)=f2(y)
x1≤x≤x2,y1≤y≤y2
求解过程如下：
对求解区域进⾏分割：将x min≤x≤x max范围内的的x轴等分成NX段，同理将y轴等分成NY段
将边界条件离散到格点上
⽤五点差分格式建⽴求解⽅程，求出各个格点的函数值
程序设计：
实现函数格式为u = peEllip5(nx, minx, maxx, ny, miny, maxy)
变量名变量作⽤
nx x⽅向上的节点数
minx求解区间x的左端
maxx求解区间x的右端
ny y⽅向的节点数
miny求解区间y的左端
maxy求解区间y的右端
u求解区间上的数值解
建⽴边界条件函数
``
{
Processing math: 100%。
偏微分方程数值解
02
常用的数值解法包括有限差分法、有限元法、 谱方法等。
03
数值解法的精度和稳定性是衡量其好坏的重要 指标。
非线性偏微分方程的有限差分法
有限差分法是一种将偏微分方程转化为差分方程的方法，通过在离散点上逼近偏导数，得到离散化的 数值解。
有限差分法的优点是简单直观，易于实现，适用于规则区域。
有限差分法的缺点是对不规则区域适应性较差，且精度较低。
波动问题
谱方法在求解波动问题中也有广泛应用，如 Helmholtz方程、Wave equation等。
固体力学问题
谱方法在求解固体力学问题中也有应用，如 Elasticity equations等。
05
非线性偏微分方程的数值解 法
非线性偏微分方程的解析解法难度
01
非线性偏微分方程的解析解法通常非常复杂，需要深
02
有限差分法的基本思想是将连 续的偏微分方程转化为离散的 差分方程，通过求解差分方程 得到偏微分方程的近似解。
03
有限差分法的精度取决于离散 点之间的间距，间距越小，精 度越高。
一阶偏微分方程的有限差分法
一阶偏微分方程的有限差分法有 多种形式，如向前差分法、向后 差分法和中心差分法等。
中心差分法是向前差分法和向后 差分法的平均值，具有二阶精度 。
通过将微分转化为差分，将原方程转化为离散的差分方程，然后求解差分方程得到近似解。
有限元法
将连续问题离散化，将微分方程转化为线性方程组，通过求解线性方程组得到近似解。
谱方法
利用函数的谱展开来求解偏微分方程，具有高精度和低数值弥散性的优点。
02
有限差分法
有限差分法的原理
01
有限差分法是一种将偏微分方 程转化为差分方程的方法，通 过在离散点上逼近偏微分方程 的解，得到数值解。
4. 偏微分方程的数值解法
§4 偏微分方程的数值解法一、 差分法差分法是常用的一种数值解法.它是在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值. 1. 网格与差商在平面 (x ,y )上的一以S 为边界的有界区域D 上考虑定解问题.为了用差分法求解，分别作平行于x 轴和y 轴的直线族.⎩⎨⎧====jh y y ihx x i i (i ,j =0,±1,±2,…,±n ) 作成一个正方形网格，这里h 为事先指定的正数，称为步长；网格的交点称为节点，简记为(i ,j ).取一些与边界S 接近的网格节点，用它们连成折线S h ，S h 所围成的区域记作D h .称D h 内的节点为内节点，位于S h 上的节点称为边界节点（图14.7）.下面都在网格D h + S h 上考虑问题：寻求各个节点上解的近似值.在边界节点上取与它最接近的边界点上的边值作为解的近似值，而在内节点上，用以下的差商代替偏导数：()()[]()()[]()()()[]()()()[]()()()[]y x u h y x u y h x u h y x u hy x u h y x u y x u h y x u h y u y h x u y x u y h x u h x u y x u h y x u hy u y x u y h x u h x u ,),(,,1,,2,1,,2,1,,1,,122222222++-+-+≈∂∂∂-+-+≈∂∂-+-+≈∂∂-+≈∂∂-+≈∂∂注意， 1︒ 式中的差商()()[]y x u y h x u h ,,1-+称为向后差商，而()()[]y h x u y x u h,,1--称为向前差商，()()[]y h x u y h x u h,,21--+称为中心差商.也可用向前差商或中心差商代替一阶偏导数.2︒ x 轴与y 轴也可分别采用不同的步长h ,l ，即用直线族⎩⎨⎧====jh y y ihx x j i (i,j =0, ±1, ±2 , ) 作一个矩形网格.2. 椭圆型方程的差分方法[五点格式] 考虑拉普拉斯方程的第一边值问题图14.7()()⎪⎪⎩⎪⎪⎨⎧=∈=∂∂+∂∂y x u D y x y ux u S ,,02222μ 式中μ(x ,y )为定义在D 的边界S 上的已知函数.采用正方形网格，记u (x i ,y j )=u ij ，在节点(i ,j )上分别用差商u u u h u u u h i j ij i j i j ij i j -+-+-+-+11211222,,,,,代替2222,yux u ∂∂∂∂，对应的差分方程为u u u h u u u hi j ij i j i j ij i j -+-+-++-+=112112220,,,, （1） 或u u u u u ij i j i j i j i j =+++-+-+141111,,,,即任一节点(i ,j )上u ij 的值等于周围相邻节点上解的值的算术平均，这种形式的差分方程称为五点格式，在边界节点上取()()()h j i ij S j i y x u ∈=,,**μ （2）式中(x i *,y j *)是与节点(i ,j )最接近的S 上的点.于是得到了以所有内节点上的u ij 值为未知量的若干个线性代数方程，由于每一个节点都可列出一个方程，所以未知量的个数与方程的个数都等于节点的总数，于是，可用通常的方法（如高斯消去法）解此线性代数方程组，但当步长不很大时，用高斯消去法将会遇到很大困难，可用下面介绍的其他方法求解.若h →0时，差分方程的解收敛于微分方程的解，则称差分方程为收敛的.在计算过程中，由于进行四则运算引起舍入误差，每一步计算的舍入误差都会影响以后的计算结果，如果这种影响所产生的计算偏差可以控制，而不至于随着计算次数的增加而无限增大，则称差分方程是稳定的.[迭代法解差分方程] 在五点格式的差分方程中，任意取一组初值{u ij }，只要求它们在边界节点(i ,j )上取以已知值μ(x i *,y j *)，然后用逐次逼近法（也称迭代法）解五点格式：()()()()()[]() ,2,1,0411,1,,1,11=+++=+-+-+n u u u u u n j i n j i n j i n j i n ij 逐次求出{u ij (n )}.当(i+1,j )，(i －1,j )，(i ,j －1)，(i ,j+1)中有一点是边界节点时，每次迭代时，都要在这一点上取最接近的边界点的值.当n →∞时，u ij (n )收敛于差分方程的解，因此n 充分大时，{u ij (n )}可作差分方程的近似解，迭代次数越多，近似解越接近差分方程的解.[用调节余数法求节点上解的近似值] 以差商代替Δu 时，用节点(i+1,j )，(i －1,j )，(i ,j+1)，(i ,j －1)上u 的近似值来表示u 在节点(i ,j )的值将产生的误差，称此误差为余数R ij ，即()()()()()ij j i j i j i j i j i R y x u h y x u h y x u y h x u y h x u =--+++-++,4,,,,设在(i ,j )上给u ij 以改变量δu ij ，从上式可见R ij 将减少4δu ij ，而其余含有u (x i ,y j )的差分方程中的余数将增加δu ij ，多次调整δu ij 的值就可将余数调整到许可的有效数字的范围内，这样可获得各节点上u (x ,y )的近似值.这种方法比较简单，特别在对称区域中计算更简捷.例 求Δu =0在内节点A ,B ,C ,D 上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4（图14.8）解 记u (A )=u A ，点A ,B ,C ,D 的余数分别为图14.8－4u A + u B + u c +5=R A u A －4 u B + u D +7=R Bu A－4 u c + u D +3=R C u B + u c －4u D +5=R D以边界节点的边值的算术平均值作为初次近似值，即u A (0)=u B (0)=u C (0)=u D (0)=2.5则相应的余数为：R A =0, R B =2, R C = －2, R D =0最大余数为±2.先用δu C =－0.5把R C 缩减为零，u C 相应地变为2，这时R A , R D 也同时缩减(－0.5)，新余数是R A =－0.5,R B =2,0=C R , R D =－0.5.类似地再变更δu B =0.5，从而 u B 变为3，则得新余数为0====D C B A R R R R .这样便可消去各节点的余数，于是u 在各节点的近似值为：u A =2.5, u B =3, u C =2, u D =2.5现将各次近似值及余数列表如下：[解重调和方程的差分方法] 在矩形D (x 0≤x ≤x 0+a ,y 0≤y ≤y 0+a )中考虑重调和方程024*******=∂∂+∂∂∂+∂∂=yu y x u x u u ∆取步长h an=，引直线族⎩⎨⎧+=+=jh y y ihx x 00 (i , j = 0, 1, 2,, n ) 作成一个正方形网格.用差商代替偏导数()()()()()[]{()()()()[]()()()()[]}h y x u h y x u y h x u y h x u h y h x u h y h x u h y h x u h y h x u h y x u h y x u y h x u y h x u y x u 2,2,,2,2,,,,2,,,,8201,-+++-++---++-+-++++--+++-++= 上式表明了以(x ,y )为中心时，u (x ,y )的函数值与周围各点函数值的关系，但对于邻近边界节点的点(x ,y )，如图14.9中的A ，就不能直接使用上式，此时将划分网格的直线族延伸，在延伸线上定出与边界距离为h 的点，称这些点为外邻边界节点，如图14.9以A 为中心时，点E ,C 为边界节点，点J ,K 为E ,C 的外邻边界节点，用下法补充定义外邻边界节点J 处函数的近似值u J ，便可应用上面的公式.1︒ 边界条件为()()()S P P x uP u SS ∈==21,μ∂∂μ图14.9时，定义u J =u A －2μ2(E )h .2︒ 边界条件为()()()S P P x uP u SS ∈=∂∂=2221,μμ时，定义u J =2μ1(E )－u A －h 2μ2(E ). [其他与Δu 有关的网格]1︒ 三角网格（图14.10(a )）取P 0(x ,y )为中心，它的周围6个邻近节点分别为：()()⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛---⎪⎪⎭⎫⎝⎛+-⎪⎪⎭⎫ ⎝⎛+++h y h x P h y h x P y h x P h y h x P h y h x P y h x P 23,2,23,2,,23,223,2,,654321 则 R u h u u u h i i +∆+∆=⎪⎭⎫⎝⎛-∑=226102161632式中u i =u (P i ), u 0=u (P 0)，R 表示余项. 2︒ 六角网格（图14.10(b )）取P 0(x ,y )为中心，它的三个邻近节点分别为()⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛++h y h x P y h x P h y h x P 23,2,,23,2321则 R u u u h i i +∆=⎪⎭⎫⎝⎛-∑=0312334.图14.103︒ 极坐标系中的网格（图14.10(c )）取P 0(r ,θ)为中心，它的四个邻近节点分别为()()()()l r P h r P l r P h r P ++--θθθθ,,,,,,4321而拉普拉斯方程01122222=∂∂+∂∂+∂∂=θ∆u r r u r ru u 的相应的差分方程为()()()011221110222134222312=⎪⎭⎫ ⎝⎛+--++++u l r hu u rh u u l r u u h 3. 抛物型方程的差分方法 考虑热传导方程的边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=><<=∂∂-∂∂0,,,,00,0,0,0,021222t t t b u t t u bx x x u t b x x u a tu μμϕ 将[0,b ]分为n 等份，每段长为∆x bn=.引两族平行线（图14.11）图14.11x =x i =i ∆x (i =0,1,2,, n )y =y j =j ∆t (j =0,1,2,, ∆t 取值见后)作成一个长方形的网格，记u (x i ,t j )为u ij ，节点(x i ,t j )为(i ,j )，在节点(i ,j )上分别用(),2,1,1,,2,1Δ2,Δ2,1,11,=-=+---++j n i x u u u t u u ji ij j i ij j i 代替22,xut u ∂∂∂∂，于是边值问题化为差分方程()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-===-==+----++ ,2,1,0,Δ,Δ1,,2,1,Δ,2,1,0,1,,2,10Δ2Δ21002,1,121,j t j u t j u n i x i u j n i x u u u a tu u nj j i j i ij j i ijj i μμϕ 记()22x ta ∆∆=λ，差分方程可写成 () ,2,1,1,,2,121,1,11,=-=+-+=-++j n i u u u u j i ij j i j i λλλ （1） 由此可按t 增加的方向逐排求解.在第0排上u i 0的值由初值ϕ(i ∆x )确定，j +1排u i ,j +1的值可由第j 排的三点(i +1,j ),(i ,j ),(i －1,j )上的值u i +1,j , u ij ,u i -1,j 确定，而u 0,j +1,u n ,j +1已由边界条件μ1((j +1)∆t )及μ2((j +1)∆t )给定，于是可逐排计算一切节点上的u ij 值.当ϕ(x ), μ1(x )和μ2(x )充分光滑，且λ≤12时，差分方程收敛而且稳定.所以利用差分方程（1）计算时，必须使λ≤12，即()22Δ21Δx a t ≤.热传导方程还可用差分方程()0Δ2Δ21,11,1,121,=+---+-++++x u u u a t u u j i j i j i ij j i 代替，此时如已知前j 排u ij 的值，为求第j +1排的u i ,j +1 必须解包含n －1个未知量u u j n j 1111,,,,+-+ 的线性代数方程组，这种差分方程称为隐式格式的差分方程，前面所提的差分方程称为显式格式差分方程.隐式格式差分方程对任意的λ都是稳定的. 4. 双曲型方程的差分方法考虑弦振动方程的第一边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=∂∂=><<=∂∂-∂∂0,,,,00),()0,(,0,0,0,02122222t t t b u t t u b x x t x u x x u t b x x u a tu μμψϕ 用矩形网格，列出对应的差分方程：()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-=∆=∆-==-==+--+--+-+ ,2,1,0,Δ,Δ1,,2,1),(,Δ,2,1,1,,2,1,0Δ2)(Δ22100102,1,1221,1,j t j u t j u n i x i t u u x i u j n i x u u u a t u u u nj j i i i j i ij j i j i ij j i μμψϕ 记ω=a tx∆∆与上段一样，利用u u n 022,和在第0排及第1排的已知数值（初始条件）u i 0 , u i 1可计算u i 2，然后用已知的u i 1 , u i 2及u u n 033,可计算u i 3,类似地可确定一切节点上的u ij 值.当ϕ(x ),ψ(x ),μ1(x )和μ2(x )充分光滑，且ω≤1时，差分方程收敛且稳定，所以要取∆∆t ax ≤1.二、 变分方法1. 自共轭边值问题将§3定义的共轭微分算子的概念推广到一般方程.设D 是n E 中的有界区域，S 为其边界，在D 上考虑2k 阶线性微分方程()x f x x u a Lu km m i i i ni m mi i n n n =∂∂≡∑∑==++201111 ∂ 的齐次边值问题()r j u l S j ,,2,10== 式中f (x )是D 内的已知函数，l j u 是线性微分算子. 将 ⎰DvLud Ω分部积分k 次得()⎰∑⎰⎪⎪⎭⎫ ⎝⎛+=Ω=S j j j D S v R u R v u vLu d ~,Λd k 1 式中Λ(u ,v )是一个D 上的积分，其被积函数包含u ,v 的k 阶导数；R j 和 R j是定义在边界S 上的两个线性微分算子.再将Λ(u ,v )分部积分k 次得()()⎰∑⎰⎪⎪⎭⎫⎝⎛-Ω=Λ=S k j j j D S u R v R v uL v u d ~d ,1***式中L*是一个2k 阶的微分算子，称为L 的共轭微分算子.若L=L*，则称L 为自共轭微分算子.从上面可推出格林公式()()⎰∑⎰=-=Ω-Skj j j jj D S u R v R v R u R v uL vLu 1***d ~~d 如从l j u |S =l j v |S =0可推出在边界S 上()∑==-kj jjjju R v R v R u R 1**0~~ 则称l j u |S =0为自共轭边界条件.如果微分算子及边界条件都是自共轭的，则称相应的边值问题为自共轭边值问题，此时有()0d ][=Ω-⎰DuLv vLu每个边值问题对应于某希尔伯特空间H （例如L 2(D )，见第九章§7）中的一个算子A ，其定义域M A 是H 中一线性稠密集合，它由足够次连续可微且满足边界条件的函数组成，在M A 上，Au 的数值与Lu 的数值相同，从而求解边值问题化为解算子方程Au f =的问题.设A 为定义在实的希尔伯特空间H 中的某线性稠密集合M A 上的线性算子.若对于M A 的任意非零元素,,v u 成立(Au ,v )=(u ,Av )则称A 为对称算子.若对任意非零元素u 成立()0,>u Au则称A 为正算子.如成立更强的不等式(Au ,u )≥r ||u ||2 (r>0)则称A 为正定算子.此处(u ,v )表示希尔伯特空间的内积，||u ||2=(u ,u ). 2. 变分原理与广义解定理 设A 是正定算子，u 是方程Au =f 在M A 上的解的充分必要条件是： u 使泛函F (u )=(Au ,u )－2(f ,u )取极小值.上述将边值问题化为等价的求泛函极值问题的方法称为能量法.在算子的定义域不够大时，泛函F (u )的极值问题可能无解.不过对于正定算子，可以开拓集合M A ，使在开拓了的集合上，泛函的极值问题有解.为开拓M A ，在M A 上引进新的内积[u ,v ]=(Au ,v )，定义模||u ||2=[u ,u ]=(Au ,u )，在模||u ||的意义下，补充极限元素，得到一个新的完备希尔伯特空间H 0，在H 0上，泛函F (u )仍然有意义，而泛函的极值问题有解.但必须注意，此时使泛函F (u )取极小的元素u 0不一定属于M A ，因此它不一定在原来的意义下满足方程Au=f 及边界条件.称u 0为广义解. 3. 极小化序列与里兹方法在处理变分问题中，极小化序列起着重要的作用.考虑泛函F (u )=(Au ,u )－2(f ,u )以d 表示泛函的极小值.设在希尔伯特空间中存在一列元素{u n } (n =1,2 ,)，使()d u F n n =∞→lim则称{u n }为极小化序列.定理 若算子A 是正定的，则F (u )的每一个极小化序列既按H 空间的模也按H 0的模收敛于使泛函F (u )取极小的元素.这个定理不但指出利用极小化序列可求问题的解，而且提供一种近似解的求法，即把极小化序列中的每一个元素当作问题的近似解.设算子A 是正定的，构造极小化序列的里兹方法的主要步骤是：(1) 在线性集合M A 中选取H 0中完备的元素序列{ϕi } ， (i =1,2 ,) 并要求对任意的n ,ϕ1,ϕ2,…,ϕn 线性无关.称这样的元素为坐标元素.(2) 令u a n k k k n==∑ϕ1 ，其中a k 为待定系数.代入泛函F (u )，得自变量a 1,a 2,…,a n 的函数()()()∑∑==-=nj jjn k j kjkj n f a A a a u F 11,,2,ϕϕϕ(3) 为使函数F (u n )取极小，必须()()n j a u F jn ,,2,10 ==∂∂，从而求出a k (k =1,2,…,n ).序列{u n }即为极小化序列，u n 可作为问题的近似解. 4. 里兹方法在特征值问题上的应用 算子方程Au －λu =0的非零解λ称为算子A 的特征值，对应的非零解u 称为λ所对应的特征函数. 对线性算子A ，若存在常数K ，使对任何M A 的元素ϕ成立(A ϕ,ϕ)≥K ||ϕ||2则称A 为下有界算子，正定算子是下有界的（此时K =0）.记(A ϕ,ϕ)/||ϕ||2的下确界为d . 定理1 设A 为下有界对称算子，若存在不为零的元素ϕ0∈M A ，使()d A =200,ϕϕϕ则d 就是A 的最小特征值，ϕ0为对应的特征函数.于是求下有界对称算子的最小特征值问题化为变分问题，即在希尔伯特空间中求使泛函(A ϕ,ϕ)/||ϕ||2取极小的元素，或在||ϕ||=1的条件下求使泛函(A ϕ,ϕ)取极小的元素.定理2 设A 是下有界对称算子，λ1≤λ2≤…≤λn 是它的前n 个特征值，ϕ1,ϕ2,…,ϕn 是对应的标准正交特征函数，如果存在不为零的元素1+n ϕ，在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0, …, (ϕ,ϕn )=0下使泛函(A ϕ,ϕ)取极小，则ϕn +1是算子A 的特征函数，对应的特征值()11,++=n n A ϕϕλ就是除λ1 ,，λn 外的最小的一个特征值.于是求第n +1个特征值就化为变分问题，即在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0 ,, (ϕ,ϕn )=0下求使泛函(A ϕ,ϕ)取极小的元素.为了利用里兹方法求特征值，在M A 中选取一列在H 0中完备的坐标元素序列{ϕi },(i =1,2 ,)， 令u a n k k k n==∑ϕ1，确定a k ，使在条件 (u n ,u n )=1下，(Au n ,u n )取极小，这个问题化为求n个变元a 1,a 2,…,a n 的函数()()∑==nm k m k k m n n a a A u Au 1,,,ϕϕ在条件()()∑===nm k m k m k n n a a u u 1,1,,ϕϕ下的极值问题，一般可用拉格朗日乘数法解（见第九章§3，t ），此时()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n A A A A A A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的最小的根即为特征值的近似值，如果将上式的根按大小排列，就依次得后面的特征值的近似值，但精确度较差. 对一般算子方程Au －λBu=0如果A 为下有界对称算子，B 为正定算子，则()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n B A B A B A B A B A B A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的根就是特征值的近似值. 5. 迦辽金方法用里兹方法解数学物理问题有很多限制，最主要的限制是要求算子正定，但很多问题不一定满足这个条件，迦辽金方法弥补了这个缺陷. 迦辽金方法的主要步骤是：(1) 在M A 中选取在空间H 中完备的元素序列{ϕi } (i =1,2 ,)，其中任意n 个元素线性无关，称{ϕi } (i =1,2,…)为坐标元素序列. (2) 把方程的近似解表示为u a n k k k n==∑ϕ1式中a k 是待定常数，把u n 代入方程Au=f 中的u ，两端与ϕj (j =1,2,…,n )求内积，得 a k 的n 个代数方程()()()n j f A a j nk j kk ,,2,1,,1==∑=ϕϕϕ(3) 求出a k ，代回u n 的表达式，便得方程的近似解u n .在自共轭边值问题中，当算子是正定时，由迦辽金方法和里兹方法得到的关于a k 的代数方程组是相同的.。
第5章偏微分方程数值解
5.1
5.2
5.3
5.4
总目录
5.2 基本离散化公式
在偏微分方程中，自变量都在两个或两个以上，应变量随两个或两 个以上的自变量变化而变化。在化工或化学动态模拟方程中，常常有一个 自变量是时间，其它的自变量为空间位置。如果只考虑一维空间，则只有 两个自变量；如果考虑两维空间，则有3个自变量。一般我们将自变量在 时间和空间以一定的间隔进行离散化，则应变量就变成了这些离散变量的 函数，以3维空间为例，我们将离散化的应变量表示成，它所表示的真正 n 含义如下 : u i , j , k = u ( t , x , y , z ) t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k ∆ z 有了以上的定义，对于一阶偏 u in, +1k − u in, j ,k ∂u j, = 导我们可以利用第四章的欧拉 ∂ t t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆t 公式直接得出向前欧拉公式： u in+1, j , k − u in, j ,k ∂u 对于时间偏导而言，有时我们 = ∂ x t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆x 常常采用向后欧拉公式，时间的 u in, j +1, k − u in, j ,k ∂u 向后欧拉公式如下： = ∂ y t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆y u in, +,1k − u in, j ,k ∂u j = u in, j ,k +1 − u in, j ,k ∂u ∂t t =( n +1) ∆t , x =i∆x , y = j∆y , z = k∆z ∆t = ∂ z t = n ∆ t , x = i ∆ x , y = j ∆ y , z = k∆ z ∆x
第8章 偏微分方程数值解
u ( x j,t k 1 ) u ( x j,t k ) a u ( x j 1 t k ) 2 u ( h x 2 j,t k ) u ( x j 1 ,t k )
[ U]kj o( h2)
(5)
最新课件
35
令
fjk f(xj,tk)
，u
k j
视为 u (xj, tk) 的近似值。
U
u t
a
2u x2
f(x,t)
a 为正常数
（3）
u(x,0)(x)
u(0,t)u(1,t)0
0x1 0tT
(4)
最新课件
34
于结点（j, k）处偏导数与差商之间有如下近似的关系：
u(xj,tk1)u(xj,tk)u tkj o()
u (xj 1tk) 2 u (h x 2j,tk)u (xj 1,tk) x 2 u 2 k j o (h 2) 利用上述表达式得到 LU 在 (j, k) 处的关系式：
在结点上采用离散化方法（数值微分、数 值积分、泰勒展开等）将微分方程的初边值 问题化成关于离散变量的相应问题，这个相 应问题的解就是方程在点xi上的数值解f(x)， 或在点（xi , ti）上的数值解U( xi , ti)。
一般来说，不同的离散化导致不同的方法。
最新课件
10
例：取一边长为1的正方形均匀薄板，上下侧面绝热， 四周保持恒温，求板内各点的稳定温定分布。
u 2u 2u 0 x2 y2
泊松方程
u2u2uf(x,y) x2 y2
最新课件
19
考虑泊松方程第一边值问题：
ux2u2 y2u2 f(x,y), (x,y)
u(x,y), (x,y)面上一有界区域，为其边界，
偏微分方程数值解的计算方法
偏微分方程数值解的计算方法偏微分方程是研究自然和社会现象的重要工具。
然而，大多数偏微分方程很难用解析方法求解，需要用数值方法求解。
本文将介绍偏微分方程数值解的计算方法，其中包括有限差分方法、有限体积法、谱方法和有限元方法。
一、有限差分方法有限差分法是偏微分方程数值解的常用方法，它将偏微分方程中的空间变量转换为网格点上的差分近似。
例如，对于一个二阶偏微分方程：$$\frac{\partial^{2}u}{\partialx^{2}}+\frac{\partial^{2}u}{\partial y^{2}}=f(x,y,u)$$可以使用中心差分方法进行近似：$$\frac{\partial^{2}u}{\partial x^{2}}\approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^{2}}$$$$\frac{\partial^{2}u}{\partial y^{2}}\approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^{2}}$$其中，$u_{i,j}$表示在第$i$行第$j$列的网格点上的函数值，$\Delta x$和$\Delta y$表示网格步长。
将差分近似代入原方程中，得到如下的差分方程：$$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Deltax)^{2}}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Deltay)^{2}}=f_{i,j,u_{i,j}}$$该方程可以用迭代法求解。
有限差分方法的优点是易于实现，但在均匀网格下准确性不高。
二、有限体积法有限体积法是将偏微分方程中的积分形式转换为求解网格单元中心值的方法。
例如，对于如下的扩散方程：$$\frac{\partial u}{\partial t}=\frac{\partial}{\partialx}\left(D(u)\frac{\partial u}{\partial x}\right)$$可以使用有限体积法进行近似。
偏微分方程的数值解方法及源程序
-240-第二十章 偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。
这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。
我们将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。
方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。
如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。
初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。
对于一个具体的问题，定解条件与泛定方程总是同时提出。
定解条件与泛定方程作为一个整体，称为定解问题。
§1 偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。
其最典型、最简单的形式是泊松(Poisson)方程),(2222y x f y ux u u =∂∂+∂∂=Δ （1）特别地，当0),(≡y x f 时，即为拉普拉斯(Laplace)方程，又称为调和方程02222=∂∂+∂∂=Δyux u u （2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。
Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222y x y x u y x y x f y uxu y x ϕ （3）其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩU 称为定解区域，),(),,(y x y x f ϕ分别为ΓΩ,上的已知连续函数。
第二类和第三类边界条件可统一表示成),(),(y x u n u y x ϕα=⎟⎠⎞⎜⎝⎛+∂∂Γ∈ （4） 其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件，0≠α时为第三类边界条件。
在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。
偏微分方程 数值解
对Richardson格式，考察其稳定性。
设 为计算 产生的误差，设 无误差，且 层前和本层中其它点的计算均无误差，而计算过程中也不再产生新的误差。取 ，则 满足：
设 时，中间点的计算有误差 ，即 ， 。则由上式，得 ， ，…… ， 等等，绝对值越来越大。故此格式不稳定。实际上对任意 ，此格式均不稳定。
y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);
u[n][1]=y[n];
for(i=n-1;i>0;i--)
u[i][1]=y[i]-beta[i]*u[i+1][1];
}
int main() //一维热传导方程的Richardson格式
{
int k,i;
printf("xj tk真实值x[i][k]近似值u[i][k]误差err[i][k]\n");
b[1]=1+r;
c[1]=-r/2;
a[n]=-r/2;
b[n]=1+r;
f[1]=r/2*u[2][0]+(1-r)*u[1][0]+r/2+r/2*u[0][1];
f[n]=r/2*u[n+1][0]+(1-r)*u[n][0]+r/2*u[n-1][0]+r/2*u[n+1][1];
0.166667 0.666667 2.300976 1.214156 1.086820
0.333333 0.666667 2.718282 3.293822 0.575541
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
4
常微分方程的数值解
大气科学中
常微分方程和偏微分方程的关系
1. 大气行星边界层(近地面具有湍流运动特性的大 气薄层,1—1.5km), 埃克曼（V.W.Ekman）（瑞典） 螺线的导出；
2. 1963年，美国气象学家Lorenz在研究热对流的 不稳定问题时，使用高截断的谱方法，由 Boussinesq流体的闭合方程组得到了一个完全确 定的三阶常微分方程组，即著名的Lorenz系统。
3. Charney, Fjortoft, and Von Neumann(1950), 借助于Princeton大学的的计算机(ENIAC)，利 用一个简单的正压涡度方程 （C.G.Rossby,1940精）选P对PT 500mb的天气形式作3 了24小时预报---成功；
The Electronic Numerical Integ精ra选toPrPaTnd Computer (ENIAC).
精选PPT
6
50 40 30 20 10
0 -10 -20 -30
0
5
10
15
20 精选P2P5T 30
35
40 45
50
7
50 40 30 20 10
0 20
0
-20 30
20精选PPT10
0
-10
-20
-30
8
精选PPT
9
精选T
10
Franceshini 将Navier-Stokes方程截断为五维的
精选PPT
12
2.差分格式求解 将积分方程通过差分方程转化为代数方程求
解，一般常用递推算法。
在常微分方程差分法中最简单的方法是 Euler方法，尽管在计算中不会使用，但从 中可领悟到建立差分格式的技术路线，下 面将对其作详细介绍:
精选PPT
13
差分方法的基本思想“就是以差商 代替微商”
考虑如下两个Taylor公式：
2. Curtis F.Gerald and Patrick O., Applied Numerical Analysis, Person Education, Inc., 2004.
3. Eugenia Kalnay, Atmospheric Modeling, Data Assimilation and Predictability, the press Syndicate of the University of Cambridge,2003.
精选PPT
5
Lorenz系统
dx / dt = a (y - x) dy / dt = x (b - z) - y dz / dt = xy - c z
其中，a=10,(Prandtl number); b=28(Rayleigh number); c=8/3; (x,y,z)_0=(0.01;0.01;1e-10)
4. Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press,1996.
5. 李荣华，冯国忱. 微分方程数值解. 北京：人民教育出版社，1980. 6. 徐长发，李红. 实用偏微分方程数值解法. 华中科技大学出版社，2003. 7. 沈桐立，田永祥等. 数值天气预报. 北京：气象出版社，2007.
偏微分方程数值解 (Numerical Solution of
Partial Differential Equations)
主讲：王曰朋
eduwyp@
精选PPT
1
参考数目
1. George J. Haltiner, Roger Terry Williams, Numerical Prediction and Dynamic 2. Meteorology(2nd Edition), the United States of America, 1979.
精选PPT
2
数值天气预报—PDE数值解
1. 挪威气象学家V.Bjerknes（1904）提出数值预 报的思想：通过求解一组方程的初值问题可以 预报将来某个时刻的天气—思想；
2. L.F.Richardson(1922)：开创了利用数值积分 进行预报天气的先例，由于一些原因（如，计 算稳定性问题“Courant,1928”）并没有取得预 期的效果—尝试；
• 有限差分法是常微分方程中数值解法中通 常有效 的方法；
• 建立差分算法的两个基本的步骤：
1. 建立差分格式，包括：a. 对解的存在域剖分； b. 采用不同的算法可得到不同的逼近误差—截断 误差（相容性）；c.数值解对真解的精度—整体 截断误差（收敛性）；d.数值解收敛于真解的速 度；e. 差分算法—舍人误差（稳定性）.
精选PPT
15
对经典的初值问题
从（1）得到：
u(ti)u(ti1)hu(ti)O(h)
精选PPT
14
从（2）得到：
u(ti)u(ti1)hu(ti)O(h)
从（1）-（2）得到：
u(ti)u(ti1)2 hu(ti1)O (h2)
从（1）+（2）得到：
u (ti)u (ti 1) 2 u h (2 ti) u (ti 1 ) O (h 2)
u ( t h ) u ( t ) u ( t ) h 1 u ( t ) h 2 1 u ( t ) h 3 L 1 u ( n ) ( t ) h n O ( h n 1 ) （1）
2 ! 3 !
n !
u ( t h ) u ( t ) u ( t ) h 2 1 ! u ( t ) h 2 3 1 ! u ( t ) h 3 L n 1 ! u ( n ) ( t ) h n O ( h n 1 ) （2）
截谱模型如下：
ìïïïïïïïïïíïïïïïïïïî
x1¢=
x
¢
2
=
x
¢
3
=
x
¢
4
=
x
¢
5
=
-
2 x1 + 4 x2 x3 + 9 x2 + 3 x1x3 5 x3 - 7 x1x2 + 5 x4 - x1x6 x5 - 3 x1x4
4 x4x5 Re
精选PPT
11
欧拉法—折线法
• 常微分方程能直接进行积分的是少数，而多数是 借助于计算机来求常微分方程的近似解；