復(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：clcclearsyms x %定義自變量xf

2、=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,2) %求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn)，也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對應(yīng)的x值for n=2:1000000 %求等分?jǐn)?shù)n rn=-(b

3、-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)fprintf(等分?jǐn)?shù) n=)disp(n) %輸出等分?jǐn)?shù)fprintf(已知值與計算值的誤差 r=)disp(r)輸出結(jié)果顯示：用復(fù)化梯形算法計算的結(jié)果 tn= 7

4、.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 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ù)

5、，輸入程序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) %求負(fù)的四階導(dǎo)數(shù)的最小值點(diǎn)，也就是求四階導(dǎo)數(shù)的最大值點(diǎn)對應(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)=8 error(為了保證newtoncotes積分的穩(wěn)定性，最多只能有9個等距節(jié)點(diǎn)！) elseif nn=2

6、 error(fun構(gòu)成應(yīng)為：第一列為x，第二列為y，并且個數(shù)為小于10的等距節(jié)點(diǎn)！) end xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin=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）fun

7、ction ck=cotescoeff(n) for i=1:n+1 k=i-1; ck(i)=(-1)(n-k)/factorial(k)/factorial(n-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á)式，這里有兩種選擇

8、 %(1)積分函數(shù)句柄，必須能夠接受矢量輸入，比如fun=(x)sin(x).*cos(x) % (2)x,y坐標(biāo)的離散點(diǎn)， 第一列為x， 第二列為y， 必須等距， 且節(jié)點(diǎn)的個數(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，即科特斯公式 % y，數(shù)值積分結(jié)果 % ck，科特斯系數(shù) % ak，求積系數(shù) % %

9、example % fun1=(x)sin(x);%必須可以接受矢量輸入 % fun2=0:0.1:0.5;sin(0:0.1:0.5);%最多8個點(diǎn)，必須等距 % 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é)點(diǎn)！) elseif nn=2 error(fun構(gòu)成應(yīng)為：第一列為x，第二列為y，并且個數(shù)為小于10的等距節(jié)點(diǎn)！) end xk=fun(1,:); fk=fun(2,:); a

10、=min(xk); b=max(xk); n=mm-1; elseif nargin=4 % 計算積分節(jié)點(diǎn)xk和節(jié)點(diǎn)函數(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)/factorial(k)/factorial(n-k)/n*quadl(t)int