偏微分中心差分格式实验报告(含matlab程序)

二阶常微分方程的中心差分求解

学校:中国石油大学（华东）理学院 姓名：张道德

一、 实验目的

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);