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

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

ID:9901113

大小:272.50 KB

页数:12页

时间:2018-05-14

利用中心差分格式数值求解导数_第1页
利用中心差分格式数值求解导数_第2页
利用中心差分格式数值求解导数_第3页
利用中心差分格式数值求解导数_第4页
利用中心差分格式数值求解导数_第5页
资源描述:

《利用中心差分格式数值求解导数》由会员上传分享,免费在线阅读,更多相关内容在学术论文-天天文库

1、利用中心差分格式数值求解导数目录一、问题描述1二、格式离散1二阶导数中心差格式离散2追赶法求解线性方程组简述3计算流程图4三、程序中主要符号和数组意义5四、计算结果与讨论5五、源程序7一、问题描述利用中心差分格式近似导数,数值求解步长分别取二、格式离散将x轴上[0,1]之间的线段按上述步长,等步长的离散为n个小段,包括端点,共n+1个网格节点,示意图如下:线段上边的数字表示x轴上的坐标值,线段下边的数字表示节点编号,从0到n编号。二阶导数中心差格式离散整理为线性方程形式其中,为空间离散步长;i=1,2,……,n-1包括边界条件的线性方程组如下:改写成矩阵形式:其中,,,系数矩阵A中仅三对角线

2、上的数值不全为0,其余位置上的数值全为0,是典型的对角占优的三对角矩阵,列向量f中,,且,作为边界条件。追赶法求解线性方程组简述对A做Crout分解,即,,其中为待定常数,由矩阵乘法可以得到下面的式子:将对角占优三对角矩阵线性方程组等价为如下两个方程组,求解对角占优三对角矩阵线性方程组的追赶法步骤:①输入数据②计算③求解方程组④求解方程组⑤输出计算流程图三、程序中主要符号和数组意义符号或数组意义A、B、C、D、filename用于自动更改dat文件名的字符串变量h离散步长n离散网格数,共n+1个网格节点p辅助变量,暂时记录网格节点上的y值数组x,y离散节点的x,y坐标子程序数组a,b,c记录

3、系数矩阵占优对角线上的值子程序数组f记录线性方程组常数向量子程序数组s,r,t,g追赶法求解线性方程组需要用到的中间辅助向量四、计算结果与讨论不同步长的数值结果函数曲线与精确解的对比从对比图中可以看到,在所取的四个计算步长下,数值计算结果与精确解都符合得相当好StepAccuracy(Max-error)0.058.7018E-0050.013.2783E-0060.0013.2867E-0080.00019.1491E-009不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量,即上表中的Max-error来度量。从上表中可以看出,随着网格节点的加密,Max-

4、error的数量级在降低,即数值解的精度提高。上述四个步长中,将线段离散成一万个网格时,数值结果的精度最高。五、源程序programmainimplicitnonecharacter(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(*,*)"输入步

5、长:"read(*,*)h!读入空间步长n=NINT(1.0/h)!nint命令,取与浮点数最接近的整数allocate(y(0:n),x(0:n))!指定可变数组的长度A="po-"!po代表problemopen(unit=11,file="help.txt")!打开一个txt文件,用于帮助更改dat文件的文件名write(11,"(f6.4)")hrewind(11)!将文件的读写位置移回到文件的最前面read(11,"(A6)")BC=".dat"filename=A//B//Ccallsubsolution(y,n,h)!调用追赶法求解线性方程组的子程序open(unit=10,f

6、ile=filename)doi=0,nx(i)=real(i)*henddowrite(10,*)'TITLE="X-YCURVE"'!写入到dat文件中的一段字符,便于用tecplot软件后处理计算数据write(10,*)'VARIABLES="X","Y"'write(10,"('ZONET=""Problem1-',f8.5,'"",I=',I6,',F=POINT')")h,n+1doi=0,nwrite(*,"(F10.8,10x,F10.8)")x(i),y(i)write(10,"(F10.8,10x,F10.8)")x(i),y(i)!将数值解数据写入dat文件enddo

7、y(0)=0y(n)=1error=0.0doi=1,n-1p=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)))!寻找各个节点误差的最大值enddowrite(10,"('ZONET=""Problem1-',f8.5,'-exact"",I=',I6,',F=POINT')")h,n+

当前文档最多预览五页,下载文档查看全文

此文档下载收益归作者所有

当前文档最多预览五页,下载文档查看全文
温馨提示:
1. 部分包含数学公式或PPT动画的文件,查看预览时可能会显示错乱或异常,文件下载后无此问题,请放心下载。
2. 本文档由用户上传,版权归属用户,天天文库负责整理代发布。如果您对本文档版权有争议请及时联系客服。
3. 下载前请仔细阅读文档内容,确认文档内容符合您的需求后进行下载,若出现内容与标题不符可向本站投诉处理。
4. 下载文档时可能由于网络波动等原因无法下载或下载错误,付费完成后未能成功下载的用户请联系客服处理。