數(shù)值分析課程設(shè)計

1、摘要實驗一拉格朗日插及數(shù)值求解1.1實驗?zāi)康牧私釲agranger差值的根本原理和方法通過實例掌握用MATLAB求插值的方法根據(jù)實際計算理論,利用Lagranger插值多項式計算1.2實驗原理設(shè)X.,XiX2,.,Xn及y.fXii=0,i,.,n,Lnx為不超過n次多項式且滿足Ln(Xi)=yi(i=0,1,.n)易知LnXv|0Xy.lnXyn其中,liX均為n次多項式,再由Xjj*i為n次多項式1iX的n個根知nli(X)=C|IX-Xjj=0.最后,由nli(Xj)=C|1(Xi-Xj)=1=jfj號1nI1(Xi-Xj)jfje,i=0,1,.,n.nlOi總之,Ln(X)=i=0n

2、X-Xj二.j=0X-Xjli(X)=j#式為n階Lagrange插值公式,其中,liXi=0,1,.n稱為n階Lagrange插值的基函數(shù)(X-X0).(X-Xi4)(X-Xi1).(X-Xn)li(X)=,i(X-M).(Xi-Xy)(Xi-Xi1).(Xi-Xn)=0,12.,n1.3實驗內(nèi)容functiony=lagranger(X0,y0,X);%UNTITLEDSummaryofthisfunctiongoeshere%Detailedexplanationgoesheren=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:nl

3、i=1.0;forj=1:nifj=kli=li*(z-x0(j)/(x0(k)-x0(j);endends=li*y0(k)+s;endy(i)=s;end1.4實驗案例及結(jié)果分析(1)輸入：x0=4,5,6;y0=10,5.25,1;x=5;y=lagranger(x0,y0,x)y-5. 2500輸入：X0=1,4,8;y0=6,3.2,4;x=4;y=lagranger(x0,y0,x)y=實驗二LU分解法解線性方程組2.1 實驗?zāi)康? .了解LU分解法解線性方程組的根本原理；2 .熟悉計算方法的技巧和過程,能用LU分解法解實際問題;3 .用matlab實現(xiàn)LU分解.2.2 實驗原理1

4、 .假設(shè)一個線性方程組系數(shù)矩陣為n階方陣A且各階順序主子式均不為0那么A的LU分解存在且唯一.將高斯消去法改寫為緊湊形式,可以直接從矩陣A的元素得到計算L,U元素的遞推公式,而不需任何中間步驟,這就是所謂直接三角分解法.一旦實現(xiàn)了矩陣A的LU分解,那么求解Ax=b的問題就等價于求解兩個三角形方程組：Ly=b,求y;Ux=y,求x.2 .在滿足1的條件下課推導(dǎo)得出以下公式j(luò)1j4uij-ajiljkukik1i4(2)yi刊-xlikykk1nxi=(yi-'Uikxk)/Uiiki1(1)lij=(aij-工likuij)/uijk13 .公式1用于求解矩陣L、U,公式2用于會帶求解V

5、、X.從公式中可以看出：L對角線上元素為1,U第一行與A第一行相同.4 .LU分解的具體過程和順序如下：(1)第一步分解：51=41(2)第二步分解：l21=a21/U11U12=a12U22-a22-l21U12(3)第三步分解：I31=a31/U11l32=(a32-I31U12)/U22U13=a13U23=a23-I21U13U33=a33-I31U13-I32U23(0)第n步分解：依次計算：In1、In2Inn,UinUnn1"U11U12.Um11211U22U2n12.3 實驗內(nèi)容編寫一個M文件fUnctionL,U,x=LU_x(A,b)n,m=size(A);if

6、n-=merror('TherowsandcoIUmnsofmatrixAmUStbeeqUaI!');retUrn;endfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)=0)error('ThematrixcannotbedividedbyLU!')return;endendAn,n=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:nL(i,i)=1;endfork=1:nforj=k:nU(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)'

7、;);endfori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)')/U(k,k);endendy(1)=d(1);fori=2:nforj=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i)；endx(n)=y(n)/U(n,n);fori=(n-1):-1:1forj=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end2.4 實驗案例及結(jié)果分析一1070-32.0999996515101反12x22一x4x3一815.9000015-11在MATLA

