向前差分格式

偏微分方程数值解法作业
题目：抛物型方程最简差分格式
学院：理学院
专业年级：2008级信息与计算科学二班*名：**
学号：********
指导老师：***
2011年6月
一 考虑问题及条件
考虑一维热传导方程：
22(),0u u
a f x t T t x
∂∂=+<≤∂∂…………（1.1） 其中a 是正常数，()f x 是给定的连续函数。

现在考虑第二类初边值问题的差分逼近：
初始条件：(,0)(),0u x x x l ϕ=<<…………（1.2） 边值条件：(0,)()u t t η=，(,)()u l t t γ=，0t T ≤≤………（1.3）
假设()f x 和()x ϕ在相应区域光滑，并且在0,x l =满足相容条件，使上述问题有惟一充分光滑的解。

用向前差分格式计算如下热传导方程的初边值问题
⎪⎪⎩
⎪⎪⎨⎧>==<<-=>=∂∂-∂∂,0,0),1(,1),0(,
10,1)0,(u )0(,022t t u t u x x x a x u
a t
u ，
是常数
已知其精确解为
u(x,t)=1-x.
二 网格剖分与差分格式
（1）.区格剖分
取空间步长l h N =和时间步长T M τ=，其中,N M 都是正整数。

用两族平行直线
(0,1,
,)
j x x jh j N ===和
(0,1,
,)k t t k k M τ===将矩形域
{}0;0G x l t T =≤≤≤≤分割成矩形网格，网格节点为(,)j k x t 。

以h G 表示网格内点
集合，即位于开矩形的网点集合；h G 表示所有位于闭矩形G 的网点集合；
h h h G G -=Γ是网格界点集合。

其次，用k j u 表示定义在网点(,)j k x t 的函数，M k N j ≤≤≤≤0,0
（2）.微分方程的离散，建立相应差分格式
将方程在节点(,)j k x t 离散化，
22()k
k
j j j
u u a f x t x ⎡⎤∂∂⎡⎤
=+⎢⎥⎢⎥∂∂⎣⎦⎣⎦，1,2,,1j N =- 1,2,,1k M =-…………（1.4）
对充分光滑的解u ，由Taylor 展式：
22312
(,)(,)(,)(,)()2j k j k j k j k u x t u x t u x t u x t O t t τττ+∂∂=+++∂∂…………（1.5） 234
23451234
(,)(,)(,)(,)(,)(,)()23!4!j k j k j k j k j k j k u x t u x t u x t u x t h h h u x t u x t h
O h x
x x x
+∂∂∂∂=+++++∂∂∂∂…………（1.6）
234
23451234
(,)(,)(,)(,)(,)(,)()23!4!j k j k j k j k j k j k u x t u x t u x t u x t h h h u x t u x t h
O h x
x x x
-∂∂∂∂=-+-++∂∂∂∂…………（1.7）
（1.5）移项得：
212
2
(,)(,)(,)(,)()2j k j k j k j k u x t u x t u x t u x t O t
t
τττ+∂-∂=-+∂∂…………（1.8） （1.6）（1.7）相加得：
24
21132
2
4(,)(,)2(,)(,)
(,)()12j k j k j k j k j k u x t x t u x t u x t u x t h O h x h x
+-∂-+∂=
-+∂∂…………（1.9） 将（1.8）（1.9）代入（1.4）得：
1112
(,)(,)
(,)2(,)(,)
()()j k j k j k j k j k k j j u x t u x t x t u x t u x t a
f x R u h
τ
++---+=++…………
（1.10）
其中，
24
2232
4
(,)(,)
()()()2
12
j k j k k j
u x t u x t h R u O O h t x ττ∂∂=
-
++∂∂
舍去()k j R u ，得到逼近（1.1）的向前格式差分方程：
111
2
2k k k k k
j j
j j j j u u u u u a
f h
τ
++---+=+，1,2,
,1j N =- 1,2,,1k M =-……（1.11）
其中，(,)k j j k u u x t =，()j j f f x =
其中j=1,2,…,N-1, k=0,1,2,…,M-1.以2/h a r τ=表示网比,则可将（1.11）
式化为：
.)21(111j k j k j k j k j f ru u r ru u τ++-+=-++
矩阵形式：
=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+-++111211k N k k u u u ⎥
⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎣
⎡---r r r r r r r 212121 +⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-k N k k u u u 121 ⎥⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎢⎣⎡++-k N N k
ru f f ru f 1201τττ .
三、截断误差
将（1.8）（1.9）代入（1.4）得：
1112
(,)(,)
(,)2(,)(,)
()()j k j k j k j k j k k j j u x t u x t x t u x t u x t a
f x R u h
τ++---+=++…………
（1.10）
其中，24
2232
4
(,)(,)
()()()2
12
j k j k k
j
u x t u x t h R u O O h t x ττ∂∂=
-
++∂∂
而得到的向前逼近差分方程是舍去()k j R u 的结果，固()k
j R u 即为向前差分方程
的局部截断误差。

四、稳定性
如果在应用差分格式（1.1）时，计算右端函数)(x f 有误差)(x g ，计算初值)
(x ϕ时有误差)(x ψ，令),(),(),(t x u t x v t x -=ε，则计算可得)(22x g x
a t =∂∂-∂∂ε
ε，)()0,(x x ψε=，由定理可知，当网比21≤r 时，有∑-=∞∞∞+≤1
||||||||||||k t l k g τψε.
上式说明当∞||||ψ和∑-=∞1
||||k t l g τ很小时，误差∞||||k ε也很小，即当网比21≤r 时，
差分格式是稳定的。

五、数值例子
实现结果：
六、参考文献
[1] 李荣华偏微分方程数值解法高等教育出版社2005。

[2]孙志忠《偏微分方程数值解法》科学出版社，2005。

[3] 汤怀民、胡建伟，《微分方程数值解法》，南开大学出版社，1990。

[4] 李庆扬《数值分析》清华大学出版社，2008 。

[5] 余德浩，汤华中编著：《微分方程数值解法》（第三版）科学出版社2004.6。

七、附件
1.流程图
程序实现流程图
2、详细设计
1.建立函数文件及执行文件
文件1：fun_ux0.m 用于处理初值条件；
文件2：fun_u0t.m 用于处理边值条件；
文件3：fun_u1t.m 用于处理边值条件；
文件4：fun_u.m 计算真值；
文件5：fun_f.m 计算f(x)
文件 6：fun.m
2.程序实现代码
2.1 输入
输入时间上限T，时间方向的节点数M；输入空间上限l，空间方向的节点数N；输入大于0的常数a.T,l大于0；M，N为正整数.
2.2 设置初值
h=l/N;v=T/M;
节点： x=0+h:h:1;t=v:v:T;
网比：r=a*v/(h*h);
u0=ones(1,N-1);% u0为初值,f是一个1*（N-1）的矩阵
for k=1:N-1
u0(k)=fun_ux0(x(k))*u0(k);
f(k)=v*fun_f(x(k));
end
f(1)=f(1)+r*fun_u0t(0);%边值条件
f(N-1)=f(N-1)+r*fun_u1t(0);
A=diag(r*ones(1,N-2),-1)+diag((1-2*r)*ones(1,N-1))+diag(r*ones(1,N-2),
1);
2.3 计算节点的值
u1=A*u0'+f';
u0=u1
u(k,:)=u1;
f(1)=f(1)-r*fun_u0t(0);
f(N-1)=f(N-1)-r*fun_u1t(0);
for k=2:M
f(1)=f(1)+r*fun_u0t(t(k-1));
f(N-1)=f(N-1)+r*fun_u1t(k-1);
u1=A*u0+f';
u0=u1;
u(k,:)=u1;
f(1)=f(1)-r*fun_u0t(t(k-1));
f(N-1)=f(N-1)-r*fun_u1t(k-1);
end。