数值微分计算方法实验

1、数值微分计算方法实验【摘要】数值微分(numerical differentiation根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟

2、合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。关 键 词 中心差分公式；理查森外推法；牛顿多项式微分；1、实验目的1.通过本次实验熟悉并掌握各种数值微分算法。2.掌握如何通过程序设计实现数值微分算法，从而更好地解决实际中的问题。2、实验原理1. 精度为的中心差分公式：设，且，则 而且存在数，满足其中项称为截断误差。2. 精度为的中心差分公式：设，且，则 而且存在数，满足其中项称为截断误差。3.理查森外推法：利用低阶公式推出高阶求解数值微分的公式，定理如下：设的两个精度为的近似值分别为和，而且它们满足和这

3、样可得到改进的近似值表达式4. 牛顿多项式微分：利用的个插值点可得近似的次牛顿多项式，其可表示为则的导数可以表示为由此可知在处的导数值3、实验内容1.P260 1用程序6.1求解下列函数在x处的倒数近似值，精度为小数点后13位。注：有必要改写程序中的max1的值和h的初始值。P260 2修改程序6.1，实现精度为O（h4）的中心差分公式（10）。用这个程序求解上题中给出的函数导数的近似值。精度为小数点后13位。P260 3使用程序6.2求解第1题中函数导数的近似值。精度为小数点后13位。注：有必要改变err，relerr和h的初始值。2.P270 1 修改程序6.3，使得可用它计算。四、实验结

4、果及分析（一）实验描述1.P260 1已知下列函数，用数值微分方法求解在处的导数近似值，精度为小数点后面的13位。（1）P260 1使用精度的中心差分公式求解（2）P260 2使用精度的中心差分公式求解（3）P260 3使用理查森外推法求解2.P270 1 计算阶牛顿插值多项式的在插值点的导数值（二）实验算法1.P260 1计算精度为导数近似值的算法.依据精度为的中心差分公式，利用步长序列。Step1：输入函数，要计算的导数的横坐标及相对误差容限toler，步长序列长度；Step2：设置初始参量，初始步长，在处的近似导数值令，计算并令近似导数值序列前两个初值，相对误差序列前两个初值。Step3

5、：for i=2 upto max1） 计算步长时的近似导数值，；2） 近似导数值序列间隔，相对误差序列。3） 若或，则停止计算循环；Step4：输出即可得到符合精度要求的近似导数值，记录。Step5：将储存到向量中列表输出。2.P260 2计算精度为导数近似值的算法.依据精度为的中心差分公式，利用步长序列。Step1：输入函数，要计算的导数的横坐标及相对误差容限toler，步长序列长度；Step2：设置初始参量，初始步长，在处的近似导数值令，计算并令近似导数值序列前两个初值，相对误差序列前两个初值。Step3：for i=1 upto max1） 计算步长时的近似导数值，；2） 近似导数值序

6、列间隔，相对误差序列。3） 若或，则停止计算循环；Step4：输出即可得到符合精度要求的近似导数值，记录。Step5：将储存到向量中输出。3.P260 3利用理查森外推法设计计算满足精度要求的近似导数值的算法，Step1：输入函数，自变量，相对误差容限toler，误差容限delta，初始步长，计算行数上限；Step2：计算初始化量1） 初始近似导数值；2） 初始误差，初始相对误差；3） ；Step3：计算出符合精度要求的近似导数值：若且且，则进行下列计算，否则进入step41）,;For k=1 upto j2 计算误差3计算相对误差4,返回step3Step4：输出 ，为符合要求的导数近似值

7、。4.P270 1 利用牛顿多项式微分求解各插值点的近似的导数值算法流程如下：Step1：输入及其对应的函数值，Step2：for k=0 upto N，计算处的近似导数值1），2令，for j=1 upto N+1for i=j+1 upto N+13）p=1，for n=2 upto NStep3：向量即处的近似导数值，（3）实验结果1.P260 （1）利用精度为的中心差分公式求解结果（2）利用精度为的中心差分公式求解结果（3）利用理查森外推法求解结果3种方法所求的结果如下，结果的精度为13位有效数字表1 三种方法求解各个函数在其相应x处近似导数值对应函数及其求解点（以字母序号代替）理查森

8、外推法（a）75.1735024617575.1734947073175b）1.2286029410281.2285974234041.228733630836（c）1.9515627659521.9515596090241.951752228326（d）2.9655290467632.9655148293492.966366834548（e）1.0152010645781.0150610589881.015008921184由上表可知，三种方法的求解结果基本相同，但是因为函数及求解点的设置的初始步长和终止条件的不同，三种方法的求解结果没有做到高度一致，因此可知求解结

9、果与初始步长及终止条件的设定有较大的关系。五、实验结论本次数值实验利用了精度为、的中心差分公式及理查森外推公式计算了各个函数在其对应点的近似导数值，三种方法的求解结果基本相同，但是因为函数及求解点的设置的初始步长和终止条件的不同，三种方法的求解结果没有做到高度一致，因此可知求解结果与初始步长及终止条件的设定有较大的关系。另外还设计了计算利用牛顿插值多项式的任意插值点导数计算算法，并且用相应程序实现，其中将插值点重新排列，使得所求点为第一个点是一个较大的创新，这样使得各个插值点的导数计算更为快捷方便。附件（代码）：1.P260 精度为近似导数值的计算子程序function df=cetr2df(

10、f,x,toler,max,h0%输入输出参数说明：%df 输出的近似导数值%f、x相应函数及其点，toler误差容忍上限，max最大计算次数，h0初始步长for i=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;%计算导数近似值相对误差序列的前两个值for i=3:max%计算导数近似值序列h=h/10;D(i=(feval(f,x+h-feval(f,x-h

11、/(2*h;E(i=abs(D(i-D(i-1;R(i=2*E(i/(abs(D(i+abs(D(i-1;if (E(i>E(i-1|(R(i<=toler%计算导数近似值序列的终止条件break;endendn=i;df=D(n;%输出满足精度要求的数值导数近似值end精度为近似导数值的计算子程序function df=cetr4df(f,x,toler,max,h0%输入输出参数说明：%df 输出的近似导数值%f、x相应函数及其点，toler误差容忍上限，max最大计算次数，h0初始步长for i=1:2%计算导数近似值序列的前两个值h=10(-(i-1*h0;D(i=(-fe

12、val(f,x+2*h+8*feval(f,x+h-8*feval(f,x-h+feval(f,x-2*h/(12*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; %输出满足精度要求的数值导数近似值for i=3:max %计算导数近似值序列h=h/10;D(i=(-feval(f,x+2*h+8*feval(f,x+h-8*feval(f,x-h+feval(f,x-2*h/(12*h;E(i=abs(D(i-D(i-1;R(i=2*E(i/(abs(D(i+abs(D(i-1;if

13、 (E(i>E(i-1&(R(i<=tolerbreak;endendn=i;df=D(n; %输出满足精度要求的数值导数近似值err=E(n;relerr=R(n;end理查森外推法计算导数近似值的子函数function df=Richarson(f,x,toler,delta,max,h0%输入输出参数说明：%df 输出的近似导数值%f、x相应函数及其点delta相对误差容忍上限，toler误差容忍上限，max最大计算次数，h0初始步长h=h0;D(1,1=(feval(f,x+h-feval(f,x-h/(2*h;err=1;relerr=1;j=1;while(er

14、r>delta&(relerr>toler&(j<=max%进行是否终止判断h=h/2;D(j+1,1=(feval(f,x+h-feval(f,x-h/(2*h;for k=1:j %理查森外推公式计算高阶精度的导数计算公式D(j+1,k+1=D(j+1,k+(D(j+1,k-D(j+1,k/(4(k-1;enderr=abs(D(j+1,j+1-D(j,j;relerr=2*err/(abs(D(j+1,j+1+abs(D(j,j;j=j+1;endn=jdf=D(n,n;%输出得到满足得到精度要求的近似导数值errrelerrEnd2.P270 1 function df=difnew(X,YN=length(X-1;A=Y;for M=1:N+1 %对X序列重新排列得到新的插值点序列if M=1t=X;else if M=N+1for i=2:N+1t(i=X(i-1;en