matlab實(shí)現(xiàn)數(shù)值分析插值及積分

1、 Matlab實(shí)現(xiàn)數(shù)值分析插值及積分摘要： 數(shù)值分析(numerical analysis)是研究分析用計(jì)算機(jī)求解數(shù)學(xué)計(jì)算問題的數(shù)值計(jì)算方法及其理論的學(xué)科，是數(shù)學(xué)的一個(gè)分支，它以數(shù)字計(jì)算機(jī)求解數(shù)學(xué)問題的理論和方法為研究對象。在實(shí)際生產(chǎn)實(shí)踐中，常常將實(shí)際問題轉(zhuǎn)化為數(shù)學(xué)模型來解決，這個(gè)過程就是數(shù)學(xué)建模。學(xué)習(xí)數(shù)值分析這門課程可以讓我們學(xué)到很多的數(shù)學(xué)建模方法。分別運(yùn)用matlab數(shù)學(xué)軟件編程來解決插值問題和數(shù)值積分問題。題目中的要求是計(jì)算差值和積分，對于問題一，可以分別利用朗格朗日插值公式，牛頓插值公式，埃特金逐次線性插值公式來進(jìn)行編程求解，具體matlab代碼見正文。編程求解出來的結(jié)果為：=+。其中

2、Aitken插值計(jì)算的結(jié)果圖如下： 對于問題二，可以分別利用復(fù)化梯形公式，復(fù)化的辛卜生公式，復(fù)化的柯特斯公式編寫程序來進(jìn)行求解，具體matlab代碼見正文。編程求解出來的結(jié)果為： 0.6932其中復(fù)化梯形公式計(jì)算的結(jié)果圖如下： 問題重述問題一：列表函數(shù)表格 101234121782257分別用拉格朗日，牛頓，埃特金插值方法計(jì)算。問題二：用復(fù)化的梯形公式，復(fù)化的辛卜生公式，復(fù)化的柯特斯公式計(jì)算積分，使精度小于5。問題解決問題一：插值方法對于問題一，用三種差值方法：拉格朗日，牛頓，埃特金差值方法來解決。一、拉格朗日插值法：拉格朗日插值多項(xiàng)式如下： 首先構(gòu)造個(gè)插值節(jié)點(diǎn)上的插值基函數(shù)，對任一

3、點(diǎn)所對應(yīng)的插值基函數(shù)，由于在所有取零值，因此有因子。又因是一個(gè)次數(shù)不超過的多項(xiàng)式，所以只可能相差一個(gè)常數(shù)因子，固可表示成：利用得：于是因此滿足 的插值多項(xiàng)式可表示為：從而次拉格朗日插值多項(xiàng)式為：matlab編程：編程思想：主要從上述朗格朗日公式入手：依靠循環(huán)，運(yùn)用poly函數(shù)和conv函數(shù)表示拉格朗日公式，其中的polyi函數(shù)表示以i作為根的多項(xiàng)式的系數(shù)，例如poly1表示x-1的系數(shù)，輸出為1 -1，而polypoly1表示x-1*x-1=x2-2*x+1的系數(shù)，輸出為1 -2 1；而conv表示多項(xiàng)式系數(shù)乘積的結(jié)果，例如convpoly1，poly1輸出為1 -2 1；所以程序最后結(jié)果為x

