中心差分格式

中心差分格式
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
121121
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 程序设计教程》刘卫国 中国水利水电出版社
附件
1、程序流程图
开始
计算各初
始条件
构造系数矩阵d
计算矩阵d的逆D
构造格式
右端矩阵B
计算数值解
u=D*B
输出u的数值解
结束
2 、程序代码
%zah.m
%u(x)=exp(x), q(x)=1+sinx, f(x)=exp(x)*sinx 0<x<1
a=input('input a:');%输入区间条件a
b=input('input b:');%输入区间条件b
N=input('input N:');%输入等分数N
h=(b-a)/N;%步长
for i=1:N-1
x(i)=a+i*h;%网格节点
q(i)=1+sin(x(i));%连续函数q网格节点上各值
f(i)=exp(x(i))*sin(x(i));%连续函数f网格节点上各值
end
c1=linspace(2+q(1)*h*h,2+q(N-1)*h*h,N-1); %系数矩阵主对角线上各值c2=linspace(-1,-1,N-2); %系数矩阵次对角线上各值
d1=diag(c1,0);
d2=diag(c2,-1);
d3=diag(c2,1);
d=d1+d2+d3 %构造系数矩阵
D=inv(d); %系数矩阵的逆
b1=linspace(h*h*f(1),h*h*f(N-1),N-1);
for i=1:N-1
if i==1
b2(i)=exp(a);
else
if i==N-1
b2(i)=exp(b);
else
b2(i)=0;
end
end
end
B1=b1+b2;%构造差分格式等式右边的向量
B=B1' %差分格式等式右边向量的转置
u=D*B %输出u（x）数值解
输出结果：
>> clear
>> zah
input a:0
input b:1
input N:10
d =
2.0110 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0119 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0127 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0136 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0144 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0153 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0161 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0170 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0178
B =
1.0011
0.0034
0.0056
0.0079
0.0102
0.0125
0.0147
0.0170
2.7375
u =
1.1117
1.2345
1.3685
1.5143
1.6727
1.8451
2.0331
2.2393
2.4664。