二阶常微分方程的中心差分求解
学校:中国石油大学（华东）理学院 姓名：张道德
一、 实验目的
1、 构造二阶常微分边值问题:
22,(),(),d u
Lu qu f a x b
dx u a u b αβ⎧=-+=<<⎪⎨⎪==⎩
其中,q f 为[,]a b 上的连续函数，0;,q αβ≥为给定常数的中心差分格式；
2、 根据中心差分格式求解出特定例题的数值解，并与该
例题的精确解进行比较。
二、 中心差分格式的构造
将区间[,]a b 分成N 等分，分点为： 0,1,2,
,i x a ih i N =+=
()/h b a N =-。
于是我们得到区间的一个网络剖分。
称为网
格的节点称为步长。
得中心差分格式为：
112
02,1,2,,1,,.
i i i h i i i i N u u u L u q u f i N h u u αβ+--+⎧=-+==-⎪⎨⎪==⎩
其中式中(),()i i i i q q x f f x ==。
三、 差分格式求解
根据中心差分格式可以构造出：
11122222222
33322
212
21121001210
12
010012
00
N 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 ⎧=-+=-<<⎪⎨⎪==⎩
其精确解为：2
x u e
=。
则我们具体求解出的解如下：
1、 当10N =时，数值解与精确解为： (1) 表1、10N =时，数值解与精确解统计表
x 的值 0.1
0.2
0.3
0.4
0.5 u 的数值解 1.011069 1.042744 1.096904 1.176896 1.28789 u 的精确解 1.01005 1.040811 1.094174 1.173511 1.284025 两者之差 0.001019 0.001934 0.002729 0.003385 0.003864 x 的值 0.6
0.7
0.8
0.9
u 的数值解 1.437443 1.636363 1.900001 2.250209 u 的精确解 1.433329 1.632316 1.896481 2.247908 两者之差
0.004114
0.004046
0.00352
0.002301
将两者绘于同一图中如下：
(2)结论：可以看出数值解与精确解之间的误差很小， 在 210-这样一个数量级上。
我们也可以求出的
|-|数值解精确解范数 来，如下：
Norm1(|-|数值解精确解)=0.0269; Norm2(|-|数值解精确解)=0.0095; Normoo(|-|数值解精确解)=0.0041;
所以，可以得出中心差分格式求解该方程效果挺好。
2、当20N =时，将数值解与精确解绘于同一图像中，如下：
0.1
0.20.30.4
0.50.60.70.80.9
x 的值
u 的数值解与精确解
图1、N=10时，数值解与精确解图像
(2)结论：可以看出数值解与精确解之间的误差很小，
在 4
10-这样一个数量级上。
我们也可以求出的
|-|数值解精确解范数 来，如下：
Norm1(|-|数值解精确解)=0.0027; Norm2(|-|数值解精确解)=4
3.018810-⨯; Normoo(|-|数值解精确解)=5
4.166910-⨯;
所以，可以得出中心差分格式求解该方程效果挺好。
五、 程序 程序1
11.21.41.61.822.22.4
2.62.8x 的值
u 的数值解与精确解
图2、当N=100时，u 的数值解与精确解
%**************************************************************
%f221.m
function [q,f]=f211(x)
%q函数，f函数
q=4*x^2;
f=-2*exp(x^2);
%***************************************************************
程序2
%******************************************************************** %追赶法
function [x]=zhuiganfa(a,b,c,d)
%对角线下方的元素，个数比A少一个
% %对角线元素
%对角线上方的元素，个数比A少一个
%d为方程常数项
%用追赶法解三对角矩阵方程
r=size(a);
m=r(2);
r=size(b);
n=r(2);
if size(a)~=size(c)|m~=n-1|size(b)~=size(d)
error('变量不匹配,检查变量输入情况!');
end
%%
%LU分解
u(1)=b(1);
for i=2:n
l(i-1)=a(i-1)/u(i-1);
u(i)=b(i)-l(i-1)*c(i-1);
v(i-1)=(b(i)-u(i))/l(i-1);
end
%追赶法实现
%%
%求解Ly=d，"追"的过程
y(1)=d(1);
for i=2:n
y(i)=d(i)-l(i-1)*y(i-1);
end
%%
%求解Ux=y,"赶"的过程
x(n)=y(n)/u(n);
for i=n-1:-1:1
x(i)=y(i)/u(i);
x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
%********************************************************************
程序3
%******************************************************************** *
%ODE2.m
function [x]=ODE2(x0,xN,u0,uN,N)
%中心差分求解
%x0,xN初始条件
%u0,uN边值条件
%N等分数
%步长
h=1/N;
%%
a(1:N-2)=-1/h^2; %对角线下方的元素，个数比A少一个
for i=1:N-1
z(i)=x0+i*h;
[q,f]=f211(z(i));
b(i)=2/h^2+q; %对角线元素
d(i)=f;
end
%对对角线元素进行调整
d(1)=d(1)+u0/h^2;
d(N-1)=d(N-1)+uN/h^2;
%%
[x]=zhuiganfa(a,b,a,d);%数值解
%******************************************************************** **
程序4
%**************************************************************
%main_chapter.m
%x0,xN初始条件
%u0,uN边值条件
%N等分数
%h步长
x0=0;xN=1;u0=1;
uN=exp(1);N=10;h=1/N;
[x]=ODE2(x0,xN,u0,uN,N);%数值解
z=x0+h:h:xN-h;
y=exp(z.^2); %真解
nr1(x-y)
plot(z,x,'og',z,y,'.r')
%***************************************************************。