版權(quán)說(shuō)明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請(qǐng)進(jìn)行舉報(bào)或認(rèn)領(lǐng)
文檔簡(jiǎn)介
1、 Matlab實(shí)現(xiàn)數(shù)值分析插值及積分摘要: (numerical analysis)是研究分析用求解計(jì)算問(wèn)題的及其理論的學(xué)科,是的一個(gè)分支,它以求解數(shù)學(xué)問(wèn)題的理論和方法為研究對(duì)象。在實(shí)際生產(chǎn)實(shí)踐中,常常將實(shí)際問(wèn)題轉(zhuǎn)化為數(shù)學(xué)模型來(lái)解決,這個(gè)過(guò)程就是數(shù)學(xué)建模。學(xué)習(xí)數(shù)值分析這門課程可以讓我們學(xué)到很多的數(shù)學(xué)建模方法。分別運(yùn)用matlab數(shù)學(xué)軟件編程來(lái)解決插值問(wèn)題和數(shù)值積分問(wèn)題。題目中的要求是計(jì)算差值和積分,對(duì)于問(wèn)題一,可以分別利用朗格朗日插值公式,牛頓插值公式,埃特金逐次線性插值公式來(lái)進(jìn)行編程求解,具體matlab代碼見正文。編程求解出來(lái)的結(jié)果為:=+。其中Aitken插值計(jì)算的結(jié)果圖如下: 對(duì)于問(wèn)題
2、二,可以分別利用復(fù)化梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式編寫程序來(lái)進(jìn)行求解,具體matlab代碼見正文。編程求解出來(lái)的結(jié)果為: 0.6932其中復(fù)化梯形公式計(jì)算的結(jié)果圖如下: 問(wèn)題重述問(wèn)題一:已知列表函數(shù)表格 101234121782257分別用拉格朗日,牛頓,埃特金插值方法計(jì)算。問(wèn)題二:用復(fù)化的梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式計(jì)算積分,使精度小于5。問(wèn)題解決問(wèn)題一:插值方法對(duì)于問(wèn)題一,用三種差值方法:拉格朗日,牛頓,埃特金差值方法來(lái)解決。一、拉格朗日插值法:拉格朗日插值多項(xiàng)式如下:首先構(gòu)造個(gè)插值節(jié)點(diǎn)上的插值基函數(shù),對(duì)任一點(diǎn)所對(duì)應(yīng)的插值基函數(shù),由于在所有取零值,因此有因子。
3、又因是一個(gè)次數(shù)不超過(guò)的多項(xiàng)式,所以只可能相差一個(gè)常數(shù)因子,固可表示成:利用得:于是因此滿足 的插值多項(xiàng)式可表示為:從而次拉格朗日插值多項(xiàng)式為:matlab編程:編程思想:主要從上述朗格朗日公式入手:依靠循環(huán),運(yùn)用poly()函數(shù)和conv()函數(shù)表示拉格朗日公式,其中的poly(i)函數(shù)表示以i作為根的多項(xiàng)式的系數(shù),例如poly(1)表示x-1的系數(shù),輸出為1 -1,而poly(poly(1)表示(x-1)*(x-1)=x2-2*x+1的系數(shù),輸出為1 -2 1;而conv()表示多項(xiàng)式系數(shù)乘積的結(jié)果,例如conv(poly(1),poly(1)輸出為1 -2 1;所以程序最后結(jié)果為xn+xn
4、-1+x2+x+1(n的值據(jù)結(jié)果的長(zhǎng)度為準(zhǔn))的對(duì)應(yīng)系數(shù)。在命令窗口輸入edit lagran來(lái)建立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.0000 -0.0000 -0.0000 0 1.0
5、000結(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ù)表示拉格朗日公式,其中的poly(i)函數(shù)表示以i作為根的多項(xiàng)式的系數(shù),例如poly(1)表示x-1的系數(shù),輸出為1 -1,而poly(poly(1)表示(x-1)*(x-1)=x2-2*x+1的系數(shù),輸出為1 -2 1;而conv()表示多項(xiàng)式系數(shù)乘積的結(jié)果,例如conv(poly(1),poly(1)輸出為1 -
6、2 1;所以程序最后結(jié)果為xn+xn-1+x2+x+1(n的值據(jù)結(jié)果的長(zhǎng)度為準(zhǔn))的對(duì)應(yīng)系數(shù)。在命令窗口輸入edit nowpoly來(lái)建立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)將公式表示出來(lái)。在命令窗口輸入edit Aitken 來(lái)建立Aitken.m文件,文件中的程序如下:function f = Aitken(x,y) syms
8、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:圖表 3問(wèn)題二:復(fù)化積分對(duì)于問(wèn)題二來(lái)說(shuō),用復(fù)化的梯形公式,復(fù)化的辛卜生公式,復(fù)化的柯特斯公式結(jié)局問(wèn)題(計(jì)
9、算積分,使精度小于5)。一 復(fù)化的梯形公式:復(fù)化梯形的迭代公式為:;matlab編程:程序1(求f(x)的n階導(dǎo)數(shù):)在命令窗口輸入edit qiudao 來(lái)建立qiudao.m文件,文件中的程序如下:function d=qiudao(x,n)syms x;f=1/x; n=input(輸入導(dǎo)數(shù)階數(shù): ); d=diff(f,x,n);輸入:qiudao(x,n)輸入所求導(dǎo)數(shù)階數(shù):2顯示:n = 2ans = 2/x3結(jié)果為:f2 =2/x3如圖表4:圖表 4程序2:在命令窗口輸入edit tixing 來(lái)建立tixing.m文件,文件中的程序如下:function y=tixing()sy
10、ms 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á)式帶引號(hào),故取負(fù)號(hào),以便求最大值e=5*10(-5); %精度要求值 a=1; %積分下限b=2; %積分上限x1=fminbnd(f3,1,2); %求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn),也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對(duì)應(yīng)的x值for n=2:1000000; %求等分?jǐn)?shù)n Rn=-(b-a)/12*(b-a)/n)2*f2(x1); %計(jì)算余項(xiàng)
11、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=)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(求f(x)的n
12、階導(dǎo)數(shù):)在命令窗口輸入edit qiudao 來(lái)建立qiudao.m文件,文件中的程序如下:function d=qiudao(x,n)syms x;f=1/x; n=input(輸入導(dǎo)數(shù)階數(shù): ); d=diff(f,x,n);輸入:qiudao(x,n)輸入所求導(dǎo)數(shù)階數(shù):4顯示:n = 4ans =24/x5結(jié)果為:f2 = 24/x5如圖表6:圖表 6程序2:在命令窗口輸入edit xinpusheng 來(lái)建立xinpusheng.m文件,文件中的程序如下:function y=xinpusheng()syms x ; f=inline(1/x,x); f2=inline(24/x5,
13、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) clear xinpusheng()運(yùn)行結(jié)果為用Simpson公式計(jì)算的結(jié)果 Sn= 0.6932等分?jǐn)?shù) n= 4結(jié)果如圖表7圖表 7三 復(fù)化的柯特斯公式:牛頓-柯特斯公式如下:matlab編程:(1)在命令窗口輸入edit NewtonCotes 來(lái)建立NewtonCotes.m文件,文件中的程序如下:function y,Ck,Ak=NewtonCotes
14、(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,并且個(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積分
15、函數(shù)的句柄,且必須能夠接受矢量輸入!) end else error(輸入?yún)?shù)錯(cuò)誤,請(qǐng)參考函數(shù)幫助!) end Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx; (2)在命令窗口輸入edit cotescoeff來(lái)建立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來(lái)建立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, 比如: fun
溫馨提示
- 1. 本站所有資源如無(wú)特殊說(shuō)明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請(qǐng)下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請(qǐng)聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁(yè)內(nèi)容里面會(huì)有圖紙預(yù)覽,若沒(méi)有圖紙預(yù)覽就沒(méi)有圖紙。
- 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
- 5. 人人文庫(kù)網(wǎng)僅提供信息存儲(chǔ)空間,僅對(duì)用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對(duì)用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對(duì)任何下載內(nèi)容負(fù)責(zé)。
- 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請(qǐng)與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時(shí)也不承擔(dān)用戶因使用這些下載資源對(duì)自己和他人造成任何形式的傷害或損失。
最新文檔
- 商務(wù)框架完整年終總結(jié)
- 2024護(hù)理核心制度培訓(xùn)
- 物業(yè)品質(zhì)管理培訓(xùn)
- 培訓(xùn)家政人的禮儀課件
- 2024-2025學(xué)年第一學(xué)期初二物理期中考試卷
- 初中八年級(jí)英語(yǔ)上學(xué)期期中考前測(cè)試卷(仁愛版)含答案解析
- 中國(guó)商業(yè)地產(chǎn)物業(yè)與資產(chǎn)管理白皮書 2024
- T-ZFDSA 09-2024 茯苓煲雞湯制作標(biāo)準(zhǔn)
- 高中物理第十一章機(jī)械振動(dòng)1簡(jiǎn)諧運(yùn)動(dòng)課件新人教版選修3-
- 語(yǔ)篇型語(yǔ)法填空題的研究
- 干部基本信息審核認(rèn)定表
- 新疆地方教材五年級(jí)可愛的中國(guó)計(jì)劃、教案
- 求職能力展示
- 《羊水異?!氛n件
- 基于PLC的熱水箱恒溫控制系統(tǒng)
- 痔瘡護(hù)理課件
- 《昆蟲病原線蟲》課件
- 船員勞務(wù)市場(chǎng)分析
- 研學(xué)旅行機(jī)構(gòu)可行性方案
- 小學(xué)勞動(dòng)教育典型案例
- 工業(yè)大數(shù)據(jù)分析與決策智能
評(píng)論
0/150
提交評(píng)論