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

/Matlab實現(xiàn)數(shù)值分析插值與積分摘要：數(shù)值分析<numericalanalysis>是研究分析用計算機求解數(shù)學計算問題的數(shù)值計算方法與其理論的學科,是數(shù)學的一個分支,它以數(shù)字計算機求解數(shù)學問題的理論和方法為研究對象.在實際生產(chǎn)實踐中,常常將實際問題轉(zhuǎn)化為數(shù)學模型來解決,這個過程就是數(shù)學建模.學習數(shù)值分析這門課程可以讓我們學到很多的數(shù)學建模方法.分別運用matlab數(shù)學軟件編程來解決插值問題和數(shù)值積分問題.題目中的要求是計算差值和積分,對于問題一,可以分別利用朗格朗日插值公式,牛頓插值公式,埃特金逐次線性插值公式來進行編程求解,具體matlab代碼見正文.編程求解出來的結(jié)果為：=+.其中Aitken插值計算的結(jié)果圖如下：對于問題二,可以分別利用復(fù)化梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式編寫程序來進行求解,具體matlab代碼見正文.編程求解出來的結(jié)果為：0.6932其中復(fù)化梯形公式計算的結(jié)果圖如下：問題重述問題一：已知列表函數(shù)表格SEQ表格\*ARABIC101234121782257分別用拉格朗日,牛頓,埃特金插值方法計算.問題二：用復(fù)化的梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式計算積分,使精度小于5.問題解決問題一：插值方法對于問題一,用三種差值方法：拉格朗日,牛頓,埃特金差值方法來解決.一、拉格朗日插值法：拉格朗日插值多項式如下：

