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

CFD实验报告一姓名： 学号：一、题目：利用中心差分格式近似导数，数值求解常微分方程 （） 步长分别取0.05, 0.01, 0.001，0.0001。二、报告要求：1）列出全部计算公式和步骤；2）表列出程序中各主要符号和数组意义；3）绘出数值计算结果的函数曲线，并与精确解比较；4）比较不同差分格式和不同网格步长计算结果的精度和代价；5）附源程序。三、相关差分格式 二阶导数的三点差分格式有向前差分、向后差分和中心差分，表达式分别如下：代入微分方程可以得到差分方程，表达式分别如下：对于三种差分格式，差分格式可以改写成的形式，其中是相同的，非齐次项不同，如下所示：系数矩阵一阶向前差分一阶向后差分二阶中心差分求解可以得到各节点的值。四、计算公式和步骤；1.关于精确解的推导：已知，对x进行两次积分，得到，再结合边界条件和得到相对应的和，确定最后精确解为：。2.关于数值求解方法：对于方程组可直接求解，也可以使用追赶法求解，下面介绍简单追赶法求解三对角方程组的过程。对于三对角方程组：系数矩阵A中仅三对角线上的数值不全为0，其余位置上的数值全为0，是典型的对角占优的三对角矩阵，利用高斯消去法，经过n-1次消元可以化为同解方程组：如上所示，求这些值的过程（即消元过程）称为追：再利用回代过程求出方程组的各变量：这一逆序求变量的过程（回代过程）称为赶。五、计算结果与分析根据题意需要分别取步长0.05, 0.01, 0.001，0.0001进行计算，因此采用MATLAB进行编程运算，然后将数值解与精确解进行比较，如下图所示：步长0.05步长0.01步长0.001步长0.0001从上面四个图可以看出，这几个格式的精确解和数值解之间符合度很好，其中中心差分格式精确度更高，并且随着步长的减小与精确解的符合度越高。下面将计算结果用表格的形式列出：追赶法求解0.050.010.0010.0001计算耗时前差0.0.0.1.后差0.0.0.1.中心差0.0.0.1.表1 不同格式不同步长计算耗时（单位s）追赶法求解0.050.010.0010.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、从时间代价来看，三种格式采用追赶法求解所需要的时间基本一样，从计算精度来看，中心差分格式的计算精度明显优于向前差分和向后差分格式。2、随着步长的减小，虽然求解的计算精度增高，但所需要的时间却大大增加。3、对于这三种格式而言，中心差分格式明显优于另两种格式。4、从上可知，提高计算精度可以通过减少步长和采用高阶格式的方法，由于减少步长需要用时间作为代价，所以采用高阶格式相对而言比较好。六、源程序及其说明符号说明：long变量取值范围dx步长n+1节点数y0始端边界值yl末端边界值A系数矩阵b非齐次项向量Y数值解向量Y_e精确解向量T_E节点截断误差Q精度源程序如下：clear; %清空long=1;%变量取值范围dx=0.05; %设置步长n=long/dx; % 此处n=节点数1y0=0;%设定始端边界条件yl=1-sin(2)/4; %设定末端边界条件A=zeros(n-1,n-1);%设置系数矩阵AA(1,1:2)=-2 1;%始端边界for i=2:n-2 A(i,i-1:i+1)=1,-2,1;endA(n-1,n-2:n-1)=1,-2;%末端边界b=zeros(n-1,1);%设置矩阵bfor i=1:n-1% b(i)=dx2*sin(i-1)*2*dx); %一阶向前差分% b(i)=dx2*sin(i+1)*2*dx); %一阶向后差分 b(i)=dx2*sin(i*2*dx); %中心差分endb(1)=b(1)-y0;%设置边界条件b(n-1)=b(n-1)-yl;x=dx*(0:n);%求解精确解Y_e=x-sin(2*x)/4;disp(SC李文建);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);endfor i=2:n-1 q(i)=(b(i)-q(i-1)/(c-u(i-1)*a);endy(n-1)=q(n-1);for i=n-2:-1:1 y(i)=q(i)-u(i)*y(i+1);endY(2:n)=y(1:n-1);%加上边界值Y(1)=y0;Y(n+1)=yl;toc; %结束计时T_E=abs(Y-Y_e);