中心差分解两点边值问题.docx
《中心差分解两点边值问题.docx》由会员分享，可在线阅读，更多相关《中心差分解两点边值问题.docx（19页珍藏版）》请在冰豆网上搜索。
中心差分解两点边值问题
一•题目
用中心差分格式计算如下两点边值问题
dxdux
esinx1xuxe2x1sinx1xxx1,1x2
dxdx
u
（1）0,u
（2）2
已知其精确解为
u（x）x（x1）
二.理论
作为模型，考虑两点边值问题：
Lu-d（p（x）du）晋q（x）uf（x）,axb（1.1）
dxdxdx
u（a）,u（b）（1.2）
1
是给定的常数。
假定pC[a,b],p（x）Pmin0,r,q,fC[a,b],
1.建立差分格式
（1）•区域网格剖分首先取N1个节点:
将区间I［a,b］分成N个小区间:
于是得到区间I的一个网格剖分。
记
Ii:
xi1xxi,i1,2,LN.
xixi1，称hmaxh为网格最大步长。
用Ih表
i
示网格内点x1,x2,L,xN』勺集合，
Ih表示内点和界点沧a,xNb的集合。
取相邻节点xi1,xi的中点x1
i一2
1.
（Xi1xj（i1,2,L,N），称为半整数点。
则由节点
2
xo
x-ix?
L
22
Xnb
又构成［a,b］的一个网格剖分，称为对偶剖分。
用差商代替微商，将方程（
（2）•微分方程的离散，建立相应差分格式
1.1）在内点Xi离散化•注意对充分光滑的u，由Taylor展式有
u（XiJu（XiJ
hihi1
hi1hi
2
[霁i
dx
O（h2）
（1.3）
P©;）
U（Xi）.U（Xi1）[P乳1
dxi7
1
2
iO（h3）
[p%.dxi
.2,3
hirdU3.
:
[p3]j1O（h）
24dx|?
（1.4）
p（X1）U（Xi1）U（Xi）
i—
2
hi1
[p^l1
dxi2
23
h1rdu】-rylPaliO（h）
24dx
（1.5）
hih1得
—2—，寸
[p（X1严J呛）
2
2dun
b（[PdX]i
由（1.5）减（1.4）,并除以
hi1
呂（P^lih
dxdx
rdu_[P]dxi
hi「d2
1）
2
p（xi1）U（Xi）hU（Xi1）]
i2hi
3[P吗O（h2）
12dx3i
[p岁]iO（h2）dx
duXh[硬（P『]i
1.6）
令P.1
i
2
满足方程:
p（xi
1）,rir（xj,qi
2
q（xj,fi
f（X）,则由
（1.3）（1.6）知，边值问题的解
u（x）
LhU（xJ
u（X
1）u（xj
hi1
u（x）
1
2
u（xi1）]
hi
1.7）
其中
亡曲1）
u（Xi1）]
qiU（Xi）
R（u）
R（u）（hi1
1
hi）（；[dx--dx
1d3u
扣％i）O（h2）
2dx
（1.8）
为差分算子Lh的截断误差，舍去Ri（u），便得逼近边值问题（1.1）（1.2）的差分方程:
LhU（Xi）土[p1U（X1）U（Xi）
hihi1i2
hi1
U（Xi）U（Xi1）
1.]
2
hi
（1.9）
[u（Xi1）u（Xi1）]qiU（Xi）
hih1
i=1,2,…,N-1,
Uo,Un
由方程（1.7）（1.9）,截断误差R（u）可表示为
（1.10）
当网格均匀，即
R（u）LhU（xJLhUiLh（u（xJuj
h（i1,2,L,N）时差分方程（1.9）简化为
（1.11）
（1.12）
x
e2x1sinx1xxx1
所以截断误差按||Rh（u）||o或||Rh（u）|C的阶为O（h）。
在本题中，p（x）ex,r（x）0,q（x）sinx1x,f（x）
因为r=0方程（1.11）的系数对角矩阵是三对角矩阵。
我们可以用消元法或迭代法求解方程
组（1.1）（1.2）
式（1.11）用方程组展开：
P1）弭
2
1f
2P1U0f1
h2
k
u
1-2
k
p
1-2
pk
/V
k
u
1-2-
pk
UN
3-2
N
q
PN
11
P3U[2（P3
h2h2
占（5Pi）q
h22
1
1h2P3
0
LL
0
0
1
h2P3
2（p5p3）Q
h22
1
严
LL
0
0
L
L
L
0
LL
0
1
h2PN|
右（Pn
3PN5）qN2
22
1
h2PN|
0
LL
0
0
1
RPn；
*Pn
1PN3）Ch
■2~2
写成矩阵形式为:
u
U2
M
UN2
UN1
fl
2•收敛性分析
根据（1.10）我们引进误差
eu（x）Ui
则误差函数eh（x）e满足下列差分方程:
Lei驾i1,2,L,N
eoeN0
的问题。
由（1.12）我们知，有
lhimllR（u）ll0
从而差分方程满足相容条件。
若引进记号
设PminC0
则可将（1.9）改写为
LhUi（p.i（uj）qqfi
I—XX
2
将差分解U|表成
U|UIUl,iIh（2.1）
其中UI满足
（ui）0,lIh,U0Uo,UNUn（2.2）
XX
而Ui满足
LhUiflLhUi,iIh,Uoun0（2.3）
先估计Uh，由
Coll（Uh）xl2（fhLhUh,Uh）|h（fh,Uh）|h（LhUh,Uh）|h（2.4）
据差分格林公式
（LhUh,Uh）ih（Ph（Uh）,（Uh））ih（qhUh,Uh）ih
XX
再利用柯西不等式，有常数c1使
|（LhUh,Uh）|h|[Cill（Uh）lloC2||Uh|0]||（Uh）ll（2.5）
XX
将不等式（2.6）用于（2.5）右端，则
~1c1C2
ll（Uh）||0-||fhll1-ll（Uh）||0-||Uh||0（2.6）
XQQXQ
解差分方程（2.2，易得）
从而
UnUi
b
U0（
（Xi
a
a）U0
1
ll（Uh）||0
i
|UnU0l
X
.（ba）
||uh||°,（ba）max|ui|、（ba）（|u°||Un|）
IIh
这样,
l|uh|1..._（ba）（ba）1（|Uo||Un|）
（2.7）
1
ll（Uh）||o-||fh||1
xCo
C1色|3|1
Co
（2.8）
因为
~
ba.
l|Uh||o
_||（uh）||o
因此
2x
利用范数||uh仏，从（2.7）推出
l|Uh|1||Uh||。
||（Uh）||o
X
〜
（1¥）ll（Uh）||0
2x
ba1c1c2
（1）[llfhlli-2||uh||i]
2CoCo
（2.9）
联结（2.1）（2.7）及（2.9）即得差分解的先验估计:
||Uh|1||uh|1||uh|1
M1||fh||1M2（|u°||Un|）
其中
（2.10）
1“ba、
M1-
（1）,
co2
M2；（ba）（ba）1[1（1呼声绍
2Co
不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1.11）关于右端及边值
稳定.
根据定理1.1:
若边值问题的解u充分光滑，差分方程按|g|R满足相容条件且关于右端
稳定，则差分解uh按|g|收敛到边值问题的解，且有和||Rh（u）|R相同的收敛阶。
所以差分
方程的解的收敛速度为O（h2）。
三．程序代码：
clc
clf
clf
symsx;
a=1;%区间界点
b=2;%区间界点
p=exp（x）;%这是p函数
q=sin（x）+1+x;%这是q函数
f=-exp（x）*（2*x+1）+（sin（x）+1+x）*x*（x-1）;%这是f函数r=0;%这是r函数.
N=10;%将区间划分的等分,这里控制！
h=（b-a）/N;%这里确定步长
value_of_f=zeros（N-1,1）;%这是fdiag_0=zeros（N-1,1）;%确定A的对角元diag_1=zeros（N-2,1）;%确定A的偏离对角的上对角元diag_2=zeros（N-2,1）;%确定A的偏离对角的下对角元X=a:
h:
b;
u_a=0;%边界条件
u_b=2;%边界条件
forj=2:
N
diag_0（j-1）=（（subs（p,{x},{（X（j+1）+X（j））/2}））+（subs（p,{x},{（X（j-1）+X（j））/2}）））/（hA2）+（subs（q,{x},{X（j）}））;
end%获取对角元素
forj=3:
Ndiag_2（j-2）=-（（subs（p,{x},{（X（j-1）+X（j））/2}）））/（hA2）-subs（r,{x},{X（j）}）/（2*h）;
end%获取A的第三条对角
forj=2:
N-1diag_1（j-1）=-（（subs（p,{x},{（X（j+1）+X（j））/2}）））/（hA2）+subs（r,{x},{X（j）}）/（2*h）;
end%获取A的第二条对角
forj=2:
N;
value_of_f（j-1）=subs（f,{x},{X（j）}）;
end%获取F值
value_of_f
（1）=value_of_f
（1）+u_a*（subs（p,{x},{（X
（2）+X
（1））/2}））/（hA2）;
value_oLf（N-1）=value_of_f（N-1）+u_b*（subs（p,{x},{（X（N）+X（N+1））/2}））/（hA2）;A=diag（diag_0）+diag（diag_1,1）+diag（diag_2,-1）;%组装系数矩阵formatlong
U=inv（A）*value_of_f%差分解
%fprintf（'%11.5f',U）
fprintf（'\n'）;
dx=X（2:
N）;
precise_value=dx.*（dx-1）%精确解
%fprintf（'%11.5f',precise_value）
deta=U-precise_value';%误差
deta_max=max（abs（deta））;%最大误差
fprintf（'最大的误差是%f\n',deta_max）plot（X（2:
N）,U,'b*',X（2:
N）,precise_value,'r--'）%差分解与精确解对比表
figure。
；
plot（X（2:
N）,deta）%误差图
结果：
的值
步长卜\、
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
最大误差
0.1
0.11015
0.24026
0.39033
0.56037
0.75038
0.96038
1.19031
1.44022
1.71012
0.000380
0.05
0.11003
0.24006
0.39008
0.56009
0.75009
0.96008
1.19007
1.44005
1.71003
0.000095
0.025
0.11000
0.24001
0.39002
0.56002
0.75002
0.96002
1.19002
1.44001
1.71001
0.000024
0.0125
0.11000
0.24000
0.39000
0.56001
0.75001
0.96000
1.19000
1.44000
1.71000
0.000006
精确解
0.11000
0.24000
0.39000
0.56000
0.75000
0.96000
1.19000
1.44000
1.71000
0
以下仅给出步长为N=20,h=0.05的精确值和差分值图与误差图，其他的只要修改程序中的步长（N）值，即可同理得出。
N=20,h=0.05（图中蓝点和紫线分另忧精确值图和差分图）
2
1.3
1.6
14
1.2
1
D8
D.6
0.4
D2
□
111121314151.61.71.61.S2
四•结论：
由本题可以总结出，通过用中心差分解两点边值所得的数值解能够较好地逼近方程的精
确解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确
解。
PNJ1]UN1
升级会员 