8数值微分.docx

1、8数值微分8.2 一阶导数的数值计算及其MATLAB程序8.2.1 差商求导及其MATLAB程序例 8.2.1 设.（1）分别利用前差公式和后差公式计算的近似值和误差，取4位小数点计算，其中步长分别取，80,.（2）将（1）中计算的的近似值分别与精确值比较.解 （1）编写计算的一阶导数计算的近似值和误差估计的MATLAB程序，并输入 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).

2、/h,wu=abs(h.*M/2), syms x,f=sin(5.*x.2-21); yx=diff(f,x)运行后屏幕显示利用前差公式和后差公式计算的近似值yq，yh和误差估计wu，取4位小数点计算，其中步长分别取，M=80，导函数yxyq =1.46596380397978 4.22848550173043 4.44250759584697 4.46320955293622yh =5.96885352366536 4.68672022108227 4.48833808130555 4.46779260847907wu =4.00000000000000 0.40000000000000

3、0.04000000000000 0.00400000000000yx =10*cos(5*x2-21)*x（2）计算的值.输入程序 x=0.79; yx =10*cos(5*x2-21)*x,wuq=abs(yq-yx), wuh=abs(yh-yx)运行后屏幕显示利用前差公式和后差公式计算的近似值与精确值的绝对误差wuq，wuh和的精确值yx如下yx = 4.46550187104484wuq =2.99953806706506 0.23701636931441 0.02299427519787 0.00229231810861wuh =1.50335165262053 0.2212183

4、5003744 0.02283621026072 0.002290737434248.2.2 中心差商公式求导及其MATLAB程序利用精度为的三点公式计算的近似值和误差估计的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-

5、1); yn2=fi(n-2);for k=2:n-1yx(1)=(-3*y1+4*y2-y3)/(2*h); 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利用精度为的三点公式计算的近似值和误差估计的MATLAB主程序function x,yxj, wuc=sandian3(h,xi,fi,M)yxj=zeros

