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

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

2、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)'%Ufminbnd()函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，以便求最大值e=5*10A(-8)%精度要求值a=1%積分下限b=2%R分上限x1=fminbnd(f3,1,2)%求負(fù)的二階導(dǎo)數(shù)的最小值

3、點，也就是求二階導(dǎo)數(shù)的最大值點對應(yīng)的x值forn=2:1000000就等分?jǐn)?shù)nRn=-(b-a)/12*(b-a)/n)A2*f2(x1)%+算余項ifabs(Rn)<e%ffl余項進(jìn)行判斷break%符合要求時結(jié)束endendh=(b-a)/n%<hTn1=0fork=1:n-1%連加和xk=a+k*hTn1=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ù)fp

4、rintf('已知值與計算值的誤差R=')disp(R)輸出結(jié)果顯示:用復(fù)化梯形算法計算的結(jié)果Tn=等分?jǐn)?shù)n=7019已知值與計算值的誤差R=2 .Simpson公式程序：程序1:(求f(x)的n階導(dǎo)數(shù))：symsx%定義函數(shù)f(x)f=x*exp(x)n=input(,輸入所求導(dǎo)數(shù)階數(shù):")f2=diff(f,x,n)%f<f(x)的n階導(dǎo)數(shù)結(jié)果1輸入n=4f2=4*exp(x)+x*exp(x)程序2:clcclearsymsx流義自變量xf=inline('x*exp(x)','x')”定義函數(shù)f(x)=x*exp(x),換

5、函數(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=fminbnd(f3,1,2)%求負(fù)的四階導(dǎo)數(shù)的最小值點，也就是求四階導(dǎo)數(shù)的最大值點對應(yīng)的x值forn=2:1000000啾等分?jǐn)?shù)nRn=-(b-a)/180*(b-a)/(2*n)A4*f2(x1

6、)府算余項ifabs(Rn)<e%ffl余項進(jìn)行判斷break%符合要求時結(jié)束endendh=(b-a)/n%<hSn1=0Sn2=0fork=0:n-1%求兩組連加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%HSn2多力口了k=0時的值，故減去f(a)z=exp(2)R=Sn-z魅已知值與計算值的差fprintf('用Simpson公式計算的結(jié)果Sn=')disp(Sn)fprintf('等分?jǐn)?shù)n=')disp(n)fpri

7、ntf('已知值與計算值的誤差R=')disp(R)輸出結(jié)果顯示：用Simpson公式計算的結(jié)果Sn=等分?jǐn)?shù)n=24已知值與計算值的誤差R=用復(fù)化梯形公式計算的結(jié)果為：，與精確解的誤差為：。等分?jǐn)?shù)n=7019用復(fù)化Simpson公式計算的結(jié)果為：，與精確解的誤差為：。等分?jǐn)?shù)n=243、柯斯特公式求積分：程序代碼：(1) functiony,Ck,Ak=NewtonCotes(fun,a,b,n)ifnargin=1mm,nn=size(fun);ifmm>=8error('為了保證NewtonCotes積分的穩(wěn)定性，最多只能有9個等距節(jié)點！,)elseifnn=2

8、error('fun構(gòu)成應(yīng)為：第一列為x,第二列為y,并且個數(shù)endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseifnargin=4xk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);elseerror('fun積分函數(shù)的句柄，且必須能夠接受矢量輸入！')endelseerror('輸入?yún)?shù)錯誤，請參考函數(shù)幫助！')endCk=cotescoeff(n);Ak=(b-a)*Ck;y=Ak*fx'(2

9、) functionCk=cotescoeff(n)fori=1:n+1k=i-1;Ck(i)=(-1)八(n-k)/factorial(k)/factorial(n-k)/n*quadl(t)intfun(t,n,k),0,n);end(3) functionf=intfun(t,n,k)f=1;fori=0:k-1,k+1:nf=f.*(t-i);end代碼解釋：functiony,Ck,Ak=NewtonCotes(fun,a,b,n)%y=NewtonCotes(fun,a,b,n)%牛頓-科特斯數(shù)值積分公式%參數(shù)說明：%fun,積分表達(dá)式，這里有兩種選擇%(1)積分函數(shù)句柄，必須能夠

10、接受矢量輸入，比如fun=(x)sin(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ù)，必須滿足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);%必

11、須可以接受矢量輸入%fun2=0:;sin(0:;%最多8個點，必須等距%y1=NewtonCotes(fun1,0,6)%y2=NewtonCotes(fun2)ifnargin=1mm,nn=size(fun);ifmm>=8error('為了保證NewtonCotes積分的穩(wěn)定性，最多只能有9個等距節(jié)點！)elseifnn=2error('fun構(gòu)成應(yīng)為：第一列為x,第二列為y,并且個數(shù)為小于10的等距節(jié)點！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseifnargin=4%計算積分節(jié)

12、點xk和節(jié)點函數(shù)值fxxk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);elseerror('fun積分函數(shù)的句柄，且必須能夠接受矢量輸入!')endelseerror('輸入?yún)?shù)錯誤，請參考函數(shù)幫助！')end%計算科特斯系數(shù)Ck=cotescoeff(n);%計算求積系數(shù)Ak=(b-a)*Ck;%求和算積分y=Ak*fx'functionCk=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ù)fori=1:n+1k=i-1;Ck(i)=(-1)八(n-