1、 数值方法实验报告 1 数值微分计算方法实验 【摘要】 数值 微分 (numerical differentiation)根据函数在一些 离散点 的 函数值 ，推算它在某点的 导数 或 高阶导数 的近似值的方法。通常用 差商 代替 微商 ，或者用一个能够近似代替该 函数 的较简单的 可微 函数（如 多项式 或 样条函数 等）的相应导数作为能 求导 数的近似值。例如一些常用的数值微分 公式 (如两点公式、三点公式等 )就是在等距步长情形下用 插值 多项式的导数作为近似值的。此外，还可以采用 待定系数法 建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的 精确度 。当函数可微性不太好时，
2、利用 样条 插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的 随机误差 ，应该用 曲线拟合代替 函数 插值，然后用拟合曲线的 导数 作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值 微分 公式还是微分方程数值解法的重要依据。 关 键 词 中心差分公式；理查森外推法；牛顿多项式微分； 一、 实验目的 1.通过本次实验熟悉并掌握各种数值微分算法。 2.掌握如何通过程序设计实现数值微分算法，从而更好地解决实际中的问题。 二、 实验原理 1. 精度为 2()Oh 的中心差分公式： 设 3 ,f C a b ， 且 , , ,x h x x h a b ， 则 2f x
3、 h f x hfx h 而且存在数 ,c c x a b，满足 数值方法实验报告 2 ,2 tr u n cf x h f x hf x E f hh 其中 32 2, 6tr u n c h f cE f h O h 项 ,truncE f h 称为截断误差。 2. 精度为 4()Oh 的中心差分公式： 设 5 ,f C a b ， 且 2 , , , , 2 ,x h x h x x h x h a b ， 则 2f x h f x hfx h 而且存在数 ,c c x a b，满足 2 8 8 2 12f x h f x h f x h f x hfx h 其中 54 4, 30tru
4、 n c h f cE f h O h 项 ,truncE f h 称为截断误差。 3.理查森外推法： 利用低阶公式推出高阶求解数值微分的公式，定理如下： 设 0fx的两个精度为 2()kOh 的近似值分别为 1kDh 和 1 2kDh ，而且它们满足 2 2 20 1 1 2 kkkf x D h c h c h 和 2 1 2 20 1 1 2 2 4 4k k k kkf x D h c h c h 这样可得到改进的近似值表达式 数值方法实验报告 3 112 2 2 201 2 ( ) | ( )41kkkkkk kD h D hf x D h O h D h O h 4. 牛顿多项式微
5、分： 利用 ft的 1N 个插值点可得近似 ft的 N 次牛顿多项式 NPt，其可表示为 0 1 0 2 0 1 0 1N N NP t a a t t a t t t t a t t t t 则 NPt的 导数可以表示为 1 11 2 0 1 00 N NN N j ijkkP t a a t t t t a t t 由此 可知在 , 0,1, ,Mx M N 处 的导数值 1 1 0 01 2 0 1 0 1 1 10 M N i iN M j i M i M i N i N ij k i MkP t a a t t t t a t t a t t a t t 三、 实验内容 1.P260
6、 1 用程序 6.1 求解下列函数在 x 处的倒数近似值，精度为小数点后 13 位。注：有必要改写程序中的 max1 的值和 h 的初始值。 4 5 3 3 5 22326 0 3 2 2 3 3 4 7 7 7 ; 1 35 sin 15ta n c o s ;13sin c o s 1 ; 1 / 215sin 7 6 8 ;2; 0 .0 0 0 1xxa f x x x x x xxb f x xxc f x x xd f x x x x xe f x x x P260 2 修改程序 6.1，实现精度为 O（ h4）的中心差分公式（ 10）。用这个程序求数值方法实验报告 4 解上题中给
7、出的函数导数的近似值。精度为小数点后 13 位。 P260 3 使用程序 6.2 求解第 1 题中函数导数的近似值。精度为小数点后 13 位。注：有必要改变 err， relerr 和 h 的初始值。 2.P270 1 修改程序 6.3，使得可用它计算 1,.2,1MxP M N），（ 。 四、实验结果及分析 （一）实验描述 1.P260 1 已知下列函数，用数值微分方法求解在 x 处的导数近似值，精度为小数点后面的 13 位。 4 5 3 3 5 22326 0 3 2 2 3 3 4 7 7 7 ; 1 35 sin 15ta n c o s ;13sin c o s 1 ; 1 / 21
8、5sin 7 6 8 ;2; 0 .0 0 0 1xxa f x x x x x xxb f x xxc f x x xd f x x x x xe f x x x （ 1） P260 1 使用精度 2()Oh 的中心差分公式求解 （ 2） P260 2 使用精度 4()Oh 的中心差分公式求解 （ 3） P260 3 使用理查森外推法 求解 2.P270 1 计算 N 阶牛顿插值多项式 NPx的在插值点的导数值数值方法实验报告 5 , 0 ,1, ,NMP x M N （二）实验算法 1.P260 1 计算 精度 为 2()Oh 导数 近似值的算法 . 依据精度为 2()Oh 的中心差分公式
9、， 利用 步长 序列 。 Step1： 输入 函数 fx，要计算 的 导数的 横坐标 x 及 相对误差容限 toler， 步长序列长度 max ； Step2：设置 初始参量，初始 步长 0h ， fx在 x 处 的近 似 导数值 00002f x h f x hD h 令 10/10hh ，计算 111 12f x h f x hD h 并令 近似 导数值 序列 前 两个初值 0 1 1 00,E E D D ， 相对 误差 序列前 两个初值 0 1 1 1 01 , 2 /R R E D D 。 Step3： for i=2 upto max 1） 计算 步长 1/10iihh 时 的近似
10、导数值， 2iiiif x h f x hD h ； 2） 近似 导数 值序列 间隔 1i i iE D D ， 相对误差序列 12/i i i iR E D D 。 3） 若 1iiEE 或 iR toler ， 则 停止 计算 循环； 数值方法实验报告 6 Step4： 输出 iD 即 可 得到 符合精度要求的近似导数值 ， 记录 ni 。 Step5： 将 , , , 0 ,1 ,k k kD E R k n 储存 到向量 ,DER 中 列表 输出 。 2.P260 2 计算 精度 为 4()Oh 导数 近似值的算法 . 依据精度为 2()Oh 的中心差分公式， 利用 步长 序列 。 S
11、tep1： 输入 函数 fx，要计算 的 导数的 横坐标 x 及 相对误差容限 toler， 步长序列长度 max ； Step2：设置 初始参量，初始 步长 0h ， fx在 x 处 的近 似 导数值 0 0 0 0002 8 8 212f x h f x h f x h f x hD h 令 10/10hh ，计算 1 1 1 11 12 8 8 212f x h f x h f x h f x hD h 并令 近似 导数值 序列 前 两个初值 0 1 1 00,E E D D ， 相对 误差 序列前 两个初值 0 1 1 1 00 , 2 /R R E D D 。 Step3： for
12、i=1 upto max 1） 计算 步长 1 /10iihh 时 的近似导数值， 2 8 8 212i i i iiif x h f x h f x h f x hD h ； 2） 近似 导数 值序列 间隔 1i i iE D D ， 相对误差序列 12/i i i iR E D D 。 3） 若 1iiEE 或 iR toler ， 则 停止 计算 循环； Step4： 输出 iD 即 可 得到 符合精度要求的近似导数值 ， 记录 ni 。 Step5： 将 , , , 0 ,1 ,k k kD E R k n 储存 到向量 ,DER 中输出 。 数值方法实验报告 7 3.P260 3 利
13、用 理查森外推 法设计计算 满足精度要求的近似导数值 的 算法 ， Step1： 输入函数 fx， 自变量 x ， 相对误差容限 toler， 误差容限 delta， 初始 步长 0h ， 计算行数上限 max ； Step2：计算 初始化 量 1） 初始 近似导数值 000 , 0 02f x h f x hD h ； 2） 初始 误差 0 0E ， 初始 相对 误差 0 1R ； 3） 1j ； Step3：计算出 符合精度要求的近似导数值： 若 1jE delta 且 1jR toler 且 1 maxj ，则进行 下 列 计算 ， 否则进入 step4 1） 1 /2jjhh , ,0
14、 2j j jjD f x h f x h h ; For k=1 upto j , 1 1 , 1, , 1 41j k j kj k j k kDDDD 2) 计算 误差 , 1, 1j j j j jE D D 3)计算相对 误差 , 1 , 12/ii j j j jR E D D 4) 1jj, 1nj返回 step3 Step4：输出 , , , nnn err relerr D， ,nnD为 符合要求的导数近似值。 4.P270 1 利用 牛顿多项式微分求 解各 插值点的近似的导数值 算法流程如下： Step1： 输入 ix 及 其对应的函数值 ifx ， 0,1, ,iN St
15、ep2： for k=0 upto N，计算 Mx 处 的近似导数值 NMPx 1） 0 Mtx ， 1 , 1 1 , , 1i i i it x i M t x M i N ， 2)令 01 ( ) , ( ) , , ( ) NA f t f t f t ， 数值方法实验报告 8 for j=1 upto N+1 for i=j+1 upto N+1 ( ) ( ( ) ( 1 ) ) / i i jA i A i A i t t 3） p=1， ( ) 2df k A ， for n=2 upto N 01( ) ( 1 )id f k d f k p t t A k Step3： 向
16、量 df 即 Mx 处 的近似导数值 NMPx， 0,1, ,MN （三） 实验结果 1.P260 （ 1）利用精度为 2()Oh 的 中心差分 公式 求解结果 （ 2）利用精度为 4()Oh 的 中心差分 公式 求解结果 （ 3）利用理查森 外推 法 求解结果 3 种 方法 所求的 结果如下 ，结果 的精度为 13 位 有效数字 表 1 三种方法求解各个函数在其相应 x 处近似导数值 对应函数及其求解点 （以 字母序号代替 ） 2()Oh 4()Oh 理查森 外推 法 （ a） 75.17350246175 75.17349470731 75.18107923463 （ b） 1.22860
17、2941028 1.228597423404 1.228733630836 （ c） 1.951562765952 1.951559609024 1.951752228326 （ d） 2.965529046763 2.965514829349 2.966366834548 （ e） 1.015201064578 1.015061058988 1.015008921184 由 上表可知，三种方法的求解结果 基本 相同， 但是 因为 函数及 求解点的设置的初始步长和 终止 条件的不同，三种方法的求解结果 没有 做到高度 一致， 因此可知求解结 果 与初始步长 及 终止 条件 的设定 有较大 的关
18、系。 五、实验结论 数值方法实验报告 9 本次 数值实验利用了精度 为 4()Oh 、 2()Oh 的中心差分公式及 理查森外推公式计算了各个函数在其 对应 点 的 近似导数值，三种方法的求解结果 基本 相同， 但是 因为 函数及 求解点的设置的初始步长和 终止 条件的不同，三种方法的求解结果没有 做到高度 一致， 因此可知求解结 果 与初始步长 及 终止 条件 的设定 有较大 的关系。 另外 还设计了计算 利用 牛顿插值多项式的任意插值点 导数 计算算法，并且用相应程序实现 ， 其中将 插值点 重新排列， 使得 所求点为第一个点是一个较大的创新，这样 使得 各个插值点 的 导数计算 更为快捷
19、 方便 。 附件（代码）： 1.P260 精度为 2()Oh 近似 导数值的计算子程序 function df=cetr2df(f,x,toler,max,h0) %输入 输出参数 说明 ： %df 输出 的 近似 导数值 %f、 x 相应函数及其点， toler 误差容忍上限， max 最大计算次数， h0 初始 步长 for i=1:2 h=10(-(i-1)*h0; D(i)=(feval(f,x+h)-feval(f,x-h)/(2*h);%计算导数近似值序列的前两个值 end E(1)=0;E(2)=abs(D(2)-D(1);%计算导数近似值间隔序列的前两个值 R(1)=1;R(2
20、)=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)/(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)E(i-1) end end n=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;
