復(fù)化梯形公式和復(fù)化Simpson公式

1、精選文檔數(shù)值計算方法上機(jī)題目3一、計算定積分的近似值：要求：（1） 若用復(fù)化梯形公式和復(fù)化Simpson公式計算，要求誤差限，分別利用他們的余項估量對每種算法做出步長的事前估量；（2） 分別利用復(fù)化梯形公式和復(fù)化Simpson公式計算定積分；（3） 將計算結(jié)果與精確解比較，并比較兩種算法的計算量。1.復(fù)化梯形公式程序：程序1（求f（x）的n階導(dǎo)數(shù)：syms xf=x*exp(x) %定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):') f2=diff(f,x,n) %求f(x)的n階導(dǎo)數(shù)結(jié)果1輸入n=2f2 =2*exp(x) + x*exp(x)程序2：clcclears

2、yms x %定義自變量xf=inline('x*exp(x)','x') %定義函數(shù)f(x)=x*exp(x)，換函數(shù)時只需換該函數(shù)表達(dá)式即可f2=inline('(2*exp(x) + x*exp(x)','x') %定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3='-(2*exp(x) + x*exp(x)' %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，以便求最大值e=5*10(-8) %精度要求值 a=1 %積分下限b=2 %積分上限x1=fminbnd(f3,1

3、,2) %求負(fù)的二階導(dǎo)數(shù)的最小值點，也就是求二階導(dǎo)數(shù)的最大值點對應(yīng)的x值for n=2:1000000 %求等分?jǐn)?shù)n Rn=-(b-a)/12*(b-a)/n)2*f2(x1) %計算余項 if abs(Rn)<e %用余項進(jìn)行推斷 break % 符合要求時結(jié)束 endendh=(b-a)/n %求hTn1=0 for k=1:n-1 %求連加和 xk=a+k*h Tn1=Tn1+f(xk)endTn=h/2*(f(a)+2*Tn1+f(b)z=exp(2)R=Tn-z %求已知值與計算值的差fprintf('用復(fù)化梯形算法計算的結(jié)果 Tn=')disp(Tn)fpri

4、ntf('等分?jǐn)?shù) n=')disp(n) %輸出等分?jǐn)?shù)fprintf('已知值與計算值的誤差 R=')disp(R)輸出結(jié)果顯示：用復(fù)化梯形算法計算的結(jié)果 Tn= 7.3891等分?jǐn)?shù) n=7019已知值與計算值的誤差 R= 2.8300e-0082. Simpson公式程序：程序1：（求f（x）的n階導(dǎo)數(shù)）：syms xf=x*exp(x) %定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):') f2=diff(f,x,n) %求f(x)的n階導(dǎo)數(shù)結(jié)果1輸入n=4f2 =4*exp(x) + x*exp(x)程序2：clcclearsyms

5、x %定義自變量xf=inline('x*exp(x)','x') %定義函數(shù)f(x)=x*exp(x)，換函數(shù)時只需換該函數(shù)表達(dá)式即可f2=inline('(4*exp(x) + x*exp(x)','x') %定義f(x)的四階導(dǎo)數(shù)，輸入程序1里求出的f2即可f3='-(4*exp(x) + x*exp(x)' %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，一邊求最大值e=5*10(-8) %精度要求值 a=1 %積分下限b=2 %積分上限x1=fminbnd(f3,1,2) %

6、求負(fù)的四階導(dǎo)數(shù)的最小值點，也就是求四階導(dǎo)數(shù)的最大值點對應(yīng)的x值for n=2:1000000 %求等分?jǐn)?shù)n Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1) %計算余項 if abs(Rn)<e %用余項進(jìn)行推斷 break % 符合要求時結(jié)束 endendh=(b-a)/n %求hSn1=0 Sn2=0for k=0:n-1 %求兩組連加和 xk=a+k*h xk1=xk+h/2 Sn1=Sn1+f(xk1) Sn2=Sn2+f(xk)end Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b) %因Sn2多加了k=0時的值，故減去f(a)z=exp

7、(2)R=Sn-z %求已知值與計算值的差fprintf('用Simpson公式計算的結(jié)果 Sn=')disp(Sn)fprintf('等分?jǐn)?shù) n=')disp(n) fprintf('已知值與計算值的誤差 R=')disp(R)輸出結(jié)果顯示：用Simpson公式計算的結(jié)果 Sn= 7.3891等分?jǐn)?shù) n=24已知值與計算值的誤差 R= 2.7284e-008用復(fù)化梯形公式計算的結(jié)果為：7.3891，與精確解的誤差為：2.8300e-008。等分?jǐn)?shù)n=7019用復(fù)化Simpson公式計算的結(jié)果為：7.3891，與精確解的誤差為：2.7284e-0

