8数值微分.docx
《8数值微分.docx》由会员分享，可在线阅读，更多相关《8数值微分.docx（30页珍藏版）》请在冰豆网上搜索。
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）./h,
wu=abs（h.*M/2）,symsx,f=sin（5.*x.^2-21）;yx=diff（f,x）
运行后屏幕显示利用前差公式和后差公式计算
的近似值yq，yh和误差估计wu，取4位小数点计算，其中步长分别取
，M=80，导函数yx
yq=
1.465963803979784.228485501730434.442507595846974.46320955293622
yh=
5.968853523665364.686720221082274.488338081305554.46779260847907
wu=
4.000000000000000.400000000000000.040000000000000.00400000000000
yx=
10*cos（5*x^2-21）*x
（2）计算
的值.输入程序
>>x=0.79;yx=10*cos（5*x^2-21）*x,
wuq=abs（yq-yx）,wuh=abs（yh-yx）
运行后屏幕显示利用前差公式和后差公式计算
的近似值与精确值的绝对误差wuq，wuh和
的精确值yx如下
yx=
4.46550187104484
wuq=
2.999538067065060.237016369314410.022994275197870.00229231810861
wuh=
1.503351652620530.221218350037440.022836210260720.00229073743424
8.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-1）;yn2=fi（n-2）;
fork=2:
n-1
yx
（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（1,3）;wuc=zeros（1,3）;x1=xi
（1）;
x2=xi
（2）;x3=xi（3）;y1=fi
（1）;y2=fi
（2）;y3=fi（3）;
fort=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;
ift==2
wuc
（2）=abs（-h.^2*M/6）;
else
wuc（1:
2:
3）=abs（h.^2*M/3）;
end
end
例8.2.3设已给出
的数据表8–5：
表8–5
x
1.00001.10001.20001.30001.40001.50001.6000
f（x）
0.25000.22680.20660.18900.17360.16000.1479
M=0.7502,试用三点公式计算下列问题：
（1）当h=0.1时，
在x=1.0000，1.1000，1.2000，1.3000，1.4000，1.5000，1.6000处的一阶导数的近似值，并估计误差；
（2）当h=0.2时，
在x=1.0000，1.2000，1.4000，1.6000处的一阶导数的近似值，并估计误差；
（3）当h=0.3时，
在x=1.0000，1.3000，1.6000处的一阶导数的近似值，并估计误差；
（4）表8–5中的数据是函数
在相应点的数值，将
（1）至（3）计算的一阶导数的近似值与
的一阶导数值比较，并求出它们的绝对误差.
解
（1）保存M文件sandian.m，sandian3.m；
（2）在MATLAB工作窗口输入如下程序
>>symsx,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.25000.22680.20660.18900.17360.16000.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.25000.20660.17360.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.3;xi=1.0000:
h:
1.6000;fi=[0.25000.18900.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）;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,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=1;xdW=1;x1=x0+h;x2=x0-h;
Dy（1,1）=（feval（fun,x1）-feval（fun,x2））/（2*h）;
while（（jdW>jdwc）&（jj;x1=x0+2^（-j）*h;x2=x0-2^（-j）*h;
Dy（j+1,1）=（feval（fun,x1）-feval（fun,x2））/（2^（1-j）*h）;
fork=1:
j
k;Dy（j+1,k+1）=Dy（j+1,k）+（Dy（j+1,k）-Dy（j,k））/（4^k-1）;
end
jdW=abs（Dy（j+1,j+1）-Dy（j+1,j））;j=j+1;
end
[n,n]=size（Dy）;jdw=abs（Dy（n,n）-Dy（n,n-1））;
dy=Dy（n,n）;
例8.2.5设
，取精度jdwc=0.0000001,用理查森外推法计算
的近似值和误差估计，将近似值与精确值4.46550187104484比较.
解取最大计算次数max1=100，输入程序
>>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=
Columns1through4
0.95036708207779000
0.874474473341400.8491769370959400
1.045433449139931.102419774406111.119302630226790
3.332918484917824.095413496843794.294946411672974.34535345582291
4.163277726770604.440064140721534.463040850313384.46570901600608
4.388729013340214.463879442196744.465467128961764.46550564132126
4.446232236290594.465399977274054.465501346279214.46550188941123
Columns5through7
000
000
000
000
4.4661809985950400
4.465504843773474.465504182820570
4.465501874697864.465501871795534.46550187123118
dy=jdw=n=wu=
4.465501871231185.643530087695581e-0107-1.863398324530863e-010
（二）精度为
的中心差商公式和误差估计及其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）;
用精度为
的中心差商公式计算
的近似值和误差估计的MATLAB程序
function[x,y,yx,wuc]=zxcs5（fun,h,x0,M）
x=zeros（1,5）;y=zeros（1,5）;
fork=1:
5
x（k）=x0+（k-3）.*h;
y（k）=feval（fun,x（k））;
end
x;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输入程序
>>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如下
x=
Columns1through3
0.590000000000000.690000000000000.79000000000000
Columns4through5
0.890000000000000.99000000000000
y=
Columns1through3
-0.398558040553540.228031972101740.82491732446828
Columns4throughColumn5
0.971513704866250.38160930304430
yx1=wuc1=
4.306405432098560.00290666666667
yx2=wuc2=
4.465484760003962.906666666666667e-007
yx3=wuc3=
4.465501869332132.906666666666667e-011
yx4=wuc4=
4.465501871034802.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=sin（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）,symsx,f=sin（5.*x.^2-21）;dy=diff（f,x）
运行后屏幕显示步长分别取
，M=872，分别根据（8.7）式和（8.24）式计算
的近似值yc2，yc4及其误差估计值wuc2，wuc4，导函数dy如下
yc2=
3.717408663822574.457602861406354.465422838576264.46550108070765
wuc2=
1.453333333333330.014533333333330.000145333333330.00000145333333
yc4=
4.306405432098564.465484760003964.465501869332134.46550187103480
wuc4=
0.002906666666670.000000290666670.000000000029070.00000000000000
dy=
10*cos（5*x^2-21）*x
（2）计算
的值.输入程序
>>x=0.79;dy=10*cos（5*x^2-21）*x,
wu2=abs（yc2-dy）,wu4=abs（yc4-dy）
运行后屏幕显示
的近似值与精确值的绝对误差wuc2，wuc4和
的精确值dy如下
dy=
4.46550187104484
wuc2=
0.748093207222270.007899009638480.000079032468570.00000079033719
wuc4=
0.159********6270.000017111040880.000000001712710.00000000001003
（三）变步长的中心差商公式及其MATLAB程序
用变步长的中心差商公式计算
的近似值和误差估计的MATLAB程序
function[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）&（kh=h/（10^（k））;H（k+1）=h;x1=x0+H（k+1）;x2=x0-H（k+1）;
Dy（k+1）=（feval（fun,x1）-feval（fun,x2））./（2*H（k+1））;
W（k+1）=abs（Dy（k+1）-Dy（k））;
k=k+1;
end
n=length（Dy（1:
k））-2;
Dy=Dy（2:
k）;W=W（2:
k）;H=H（2:
k）;
例8.2.7设
.取初始步长h0=0.2，精度wu=0.00001,用变步长的中心差商公式计算
的近似值和误差估计，与精确值4.46550187104484比较.
解输入程序
>>x0=0.79;wu=0.00001;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=
2
H=
.020*********
Dy=
4.433957165613544.465498709726184.46550187410688
W=
2.483538806618950.031541544112640.00000316438070
jdwc=
.0315********
8.2.4牛顿（Newton）多项式求导及其MATLAB程序
利用牛顿插值多项式求导的MATLAB主程序
function[df,A,P]=diffnew（X,Y）
n=length（X）;
A=Y;
forj=2:
n
fori=n:
-1:
j
A（i）=（A（i）-A（i-1））/（X（i）-X（i-j+1））;
end
end
x0=X
（1）;df=A
（2）;chsh=1;m=length（A）-1;
fork=2:
m
chsh=chsh*（x0-X（k））;df=df+chsh*（A（k+1））;
end
P=poly2sym（A）;
例8.2.9根据下表给定的一组数据（X,Y）写出（8.30）式的具体形式及其精度和名称，并用它计算
的近似值，取4位小数计算.
X
3.13523.33523.53523.73523.93524.13524.33524.5352
Y
0.1266-0.0602-0.6032-0.9980-0.11940.9953-0.65420.1581
解因为表中所给数据（X,Y）是等差数列，公差为h=0.2,即
x0=3.1352,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.7352,3.9352,
4.1352,4.3352,4.5352];
Y=[0.1266,-0.0602,-0.6032,-0.9980,-0.1194,
0.9953,-0.65420.1581];
[df,A,P]=diffnew（X,Y）
运行后屏幕显示
的近似值df和
阶牛顿多项式P及其系数向量A如下
df=-0.2428
A=
0.1266-0.9340-4.452510.508316.1667-72.481864.7309108.6155
P=
633/5000*x^7-467/500*x^6-1781/400*x^5+1478916440133901/140737488355328*x^4+4550512123488957/281474976710656*x^3-27833/384*x^2+4
升级会员 