中科大计算流体力学CFD之大作业一

CFD 实验报告一
姓名： 学号：
一、题目：
利用中心差分格式近似导数22/dx y d ，数值求解常微分方程
x dx y
d 2sin 2
2= （10≤≤x ） 00==x y 4
2
sin 11-
==x y 步长分别取x ∆＝0.05, 0.01, 0.001，0.0001。

二、报告要求：
1）列出全部计算公式和步骤；
2）表列出程序中各主要符号和数组意义；
3）绘出数值计算结果的函数曲线，并与精确解比较；
4）比较不同差分格式和不同网格步长计算结果的精度和代价； 5）附源程序。

三、相关差分格式
二阶导数22/dx y d 的三点差分格式有向前差分、向后差分和中心差分，表达
式分别如下：
()()()22122
22122
211
222
222j j j
j j j
j j j u u u u O x x x
u u u u O x x x u u u u O x x x
++--+--+∂=+∆∂∆-+∂=+∆∂∆-+∂=+∆∂∆一阶向前差分：一阶向后差分：二阶中心差分：
代入微分方程可以得到差分方程，表达式分别如下：
212
212
11
22=sin 22=sin 22=sin 2j j j j j j j j j j j j
u u u x x
u u u x x u u u x x
++--+--+∆-+∆-+∆一阶向前差分：一阶向后差分：二阶中心差分：
对于三种差分格式，差分格式可以改写成AY b =的形式，其中A 是相同的，
非齐次项b 不同，如下所示：
2112112112A -⎡⎤⎢⎥-⎢⎥
⎢
⎥=⎢⎥
-⎢⎥
⎢⎥-⎣⎦
系数矩阵 ()()()02112
3221sin 2sin 2sin 2k k y x x b x x x x y ---⎡⎤
⎢⎥
∆⎢⎥
⎢
⎥=⎢⎥∆⎢⎥
⎢⎥
∆-⎣⎦
一阶向前差分 ()()()()2202
322121sin 2sin 2sin 2sin 2k k x x y x x b x x x x y -⎡⎤∆-⎢⎥∆⎢⎥
⎢
⎥=⎢⎥∆⎢⎥
⎢⎥∆-⎣⎦一阶向后差分 ()()()()2102
232
2211sin 2sin 2sin 2sin 2k k x x y x x b x x x x y --⎡⎤
∆-⎢⎥∆⎢⎥
⎢
⎥=⎢⎥∆⎢⎥⎢⎥
∆-⎣
⎦二阶中心差分 求解AY b =可以得到各节点y 的值[]T
1
22
1k k Y y y y y --=。

四、计算公式和步骤；
1.关于精确解的推导：
已知
22
sin 2d y
x dx =，对x 进行两次积分，得到121
sin 24
y x C x C =-++，再结合
边界条件00
==x y 和4
2
sin 11-==x y 得到相对应的1C 和2C ，确定最后精确解为：
1
sin 24y x x =-+。

2.关于数值求解方法：
对于方程组AY b =可直接求解，也可以使用追赶法求解，下面介绍简单追赶法求解三对角方程组的过程。

对于三对角方程组：
11112222211111=n n n n n n n n n y b d c y b a d c y b a d c a d y b -----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
系数矩阵A 中仅三对角线上的数值不全为0，其余位置上的数值全为0，是典型的对角占优的三对角矩阵，利用高斯消去法，经过n-1次消元可以化为同解方程组：
11122211111=11n n n n n y q u y q u y q u y q ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
如上所示，求这些值的过程（即消元过程）称为追：
()
()
()()()111111111/,/,
/2,...,1/2,...,i i i i i i i i i i i i u c d q b d u c d u a i n q b q a d u a i n ---===-=-=--=
再利用回代过程求出方程组的各变量：
()1
,1,2,...,1n n i i i i y q y q u y i n n +==-=--
这一逆序求变量的过程（回代过程）称为赶。

五、计算结果与分析
根据题意需要分别取步长x ∆＝0.05, 0.01, 0.001，0.0001进行计算，因此采用MATLAB 进行编程运算，然后将数值解与精确解进行比较，如下图所示：
∆＝0.05 步长x
∆＝0.01步长x
∆＝0.001
步长x
∆＝0.0001
步长x
从上面四个图可以看出，这几个格式的精确解和数值解之间符合度很好，其中中心差分格式精确度更高，并且随着步长的减小与精确解的符合度越高。

下面将计算结果用表格的形式列出：
）
表1 不同格式不同步长计算耗时（单位s
由不同格式不同步长计算耗时和计算精度的结果对比可知：
1、从时间代价来看，三种格式采用追赶法求解所需要的时间基本一样，从计算精度来看，中心差分格式的计算精度明显优于向前差分和向后差分格式。

2、随着步长的减小，虽然求解的计算精度增高，但所需要的时间却大大增加。

3、对于这三种格式而言，中心差分格式明显优于另两种格式。

4、从上可知，提高计算精度可以通过减少步长和采用高阶格式的方法，由于减
少步长需要用时间作为代价，所以采用高阶格式相对而言比较好。

六、源程序及其说明
符号说明：
源程序如下：
clear; %清空
long=1;%变量取值范围
dx=0.05; %设置步长
n=long/dx; % 此处n=节点数—1
y0=0;%设定始端边界条件
yl=1-sin(2)/4; %设定末端边界条件
A=zeros(n-1,n-1);%设置系数矩阵A
A(1,1:2)=[-2 1];%始端边界
for i=2:n-2
A(i,i-1:i+1)=[1,-2,1];
end
A(n-1,n-2:n-1)=[1,-2];%末端边界
b=zeros(n-1,1);%设置矩阵b
for i=1:n-1
% b(i)=dx^2*sin((i-1)*2*dx); %一阶向前差分
% b(i)=dx^2*sin((i+1)*2*dx); %一阶向后差分
b(i)=dx^2*sin(i*2*dx); %中心差分
end
b(1)=b(1)-y0;%设置边界条件
b(n-1)=b(n-1)-yl;
x=dx*(0:n);%求解精确解
Y_e=x-sin(2*x)/4;
disp('~~~~~~~~~~~~~~~~~~~SC16005041李文建~~~~~~~~~~~~~~~~~~~~~~~~'); disp(' ');
tic; % 追赶法求解并开始计时
a=1;
c=-2;
u(1)=a/c;
q(1)=b(1)/c;
for i=2:n-2
u(i)=a/(c-u(i-1)*a);
end
for i=2:n-1
q(i)=(b(i)-q(i-1))/(c-u(i-1)*a);
end
y(n-1)=q(n-1);
for i=n-2:-1:1
y(i)=q(i)-u(i)*y(i+1);
end
Y(2:n)=y(1:n-1);%加上边界值
Y(1)=y0;
Y(n+1)=yl;
toc; %结束计时
T_E=abs(Y-Y_e); %计算截断误差Q=sum(T_E.^2); %求平方和为精度disp('所求计算精度为：');
Q
plot(x,Y,'m') %画出差分近似解hold on%将两个曲线绘在同一图上plot(x,Y_e,'g') %画出精确解。