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

1、數(shù)值計算方法上機(jī)題目 3一、計算定積分的近似值：22xe二xe dxi要求：(1)若用復(fù)化梯形公式和復(fù)化 Simpson公式計算，要求誤差限10，分別利用他們的余項估計對每種算法做出步長的事前估2計；(2) 分別利用復(fù)化梯形公式和復(fù)化 Simps on 公式計算定積分；(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、序 2 :clcclearsyms 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 即可。%因 fminbnd ()函數(shù)求的是表達(dá)式的最小值，且要求表 達(dá)式帶引號，故取負(fù)號，以便求最大值e=5*10A(-8)%精度要求值a=1%積分下限b=2%積分上限x1=fmi nbn d(f3,1,2)%求負(fù)的二階導(dǎo)數(shù)的最小值點，也就是求二階導(dǎo)數(shù)的最大值點對應(yīng)的 x 值f3= -(2*exp(x

3、) + x*exp(x)for n=2:10000009 求等分?jǐn)?shù) nTn 1=T n1+f(xk)endTn=h/2*(f(a)+2*T n1+f(b)z=exp(2)R=T n-z%求已知值與計算值的差fprintf(用復(fù)化梯形算法計算的結(jié)果Tn=)disp(T n)fprintf(等分?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-008Rn=-(b-a)/12*(b-a)/n)A2*f2(x1) if abs(Rn)ve

4、 breakendendh=(b-a)/nTn 仁 0for k=1:n-1%求連加和xk=a+k*h%計算余項%用余項進(jìn)行判斷%符合要求時結(jié)束%求 h2. Simps on 公式程序：程序 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=4h=(b-a)/n%求 hf2 = 4*exp(x) + x*exp(x)程序 2 :clcclearsyms x%定義自變量 xf=inline(x*exp(x) , x ) %定義函數(shù) f(x)=

5、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*10A(-8)a=1%積分下限b=2%積分上限x1=fmi nbn d(f3,1,2)最大值點對應(yīng)的 x 值for n=2:1000000endendSn 仁 0Sn2=0%精度要求值%求負(fù)的四階導(dǎo)數(shù)的最小值點，也就是求四階導(dǎo)數(shù)的Rn=-(b-a)/180*(b-a

6、)/(2* n)F4*f2(x1)%計算余項if abs(Rn)ve%用余項進(jìn)行判斷break%符合要求時結(jié)束%求等分?jǐn)?shù) nfor k=0:n-1%求兩組連加和xk=a+k*hxk 仁 xk+h/2Sn 1=S n1+f(xk1)Sn 2=S n2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%因 Sn2 多加了 k=0 時的值，故減去 f(a)z=exp(2)R=S n-z%求已知值與計算值的差fprintf(用 Simpson 公式計算的結(jié)果 Sn=)disp(S n)fprintf(等分?jǐn)?shù) n=)disp (n)fprintf(已知值與計算值的誤

7、差 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-008。等分?jǐn)?shù) n=243、柯斯特公式求積分:程序代碼：(1) function y,Ck,Ak二NewtonCotes(fun,a,b,n)if n arg in=1mm, nn =size(fu n);if mm=8error(為了保證 Newt

8、onCotes 積分的穩(wěn)定性，最多只能有 9個等距節(jié)點！)elseif nn=2error(fun 構(gòu)成應(yīng)為：第一列為 x，第二列為 y,并且個數(shù)為小于 10 的等距節(jié)點！)endxk=fu n(1,:);fk=fu n( 2,:);a=mi n(xk);b=max(xk);n=m m-1;elseif nargin=4xk=li nspace(a,b, n+1);if isa(fu n,fun cti on_han dle)fx=fu n(xk);elseerror（fun 積分函數(shù)的句柄，且必須能夠接受矢量輸入！）endelseerror(輸入?yún)?shù)錯誤，請參考函數(shù)幫助！)endCk=cot

9、escoeff( n);Ak=(b-a)*Ck;y=Ak*fx：(2) function Ck=cotescoeff(n)for i=1: n+1k=i-1;Ck(i)=(-1)A( n-k)/factorial(k)/factorial( n-k)/n*quadl(t)i ntfun (t, n,k),0,n);end(3) function f二intfun(t,n,k)f=1;for i=0:k-1,k+1: nf=f.*(t-i);end代碼解釋:fun ctio n y,Ck,Ak二Newto nCotes(fu n,a,b ,n) % y=Newt on Cotes(f un,a,

10、b ,n)%牛頓-科特斯數(shù)值積分公式%參數(shù)說明：% fun,積分表達(dá)式，這里有兩種選擇%(1)積分函數(shù)句柄，必須能夠接受矢量輸入，比如fun二(x)si n( x).*cos(x)% (2)x,y 坐標(biāo)的離散點，第一列為 x，第二列為 y，必須等距，且節(jié)點的個數(shù)小于 9， 比如：fun二1:8;sin(1:8)%如果 fun 的表采用第二種方式，那么只需要輸入第一個參數(shù)即可， 否則還要輸入 a,b,n 三個參數(shù)%a，積分下限% b，積分上限% n，牛頓-科特斯數(shù)公式的階數(shù)，必須滿足 1n=8 時不能保證公式的穩(wěn)定性%(1) n=1，即梯形公式%(2) n=2，即辛普森公式%(3) n=4，即科

11、特斯公式%y，數(shù)值積分結(jié)果% Ck，科特斯系數(shù)% Ak，求積系數(shù)% Example% fun 仁(x)sin(x);%必須可以接受矢量輸入% fun2 二0:0.1:0.5;sin(0:0.1:0.5);% 最多 8 個點，必須等距% y 仁 N ewto nCotes(fu n1,0,0.5,6)% y2=Newto nCotes(fu n2)if n arg in=1mm, nn =size(f un);if mm=8error(為了保證 NewtonCotes 積分的穩(wěn)定性，最多只能有 9 個等距節(jié)點！)elseif nn=2error(fun 構(gòu)成應(yīng)為：第一列為 x,第二列為 y，并且

12、個數(shù)為小于 10 的等距節(jié)點！)endxk=fu n(1,:);fk=fu n( 2,:);a=mi n(xk);b=max(xk);n=m m-1;elseif nargin=4%計算積分節(jié)點 xk 和節(jié)點函數(shù)值 fx xk=li nspace(a,b ,n+1);if isa(fu n,fu ncti on_han die)fx=fu n(xk);elseerror(fun 積分函數(shù)的句柄，且必須能夠接受矢量輸入！)endelseerror(輸入?yún)?shù)錯誤，請參考函數(shù)幫助！)end%計算科特斯系數(shù)Ck=cotescoeff( n);%計算求積系數(shù)Ak=(b-a)*Ck;%求和算積分y=Ak*

13、fx;fun cti on Ck=cotescoeff( n)%由于科特斯系數(shù)最多 7 階，為了方便我們可以直接使用，省 得每次都計算% A 仁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+1k=i-1;Ck(i)=(-1F( n-k)/factorial(k)/f