6、(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:3s(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=2wuc(2)=abs(-h.2*M/6);elsewuc(1:2:3)=abs(h.2*M/3);endend例 8.2.3 设已给出的数据表85：表85x1.000 0 1.100 0 1.200 0 1.300 0 1.400 0 1.500 0 1.600

7、0 f(x)0.250 0 0.226 8 0.206 6 0.189 0 0.173 6 0.160 0 0.147 9 M= 0.750 2,试用三点公式计算下列问题：(1)当h=0.1时，在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=1.000 0，1.200 0， 1.400 0，1.600 0处的一阶导数的近似值，并估计误差；(3)当h=0.3时，在x=1.000 0，1.300 0 ，1.600 0处的一阶导数的近似值，并估计误差；(4) 表85中的数

8、据是函数在相应点的数值，将(1)至(3)计算的一阶导数的近似值与的一阶导数值比较，并求出它们的绝对误差.解 （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

9、.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);n2,x2,yx2,wuc2=sandian(h,xi,fi,M)yxj2=-2./(1+xi).3,wuyxj2=abs(yxj2- yx2)h=0

10、.3; xi=1.0000:h:1.6000; fi=0.2500 0.1890 0.1479;x=1:0.001:1.6; yx3 =-24./(1+x).5; M=max(abs(yx3);x3,yx3, wuc3=sandian3(h,xi,fi,M)yxj3=-2./(1+xi).3,wuyxj3=abs(yxj3- yx3)或 h1=0.1,x=1.0000,1.1000,1.2000,1.3000,1.4000,1.5000,1.6000;f=0.2500,0.2268,0.2066,0.1890,0.1736,0.1600,0.1479;xi=x(1:3);f11=f(1:3);

11、 M= 0.7502;x11,yxj11,wuc11=sandian3(h1,xi,f11,M), xi= x(4:6);f12=f(4:6); x12,yxj12,wuc12=sandian3(h1,xi,f12,M), xi=x(5:7);f13=f(5:7); x13,yxj13,wuc13=sandian3(h1,xi,f13,M), h2=0.2, xi= x(1:2:5);f21= f(1:2:5);x21,yxj21,wuc21=sandian3(h2,xi,f21,M),xi= x(2:2:6);f22=f(2:2:6); x22,yxj22,wuc22=sandian3(h2

12、,xi,f22,M),xi= x(3:2:7);f23=f(3:2:7); x23,yxj23,wuc23=sandian3(h2,xi,f23,M), h3=0.3, xi= x(1:3:7);f31= f(1:3:7); x31,yxj31,wuc31=sandian3(h3,xi,f31,M),将运行的结果（略）.8.2.3 理查森外推法求导及其MATLAB程序（一）一般形式的理查森外推法及其MATLAB程序利用理查森外推法计算的近似值和误差估计的MATLAB程序function Dy,dy,jdw,n=diffext1(fun,x0,jdwc,max1)h=1;j=1; n=1;jdW

13、=1;xdW=1; x1=x0+h;x2=x0-h;Dy(1,1)=(feval(fun,x1)- feval(fun,x2)/(2*h); while(jdWjdwc)&(j x0=0.79;jdwc=0.0000001,max1=100;Dy,dy,jdw,n=diffext1(fun,x0,jdwc,max1),wu=4.46550187104484-dy运行后屏幕显示的近似值dy，dy与精确值的绝对误差wu，导数近似值的迭代矩阵Dy，jdW=|Dy(n,n)-Dy(n,n-1)|的值,最佳近似值dy的坐标n如下Dy = Columns 1 through 4 0.95036708207

14、779 0 0 00.87447447334140 0.84917693709594 0 01.04543344913993 1.10241977440611 1.11930263022679 03.33291848491782 4.09541349684379 4.29494641167297 4.345353455822914.16327772677060 4.44006414072153 4.46304085031338 4.465709016006084.38872901334021 4.46387944219674 4.46546712896176 4.465505641321264

15、.44623223629059 4.46539997727405 4.46550134627921 4.46550188941123 Columns 5 through 7 0 0 0 0 0 0 0 0 0 0 0 0 4.46618099859504 0 0 4.46550484377347 4.46550418282057 0 4.46550187469786 4.46550187179553 4.46550187123118dy = jdw = n= wu = 4.46550187123118 5.643530087695581e-010 7 -1.863398324530863e-0

16、10（二）精度为的中心差商公式和误差估计及其MATLAB程序用精度为的中心差商公式计算的近似值和误差估计的MATLAB程序function x0,yx,wuc=zxcs4(h,x0,fi,M)xi=x0-2*h,x0-h,x0,x0+h,x0+2*h;x1= xi(1); x2= xi(2); x3= xi(3); x4= xi(4);x5= xi(5);y1=fi(1); y2=fi(2); y3=fi(3);y4=fi(4); y5=fi(5);yx=(8*y4-8*y2-y5+y1)/(12*h); wuc=abs(h.4*M/30);用精度为的中心差商公式计算的近似值和误差估计的MAT

17、LAB程序function x,y,yx,wuc=zxcs5(fun,h,x0, M)x=zeros(1,5);y= zeros(1,5);for k=1:5x(k)=x0+(k-3).*h;y(k)= feval(fun, x(k);endx;y;yx=(8*y(4)-8*y(2)-y(5)+y(1)/(12*h); wuc=abs(h.4*M/30);例 8.2.6 设.（1）分别根据(8.7) 式和(8.24)式计算的近似值，并估计误差，取小数点后14位计算，其中步长分别取，M=872,.（2）将（1）中计算的的近似值分别与精确值比较.解 （1）计算的一阶导数的近似值和误差估计.方法1

18、输入程序x0=0.79;h1=0.1,M=872;x,y,yx1,wuc1=zxcs5(fun,h1,x0,M)h2=0.01, x,y,yx2,wuc2=zxcs5(fun,h2,x0, M), h3=0.001, x,y,yx3,wuc3=zxcs5(fun,h3,x0, M), h4=0.0001, x,y,yx4,wuc4=zxcs5(fun,h4,x0, M),运行后屏幕显示x=(-2h, -h, ,+h, +2h), 在x处函数值向量y，步长分别取时，根据 (8.24)式计算函数在处的导数的近似值yx1, yx2, yx3, yx4和误差估计wuc1，wuc2，wuc3，wuc4如

19、下x = Columns 1 through 3 0.59000000000000 0.69000000000000 0.79000000000000 Columns 4 through 5 0.89000000000000 0.99000000000000y = Columns 1 through 3 -0.39855804055354 0.22803197210174 0.82491732446828 Columns 4 through Column 5 0.97151370486625 0.38160930304430yx1 = wuc1 = 4.30640543209856 0.002

20、90666666667yx2 = wuc2 = 4.46548476000396 2.906666666666667e-007yx3 = wuc3 = 4.46550186933213 2.906666666666667e-011yx4 = wuc4 = 4.46550187103480 2.906666666666667e-015方法2编写计算的一阶导数的近似值和误差估计的MATLAB程序，并输入此程序 x=0.79;h=0.1,0.01,0.001,0.0001; M=872;x1=x+h;x2=x-h;x3=x+2*h;x4=x-2*h; y1=sin(5.*x1.2-21);y2=si

21、n(5.*x2.2-21);y3=sin(5.*x3.2-21); y4=sin(5.*x4.2-21); yc2=(y1-y2)./(2*h),wuc2=abs(h.2*M/6),yc4=(8*y1-8*y2-y3+y4)./ (12*h), wuc4=abs(h.4*M/30),syms x,f=sin(5.*x.2-21); dy=diff(f,x)运行后屏幕显示步长分别取，M= 872，分别根据(8.7) 式和(8.24)式计算的近似值yc2，yc4及其误差估计值wuc2，wuc4，导函数dy如下yc2 =3.71740866382257 4.45760286140635 4.4654

22、2283857626 4.46550108070765wuc2 =1.45333333333333 0.01453333333333 0.00014533333333 0.00000145333333yc4 =4.30640543209856 4.46548476000396 4.46550186933213 4.46550187103480wuc4 =0.00290666666667 0.00000029066667 0.00000000002907 0.00000000000000dy =10*cos(5*x2-21)*x（2）计算的值.输入程序 x=0.79; dy =10*cos(5*

23、x2-21)*x, wu2=abs(yc2-dy), wu4=abs(yc4-dy)运行后屏幕显示的近似值与精确值的绝对误差wuc2，wuc4和的精确值dy如下dy =4.46550187104484wuc2=0.74809320722227 0.00789900963848 0.00007903246857 0.00000079033719wuc4=0.159*627 0.00001711104088 0.00000000171271 0.00000000001003（三）变步长的中心差商公式及其MATLAB程序用变步长的中心差商公式计算的近似值和误差估计的MATLAB程序function

24、n,H,Dy,W=difflim(fun,x0,h0,wu,max1)Dy=zeros(1,max1); W=zeros(1, max1); H=zeros(1, max1);h=h0;H(1)=h; E(1)=0; x1=x0+h;x2=x0-h; h1=h/10;x3=x0+h1;x4=x0-h1;Dy(1)=(feval(fun,x1)- feval(fun,x2)./(2*h);Dy(2)=(feval(fun,x3)- feval(fun,x4)./(2*h1);W(1)=abs(Dy(2)- Dy(1); k=1; while( (W(k)wu)&(k x0=0.79;wu=0.0

25、0001;max1=100;h0=0.2;n,H,Dy,W= difflim(fun,x0,h0,wu,max1)jdwc=4.46550187104484- Dy运行后屏幕显示迭代次数n, 变步长数组H，用变步长的中心差商公式计算的近似值Dy，误差估计值W，jdwc，精确值与近似值Dy的差如下n = 2H =.020*Dy = 4.43395716561354 4.46549870972618 4.46550187410688W = 2.48353880661895 0.03154154411264 0.00000316438070jdwc =.0315*8.2.4 牛顿（Newton）多项

26、式求导及其MATLAB程序利用牛顿插值多项式求导的MATLAB主程序function df,A,P=diffnew(X,Y)n=length(X);A=Y;for j=2:n for i=n:-1:j A(i)=(A(i)- A(i-1)/(X(i)-X(i-j+1); end endx0=X(1);df=A(2);chsh=1;m=length(A)-1; for k=2:m chsh=chsh*(x0-X(k); df=df+chsh*(A(k+1);end P=poly2sym(A);例8.2.9 根据下表给定的一组数据（X ,Y）写出 (8.30)式的具体形式及其精度和名称，并用它计算

27、的近似值，取4位小数计算.X3.135 2 3.335 2 3.535 2 3.735 2 3.935 2 4.135 2 4.335 2 4.535 2 Y0.126 6 -0.060 2 -0.603 2 -0.998 0 -0.119 4 0.995 3 -0.654 2 0.158 1 解 因为表中所给数据（X ,Y）是等差数列，公差为h=0.2, 即x0=3.135 2, x1= x0+h, x2= x0+2h, x3= x0+3h, x4= x0+4h, x5= x0+5h, x6= x0+6h, x7= x0+7h.输入程序X=3.1352,3.3352,3.5352,3.735

28、2,3.9352,4.1352,4.3352,4.5352 ;Y=0.1266,-0.0602,-0.6032,-0.9980,-0.1194,0.9953,-0.6542 0.1581;df,A,P=diffnew(X,Y)运行后屏幕显示的近似值df和阶牛顿多项式P及其系数向量A如下df = -0.2428A =0.1266 -0.9340 -4.4525 10.5083 16.1667 -72.4818 64.7309 108.6155P=633/5000*x7-467/500*x6-1781/400*x5+1478916440133901/140737488355328*x4+4550512123488957/281474976710656*x3-27833/384*x2+4