第2章椭圆型方程的有限差分法
§1 差分逼近的基本概念
.
,;
0,],[,)2.1()(,)()1.1(,22
为给定常数上的连续函数为其中边值问题考虑二阶常微分方程的βαβ
α≥?????==<<=+-=q b a f q b u a u b x a f qu dx u
d Lu
区间剖分
.
],[.
/)(,,2,1,0],[称为步长，间距称为网格结点（节点）的一个网格剖分，于是我们得到区间等分，分点为分成将区间h x b a I N a b h N i ih a x N b a i i =-==+=
微分方程离散(差分方程）
.][)3.1(),(])([12])([)
()(2)()1.1(34
4222211点取值表示括号内函数在其中展式可得
，由的解离散化，对充分光滑在节点现在将方程i i i i i i i i x h O dx
x u d h dx x u d h
x u x u x u Taylor u x ++=+--+)5.1(),(])([12)()4.1(),()()()()()(2)()1.1(3
4
4
22
11h O dx
x u d h u R u R x f x u x q h
x u x u x u x i i i i i i i i i i +-=+=++---+其中写成
可将方程于是在
.
)6.1()()(][).(),()6.1(,2)1.1()(.)(2
1
1的截断误差为差分方程称，记式中的差分方程：，则得逼近方程若舍去的二阶无穷小量是足够小，当u R x f Lu x f f x q q f u q h
u u u u L u R h u R h i i i i i i i i i i i i i i h i i ====++--=-+.
)()7.1(][)()(截断误差所引起的代替微分算子是用差分算子，截断误差L L u R Lu x u L u R h i i i h i -=
.
.)2.1(),1.1()9.1(),8.1(.)()
9.1(.,)8.1(,1,2,1,2,1,2,1)6.1(02
1
1式此格式称为中心差分格的差分方程或差分格式为逼近称的近似于是它的解的线性代数方程组：于加上边值条件就得到关时成立，
当差分方程i i N i i i i i i i h i x x x u u u u N i f u q h
u u u u L u N i ===-==++--=-=-+βα .
1,,,,)8.1(:121阶方程组因此它是个数的的个数等于网格内点方程注意--N x x x N
§2 一维差分格式
.
,],
,[,,,0)(],,[)2.2()(,)()1.2(,)(min 1为给定常数其中考虑两点边值问题：
βαβ
αb a C f q r p x p b a C p b u a u b x a f qu dx du r dx du p dx d Lu ∈>≥∈?????
==<<=++-=积法
直接差分化法、有限体两种方法：我们将介绍差分格式的
直接差分化
.
,2,1,:],[,1110N i x x x I N b a I b x x x x a N i i i N i =≤≤==<<<<<=+-个小区间：
分成将区间个节点：
首先取.
,,,,max ,01211的集合表示内点和界点的集合表示网格内点
为最大网格步长。用称的一个网格剖分，记于是得到区间b x a x I x x x I h h x x h I N h N h i i
i i i ===-=--
.
],[,
),,2,1)((2
1
,2
12
12
32
1012
11对偶剖分的一个网格剖分，称为又构成点
称为半整数点，则由节，的中点取相邻节点b a b x x
x x x x a N i x x x
x x N N i i i i i i =<<<<<<<==+=-
---
- 点取值。
表示括号内函数其中展式可得
，由为此，对充分光滑的解离散化，在节点方程其次用差商代替微商将i i i i i i i i i i i x h o dx
u d h h dx du
h h x u x u Taylor u x ][)3.2(),(][2][)
()()1.2(222
11
11+-+=+-++-+
)
4.2(),(][24][),
(][24][)()()(3
33221321332
2112
1h o dx
u d p h dx du p h o dx u
d p h dx du p h x u x u x p i i i i i i i i i i ++=++=------)
5.2(),(][24][)()()(3331221112
1h o dx
u
d p h dx du p h x u x u x p i i i i i i i ++=-+++++
)6.2(),()][12)]([4)]([),(][12)][]([2]
)()()()()()([2
2
)4.2()5.2(233
122
1233
121211*********
h o dx
u d p h h dx du p dx d h h dx du p dx d h o dx
u d p h h dx du
p dx du p h h h x u x u x p h x u x u x p h h h h i i i i i i i i i i i i i i i i i i i i i i i i i i +-+-+=+-+-+=---+++++-++--+++++，得
，并除以减由
满足方程：
边值问题的解知则由令)(,)6.2(),3.2(),
(),(),(),(2
12
1x u x f f x q q x r r x
p p
i i i i i i i i ====--)
7.2(),()()]()([]
)
()()()([2
)(11112
1
1
121
1u R f x u q x u x u h h r h x u x u p
h x u x u p h h x u L i i i i i i i i i
i
i i i i i i i i i i h +=+-++---+-=-++--++++
.
)2.2(),1.2(),()
8.2(),()][21][121)]([41)(()(222
33221的差分方程便得逼近边值问题的截断误差，舍去为差分算子其中u R L h o dx
u
d r dx u d p dx du
p dx d h h u R i h i i i
i i i +-+--=+)
10.2(,
,)
9.2(,1,,1,][]
[2011112
111211βα==-==+-++---+-=-++--++++N i i i i i i i i
i i i i i i i i i i i h u u N i f u q u u h h r h u u p h u u p h h u L
有限体积法
恒律具有形式
上的热量守内任一小区间程，则在方一根杆上的稳定温度场如果把它看作是分布在，考虑守恒型微分方程：
],[],[)
13.2()()()()
2()1(x x b a x f u x q dx
du
p dx d Lu =+-=)
15.2()()()
14.2(,)()(,
)()
2()1()
2()1()
2()
1()2()1()2()1()
2()
1(，其中或dx
du
x p x W fdx qudx x W x W fdx qudx dx dx du
p dx d x x
x x
x x x x x x ==+-=+-?????把微分方程写成积分守恒型后，最高阶微商由二阶降到一阶，从而可减弱对函数光滑性的要求。
,
)()15.2()()14.2(,)()(],
,[],[)14.2(2
12
12
12
12
12
12
12
1)2()1(恒连续流量”化是不合适的。但“热进一步差分允许有间断点，由考虑到则对偶单元取特别于x W x p fdx qudx x
W x
W x
x x x i i i i x
x
x
x
i i i i '
=+-?
?
+-
+
-
+-+-,
)
()
(],[,)()()15.2(111?-=-=--i
i x x i i i i dx x p x W u u x x x p x W dx du 积分，得
再沿改写成
故将
)
18.2(,
2)17.2(.])
(1[
)16.2(,
1
1
12
12
12
11
i i i i x
x x x i
i i
i i i i u d h h qudx x p dx h a h u u a W
i i i
i +---+≈=-≈??
+--又
利用中矩形公式，得)
21.2(.)(2)20.2(,)(2
1
,
)(2
1][)14.2()18.2(),16.2(212
1
1111111?+
-+++-++++=+=++----'i i x x i i i i i i i i i i i i i i i i i i dx x f h h h h u d h h h u u a h u u a ??，即得守恒型差分方程代到将
)
22.2(),
(),(),
()21.2()19.2(),17.2(,2121???
?
???======-
-i i i i i i i i i x f f x q q d x p p a f q p ?，从而和计算式光滑，则可用中矩形公及右端如果系数)
23.2(),
(21),
(21,22121212111??
?
?
?
?
???+=+=+=
+-+---i i i i i i i
i i
i i f f f q q d p p p p a 也可用梯形公式，此时
§3矩形网的差分格式
条件之一：
在边界上满足下列边值，其边界为分段光滑曲线是平面上一有界区域，方程
考虑G G y x y x f u Poisson )
1.3(),(),,(∈=?-32
1)1.3()(),()1.3()(),()1.3()(),(第三边值问题第二边值问题第一边值问题y x ku n u
y x n u
y x u γβα=+??=??=Γ
Γ
Γ连续函数。
都是及其中0),(),(),,(),,(),,(≥y x k y x y x y x y x f γβα
3.1五点差分格式
,1,0,,1,0,.)(,212122
2
1
21±==±==+=j jh y i ih x h h h h h y x 直线：
作两族与坐标轴平行的和轴方向的步长轴和取定沿.11,),(),().
,(),(),(2
121='-+'-=-+-'
'''j j i i h y y h x x y x y x j i y x jh ih i i i i j i j i j i 或如果
是相邻的和说两个节点或记为称为网点或节点，两族直线的交点
{}否则称为非正则内点。
就称为正则内点相邻点都属于的四个若内点的网点的集合是代替域就则令点为界点交点集合，并称如此的的与或表示网线以并称如此节点为内点内部的节点集合，表示所有属于以;,),(.,..),(h j i h h h h j i h j i h G y x G G G G G y y x x G G y x G Γ=Γ=Γ==Γ∈= )2.3(,
)2.3(),,(),(,),(,),()
2.3(,]22[
,,),(22
1
,1,2
1
,1,1'
=?-====+-+
+--=?--+-+h h h j i ij j i h ij j i h h h ij ij j i ij j i j
i ij j i ij h yy xx j i f u y x f f y x f u y x u f u j i u f h
u u u h
u u u u u u y x y x 可简写成
则差分方程表示网格函数，上的网函数。若以表示节点式中，则得用二阶中心差商代替方向分别为正则内点，沿现在假定
用差分法解椭圆型偏微分方程 -（Uxx+Uyy)=(pi*pi-1)e^xsin(pi*y) 0 抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a 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 = 为空间步长，M T = τ为时间步长，其中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 ,处的微商和差商之间的下列关系 （(,)k j k j u u x t t t ????≡ ? ????）： ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式 数值实验报告I 实验名称Poisson方程九点差分格式实验时间2016年 4 月 15 日姓名米瑞琪班级信息1303学号04成绩 一、实验目的，内容 1、理解Poisson方程九点差分格式的构造原理； 2、理解因为网格点的不同排序方式造成的系数矩阵格式的差异； 3、学会利用matlab的spdiags(),kron()函数生成系数矩阵； 二、算法描述 针对一个Poisson方程问题： 在Poisson方程五点差分格式的基础上，采用Taylor展开分析五点差分算子的截断误差，可以得到： 为了提高算子截断误差的精度，在(1)式中配凑出了差分算子的形式，将原Poisson方程代入（1）式有： 考虑，有： 将（3）代回（2）可得 得到Poisson方程的九点差分格式： 在计算机上实现（4）式，需要在五点差分格式 的基础上在等式两端分别增加一部分，将等式左侧新增的部分写成紧凑格式，有： 对于该矩阵，可以看成是两个矩阵的组合： 以及 则生成这两个矩阵可以采用Kroncker生成，方法类似于五点差分格式。 对于右端添加的关于f(x,y)的二阶导数，可以采用中心差分格式进行近似代替，即： 写成相应的紧凑格式有： 该式中的矩阵又可以分解为两个矩阵的和： %计算误差 u_real=@(x,y)exp(pi*(x+y))*sin(pi*x).*sin(pi*y); for i=1:N1-1 u_m((i-1)*(N2-1)+1:i*(N2-1))=u_real(x(i),y); end u_v=u_m'; err_d=max(abs(u_d-u_v)); sol=reshape(u_d,N2-1,N1-1); mesh(X,Y,sol) 四. 数值结果 针对课本P93给出的问题，分别采用步长，将计算出的误差列表如下： 步长五点差分格式误差九点差分格式误差 可见采用九点差分格式可以进一步缩小误差，达到更高阶的精度。 五. 计算中出现的问题，解决方法及体会 在生成九点差分格式的时候，等号右端涉及到了对f的二阶偏导，我最初利用符号函数定义了f,随后求出其二阶偏导（仍然是符号函数）之后带入网格点，求f二阶偏导的精确解，但是代入过程相当繁琐，运行速度非常慢，最终我改变策略，选用f关于x,y的二阶中心差分格式替代精确值，最终得到了相对满意的结果。 教 师 评 语 指导教师：年月日 在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。 偏微分方程基本形式 而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分 5) 有限元求解 6) 解的图形表达 以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适 当的模式进行建模和分析。 在Options菜单的Application菜单项下可以做选择，如下 或者直接在工具栏上选择，如下 列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。 ④ Structural Mech.，Plane Strain：结构力学平面应变。 ⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。 可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外，对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。 另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系范围 l Axes Equal：同Matlab的命令axes equal命令 建立几何模型 使用菜单Draw的命令或使用工具箱命令可以实现简单几何模型的建立，如下 各项代表的意义分别为 用五点有限差分格式求解椭圆型方程（偏微分方程）程序2 2010-04-29 10:33 function varargout=liu(varargin) a=0;b=2;c=0;d=1;h1=1/16;h2=1/16; f=inline('(pi^2-1)*exp(x)*sin(pi*y)','x','y'); g1x=inline('0'); g2x=inline('0'); g1y=inline('sin(pi*y)'); g2y=inline('exp(2)*sin(pi*y)'); [X,Y,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h1,h2); mesh(X,Y,Z); shading flat; xlabel('X','FontSize',14); ylabel('Y','FontSize',14); zlabel('error','FontSize',14); title('误差图'); function [X,T,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h1,h2) %求解下问题 %-(u_xx+u_yy)=f(x,y) x,y 在区域内x in a %h2离散y方向的步长 N=10000; x=a:h1:b; y=c:h2:d; m=length(x); n=length(y); ee=0.00001; [X,T]=meshgrid(x,y); Z=zeros(n,m); U=zeros(n,m); for i=2:m-1 U(1,i)=feval(g1x,x(i)); U(n,i)=feval(g2x,x(i)); end for j=1:n U(j,1)=feval(g1y,y(j)); U(j,m)=feval(g2y,y(j)); end %while true %下为高斯赛德尔迭代法 %---------------------------------------------------------------------- for k=1:N 对流扩散方程有限差分方法 求解对流扩散方程的差分格式有很多种，在本节中将介绍以下3种有限差分格式：中心差分格式、Samarskii 格式、Crank-Nicolson 型隐式差分格式。 3.1 中心差分格式 时间导数用向前差商、空间导数用中心差商来逼近，那么就得到了（1）式的中心差分格式]6[ 2 1 11 1122h u u u v h u u a u u n j n j n j n j n j n j n j -+-+++-=-+-τ （3） 若令 h a τ λ=，2h v τ μ=，则（3）式可改写为 )2()(2 111111 n j n j n j n j n j n j n j u u u u u u u -+-+++-+--=μλ （4） 从上式我们看到，在新的时间层1+n 上只包含了一个未知量1 +n j u ，它可以由时间层n 上的值n j u 1-，n j u ，n j u 1+直接计算出来。因此，中心差分格式是求解对 流扩散方程的显示格式。 假定),(t x u 是定解问题的充分光滑的解，将1 +n j u ，n j u 1+，n j u 1-分别在),(n j t x 处 进行Taylor 展开： )(),(),(211ττO t u t x u t x u u n j n j n j n j +??? ?????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +????????+????????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +????????+????????-==-- 代入(4)式，有 2 111 1122),(h u u u v h u u a u u t x T n j n j n j n j n j n j n j n j -+-+++---+-= τ )()()(2222 h O v x u v h O a x u a O t u n j n j n j ?-????????-?+????????++????????=τ )()()(222h O v a O x u v x u a t u n j n j n j ?-++????????-??? ?????+????????=τ 课程名称（中文）：偏微分方程数值解专题 课程名称（英文）：Some topics on numerical solutions of partial differential equations 一）课程目的和任务：有限差分方法是微分方程定解问题的最广泛的数值方法之一，其基本思想是用差商近似代替导数，用有限个未知量的差分方程组的解作为微分方程定解问题的解。本课程旨在介绍非线性抛物和椭圆型方程的有限差分方法的最新进展，并简单介绍在实际模型问题中的应用。 本课程是为计算数学专业二年级研究生开设的一门专业选修课程，也可作为其它专业研究生的选修课程。其主要目的是使学生了解不同非线性定解问题的有限差分格式的构造方法，掌握有限差分方法的基本理论和一些数值分析技巧，使他们对这一方法能够有清晰和全面的了解，并能应用于实际问题的数值模拟和数值计算。通过对本课程的学习，进一步丰富微分方程定解问题的数值方法，掌握有限差分方法的最新进展，也为学习进一步的专业知识打下坚实的基础。 主要内容：本课程以非线性抛物型和椭圆型方程的定解问题为例，介绍有限差分格式的构造方法，同时介绍相应数值分析的基本方法和技巧，也简单讨论有限差分解的一些性质。所有内容基本都是本课题的最新研究进展。具体内容如下： 第一章：非线性椭圆边值问题的有限差分方法。本章主要介绍非线性椭圆边值问题的差分解法，给出一些定性分析，包括解的存在唯一性、收敛性和误差估计，同时介绍分析非线性差分方程组的上下解方法以及求解的单调迭代算法； 第二章：非线性抛物初边值问题的有限差分方法。本章主要介绍非线性抛物边值问题的差分解法，包括一些定性分析，如解的存在唯一性、收敛性和误差估计，同时介绍分析非线性差分方程组的上下解方法以及求解的单调迭代算法； 第三章：依赖于时间数值解的渐近性。本章着重介绍非线性抛物初边值问题有限差分解的长时间渐近性，稳定性及在实际模型问题中的应用； 第四章：非线性椭圆方程组的有限差分方法。本章主要介绍非线性椭圆型方程组的有限差分解法，着重介绍分析非线性差分方程组的上下解方法以及求解的单调迭代算法。 第五章：非线性抛物方程组的有限差分方法。本章主要介绍非线性抛物型方程组的有限差分解法，介绍分析非线性差分方程组的上下解方法以及求解的单调迭代算法。除此之外，还介绍特别介绍有限差分解的长时间渐近性，稳定性及在实际模型问题中的应用。 二）预备知识：数值代数、数值逼近、微分方程数值解、泛函分析。 三）教材及参考书目： 教材：自编讲义，2012. 参考书目：1) C.V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, New York, 1992. 2) A. Berman, R. Plemmons, Nonnegative Matrix in the Mathematical Science, Academic Press, New York, 1979. 四）讲授大纲（中英文） 第一章非线性椭圆边值问题的有限差分方法 1)有限差分方程 2)上下解方法 3)单调迭代算法 4)上下解的构造方法 利用中心差分格式数值求解导数 目录 一、问题描述 (2) 二、格式离散 (2) 二阶导数中心差格式离散 (2) 追赶法求解线性方程组简述 (3) 计算流程图 (5) 三、程序中主要符号和数组意义 (5) 四、计算结果与讨论 (6) 五、源程序 (9) 一、问题描述 利用中心差分格式近似导数22/dx y d ，数值求解 ()x dx y d 2sin 22= ()10≤≤x 1 /,0/10====x x y y 步长分别取 0001.0,001.0,01.0, 05.0=?x 二、格式离散 将x 轴上[0,1]之间的线段按上述步长，等步长的离散为n 个小段，包括端点，共n+1个网格节点，示意图如下： 线段上边的数字表示x 轴上的坐标值，线段下边的数字表示节点编号，从0到n 编号。 二阶导数中心差格式离散 211222)2sin(x y y y dx y d x i i i ?+-==+- 整理为线性方程形式 )2sin(2211x x y y y i i i ?=+-+- 其中，x ? 为空间离散步长；i=1,2，……，n-1 包括边界条件的线性方程组如下： 边界条件 边界条件0 ) *)1(*2sin(2......... ..........) **2sin(2..................) *1*2sin(20 21221122100=?-?=+-??=+-??=+-=--+-n n n n i i i y x n x y y y x i x y y y x x y y y y 改写成矩阵形式： f Ay = 其中，?????? ????????????????????----=1012112112112101 A ，??????????????????????=-n n i y y y y y y 110 ，??????????????????????=-n n i f f f f f f 110 系数矩阵A 中仅三对角线上的数值不全为0，其余位置上的数值全为0，是 典型的对角占优的三对角矩阵，列向量f 中，)2sin(2x i x f i ??=，且10==n f f ，作为边界条件。 追赶法求解线性方程组简述 ????? ?????????????????=??????????????????????????----=---n n n n n i i i b a c b a c b a c b a c b A 1111110 01012112112112101 两点边值问题有限元法（必做） 从Galerkin 原理出发用线性元解两点边值问题： "2,01(0)(1)0u u x x u u ?-+=< 有限差分法的Matlab程序（椭圆型方程） function FD_PDE(fun,gun,a,b,c,d) % 用有限差分法求解矩形域上的Poisson方程 tol=10^(-6); % 误差界 N=1000; % 最大迭代次数 n=20; % x轴方向的网格数 m=20; % y轴方向的网格数 h=(b-a)/n; % x轴方向的步长 l=(d-c)/m; % y轴方向的步长 for i=1:n-1 x(i)=a+i*h; end % 定义网格点坐标 for j=1:m-1 y(j)=c+j*l; end % 定义网格点坐标 u=zeros(n-1,m-1); %对u赋初值 % 下面定义几个参数 r=h^2/l^2; s=2*(1+r); k=1; % 应用Gauss-Seidel法求解差分方程 while k<=N % 对靠近上边界的网格点进行处理 % 对左上角的网格点进行处理 z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s; norm=abs(z-u(1,m-1)); u(1,m-1)=z; % 对靠近上边界的除第一点和最后点外网格点进行处理 for i=2:n-2 z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s; if abs(u(i,m-1)-z)>norm; norm=abs(u(i,m-1)-z); end u(i,m-1)=z; end % 对右上角的网格点进行处理 z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s; if abs(u(n-1,m-1)-z)>norm norm=abs(u(n-1,m-1)-z); end u(n-1,m-1)=z; % 对不靠近上下边界的网格点进行处理 for j=m-2:-1:2 % 对靠近左边界的网格点进行处理 椭圆方程差分法 1 矩形网上差分方程 考虑二阶椭圆型偏微分方程的第一边值问题 （1.1） ()()()?????=∈=+++--Γy x y x u y x F Eu Du Cu u u y x yy xx ,,,αG 其中C ,E D ,是常数；0≥E ；()()G C 0,∈=y x F F ；(,)x y α是给定的光滑函数。假设（5.1）存在光滑的唯一解。 为简单起见，假设G 是矩形区域，其四个边与相应坐标轴平行。考虑矩形网格：1h 和2h 分别为x 和y 方向的步长，h G 为网格内点节点集合，h Γ为网格边界点集合，=h G h G h Γ。 对于内点()j i y x ,h G ∈用如下的差分方程逼近（1.1） (1.2) 21 ,1,12h u u u j i ij j i -++---221,1,2h u u u j i ij j i -++-+1,1,12h u u C j i j i -+-+21,1,2h u u D j i j i -+-+ij Eu =ij F 其中),(j i ij y x F F =。（1.2）通常称为五点差分格式。 用（1.1）的真解(,)u x y 在网点上的值(,)i j u x y 、1(,)i j u x y -等等分别替换（1.2）中的ij u 、1,i j u -等等，然后在(,)i j x y 点处作Tailor 展开，便知（1.2）逼近（1.1） 的截断误差阶为() 2221h h O +。 方程（1.2）可以改写为 （1.3） j i a ,1-j i u ,1-+j i a ,1+j i u ,1++1,-j i a 1,-j i u +1,+j i a 1,+j i u +j i a ,j i u ,ij F = 对每一内点都可以列出这样一个方程。遇到边界点时，因为边界点u 的函数值已知，将相应的项挪到右端去。最后，得到一个以u 的内点近似值为未知数的线性方程组。这个方程组是稀疏的，并且当1h 和2h 足够小时是对角占优的。 可以证明，五点差分格式关于右端和初值都是稳定的，收敛阶为2212()O h h +。 用五点有限差分格式求解椭圆型方程（偏微分方程）程序3 2010-06-16 06:55 function varargout=liu(varargin) a=0;b=2;c=0;d=1;h1=1/32;h2=1/32; x=a:h1:b; y=c:h2:d; m=length(x); n=length(y); C=tri1(-1/h1^2,2*(1/h1^2+1/h2^2),-1/h1^2,n-2); D=-1/h1^2*eye(n-2); AA=tri2(D,C,D,m-2); BB=fc1t(a,b,c,d,h1,h2); XX=AA\BB; UU=zeros(m,n); for r=1:m for j=1:n UU(r,j)=fU(x(r),y(j)); end end UL=UU; for r=2:m-1 UL(r,2:n-1)=XX((r-2)*(n-2)+1:(r-1)*(n-2)); end [Y,X]=meshgrid(y,x); Z=abs(UL-UU); mesh(Y,X,Z); shading flat; xlabel('X','FontSize',14); ylabel('Y','FontSize',14); zlabel('error','FontSize',14); title('误差图') %求解下问题 %-(u_xx+u_yy)=f(x,y) x,y 在区域内x in a 第4章 椭圆型方程的有限差分法 §2 一维差分格式 1、用积分插值法导出逼近微分方程的差分格式。 d du du Lu=-(p )+r +qu=f,a 中心差分格式的程序实现 数学10-1班 余帆 10072121 1、考虑问题 考虑二阶常微分方程边值问题： f qu dx u d Lu =+-=22 （1） βα==)(,)(b u a u 其中q ，f 为[a,b]上的连续函数，βα,为常数。 2、网格剖分与差分格式 将区间[a,b]分成N 等分,分点为 N i ih a x i ,,1,0,???=+= , h=(b-a)/N,于是我们得到区间I=[a,b]的网格剖分，i x 为网格节点，h 为步长。 差分格式为： . ,,1,,2,1202 1 1βα==-???==++--=-+N i i i i i i i h u u N i f u q h u u u u L 3、截断误差 将方程（1）在节点离散化，由泰勒公式展开得 )()(12)()()(2)(344 2222 1 1h dx x u d h dx x u d h x u x u x u i i i i i O +??????+??????=+--+ 所以截断误差为 )()(12)(3 44 2h dx x u d h u R i i O +? ?????-= 4、数值例子 x x q e x u x sin 1)()(+== x e x f x sin )(= 其中[]1 ,0∈x 5、求解 由f qu dx u d Lu =+-=22， x e x f x sin )(= x x q e x u x sin 1)()(+== 将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为 ????????????? ?+--+--+-212 22 12112112h q h q h q N ????????????-121N u u u = ?????? ? ????? ??++-βα 12 2 212N f h f h f h 系数矩阵A=??? ? ? ? ? ??? ? ???+--+--+-212 22 12112112h q h q h q N ，我们可以利用高斯消去 法求得u （x ）的数值解。 6、实验结果 程序输出结果： 取N=10； 逼近解u1 真解u 1.10521961652189 1.10517091807565 1.22149147782632 1.22140275816017 二阶常微分方程的中心差分求解 学校:中国石油大学（华东）理学院 姓名：张道德 一、 实验目的 1、 构造二阶常微分边值问题: 22,(),(), d u Lu qu f a x b dx u a u b αβ?=-+=< 11122 222222333222122112 100121012010012 00N N N u f q h h u f q h h h u f q h h h q u f h h ---???? ??+-???? ??? ???? ???????-+-? ?????? ???????????=-+? ?????? ???????????-???? ????????-+????? ?? ????? 可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。 四、 举例求解 我们选取的二阶常微分方程边值问题为： 2 22242,01 (0)1,(1), x d u Lu x u e x dx u u e ?=-+=-< j. 椭圆方程数值解法 本章考虑椭圆微分方程数值解法。首先以二维二阶椭圆方程为例，给出矩形网和三角网上的差分法。然后以一维二阶椭圆方程为例，简要描述有限元法的基本思想。 J.1 矩形网上差分方程 考虑二维区域（区域=连通的开集）G 上的二阶椭圆型偏微分方程第一边值问题 （j.1） ()()() ,,,xx yy x y u u Cu Du Eu F x y u x y x y αΓ?--+++=∈?? =??G 其中C ,E D ,是常数；0≥E ；()()G C 0,∈=y x F F ；(,)x y α是给定的光滑函数； Γ是G 的边界；G =ΓG 。假设（J.1）存在光滑的唯一解。 考虑一种简单情形，即求解区域G 是矩形区域，并且其四个边与相应坐标轴平行。令1h 和2h 分别为x 和y 方向的步长，用平行于坐标轴的直线段分割区域 G ，构造矩形网格： h G 为网格内点节点集合，h Γ为网格边界节点集合，=h G h G h Γ。 对于内点()j i y x ,h G ∈，用如下的差分方程逼近微分方程（J.1）： (J.2) 1,1,,1,1 1,1,,1,1 2 2 212 1 2 2222i j ij i j i j ij i j i j i j i j i j ij ij u u u u u u u u u u C D Eu F h h h h +++-+-+--+-+--- - +++=其中),(j i ij y x F F =。（J.2）通常称为五点差分格式。 方程（J.2）可以整理改写为 （J.3） j i a ,1-j i u ,1-+j i a ,1+j i u ,1++1,-j i a 1,-j i u +1,+j i a 1,+j i u +j i a ,j i u ,ij F = 对每一内点()j i y x ,都可以列出这样一个方程。方程中遇到边界点时，注意到边界点上函数值u 已知，将相应的项挪到右端去。最后得到以u 的内点近似值为未知数的线性方程组。这个方程组是稀疏的，并且当1h 和2h 足够小时是对角占优的。 中心差分格式 1、考虑问题 考虑二阶常微分方程边值问题： f qu dx u d Lu =+-=22 （1） βα==)(,)(b u a u 其中q ，f 为[a,b]上的连续函数，βα,为常数。 2、网格剖分与差分格式 将区间[a,b]分成N 等分,分点为 N i ih a x i ,,1,0,???=+=, h=(b-a)/N,于是我们得到区间I=[a,b]的网格剖分，i x 为网格节点，h 为步 长。 差分格式为： .,,1,,2,120211βα==-???==++-- =-+N i i i i i i i h u u N i f u q h u u u u L 3、截断误差 将方程（1）在节点离散化，由泰勒公式展开得 )()(12)()()(2)(344222211h dx x u d h dx x u d h x u x u x u i i i i i O +??????+??????=+--+ 所以截断误差为 )()(12)(3442h dx x u d h u R i i O +??????-= 4、数值例子 x x q e x u x sin 1)()(+== 其中[]1,0∈x 5、求解 由f qu dx u d Lu =+-=22，且已知 x x q e x u x sin 1)()(+== 可得x e x f x sin )(= 将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为 ??????????????+--+--+-21222 1211 21 12h q h q h q N ????????????-121N u u u =??????????????++-βα122212N f h f h f h 系数矩阵A=?????? ????????+--+--+-212 221211211 2h q h q h q N ，我们求出矩阵A 极其逆便可求得u （x ）的数值解。 6、参考文献 《偏微分方程数值解法》李荣华 高等教育出版社 《科学计算中的有限差分法》 《MATLAB 程序设计教程》刘卫国 中国水利水电出版社 有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛(完整版)大连理工大学高等数值分析抛物型方程有限差分法
Poisson方程九点差分格式_米瑞琪
Matlab PDE工具箱有限元法求解偏微分方程
用五点有限差分格式求解椭圆型方程(偏微分方程) 程序2
对流扩散方程有限差分方法.
课程名称(中文)偏微分方程数值解专题
利用中心差分格式数值求解导数
椭圆型方程的有限元法
有限差分法的Matlab程序(椭圆型方程)
大连理工大学 高等数值分析 椭圆方程差分法
用五点有限差分格式求解椭圆型方程(偏微分方程) 程序3
椭圆型方程的有限差分法
一阶中心差分格式
偏微分中心差分格式实验报告(含matlab程序)
椭圆方程数值解
中心差分格式
有限差分法
相关主题
文本预览