8、08。等分?jǐn)?shù)n=243、柯斯特公式求積分：程序代碼：（1）function y,Ck,Ak=NewtonCotes(fun,a,b,n) if nargin=1 mm,nn=size(fun); if mm>=8 error('為了保證NewtonCotes積分的穩(wěn)定性，最多只能有9個等距節(jié)點！') elseif nn=2 error('fun構(gòu)成應(yīng)為：第一列為x，其次列為y，并且個數(shù)為小于10的等距節(jié)點！') end xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin

9、=4 xk=linspace(a,b,n+1); if isa(fun,'function_handle') fx=fun(xk); else error('fun積分函數(shù)的句柄，且必需能夠接受矢量輸入！') end else error('輸入?yún)?shù)錯誤，請參考函數(shù)掛念！') end Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx' （2）function Ck=cotescoeff(n) for i=1:n+1 k=i-1; Ck(i)=(-1)(n-k)/factorial(k)/factorial(n-

10、k)/n*quadl(t)intfun(t,n,k),0,n); end （3）function f=intfun(t,n,k) f=1; for i=0:k-1,k+1:n f=f.*(t-i); end代碼解釋：function y,Ck,Ak=NewtonCotes(fun,a,b,n) % y=NewtonCotes(fun,a,b,n) % 牛頓-科特斯數(shù)值積分公式 % 參數(shù)說明： % fun，積分表達(dá)式，這里有兩種選擇 %(1)積分函數(shù)句柄，必需能夠接受矢量輸入，比如fun=(x)sin(x).*cos(x) % (2)x,y坐標(biāo)的離散點， 第一列為x， 其次列為y， 必需等距，

11、且節(jié)點的個數(shù)小于9， 比如： fun=1:8;sin(1:8)' % 假如fun的表接受其次種方式，那么只需要輸入第一個參數(shù)即可，否則還要輸入a,b,n三個參數(shù) % a，積分下限 % b，積分上限 % n，牛頓-科特斯數(shù)公式的階數(shù)，必需滿足1<n<7，由于n>=8時不能保證公式的穩(wěn)定性 % (1)n=1，即梯形公式 % (2)n=2，即辛普森公式% (3)n=4，即科特斯公式 % y，數(shù)值積分結(jié)果 % Ck，科特斯系數(shù) % Ak，求積系數(shù) % % Example % fun1=(x)sin(x);%必需可以接受矢量輸入 % fun2=0:0.1:0.5;sin(0:0

12、.1:0.5);%最多8個點，必需等距 % y1=NewtonCotes(fun1,0,0.5,6) % y2=NewtonCotes(fun2) if nargin=1 mm,nn=size(fun); if mm>=8 error('為了保證NewtonCotes積分的穩(wěn)定性，最多只能有9個等距節(jié)點！') elseif nn=2 error('fun構(gòu)成應(yīng)為：第一列為x，其次列為y，并且個數(shù)為小于10的等距節(jié)點！') end xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif n

13、argin=4 % 計算積分節(jié)點xk和節(jié)點函數(shù)值fx xk=linspace(a,b,n+1); if isa(fun,'function_handle') fx=fun(xk); else error('fun積分函數(shù)的句柄，且必需能夠接受矢量輸入！') end else error('輸入?yún)?shù)錯誤，請參考函數(shù)掛念！') end % 計算科特斯系數(shù) Ck=cotescoeff(n); % 計算求積系數(shù) Ak=(b-a)*Ck; % 求和算積分 y=Ak*fx' function Ck=cotescoeff(n) % 由于科特斯系數(shù)最多7階，為了便利我們可以直接使用，省得每次都計算 % A1=1,1/2% A2=1,4,1/6 % A3=1,3,3,1/8 % A4=7,32,12,32,1/90 % A5=19,75,50,50,75,19/288 % A6=41,216,27,272,27,216,41/840 % A7=751,3577,1323,2989,2989,1323,3577,751/17280 % 當(dāng)時為了體現(xiàn)公式，我們使用程序計算n階科特斯系數(shù) for i=1:n+1 k=i-1; Ck(i)=(-1)(n-k