第八章 数值微分

8.2 一阶导数的数值计算及其MATLAB 程序

8.2.1 差商求导及其MATLAB 程序

例 8.2.1 设)215sin()(2-=x x f .

（1）分别利用前差公式和后差公式计算)79.0('f 的近似值和误差，取4位小数点计

算，其中步长分别取1000.0,001.0,01

.0,1.0=h ，≤)("x f 80,]1,0[∈x . （2）将（1）中计算的)79.0('f 的近似值分别与精确值比较.

解 （1）编写计算)(x f y =的一阶导数计算)('x f 的近似值和误差估计的MA TLAB 程序，并输入

>> x=0.79;h=[0.1,0.01,0.001,0.0001];

M=80;x1=x+h;x2=x-h; y=sin(5.*x.^2-21);

y1=sin(5.*x1.^2-21); y2=sin(5.*x2.^2-21); yq=(y1-y)./h,

yh=(y-y2)./h,

wu=abs(h.*M/2), syms x,f=sin(5.*x.^2-21); yx=diff(f,x)

运行后屏幕显示利用前差公式和后差公式计算)79.0('f 的近似值yq ，yh 和误差估计wu ，

取4位小数点计算，其中步长分别取1000.0,001.0,01

.0,1.0=h ，M =80，导函数yx yq =

1.46596380397978 4.22848550173043 4.44250759584697 4.46320955293622

yh =

5.96885352366536 4.68672022108227 4.48833808130555 4.46779260847907

wu =

4.00000000000000 0.40000000000000 0.04000000000000 0.00400000000000

yx =

10*cos(5*x^2-21)*x

（2）计算)79.0(f '的值.输入程序

>> x=0.79; yx =10*cos(5*x^2-21)*x,

wuq=abs(yq-yx), wuh=abs(yh-yx)

运行后屏幕显示利用前差公式和后差公式计算)79.0(f '的近似值与精确值的绝对误差wuq ，wuh 和)79.0(f '的精确值yx 如下

yx =

4.46550187104484

wuq =

2.99953806706506 0.23701636931441 0.02299427519787 0.00229231810861

wuh =

1.50335165262053 0.22121835003744 0.02283621026072 0.00229073743424

8.2.2 中心差商公式求导及其MATLAB 程序

利用精度为)(2h O 的三点公式计算)(x f '的近似值和误差估计的MATLAB 主程序

function [n,xi,yx,wuc]=sandian(h,xi,fi,M)

n=length(fi); yx=zeros(1,n); wuc=zeros(1,n); x1= xi(1);

x2= xi(2); x3= xi(3);

y1=fi(1); y2=fi(2); y3=fi(3); xn= xi(n); xn1= xi(n-1);

xn2= xi(n-2); yn=fi(n); yn1=fi(n-1); yn2=fi(n-2);

for k=2:n-1

yx(n)=(yn2-4*yn1+3*yn)/(2*h);

yx(2)=( fi(3)- fi(1))/(2*h); yx(k)=( fi(k+1)-

fi(k-1))./(2*h);

wuc(1)=abs(h.^2.*M./3); wuc(n)=abs(h.^2.*M./3);

wuc(2:n-1)=abs(-h.^2.*M./6);

end

利用精度为)(2h O 的三点公式计算)('x f 的近似值和误差估计的MATLAB 主程序

function [x,yxj, wuc]=sandian3(h,xi,fi,M)

yxj=zeros(1,3); wuc=zeros(1,3); x1= xi(1);

x2= xi(2); x3= xi(3); y1=fi(1); y2=fi(2); y3=fi(3);

for t=1:3

s(t)=((2*t-5)*y1-4*(t-2)*y2+(2*t-3)*y3)/(2*h); x=xi;

y=s(t); yxj(t)=y;

if t==2

wuc(2)=abs(-h.^2*M/6);

else

wuc(1:2:3)=abs(h.^2*M/3);

end

end

例 8.2.3 设已给出)(x f y =的数据表8–5：

M (1)当h =0.1时，)(x f y =在x =1.000 0，1.100 0，1.200 0，1.300 0，1.400 0，1.500 0，1.600 0处的一阶导数的近似值，并估计误差；

(2)当h =0.2时，)(x f y =在x =1.000 0，1.200 0， 1.400 0，1.600 0处的一阶导数的近似值，并估计误差；

(3)当h =0.3时，)(x f y =在x =1.000 0，1.300 0 ，1.600 0处的一阶导数的近似值，并估计误差；

(4) 表8–5中的数据是函数2)1(1)(x x f +=

在相应点的数值，将(1)至(3)计算的一阶导数的近似值与2)

1(1)(x x f +=的一阶导数值比较，并求出它们的绝对误差.

解 （1）保存M 文件sandian.m ，sandian3.m ；

（2）在MATLAB 工作窗口输入如下程序

>> syms x,y=1/((1+x)^2); yx=diff(y,x,1),

yx3=diff(y,x,3),

运行后将屏幕显示的结果为

yx = yx3 =

-2/(1+x)^3 -24/(1+x)^5

（3）在MATLAB 工作窗口输入如下程序

>>h=0.1; xi=1.0000:h:1.6000;

fi=[0.2500 0.2268 0.2066 0.1890 0.1736 0.1600 0.1479];

x=1:0.001:1.6; yx3 =-24./(1+x).^5; M= max(abs(yx3));

[n1,x1,yx1,wuc1]=sandian(h,xi,fi,M)

yxj1=-2./(1+xi).^3,wuyxj1=abs(yxj1- yx1)

h=0.2; xi=1.0000:h:1.6000; fi=[0.2500 0.2066 0.1736 0.1479];

x=1:0.001:1.6; yx3 =-24./(1+x).^5; M= max(abs(yx3));