首先構(gòu)造個插值節(jié)點上的插值基函數(shù),對任一點所對應(yīng)的插值基函數(shù),由于在所有取零值,因此有因子.又因是一個次數(shù)不超過的多項式,所以只可能相差一個常數(shù)因子,固可表示成：利用得：于是因此滿足的插值多項式可表示為：從而次拉格朗日插值多項式為：matlab編程：編程思想：主要從上述朗格朗日公式入手：依靠循環(huán),運用poly〔〕函數(shù)和conv〔〕函數(shù)表示拉格朗日公式,其中的poly〔i〕函數(shù)表示以i作為根的多項式的系數(shù),例如poly〔1〕表示x-1的系數(shù),輸出為1-1,而poly〔poly〔1〕〕表示〔x-1〕*〔x-1〕=x^2-2*x+1的系數(shù),輸出為1-21；而conv〔〕表示多項式系數(shù)乘積的結(jié)果,例如conv〔poly〔1〕,poly〔1〕〕輸出為1-21；所以程序最后結(jié)果為x^n+x^n-1+……+x^2+x+1〔n的值據(jù)結(jié)果的長度為準〕的對應(yīng)系數(shù).在命令窗口輸入editlagran來建立lagran.m文件,文件中的程序如下：function[c,l]=lagran<x,y>w=length<x>;n=w-1;l=zeros<w,w>;fork=1:n+1v=1;forj=1:n+1ifk~=jv=conv<v,poly<x<j>>>/<x<k>-x<j>>;endendl<k,:>=v;endc=y*l;輸入：>>x=[01234];>>y=[121782257];>>lagran<x,y>運行結(jié)果為ans=1.0000-0.0000-0.000001.0000結(jié)果為：=+.如圖表1：圖表SEQ圖表\*ARABIC1二．牛頓插值法newton插值多項式的表達式如下：其中每一項的系數(shù)ci的表達式如下：即為f<x>在點處的i階差商,〔,〕,由差商的性質(zhì)可知：matlab編程：編程思想：主要從上述牛頓插值公式入手：依靠循環(huán),運用poly〔〕函數(shù)和conv〔〕函數(shù)表示拉格朗日公式,其中的poly〔i〕函數(shù)表示以i作為根的多項式的系數(shù),例如poly〔1〕表示x-1的系數(shù),輸出為1-1,而poly〔poly〔1〕〕表示〔x-1〕*〔x-1〕=x^2-2*x+1的系數(shù),輸出為1-21；而conv〔〕表示多項式系數(shù)乘積的結(jié)果,例如conv〔poly〔1〕,poly〔1〕〕輸出為1-21；所以程序最后結(jié)果為x^n+x^n-1+……+x^2+x+1〔n的值據(jù)結(jié)果的長度為準〕的對應(yīng)系數(shù).在命令窗口輸入editnowpoly來建立newpoly.m文件,文件中的程序如下：function[c,d]=newpoly<x,y>n=length<x>;d=zeros<n,n>;d<:,1>=y';forj=2:nfork=j:nd<k,j>=<d<k,j-1>-d<k-1,j-1>>/<x<k>-x<k-j+1>>;endendc=d<n,n>;fork=<n-1>:-1:1c=conv<c,poly<x<k>>>;m=length<c>;c<m>=c<m>+d<k,k>;end輸入：>>x=[01234];>>y=[121782257];>>newpoly<x,y>運行結(jié)果為ans=10001所以=+.如圖表2：圖表SEQ圖表\*ARABIC2三．埃特金插值法：Aitken插值公式如下：遞推表達式為：=+當n=1時,=+當n=2時,=+其中的帶入遞推表達式求得.由此遞推下去,最終得到的結(jié)果.matlab編程：編程思想：埃特金插值多項式又稱作Aitken逐次線性插值多項式,根據(jù)公式的特點,可以利用2次嵌套循環(huán)將公式表示出來.在命令窗口輸入editAitken來建立Aitken.m文件,文件中的程序如下：functionf=Aitken<x,y>symsz;n=length<x>;y1<1:n>=z;fori=1:n-1forj=i+1:ny1<j>=y<j>*<z-x<i>>/<x<j>-x<i>>+y<i>*<z-x<j>>/<x<i>-x<j>>;endy=y1;simplify<y1>;endsimplify<y1<n>>;f=collect<y1<n>>;輸入：>>x=[01234];>>y=[121782257];>>Aitken<x,y>運行結(jié)果為ans=z^4+1所以=+.如圖表3：圖表SEQ圖表\*ARABIC3問題二：復(fù)化積分對于問題二來說,用復(fù)化的梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式結(jié)局問題〔計算積分,使精度小于5〕.一復(fù)化的梯形公式：復(fù)化梯形的迭代公式為：;matlab編程：程序1〔求f〔x〕的n階導數(shù)：〕在命令窗口輸入editqiudao來建立qiudao.m文件,文件中的程序如下：functiond=qiudao<x,n>symsx;f=1/x;n=input<'輸入導數(shù)階數(shù)：'>;d=diff<f,x,n>;輸入:qiudao〔x,n〕輸入所求導數(shù)階數(shù):2顯示：n=2ans=2/x^3結(jié)果為：f2=2/x^3如圖表4：圖表SEQ圖表\*ARABIC4程序2：在命令窗口輸入edittixing來建立tixing.m文件,文件中的程序如下：functiony=tixing<>symsx；%定義自變量xf=inline<'1/x','x'>；%定義函數(shù)f<x>=1/xf2=inline<'2/x^3','x'>；%定義f<x>的二階導數(shù),輸入程序1里求出的f2即可.f3='-2/x^3'；%因fminbnd〔〕函數(shù)求的是表達式的最小值,且要求表達式帶引號,故取負號,以便求最大值e=5*10^<-5>；%精度要求值a=1；%積分下限b=2；%積分上限x1=fminbnd<f3,1,2>；%求負的二階導數(shù)的最小值點,也就是求二階導數(shù)的最大值點對應(yīng)的x值forn=2:1000000；%求等分數(shù)nRn=-<b-a>/12*<<b-a>/n>^2*f2<x1>；%計算余項ifabs<Rn><e%用余項進行判斷break%符合要求時結(jié)束endendh=<b-a>/n；%求hTn1=0；fork=1:n-1%求連加和xk=a+k*h；Tn1=Tn1+f<xk>；endTn=h/2*<<f<a>+2*Tn1+f<b>>>；fprintf<'用復(fù)化梯形算法計算的結(jié)果Tn='>disp<Tn>fprintf<'等分數(shù)n='>disp<n>%輸出等分數(shù)輸入：tixing<>運行結(jié)果為:用復(fù)化梯形計算的結(jié)果Tn=0.6932等分數(shù)n=58結(jié)果如圖表:5圖表SEQ圖表\*ARABIC5二復(fù)化的辛卜生公式：復(fù)化simpson迭代公式為：;matlab編程：程序1〔求f〔x〕的n階導數(shù)：〕在命令窗口輸入editqiudao來建立qiudao.m文件,文件中的程序如下：functiond=qiudao<x,n>symsx;f=1/x;n=input<'輸入導數(shù)階數(shù)：'>;d=diff<f,x,n>;輸入:qiudao〔x,n〕輸入所求導數(shù)階數(shù):4顯示：n=4ans=24/x^5結(jié)果為：f2=24/x^5如圖表6：圖表SEQ圖表\*ARABIC6程序2：在命令窗口輸入editxinpusheng來建立xinpusheng.m文件,文件中的程序如下：functiony=xinpusheng<>symsx;f=inline<'1/x','x'>;f2=inline<'24/x^5','x'>;f3='-24/x^5';e=5*10^<-5>;a=1;b=2;x1=fminbnd<f3,1,2>;forn=2:1000000Rn=-<b-a>/180*<<b-a>/<2*n>>^4*f2<x1>;ifabs<Rn><ebreakendendh=<b-a>/n;Sn1=0;Sn2=0;fork=0:n-1xk=a+k*h;xk1=xk+h/2;Sn1=Sn1+f<xk1>;Sn2=Sn2+f<xk>;endSn=h/6*<f<a>+4*Sn1+2*<Sn2-f<a>>+f<b>>;fprintf<'用Simpson公式計算的結(jié)果Sn='>disp<Sn>fprintf<'等分數(shù)n='>disp<n>調(diào)用xinpusheng.m中xinpusheng函數(shù)：輸入：>>clear>>xinpusheng<>運行結(jié)果為用Simpson公式計算的結(jié)果Sn=0.6932等分數(shù)n=4結(jié)果如圖表7圖表SEQ圖表\*ARABIC7三復(fù)化的柯特斯公式：牛頓-柯特斯公式如下：matlab編程：<1>在命令窗口輸入editNewtonCotes來建立NewtonCotes.m文件,文件中的程序如下：function[y,Ck,Ak]=NewtonCotes<fun,a,b,n>ifnargin==1[mm,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==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〕在命令窗口輸入editcotescoeff來建立cotescoeff.m文件,文件中的程序如下：functionCk=cotescoeff<n>f