8、Bt令窗體輸入:>>A=10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2;>>b=8,5.900001,5,1;>>L,U,x=Lu_x(A,b)得到結(jié)果如下：10,0000-7.000001.0000-3.00002.10006.OOCO2.00005.DOOO-1.00005.0000700002.00001.00000之00001.0e+006*0.0000000一口,0000口.0000000.0000-2.50000.000000.0000-2.4C000.OOCO0.0000一L0e+007*0,0000-0

9、.000000.00000-a.ocoo0.OOCO0.CIOOO001.50000.5750a00.aaoo-0,0098-1.01181.01061.0157實驗三龍貝格求積公式求數(shù)值積分3.1實驗?zāi)康? .熟練掌握龍貝格求積的根本思路和步驟;2 .培養(yǎng)編程與上機調(diào)試水平；3 .利用龍貝格求積方法求解積分.3.2 實驗原理(1),b-ah=置n=1,精度要求3n(2)計算TT%).fMh2n=二hn(3)置2,并計算(0)1(0)n12?=Tn()h2nXf(a(2k-1')232n9km(2Tk49k(mT-計5mm二一(6)假設(shè)m=14專(7)；否貝U,置2,k=k+1,轉(zhuǎn)(5

10、).T(k)_T(kJ1)(k)(7)假設(shè)11,那么停止計算(輸出T1),否那么轉(zhuǎn)(3)3.3 實驗內(nèi)容functionR=romberg(f,a,b,e)%參數(shù)介紹：%f-被積函數(shù)f(x)%a-x的左區(qū)間.%b-x的右區(qū)間.%e-誤差限.%結(jié)果：%R-返回Romberg表.n=1;%區(qū)間二分次數(shù)while1%在此僅代表屢次二分,在后面判斷循環(huán)終止R=zeros(n+1,n+1);%生成(n+1)*(n+1)的0矩陣R(0+1,0+1)=(b-a)/2*(feval(f,a)+feval(f,b);%始值2點梯形公式.fori=1:n%根據(jù)公式計算Romberg表的第一列.h=(b-a)/2A

11、i;s=0;fork=1:2A(i-1)s=s+feval(f,a+(2*k-1)*h);endR(i+1,0+1)=R(i-1+1,0+1)/2+h*s;endforj=1:n%計算Romberg表的其他列.fac=1/(4Aj-1);form=j:nR(m+1,j+1)=R(m+1,j-1+1)+fac*(R(m+1,j-1+1)-R(m-1+1,j-1+1);endendifabs(R(n,n)-R(n+1,n+1)<e%當(dāng)精度滿足設(shè)定要求時退出break;endn=n+1;%未到達指定精度繼續(xù)二分endt=R(i,j)3.4 實驗案例及結(jié)果分析15.dx用Romberg求積法計算

12、積分x+3,取精度要求=10-50在Matlab命令窗口輸入:fun=inline('5./(3+x)','x');romberg(fun,-1,1,1e-6)輸出結(jié)果如下：3.4657ans=3.75000000D354173.472200003.48513.46533.46590003.47063.46583.46573.4657003.46703,46573.46573.46573.465703.46603,46573.46"3.46573.46573.4657由輸出結(jié)果可知最終積分為3.4657.實驗四Runge-Kutta方法求常微分方程數(shù)值

13、解4.1 實驗?zāi)康? .熟悉Runge-Kutta常微分方程初值問題的根本原理2 .了解Runge-Kutta常微分方程初值問題的計算流程3 .能編程實現(xiàn)Runge-Kutta常微分方程初值問題4.2 實驗原理在歐才立法中(Xn,yn,h)=f(Xn,yn)的根底上增加了計算一個右函數(shù)f的值,可望p=2.假設(shè)要使得到的公式階數(shù)p更大,中就必須包含更多的f值.實際上從與方程等價的積分形式即X4n1y(Xn1-yn)=xf(X,y(x)dxxn假設(shè)要使公式階數(shù)提升,就必須使右端積分的數(shù)值求積公式精度提升,它必須要增加求積節(jié)點,為此可將上述公式的右端項表示為Xn1rXf(x,y(x)dx之h

14、3;Cif(Xn+九由,y(Xn+兒h)XnT一般來說,點數(shù)r越多,精度越高,上式右端相當(dāng)于增量函數(shù)中(X,y,h),為得到便于計算的顯示方法,可類似改良的歐拉法,將公式表示為yn1=ynh(Xn,yn,h)r(Xn,yn,h)=.cKii1iKi=f(Xnih,ynh'ijKj),i=2,.,r這里G,A,均為常數(shù),上式稱為顯示龍格一庫塔(Ruuge-Kutta)法,簡稱R-K方法.當(dāng)r=1時,就是歐拉法,當(dāng)r=2時,改良的歐拉法就是其中的一種.yn1=ynh*(K1七)/2Ki=f(Xn,yn)K2=f(Xnh,ynh*K1)依次類推,如果在區(qū)間(為兇2內(nèi)多預(yù)估幾個點上的斜率值KK

15、2、.Km,并用他們的加權(quán)平均數(shù)作為平均斜率K"的近似值,顯然能構(gòu)造出具有很高精度的高階計算公式.經(jīng)數(shù)學(xué)推導(dǎo)、求解,可以得出四階龍格-庫塔公式,也就是在工程中應(yīng)用廣泛的經(jīng)典龍格-庫塔算法：hyni=yn(Ki2K22K3K4)6ynyy-1Kh-2+h-2+2KJI2Xnhhyn2K2(xnh)ynhK3hK3=ynK22K4=ynhK34.3實驗內(nèi)容四階龍格一庫塔法的計算公式為：Ki=g(Xi,yi)hhK2=g(Xi二,yi-Ki)22hhK3=g(Xi,yiK2)22K4=g(Xih,yhK3)hyii=yiKi2K22K3K4)四階龍格一庫塔公式的Matlab程序代碼：fun

16、ctiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(i)=y0;x=a:h:b;var=findsym(f);fori=2:N+iK1=Funval(f,varvec,x(i-1)y(i-1);K2=Funval(f,varvec,x(i-1)+h/2y(i-1)+K1*h/2);K3=Funval(f,varvec,x(i-1)+h/2y(i-1)+K2*h/2);K4=Funval(f,varvec,x(i-1)+hy(i-1)+h*K3);y(i)=y(i-1)+h*(K1+2*K

17、2+2*K3+K4)/6;end函數(shù)運行時需要調(diào)用以下函數(shù)：functionfv=Funval(f,varvec,varval)var=findsym(f);iflength(var)<4ifvar(1)=varvec(1)fv=subs(f,varvec(1),varval(1);elsefv=subs(f,varvec(2),varval(2);endelsefv=subs(f,varvec,varval);end4.4實驗案例及結(jié)果分析四階龍格一庫塔求解一階常微分方程應(yīng)用實例,用四階龍格一庫塔法求下面常微分方程的數(shù)值解.J"y3y=10MxM1在MATLABt令窗口輸入:

18、symsxy;z=4*x+y+3yy=DELGKT3_kuta(z,0.1,0,1,1,xy)輸出結(jié)果如下：yy=L00000000000000L44133333333333L97114688B888892,598745836703703,334413940530384.18951S473276155.176309587856725137.603779512990729.0761436584402510.74541809935288實驗五用插值方法求解帶式錄音機的播放時間問題5.1 實驗?zāi)康倪M一步穩(wěn)固、增強插值模型的建模、求解水平.初步研究插值法的穩(wěn)定性和數(shù)值解法.學(xué)

19、習(xí)掌握用MATLAB軟件求解插值的相關(guān)命令.5.2 實驗原理1 .用interp1可以繪制出曲線；2 .用8次多項式做插值,求解8個系數(shù),從而可以顯示給出一個函數(shù)5.3 實驗內(nèi)容>>猿項式插值>>t=205,430,677,945,1233,1542,1872,2224;>>c=100,200,300,400,500,600,700,800;>>c1=0:10:1000;>>t1=interp1(c,t,c1,'nearest');>>t2=interp1(c,t,c1,'linear');

20、>>t3=interp1(c,t,c1,'spline');>>t4=interp1(c,t,c1,'cubic');>>plot(c,t,'*',c1,t1,':b',c1,t2,'-r',c1,t3,'-g',c1,t4,'.-r')>>legend(原始數(shù)據(jù)','最近點插值,線性插值','樣條插值','立方插值')輸出結(jié)果如下：>>核項式擬合>>t=2

21、05,430,677,945,1233,1542,1872,2224;>>c=100,200,300,400,500,600,700,800;> >n=8;%以合度> >p=polyfit(c,t,n)Warning:Polynomialisnotunique;degree>=numberofdatapoints.> Inpolyfitat720.0000-0.00000.0000-0.00480.8068000> >c1=linspace(0,10,1000);%強由范圍> >t1=polyval(p,c1);%5+算在c1數(shù)據(jù)點的多項式值> >plot(c,t,'*',c,t,c1,t1,':b')>