利用中心差分格式数值求解导数.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追赶法求解线性方程组需要用到的中间辅助向量四、计算结果与讨论不同步长的数值结果函数曲线与精确解的对比 从对比图中可以看到，在所取的四个计算步长下，数值计算结果与精确解都符合得相当好StepAccuracy（Max-error）0.058.152404966677018E-0050.013.264604683472783E-0060.0013.817221055912867E-0080.00019.862256566961491E-009不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量，即上表中的Max-error来度量。

从上表中可以看出，随着网格节点的加密，Max-error的数量级在降低，即数值解的精度提高上述四个步长中，将线段离散成一万个网格时，数值结果的精度最高五、源程序program mainimplicit nonecharacter(13) filename !定义了五个字符串变量，用于按输出数据的需要自动更改dat文件的文件名character(3) Acharacter(6) Bcharacter(4) Ccharacter(3) Dinteger :: n,ireal(8) :: h,error,preal(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 = "X", "Y"' write(10,"('ZONE T=""Problem1-',f8.5,'"", I=',I6,', F=POINT')")h,n+1 do i=0,n write(*,"(F10.8,10x,F10.8)")x(i),y(i) write(10,"(F10.8,10x,F10.8)")x(i),y(i) !将数值解数据写入dat文件 end do y(0)=0 y(n)=1 error=0.0 do i=1,n-1 p=y(i) y(i)=(1+0.25*sin(2.0))*(i*h)-0.25*sin(2.0*i*h) !求解节点精确值 error=max(error,abs(p-y(i))) !寻找各个节点误差的最大值 end do write(10,"('ZONE T=""Problem1-',f8.5,'-exact"", I=',I6,', F=POINT')")h,n+1 do i=0,n write(*,"(F10.8,10X,F10.8)")x(i),y(i) write(10,"(F10.8,10X,F10.8)")x(i),y(i) !将精确解数据写入dat文件 end do D="er-" !er 代表error，这里指精确值和数值计算值之间的差别 filename=D//B//C open(unit=12,file=filename) write(12,*)"max-error=",error !将误差最大值写入dat文件 write(*,*)n,error write(*,*)filename stop end !主程序结束 subroutine subsolution(y,n,h) !子程序 implicit none integer ::n,i real(8) :: h real(8) :: y(0:n) !数组和变量的声明不能同时进行 real(8),allocatable :: a(:),b(:),c(:),s(:),t(:),f(:),g(:) !声明可变长度数组 allocate(a(1:n),b(0:n),c(0:n-1),s(0:n),t(0:n-1),f(0:n),g(0:n)) !指定可变数组的长度 a=1 !对数组a，b，c赋值 b=-2 c=1 a(n)=0 b(0)=1 b(n)=1 c(0)=0 f(0)=0 !对数组f 赋值 f(n)=1 do i=1,n-1 f(i)=h*h*sin(2.0*real(i)*h) end do s(0)=b(0) !计算s，t t(0)=c(0)/b(0) do i=1,n-1 s(i)=b(i)-a(i)*t(i-1) t(i)=c(i)/s(i) end do s(n)=b(n)-a(n)*t(n-1) g(0)=f(0)/b(0) !追过程 do i=1,n g(i)=(f(i)-a(i)*g(i-1))/s(i) end do y(n)=g(n) !赶过程 do i=n-1,1,-1 y(i)=g(i)-t(i)*y(i+1) end do y(0)=0 y(n)=1 return end subroutine subsolution。