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

1、CFD实验报告姓名：学号、题目：利用中心差分格式近似导数d2y/dx2，数值求解常微分方程(0 x1)d2y=sin2xdx2sinT步长分别取Ax=0.05,0.01,0.001,0.0001。二、报告要求：列出全部计算公式和步骤；表列出程序中各主要符号和数组意义；绘出数值计算结果的函数曲线，并与精确解比较；比较不同差分格式和不同网格步长计算结果的精度和代价附源程序。三、相关差分格式二阶导数d2y/dx2的三点差分格式有向前差分、向后差分和中心差分，表达式分别如下：TOC o 1-5 h z一阶向前差分：色上j+i*j+O(Ax)dx2Ax2一阶向后差分：凹=红jT+Uj+O(Ax)dx2A

2、x2二阶中心差分：旦=U+厂2uj+j+O(Ax2)dx2Ax2代入微分方程可以得到差分方程，表达式分别如下：一阶向前差分一阶向后差分二阶中心差分u2u+u.j2/1j=sin2xAx2ju2u+u-j2jj=sin2x HYPERLINK l bookmark10 o Current Document Ax2ju2u+uj+1jj1=sin2x HYPERLINK l bookmark14 o Current Document Ax2j对于三种差分格式，差分格式可以改写成AY=b的形式，其中A是相同的,非齐次项b不同，如下所示：-211-21系数矩阵1-211-2-y0sin(2xZAx21

3、sin(2x)2hx2sin(2x呱2-yk-21sin(2x)2hx2-ysin(2x)2Ax23sin(2x)2Ax2sin(2xllhx2-yk1sin(2x)2Ax2-ysin(2x)2Ax22sin(2x)2hx2sin(2x良2-yk-11一阶向前差分一阶向后差分二阶中心差分求解AY=b可以得到各节点y的值Y=yy12四、计算公式和步骤；1.关于精确解的推导：已知山=sin2x，对X进行两次积分，得到y二dx21sin2x+Cx+C，再结合412边界条件yx=0二0和yx广1-晋得到相对应的Ci和C2，确定最后精确解为:2.关于数值求解方法：对于方程组AY=b可直接求解，也可以使用

4、追赶法求解，下面介绍简单追赶法求解三对角方程组的过程。对于三对角方程组：dc11adc222adn1n1any1y22cybn1n1n1dybnnn系数矩阵A中仅三对角线上的数值不全为0,其余位置上的数值全为0,是典型的对角占优的三对角矩阵，利用高斯消去法，经过n-1次消元可以化为同解方程组：1u1y1y2q1q21u2-1uyqn1n1n11_y_q_nn如上所示，求这些值的过程（即消元过程）称为追：u1uiqi=c/d,q=b/d,=c/(d一ua)(i=2,.,n=(b-qa)/(d-ua)(iii一1iii一1i一1)=2,.,n)再利用回代过程求出方程组的各变量：y=q,nny=q-

5、uyiiii+1(i这一逆序求变量的过程（回代过程）称为赶。五、计算结果与分析根据题意需要分别取步长Ax=0.05,0.01,0.001,0.0001进行计算，因此采用MATLAB进行编程运算，然后将数值解与精确解进行比较，如下图所示：从上面四个图可以看出，这几个格式的精确解和数值解之间符合度很好，其中中心差分格式精确度更高，并且随着步长的减小与精确解的符合度越高。下面将计算结果用表格的形式列出：追赶法求解Ax=0.05Ax=0.01Ax=0.001Ax=0.0001计算耗时、/.、八前差0.0068530.0148030.0489311.111667后差0.0066450.0141160.0

6、486621.103086中心差0.0055260.0158000.0596941.072931表1不同格式不同步长计算耗时（单位s）追赶法求解Ax=0.05Ax=0.01Ax=0.001Ax=0.0001计算精度、/.、八前差4.6976e-048.5578e-058.3727e-068.3542e-07后差3.6746e-048.1482e-058.3317e-068.3501e-07中心差6.8119e-085.4443e-105.4441e-135.3257e-16表2不同格式不同步长计算精度由不同格式不同步长计算耗时和计算精度的结果对比可知：1、从时间代价来看，三种格式采用追赶法求解

7、所需要的时间基本一样，从计算精度来看，中心差分格式的计算精度明显优于向前差分和向后差分格式。2、随着步长的减小，虽然求解的计算精度增高，但所需要的时间却大大增加。3、对于这三种格式而言，中心差分格式明显优于另两种格式。4、从上可知，提高计算精度可以通过减少步长和采用高阶格式的方法，由于减少步长需要用时间作为代价，所以采用高阶格式相对而言比较好。六、源程序及其说明符号说明：long变量取值范围dx步长n+1节点数y0始端边界值yl末端边界值A系数矩阵b非齐次项向量Y数值解向量Ye精确解向量TE节点截断误差Q精度源程序如下：clear;%清空long=l;%变量取值范围dx=0.05;%设置步长n

8、=long/dx;%此处n=节点数一1y0=0;%设定始端边界条件yl=1-sin(2)/4;%设定末端边界条件A=zeros(n-1,n-1);%设置系数矩阵AA(1,1:2)=-21;%始端边界fori=2:n-2A(i,i-1:i+1)=1,-2,1;endA(n-1,n-2:n-1)=1,-2;%末端边界b=zeros(n-1,1);%设置矩阵bfori=1:n-1%b(i)=dxA2*sin(i-1)*2*dx);%阶向前差分%b(i)=dxA2*sin(i+1)*2*dx);%阶向后差分b(i)=dxA2*sin(i*2*dx);%中心差分endb(1)=b-yO;%设置边界条件b(n-1)=b(n-1)-yl;x=dx*(0:n);%求解精确解Y_e=x-sin(2*x)/4;);disp(SC16OO5O41李文建disp();tic;%追赶法求解并开始计时a=1;c=-2;u(1)=a/c;q(1)=b(1)/c;fori=2:n-2u(i)=a/(c-u(i-1)*a);endfori=2:n-1q(i)=(b(i)-q(i-1)/(c-u(i-1)*a);endy(n-1)=q(n-1);fori=n-2:-1:1y(i)=q(i)-u(i)*y(i+1);endY(2:n)=y(l:n-1);%加上边界值Y(1