数值微分计算方法实验

《数值方法》实验报告1数值微分计算方法实验郑发进2012042020022【摘要】数值微分(numericaldifferentiation)根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等)就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。关键词中心差分公式；理查森外推法；牛顿多项式微分；一、实验目的1.通过本次实验熟悉并掌握各种数值微分算法。2.掌握如何通过程序设计实现数值微分算法，从而更好地解决实际中的问题。二、实验原理1.精度为2()Oh的中心差分公式：设3,fCab，且,,,xhxxhab，则'2fxhfxhfxh而且存在数,ccxab，满足《数值方法》实验报告2',2truncfxhfxhfxEfhh其中322,6trunchfcEfhOh项,truncEfh称为截断误差。2.精度为4()Oh的中心差分公式：设5,fCab，且2,,,,2,xhxhxxhxhab，则'2fxhfxhfxh而且存在数,ccxab，满足2882'12fxhfxhfxhfxhfxh其中544,30trunchfcEfhOh项,truncEfh称为截断误差。3.理查森外推法：利用低阶公式推出高阶求解数值微分的公式，定理如下：设0'fx的两个精度为2()kOh的近似值分别为1kDh和12kDh，而且它们满足2220112'kkkfxDhchch和21220112'244kkkkkfxDhchch这样可得到改进的近似值表达式《数值方法》实验报告3112222012'()|()41kkkkkkkDhDhfxDhOhDhOh4.牛顿多项式微分：利用ft的1N个插值点可得近似ft的N次牛顿多项式NPt，其可表示为01020101NNNPtaattattttatttt则NPt的导数可以表示为11120100'NNNNjijkkPtaattttatt由此可知在,0,1,,MxMN处的导数值1100120101110'MNiiNMjiMiMiNiNijkiMkPtaattttattattatt三、实验内容1.P2601用程序6.1求解下列函数在x处的倒数近似值，精度为小数点后13位。注：有必要改写程序中的max1的值和h的初始值。45335223260322334777;135sin15tancos;13sincos1;1/215sin768;2;0.0001xxafxxxxxxxbfxxxcfxxxdfxxxxxefxxxP2602修改程序6.1，实现精度为O（h4）的中心差分公式（10）。用这个程序求《数值方法》实验报告4解上题中给出的函数导数的近似值。精度为小数点后13位。P2603使用程序6.2求解第1题中函数导数的近似值。精度为小数点后13位。注：有必要改变err，relerr和h的初始值。2.P2701修改程序6.3，使得可用它计算1,...2,1MxPMN），（。四、实验结果及分析（一）实验描述1.P2601已知下列函数，用数值微分方法求解在x处的导数近似值，精度为小数点后面的13位。45335223260322334777;135sin15tancos;13sincos1;1/215sin768;2;0.0001xxafxxxxxxxbfxxxcfxxxdfxxxxxefxxx（1）P2601使用精度2()Oh的中心差分公式求解（2）P2602使用精度4()Oh的中心差分公式求解（3）P2603使用理查森外推法求解2.P2701计算N阶牛顿插值多项式NPx的在插值点的导数值《数值方法》实验报告5',0,1,,NMPxMN（二）实验算法1.P2601计算精度为2()Oh导数近似值的算法.依据精度为2()Oh的中心差分公式，利用步长序列。Step1：输入函数fx，要计算的导数的横坐标x及相对误差容限toler，步长序列长度max；Step2：设置初始参量，初始步长0h，fx在x处的近似导数值00002fxhfxhDh令10/10hh，计算11112fxhfxhDh并令近似导数值序列前两个初值01100,EEDD，相对误差序列前两个初值011101,2/RREDD。Step3：fori=2uptomax1）计算步长1/10iihh时的近似导数值，2iiiifxhfxhDh；2）近似导数值序列间隔1iiiEDD，相对误差序列12/iiiiREDD。3）若1iiEE或iRtoler，则停止计算循环；《数值方法》实验报告6Step4：输出iD即可得到符合精度要求的近似导数值，记录ni。Step5：将,,,0,1,kkkDERkn储存到向量,,DER中列表输出。2.P2602计算精度为4()Oh导数近似值的算法.依据精度为2()Oh的中心差分公式，利用步长序列。Step1：输入函数fx，要计算的导数的横坐标x及相对误差容限toler，步长序列长度max；Step2：设置初始参量，初始步长0h，fx在x处的近似导数值000000288212fxhfxhfxhfxhDh令10/10hh，计算111111288212fxhfxhfxhfxhDh并令近似导数值序列前两个初值01100,EEDD，相对误差序列前两个初值011100,2/RREDD。Step3：fori=1uptomax1）计算步长1/10iihh时的近似导数值，288212iiiiiifxhfxhfxhfxhDh；2）近似导数值序列间隔1iiiEDD，相对误差序列12/iiiiREDD。3）若1iiEE或iRtoler，则停止计算循环；Step4：输出iD即可得到符合精度要求的近似导数值，记录ni。Step5：将,,,0,1,kkkDERkn储存到向量,,DER中输出。《数值方法》实验报告73.P2603利用理查森外推法设计计算满足精度要求的近似导数值的算法，Step1：输入函数fx，自变量x，相对误差容限toler，误差容限delta，初始步长0h，计算行数上限max；Step2：计算初始化量1）初始近似导数值000,002fxhfxhDh；2）初始误差00E，初始相对误差01R；3）1j；Step3：计算出符合精度要求的近似导数值：若1jEdelta且1jRtoler且1maxj，则进行下列计算，否则进入step41）1/2jjhh,,02jjjjDfxhfxhh;Fork=1uptoj,11,1,,141jkjkjkjkkDDDD2)计算误差,1,1jjjjjEDD3)计算相对误差,1,12/iijjjjREDD4)1jj,1nj返回step3Step4：输出,,,,nnnerrrelerrD，,nnD为符合要求的导数近似值。4.P2701利用牛顿多项式微分求解各插值点的近似的导数值算法流程如下：Step1：输入ix及其对应的函数值ifx，0,1,,iNStep2：fork=0uptoN，计算Mx处的近似导数值'NMPx1）0Mtx，1,11,,1iiiitxiMtxMiN，2)令01[(),(),,()]NAftftft，《数值方法》实验报告8forj=1uptoN+1fori=j+1uptoN+1()(()(1))/iijAiAiAitt3）p=1，()2dfkA，forn=2uptoN01()(1)idfkdfkpttAkStep3：向量df即Mx处的近似导数值'NMPx，0,1,,MN（三）实验结果1.P260（1）利用精度为2()Oh的中心差分公式求解结果（2）利用精度为4()Oh的中心差分公式求解结果（3）利用理查森外推法求解结果3种方法所求的结果如下，结果的精度为13位有效数字表1三种方法求解各个函数在其相应x处近似导数值对应函数及其求解点（以字母序号代替）2()Oh4()Oh理查森外推法（a）75.1735024617575.1734947073175.18107923463（b）1.2286029410281.2285974234041.228733630836（c）1.9515627659521.9515596090241.951752228326（d）2.9655290467632.9655148293492.966366834548（e）1.0152010645781.0150610589881.015008921184由上表可知，三种方法的求解结果基本相同，但是因为函数及求解点的设置的初始步长和终止条件的不同，三种方法的求解结果没有做到高度一致，因此可知求解结果与初始步长及终止条件的设定有较大的关系。五、实验结论《数值方法》实验报告9本次数值实验利用了精度为4()Oh、2()Oh的中心差分公式及理查森外推公式计算了各个函数在其对应点的近似导数值，三种方法的求解结果基本相同，但是因为函数及求解点的设置的初始步长和终止条件的不同，三种方法的求解结果没有做到高度一致，因此可知求解结果与初始步长及终止条件的设定有较大的关系。另外还设计了计算利用牛顿插值多项式的任意插值点导数计算算法，并且用相应程序实现，其中将插值点重新排列，使得所求点为第一个点是一个较大的创新，这样使得各个插值点的导数计算更为快捷方便。附件（代码）：1.P260精度为2()Oh近似导数值的计算子程序functiondf=cetr2df(f,x,toler,max,h0)%输入输出参数说明：%df输出的近似导数值%f、x相应函数及其点，toler误差容忍上限，max最大计算次数，h0初始步长fori=1:2h=10^(-(i-1))*h0;D(i)=(feval(f,x+h)-feval(f,x-h))/(2*h);%计算导数近似值序列的前两个值endE(1)=0;E(2)=abs(D(2)-D(1));%计算导数近似值间隔序列的前两个值R(1)=1;R(2)=2*E(2)/(abs(D(1))+abs(D(2)));%计算导数近似值相对误差序列的前两个值fori=3:max%计算导