數(shù)值積分與數(shù)值微分

1、8 數(shù)值積分與數(shù)值微分8.1 例題解答例 8.1 給定積分,分別用梯形公式、公式、公式作近似計算.解:先輸入主要初始參數(shù)a=0.5;b=1;f=inline(x(1/2);%梯形公式I1=(b-a)/2*(feval(f,a)+feval(f,b)I1 = 0.426776695296637%simpson公式I2=(b-a)/6*(feval(f,a)+4*feval(f,(a+b)/2)+feval(f,b)I2 = 0.430934033027025%Cotes公式(n=4)tc=0;C0=7 32 12 32 7;for i=0:4 tc=tc+C0(i+1)*feval(f,a+i*

2、(b-a)/4);endI3=(b-a)/90*tc I3 = 0.430964070495876%準(zhǔn)確值I=int(char(f),a,b)vpa(I) I =-1/6*2(1/2)+2/3ans =0.43096440627115082519971854596505例 8.2 對積分,為使其精度達到.若用復(fù)化梯形公式,應(yīng)將0,1多少等分?若用復(fù)化公式,應(yīng)將0,1多少等分?解:直接按余項計算即可.復(fù)化梯形公式的余項為:復(fù)化公式余項為:對于,在課本中我們已證得以下不等式成立:直接利用上述不等式關(guān)系解答本題.先輸入誤差精度:Eps=1E-4 Eps = 1.000000000000000e-00

3、4(1) 復(fù)化梯形公式h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1) %先求出步長h1 = 0.060000000000000 N1=ceil(1/h1) %向上取整,得到等分區(qū)間數(shù)N1 = 17故可將區(qū)間17等分即可達到所要求的精度.(2) 復(fù)化公式h2=power(Eps/abs(-(1-0)/180*1/(1+4),1/4) %先求出步長h2 = 0.547722557505166 N2=ceil(1/h2) %向上取整,得到等分區(qū)間數(shù)N2 = 2 故可將區(qū)間2等分即可達到所要求的精度. 擴展:1) Matlab中復(fù)化梯形公式命令為I=trapz(x,y),復(fù)化公式命

4、令為quad().2) Matlab中有四個取整函數(shù),分別為ceil(),floor(),fix(),round(),分別表示向正無窮大方向取整、向負無窮大方向取整、向靠近零方向舍入和四舍五入.例 8.3 對積分,利用變步長方法求其近似值,使其精度達到.解：利用變步長法前先建立三種變步長復(fù)化積分公式的函數(shù). 注意在Matlab中直接用sin(0)/0得不到1，因此解此題時我們改用求極限的方法得到函數(shù)值,此函數(shù)名為limit().先建立三種復(fù)化公式的函數(shù)文件,它們分別為復(fù)化梯形公式trap.m、復(fù)化公式為simpson.m、公式為cotes.m,三個函數(shù)的源程序如下:(1) 復(fù)化梯形公式trap

5、.mfunction T=trap(f,a,b,n)%trap.m%復(fù)化梯形公式求積分值%f為積分函數(shù)%a,b為積分區(qū)間%n是等分區(qū)間份數(shù)h=(b-a)/n;%步長T=0;for k=1:(n-1) x0=a+h*k; T=T+limit(f,x0);endT=h*(limit(f,a)+limit(f,b)/2+h*T;T=double(T);(2) 復(fù)化公式simpson.m:function S=simpson(f,a,b,n)%simpson.m%Simpson公式求積分值%f為積分函數(shù)%a,b為積分區(qū)間%n是等分區(qū)間份數(shù)h=(b-a)/(2*n);%步長s1=0;s2=0;for k

6、=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0);endfor k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0);endS=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3;S=double(S);(3) 復(fù)化公式cotes.m:function C=cote(f,a,b,n)%cote.m%復(fù)化cotes公式求積分值%f為積分函數(shù)%a,b為積分區(qū)間%n是等分區(qū)間份數(shù)h=(b-a)/n;%步長C=0;for i=1:(n-1) x0=a+i*h; C=C+14*limit(f,x0);endfor k=0:(n

7、-1) x0=a+h*k; s=32*limit(f,x0+h*1/4)+12*limit(f,x0+h*1/2)+32*limit(f,x0+h*3/4); C=C+s;endC=C+7*(limit(f,a)+limit(f,b);C=C*h/90;C=double(C);再編寫主程序調(diào)用這三個函數(shù),主程序名為ex8_3.m,源程序如下:%ex8_3.mclc;syms x;f=sym(sin(x)/x);a=0;b=1;%積分上下限n=20;%作1,2,3,20次區(qū)間等分%復(fù)化梯形公式T=zeros(n,1);for i=1:n T(i)=trap(f,a,b,i);end%復(fù)化Simp

8、son公式;S=zeros(n,1);for i=1:n S(i)=simpson(f,a,b,i);end%復(fù)化Cotes公式C=zeros(n,1);for i=1:n C(i)=cote(f,a,b,i);end%準(zhǔn)確值I=int(f,a,b);I=double(I);%畫圖作出直觀觀察x=;x=1:n;figure;plot(x,ones(1,n)*I,-);hold on;plot(x,T,r-,LineWidth,2);plot(x,S,m.-,LineWidth,1);plot(x,C,c:,LineWidth,1.5);grid on;title(三種復(fù)化公式積分效果對比圖);

9、legend(準(zhǔn)確值曲線,復(fù)化梯形公式,復(fù)化Simpson公式,復(fù)化Cotes公式);hold off;disp( 復(fù)化梯形公式 復(fù)化Simpson公式 復(fù)化Cotes公式);disp(T S C);在Matlab命令窗口中輸入ex8_3,得到如下積分結(jié)果:復(fù)化梯形公式 復(fù)化Simpson公式 復(fù)化Cotes公式 0.920735492403948 0.946145882273587 0.946083004063674 0.939793284806177 0.946086933951794 0.946083069350917 0.943291429132337 0.94608383131169

10、9 0.946083070278278 0.944513521665390 0.946083310888472 0.946083070351379 0.945078780953402 0.946083168838073 0.946083070363043 0.945385730766859 0.946083117842867 0.946083070365797 0.945570776256246 0.946083095989403 0.946083070366633 0.945690863582701 0.946083085384948 0.946083070366936 0.94577318

11、8549752 0.946083079742053 0.946083070367061 0.945832071866905 0.946083076517732 0.946083070367118 0.945875637119198 0.946083074567937 0.946083070367147 0.945908771073865 0.946083073333114 0.946083070367161 0.945934556488475 0.946083072520476 0.946083070367170 0.945955016056114 0.946083071968057 0.94

12、6083070367174 0.945971521552985 0.946083071581964 0.946083070367177 0.945985029934386 0.946083071305562 0.946083070367179 0.945996225242376 0.946083071103489 0.946083070367180 0.946005606943978 0.946083070952998 0.946083070367181 0.946013546623966 0.946083070839065 0.946083070367182 0.94602032535502

13、5 0.946083070751532 0.946083070367182由以上結(jié)果可看到復(fù)化梯形公式有一個上升接近準(zhǔn)確值過程,而復(fù)化公式積分結(jié)果和復(fù)化公式積分的結(jié)果基本上和準(zhǔn)確值的曲線重疊在一塊,可見它們的精度是相當(dāng)高的.例 8.4 用積分法計算,精度.解：編寫積分法的函數(shù)M文件,源程序如下(romberg.m):function I,T=romberg(f,a,b,n,Eps)%Romberg積分計算%f為積分函數(shù)%a,b為積分區(qū)間%n+1是T數(shù)表的列數(shù)目%Eps為迭代精度%返回值中I為積分結(jié)果,T是積分表if narginEps) & (jn)| (jsyms x;%創(chuàng)建符號變量f=sy

14、m(sin(x)/x) %符號函數(shù)f =sin(x)/x I,T=romberg(f,0,1,3,1E-6) %積分計算 I = 0.9461T = 0.9207 0 0 0 0 0.9398 0.9461 0 0 0 0.9445 0.9461 0.9461 0 0 0.9457 0.9461 0.9461 0.9461 0 0.9460 0.9461 0.9461 0.9461 0.9461 其中T為積分表,由輸出結(jié)果可知計算結(jié)果為.例 8.5 利用積分法求.解:直接利用例8.4的函數(shù)即可.syms x;f=sym(4/(1+x2) f =4/(1+x2) I,T=romberg(f,0,

15、1,5,1E-6) I = 3.1416T = 3.0000 0 0 0 0 0 3.1000 3.1333 0 0 0 0 3.1312 3.1416 3.1421 0 0 0 3.1390 3.1416 3.1416 3.1416 0 0 3.1409 3.1416 3.1416 3.1416 3.1416 0 3.1414 3.1416 3.1416 3.1416 3.1416 3.1416 故結(jié)果為.例 8.6 對積分構(gòu)造其型求積公式.解:取,得到方程組:這是一個抽象方程組,可以利用Matlab的符號法來解之,該函數(shù)名為solve():%直接輸入解之:x=solve(A0+A1=2,A

16、0*x0+A1*x1=0,A0*x02+A1*x12=2/3,A0*x03+A1*x13=0,A0,A1,x0,x1) x = A0: 2x1 sym A1: 2x1 sym x0: 2x1 sym x1: 2x1 sym%顯示結(jié)果x.A0,x.A1,x.x0,x.x1 ans = 1 1ans = 1 1ans = 1/3*3(1/2) -1/3*3(1/2)ans = -1/3*3(1/2) 1/3*3(1/2)因此得到兩組解為:或: 求積公式為例 8.13 用中心差商公式計算在處的一階導(dǎo)數(shù).解：應(yīng)用中心差商公式可得到：編制如下程序計算：%ex8_13.mclc;syms x;%符號變量f

17、=sym(sqrt(x)%符號函數(shù)x0=2;%求值點h=2;%第一次初始步長N=10;%進行10次迭代計算dF=diff(f,x,1)%求導(dǎo)dFx=subs(dF,x0)%精確值df=zeros(N,6);for i=1:N df(i,1)=h; temp=(subs(f,x,x0+h/2)-subs(f,x,x0-h/2)/h;%按中心差商公式計算 df(i,2)=double(temp); h=h/10;%步長縮小10倍endtemp=dFx*ones(N,1)-df(:,2);df(:,3)=temp;h=1;%第二次取初始步長為1for i=1:N df(i,4)=h; temp=(s

18、ubs(f,x,x0+h/2)-subs(f,x,x0-h/2)/h; df(i,5)=double(temp); h=h/10;endtemp=dFx*ones(N,1)-df(:,5);df(:,6)=temp;disp( 步長h 近似值 誤差 步長h 近似值 誤差 );disp(df);運行結(jié)果如下：dFx = 0.353553390593274 步長h 近似值 誤差 步長h 近似值 誤差 Columns 1 through 4 2.000000000000000 0.366025403784439 -0.012472013191165 1.000000000000000 0.20000

19、0000000000 0.353663997049609 -0.000110606456335 0.100000000000000 0.020000000000000 0.353554495459696 -0.000001104866422 0.010000000000000 0.002000000000000 0.353553401641782 -0.000000011048508 0.001000000000000 0.000200000000000 0.353553390703976 -0.000000000110702 0.000100000000000 0.0000200000000

20、00 0.353553390597394 -0.000000000004120 0.000010000000000 0.000002000000000 0.353553390564087 0.000000000029186 0.000001000000000 0.000000200000000 0.353553389897954 0.000000000695320 0.000000100000000 0.000000020000000 0.353553397669515 -0.000000007076241 0.000000010000000 0.000000002000000 0.35355

21、3408771745 -0.000000018178471 0.000000001000000 Columns 5 through 6 0.356393958692601 -0.002840568099327 0.353581019507412 -0.000027628914138 0.353553666807604 -0.000000276214331 0.353553393355410 -0.000000002762136 0.353553390621819 -0.000000000028545 0.353553390586292 0.000000000006982 0.353553390564087 0.000000000029186 0.353553391008177 -0.000000000414903 0.353553386567285 0.000000004025989 0.353553408771745 -0.0000000181784718.2 Matlab數(shù)值積分函數(shù)介紹1. 求積分的符號運算intI=int(fun,v,a,b)fun為被積函數(shù)和符號表達式,可以為函數(shù)向量或函數(shù)矩陣,v為積分變量,積分函數(shù)中只有一個變量時可省略,a和b為積分上下限,即在a,b上計