数值微分文档格式.docx
《数值微分文档格式.docx》由会员分享，可在线阅读，更多相关《数值微分文档格式.docx（30页珍藏版）》请在冰点文库上搜索。
的近似值与精确值的绝对误差wuq，wuh和
的精确值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
function[x,yxj,wuc]=sandian3（h,xi,fi,M）
yxj=zeros（1,3）;
wuc=zeros（1,3）;
x2=xi
（2）;
y1=fi
（1）;
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）;
例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;
fi=[0.25000.20660.17360.1479];
[n2,x2,yx2,wuc2]=sandian（h,xi,fi,M）
yxj2=-2./（1+xi）.^3,wuyxj2=abs（（yxj2-yx2））
h=0.3;
fi=[0.25000.18900.1479];
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:
M=0.7502;
[x11,yxj11,wuc11]=sandian3（h1,xi,f11,M）,
xi=x（4:
6）;
f12=f（4:
[x12,yxj12,wuc12]=sandian3（h1,xi,f12,M）,
xi=x（5:
7）;
f13=f（5:
[x13,yxj13,wuc13]=sandian3（h1,xi,f13,M）,
h2=0.2,xi=x（1:
5）;
f21=f（1:
[x21,yxj21,wuc21]=sandian3（h2,xi,f21,M）,
xi=x（2:
f22=f（2:
[x22,yxj22,wuc22]=sandian3（h2,xi,f22,M）,
xi=x（3:
f23=f（3:
[x23,yxj23,wuc23]=sandian3（h2,xi,f23,M）,
h3=0.3,xi=x（1:
3:
f31=f（1:
[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）&
（j<
max1））
j;
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）;
jdW=abs（Dy（j+1,j+1）-Dy（j+1,j））;
j=j+1;
[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
4.4661809985950400
4.465504843773474.465504182820570
4.465501874697864.465501871795534.46550187123118
dy=jdw=n=wu=
4.465501871231185.643530087695581e-0107-1.863398324530863e-010
（二）精度为
的中心差商公式和误差估计及其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）;
x4=xi（4）;
x5=xi（5）;
y4=fi（4）;
y5=fi（5）;
yx=（8*y4-8*y2-y5+y1）/（12*h）;
wuc=abs（h.^4*M/30）;
function[x,y,yx,wuc]=zxcs5（fun,h,x0,M）
x=zeros（1,5）;
y=zeros（1,5）;
5
x（k）=x0+（k-3）.*h;
y（k）=feval（fun,x（k））;
x;
y;
yx=（8*y（4）-8*y
（2）-y（5）+y
（1））/（12*h）;
例8.2.6设
（1）分别根据（8.7）式和（8.24）式计算
的近似值，并估计误差，取小数点后14位计算，其中步长分别取
M=872,
解
（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=
-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程序，并输入此程序
M=872;
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=
dy=10*cos（5*x^2-21）*x,
wu2=abs（yc2-dy）,wu4=abs（yc4-dy）
的近似值与精确值的绝对误差wuc2，wuc4和
的精确值dy如下
4.46550187104484
wuc2=
0.748093207222270.007899009638480.000079032468570.00000079033719
wuc4=
0.159********6270.000017111040880.000000001712710.00000000001003
（三）变步长的中心差商公式及其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;
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<
max1））
h=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;
n=length（Dy（1:
k））-2;
Dy=Dy（2:
k）;
W=W（2:
H=H（2:
例8.2.7设
.取初始步长h0=0.2，精度wu=0.00001,用变步长的中心差商公式计算
的近似值和误差估计，与精确值4.46550187104484比较.
解输入程序
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*********
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））;
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+4555
升级会员 