利用中心差分格式数值求解导数.doc

利用中心差分格式数值求解导数 目录 一、问题描述 1 二、格式离散 1 二阶导数中心差格式离散 2 追赶法求解线性方程组简述 3 计算流程图 4 三、程序中主要符号和数组意义 5 四、计算结果与讨论 5 五、源程序 7 一、问题描述 利用中心差分格式近似导数，数值求解 步长分别取 二、格式离散 将x轴上[0,1]之间的线段按上述步长，等步长的离散为n个小段，包括端点，共n+1个网格节点，示意图如下： 线段上边的数字表示x轴上的坐标值，线段下边的数字表示节点编号，从0到n编号。 二阶导数中心差格式离散 整理为线性方程形式 其中， 为空间离散步长；i=1,2，……，n-1 包括边界条件的线性方程组如下： 改写成矩阵形式： 其中，， ， 系数矩阵A中仅三对角线上的数值不全为0，其余位置上的数值全为0，是典型的对角占优的三对角矩阵，列向量f中，，且，作为边界条件。 追赶法求解线性方程组简述 对A做Crout分解，即 ,， 其中为待定常数，由矩阵乘法可以得到下面的式子： 将对角占优三对角矩阵线性方程组等价为如下两个方程组 ， 求解对角占优三对角矩阵线性方程组的追赶法步骤： ①输入数据 ②计算 ③求解方程组 ④求解方程组 ⑤输出 计算流程图 三、程序中主要符号和数组意义 符号或数组 意义 A、B、C、D、filename 用于自动更改dat文件名的字符串变量 h 离散步长 n 离散网格数，共n+1个网格节点 p 辅助变量，暂时记录网格节点上的y值 数组x，y 离散节点的x，y坐标 子程序数组a,b,c 记录系数矩阵占优对角线上的值 子程序数组f 记录线性方程组常数向量 子程序数组s,r,t,g 追赶法求解线性方程组需要用到的中间辅助向量 四、计算结果与讨论 不同步长的数值结果函数曲线与精确解的对比 从对比图中可以看到，在所取的四个计算步长下，数值计算结果与精确解都符合得相当好 Step Accuracy（Max-error） 0.05 8.152404966677018E-005 0.01 3.264604683472783E-006 0.001 3.817221055912867E-008 0.0001 9.862256566961491E-009 不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量，即上表中的Max-error来度量。从上表中可以看出，随着网格节点的加密，Max-error的数量级在降低，即数值解的精度提高。上述四个步长中，将线段离散成一万个网格时，数值结果的精度最高。 五、源程序 program main implicit none character(13) filename !定义了五个字符串变量，用于按输出数据的需要自动更改dat文件的文件名 character(3) A character(6) B character(4) C character(3) D integer :: n,i real(8) :: h,error,p real(8),allocatable :: y(:),x(:) !声明可变长度数组，x、y轴坐标定义为可变数组，数组长度按需要自动更改 write(*,*)输入步长： read (*,*)h !读入空间步长 n=NINT(1.0/h) !nint命令，取与浮点数最接近的整数 allocate (y(0:n),x(0:n)) !指定可变数组的长度 A=po- !po 代表problem open(unit=11,file=help.txt) !打开一个txt文件，用于帮助更改dat文件的文件名 write(11,(f6.4))h rewind(11) ! 将文件的读写位置移回到文件的最前面 read(11,(A6))B C=.dat filename=A//B//C call subsolution(y,n,h) !调用追赶法求解线性方程组的子程序 open(unit=10,file=filename) do i=0,n x(i)=real(i)*h end do write(10,*)TITLE = X - Y CURVE !写入到dat文件中的一段字符，便于用tecplot软件后处理计算数据 write(10,*)VARIABLES =