第八章 数值微分.docx
- 文档编号:8890107
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:25
- 大小:126.83KB
第八章 数值微分.docx
《第八章 数值微分.docx》由会员分享,可在线阅读,更多相关《第八章 数值微分.docx(25页珍藏版)》请在冰豆网上搜索。
第八章数值微分
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)&(j 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); 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)&(k 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; 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第八章 数值微分 第八 数值 微分