4、n+xn-1+x2+x+1n的值據(jù)結(jié)果的長度為準(zhǔn)的對應(yīng)系數(shù)。在命令窗口輸入edit lagran來建立lagran.m文件，文件中的程序如下：function c,l=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1; for j=1:n+1 if k=j v=conv(v,poly(x(j)/(x(k)-x(j); end end l(k,:)=v;endc=y*l;輸入：>> x=0 1 2 3 4;>> y=1 2 17 82 257;>> lagran(x,y)運(yùn)行結(jié)果為ans =1.

5、0000 -0.0000 -0.0000 0 1.0000結(jié)果為：=+。如圖表1：圖表 1二牛頓插值法newton插值多項(xiàng)式的表達(dá)式如下：其中每一項(xiàng)的系數(shù)ci的表達(dá)式如下：即為f (x)在點(diǎn)處的i階差商，由差商的性質(zhì)可知：matlab編程：編程思想：主要從上述牛頓插值公式入手：依靠循環(huán)，運(yùn)用poly函數(shù)和conv函數(shù)表示拉格朗日公式，其中的polyi函數(shù)表示以i作為根的多項(xiàng)式的系數(shù)，例如poly1表示x-1的系數(shù)，輸出為1 -1，而polypoly1表示x-1*x-1=x2-2*x+1的系數(shù)，輸出為1 -2 1；而conv表示多項(xiàng)式系數(shù)乘積的結(jié)果，例如convpoly1，poly1輸出為1 -

6、2 1；所以程序最后結(jié)果為xn+xn-1+x2+x+1n的值據(jù)結(jié)果的長度為準(zhǔn)的對應(yīng)系數(shù)。在命令窗口輸入edit nowpoly來建立newpoly.m文件，文件中的程序如下：function c,d=newpoly(x,y)n=length(x);d=zeros(n,n);d(:,1)=y'for j=2:n for k=j:n d(k,j)=(d(k,j-1)-d(k-1,j-1)/(x(k)-x(k-j+1); endendc=d(n,n);for k=(n-1):-1:1 c=conv(c,poly(x(k); m=length(c); c(m)=c(m)+d(k,k);end輸

7、入：>> x=0 1 2 3 4;>> y=1 2 17 82 257;>> newpoly(x,y)運(yùn)行結(jié)果為ans =1 0 0 0 1所以=+。如圖表2：圖表 2三埃特金插值法：Aitken插值公式如下： 遞推表達(dá)式為： = + 當(dāng)n=1時(shí)， = + 當(dāng)n=2時(shí)， = + 其中的帶入遞推表達(dá)式求得。由此遞推下去，最終得到的結(jié)果。matlab編程：編程思想：埃特金插值多項(xiàng)式又稱作Aitken逐次線性插值多項(xiàng)式， 根據(jù)公式的特點(diǎn)，可以利用2次嵌套循環(huán)將公式表示出來。在命令窗口輸入edit Aitken 來建立Aitken.m文件，文件中的程序如下：func

8、tion f = Aitken(x,y) syms z; n = length(x); y1(1:n) = z; for i=1:n-1 for j=i+1:n y1(j) = y(j)*(z-x(i)/(x(j)-x(i)+y(i)*(z-x(j)/(x(i)-x(j); end y = y1; simplify(y1);end simplify(y1(n); f = collect(y1(n); 輸入：>> x=0 1 2 3 4;>> y=1 2 17 82 257;>> Aitken(x,y)運(yùn)行結(jié)果為ans = z4 + 1所以=+。如圖表3：圖表

9、 3問題二：復(fù)化積分對于問題二來說，用復(fù)化的梯形公式，復(fù)化的辛卜生公式，復(fù)化的柯特斯公式結(jié)局問題計(jì)算積分，使精度小于5。一 復(fù)化的梯形公式：復(fù)化梯形的迭代公式為：;matlab編程：程序1求fx的n階導(dǎo)數(shù)：在命令窗口輸入edit qiudao 來建立qiudao.m文件，文件中的程序如下：function d=qiudao(x,n)syms x;f=1/x; n=input('輸入導(dǎo)數(shù)階數(shù)： '); d=diff(f,x,n);輸入:qiudaox,n輸入所求導(dǎo)數(shù)階數(shù):2顯示：n = 2ans = 2/x3結(jié)果為：f2 =2/x3如圖表4：圖表 4程序2：在命令窗口輸入edit

10、 tixing 來建立tixing.m文件，文件中的程序如下：function y=tixing()syms x ； %定義自變量xf=inline('1/x','x')； %定義函數(shù)f(x)= 1/x f2=inline('2/x3','x') ； %定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3='-2/x3'； %因fminbnd函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，以便求最大值e=5*10(-5)； %精度要求值 a=1； %積分下限b=2； %積分上限x1=fminbnd(f3

11、,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)/12*(b-a)/n)2*f2(x1)； %計(jì)算余項(xiàng) if abs(Rn)<e %用余項(xiàng)進(jìn)行判斷 break % 符合要求時(shí)結(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)；fprintf('用復(fù)化梯形算法計(jì)算的結(jié)果 Tn=')disp(Tn)fprintf('等分?jǐn)?shù) n=&

12、#39;)disp(n) %輸出等分?jǐn)?shù)輸入：tixing()運(yùn)行結(jié)果為:用復(fù)化梯形計(jì)算的結(jié)果 Tn= 0.6932等分?jǐn)?shù) n= 58結(jié)果如圖表:5圖表 5二 復(fù)化的辛卜生公式：復(fù)化simpson迭代公式為：;matlab編程：程序1求fx的n階導(dǎo)數(shù)：在命令窗口輸入edit qiudao 來建立qiudao.m文件，文件中的程序如下：function d=qiudao(x,n)syms x;f=1/x; n=input('輸入導(dǎo)數(shù)階數(shù)： '); d=diff(f,x,n);輸入:qiudaox,n輸入所求導(dǎo)數(shù)階數(shù):4顯示：n = 4ans =24/x5結(jié)果為：f2 = 24/x5

13、如圖表6：圖表 6程序2：在命令窗口輸入edit xinpusheng 來建立xinpusheng.m文件，文件中的程序如下：function y=xinpusheng()syms x ; f=inline('1/x','x'); f2=inline('24/x5','x'); f3='-24/x5' e=5*10(-5); a=1; b=2; x1=fminbnd(f3,1,2); for n=2:1000000 Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1); if abs(Rn)<e

14、 break endendh=(b-a)/n; Sn1=0; Sn2=0;for 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); fprintf('用Simpson公式計(jì)算的結(jié)果 Sn=')disp(Sn)fprintf('等分?jǐn)?shù) n=')disp(n) 調(diào)用xinpusheng.m中xinpusheng函數(shù)：輸入：>> clear>> xinpusheng()運(yùn)行結(jié)果為用Simpson

15、公式計(jì)算的結(jié)果 Sn= 0.6932等分?jǐn)?shù) n= 4結(jié)果如圖表7圖表 7三 復(fù)化的柯特斯公式：牛頓-柯特斯公式如下：matlab編程：(1)在命令窗口輸入edit NewtonCotes 來建立NewtonCotes.m文件，文件中的程序如下：function y,Ck,Ak=NewtonCotes(fun,a,b,n) if nargin=1 mm,nn=size(fun); if mm>=8 error('為了保證NewtonCotes積分的穩(wěn)定性，最多只能有9個(gè)等距節(jié)點(diǎn)！') elseif nn=2 error('fun構(gòu)成應(yīng)為：第一列為x，第二列為y，并且

16、個(gè)數(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ù)錯(cuò)誤，請參考函數(shù)幫助！') end Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx&#

17、39; 2在命令窗口輸入edit cotescoeff來建立cotescoeff.m文件，文件中的程序如下：function 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在命令窗口輸入edit intfun來建立intfun.m文件，文件中的程序如下：function f=intfun(t,n,k) f=1; for i=0:k-1,k+1:n f=f.*(t-i); end% fun，積分表達(dá)式，這里有兩種選擇 %(1)積分函數(shù)句柄，必須能夠接受矢量輸入，比方fun=(x)1./x % (2)x,y坐標(biāo)的離散點(diǎn)， 第一列為x， 第二列為y， 必須等距， 且節(jié)點(diǎn)的個(gè)數(shù)小于9，