中心差分解两点边值问题.doc
《中心差分解两点边值问题.doc》由会员分享，可在线阅读，更多相关《中心差分解两点边值问题.doc（9页珍藏版）》请在冰点文库上搜索。
一．题目
用中心差分格式计算如下两点边值问题
已知其精确解为
二．理论
作为模型，考虑两点边值问题：
…………（1.1）
…………（1.2）
假定是给定的常数。
1.建立差分格式
（1）.区域网格剖分
首先取个节点：
将区间分成个小区间：
于是得到区间的一个网格剖分。
记，称为网格最大步长。
用表示网格内点，，，的集合，表示内点和界点的集合。
取相邻节点的中点，称为半整数点。
则由节点
又构成的一个网格剖分，称为对偶剖分。
（2）.微分方程的离散，建立相应差分格式
用差商代替微商，将方程（1.1）在内点离散化.注意对充分光滑的，由Taylor展式有
………（1.3）
………（1.5）
由（1.5）减（1.4），并除以，得
…………（1.6）
令则由（1.3）（1.6）知，边值问题的解满足方程：
…………（1.7）
其中
…………（1.8）
为差分算子的截断误差，舍去，便得逼近边值问题（1.1）（1.2）的差分方程：
…………（1.9）
i=1,2,…,N-1，
由方程（1.7）（1.9），截断误差可表示为
…………（1.10）
当网格均匀，即时差分方程（1.9）简化为
…………（1.11）
这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。
这个方程就是中心差分格式。
截断误差为：
…………（1.12）
所以截断误差按或的阶为。
在本题中，，，，
，
因为r=0方程（1.11）的系数对角矩阵是三对角矩阵。
我们可以用消元法或迭代法求解方程组（1.1）（1.2）
式（1.11）用方程组展开：
写成矩阵形式为：
2.收敛性分析
根据（1.10）我们引进误差
则误差函数满足下列差分方程：
于是收敛性及收敛速度的估计问题，就归结到通过右端（截断误差）估计误差函数的问题。
由（1.12）我们知，有
从而差分方程满足相容条件。
若引进记号
，，，
，，
设
则可将（1.9）改写为
将差分解表成
…………（2.1）
其中满足
…………（2.2）
而满足
…………（2.3）
先估计，由
…………（2.4）
据差分格林公式
再利用柯西不等式，有常数使
…………（2.5）
将不等式（2.6）用于（2.5）右端，则
…………（2.6）
解差分方程（2.2，易得）
从而
这样，
…………（2.7）
利用范数，从（2.7）推出
…………（2.8）
因为
因此
…………（2.9）
联结（2.1）（2.7）及（2.9）即得差分解的先验估计：
…………（2.10）
其中
不等式（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1.11）关于右端及边值稳定.
根据定理1.1：
若边值问题的解u充分光滑，差分方程按满足相容条件且关于右端稳定，则差分解按收敛到边值问题的解，且有和相同的收敛阶。
所以差分方程的解的收敛速度为。
三．程序代码：
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）;%这是f
diag_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}）））/（h^2）+（subs（q,{x},{X（j）}））;
end%获取对角元素
forj=3:
N
diag_2（j-2）=-（（subs（p,{x},{（X（j-1）+X（j））/2}）））/（h^2）-subs（r,{x},{X（j）}）/（2*h）;
end%获取A的第三条对角
forj=2:
N-1
diag_1（j-1）=-（（subs（p,{x},{（X（j+1）+X（j））/2}）））/（h^2）+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}））/（h^2）;
value_of_f（N-1）=value_of_f（N-1）+u_b*（subs（p,{x},{（X（N）+X（N+1））/2}））/（h^2）;
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）%误差图
结果：
X的值
步长h
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（图中蓝点和紫线分别是精确值图和差分图）
四．结论：
由本题可以总结出，通过用中心差分解两点边值所得的数值解能够较好地逼近方程的精确解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确解。
升级会员 