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

1、計(jì)目 設(shè)題共15題如下成績(jī)課程設(shè)計(jì)我再也回不到大二了，大學(xué)是那么短暫設(shè)計(jì)題目 數(shù)值分析學(xué) 生姓名李飛吾學(xué)號(hào) XXXXXXXX專業(yè)班級(jí) 信息計(jì)XXXXX班指導(dǎo)教師課 程 設(shè) 計(jì) 主 要 內(nèi) 容設(shè)計(jì)目的： 通過不同題目的理解，進(jìn)行算法分析。 通過MATLAB軟件進(jìn)行編程對(duì)題目進(jìn)行解決。 個(gè)別題目設(shè)計(jì)驗(yàn)證，加深對(duì)數(shù)值分析的理解。 函數(shù)的圖像繪制的運(yùn)用 設(shè)計(jì)題目： 題1. 1利用逆向遞推的方法求解問題，通過條件終止地推 題1. 2從某個(gè)初始值開始，利用遞推公式進(jìn)行積分估值 題1. 3 繪制Koch分形曲線，節(jié)點(diǎn)之間的關(guān)系與坐標(biāo)變換題2. 1用高斯消元法的消元過程作矩陣分解，LU分解題2. 2矩陣分解方

2、法求上題中 A的逆矩陣，針對(duì)不冋的b,而重復(fù)利用已知的 LU題2. 3驗(yàn)證希爾伯特矩陣的病態(tài)性，矩陣基本運(yùn)算 題3. 1用泰勒級(jí)數(shù)的有限項(xiàng)逼近正弦函數(shù)，由圖像觀察逼近效果 題3. 2繪制飛機(jī)的降落曲線，線性方程組求解，與繪圖題4. 1線性擬合的函數(shù)表達(dá)式的推導(dǎo)，使用了兩種代碼方法題5. 1用幾種不冋的方法求積分，觀察數(shù)值積分的逼近效果題5. 5 求空間曲線L弧長(zhǎng)。求導(dǎo)后使用符號(hào)函數(shù)積分計(jì)算題6. 1用歐拉公式和四階龍格-庫(kù)塔法分別求解下列初值問題，代碼搜索內(nèi)容。題6. 4常微分方程的解，dsolve()函數(shù)使用題8. 2差分法解常微分方程邊值問題，ode函數(shù)無能為力，Matlab中提供bvp解

3、算器。soli nit=bvp in it(x, yin it,params)sol= bvpsolver(odefu n,bcfu n, soli nit,opti ons)題& 3求解圓的半徑，圓心。線性方程組解參數(shù)設(shè)計(jì)總結(jié)：(1) 算法是題目的解題核心，好的算法可以使計(jì)算更加精確(題5.1)(2) 圖形繪制在今后的課程設(shè)計(jì)，或者是論文中可以用到。(3) 無法解決的問題，需要請(qǐng)教室友，或者上網(wǎng)查閱。(4) MATLAB是一個(gè)很強(qiáng)大的軟件，提供了很多內(nèi)置的數(shù)學(xué)函數(shù)，直接進(jìn)行解題。查閱資料時(shí)了解到一些MATLAB論壇。通過帖子閱讀，了解到了MATLAB在科學(xué)計(jì)算方面，模擬，圖形，視頻等起到的作

4、用。增加了對(duì)“計(jì)算科學(xué)“的理解。指 導(dǎo) 老 師 評(píng) 語(yǔ)建議：從學(xué)生的工作態(tài)度、工作量、設(shè)計(jì) (論文的)創(chuàng)造性、學(xué)術(shù)性、使用性及書面表達(dá)能力等方面給出評(píng)價(jià)。 簽名：20年月日數(shù)值分析課程設(shè)計(jì)1.1水手、猴子和椰子問題：五個(gè)水手帶了一只猴子來到南太平洋的 一個(gè)荒島上，發(fā)現(xiàn)那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊， 很快就入睡了。第一個(gè)水手醒來后，把椰子平分成五堆，將多余的一只 給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個(gè)水 手也陸續(xù)起來，和第一個(gè)水手一樣，把椰子分成五堆，恰多一只猴子， 私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆， 每人分一堆，正好余一只再給

5、猴子，試問原先共有幾只椰子？ (15621)試分析椰子數(shù)目的變化規(guī)律，利用逆向遞推的方法求解這一問題 解：算法分析：解該問題主要使用遞推算法，關(guān)于椰子數(shù)目的變化規(guī)律 可以設(shè)起初的椰子數(shù)為po，第一至五次猴子在夜里藏椰子后，椰子的數(shù) 目分別為Pb,p1,p2,P3,P4再設(shè)最后每個(gè)人分得x個(gè)椰子，由題：4 1pk 1 (pk n1. 2 設(shè)，In dx05 x(1) 從I。盡可能精確的近似值出發(fā)，利用遞推公式：1 In 5In 1-(n 1,2,L 20)n計(jì)算機(jī)從h到I 20的近似值；(2) 從較粗糙的估計(jì)值出發(fā)，用遞推公式：( k=0,1,2,3,4 ) x -( p5 1)5 5所以P5

6、5x 1, p p if i=0.1 &i0&i jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U (k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end RA,hl,U ,L end end以上上代碼保存為M文件，并在命令窗口輸入A=20 2 3;1 8 1; 2 -3 15;b=0 0 0;h=zhjLU(A)程序運(yùn)行結(jié)果:L =U=1.00000 020.00002.00003.00000

7、.05001.0000 007.90000.85000.1000-0.40511.00000015.0443實(shí)驗(yàn)驗(yàn)證：可以直接使用MATLAB內(nèi)置LU分解A=20 2 3;1 8 1; 2 -3 15;L U=lu(A)輸出結(jié)果與上程序輸出結(jié)果一致。2. 2用矩陣分解方法求上題中A的逆矩陣。記100b|0 ,b21 ,b30001分別求解方程組AX b1, AX b,AX由于三個(gè)方程組系數(shù)矩陣相同，可以將分解后的矩陣重復(fù)使用。對(duì)第一個(gè)方程組，由于A=LU，所以先求解下三角方程組LY b,再求解上三角 方程組UX Y，則可得逆矩陣的第一列列向量；類似可解第二、第三方 程組，得逆矩陣的第二列列向量

8、的第三列列向量。由三個(gè)列向量拼裝可 得逆矩陣A1。解:MATLAB代碼如下： b1=1;0;0; b2=0;1;0; b3=0;1;1; A=20,2,3;1,8,1;2,-3,15;L二1,0,0;0.05,1,0;0.1,-0.4051,1;U=20 2 3;0 7.9 0.85;0 0 15.0443; Y1=Lb1X1=UY1Y2=Lb2X2=UY2Y3=Lb3X3=UY3Y1 =1.0000-0.0500-0.1203X1 =0.0517-0.0055-0.0080Y2 =0實(shí)驗(yàn)驗(yàn)證：X1 X2 X3 ans =0.0517-0.0164-0.0257-0.00550.12370.1

9、165-0.00800.02690.0934而：inv （ A）=X1 X2 X3得證1.00000.4051X2 =-0.01640.12370.0269Y3 =01.00001.4051X3 =-0.02570.11650.0934-9 -ii2. 3驗(yàn)證希爾伯特矩陣的病態(tài)性：對(duì)于三階矩陣11/21/3H1/21/31/41/31/41/5取右端向量b 11/6 13/12 1用泰勒級(jí)數(shù)的有限項(xiàng)逼近正弦函數(shù)y（x） sin x,x 0,y1（x） x,x 0,/23y2（x） x x /6,x 0,/235 y3（x） x x /6 x /120, x 0,7/60T，驗(yàn)證：(1) 向量X

10、 X1 X2 X3T 1 1 1T是方程組HX b的準(zhǔn)確解；(2) 取右端向量b的三位有效數(shù)字得b 183 M8 0.783t，求方程組 的準(zhǔn)確解X，并與X的數(shù)據(jù)1,1,1作比較。說明矩陣的病態(tài)性。解: (1)H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;X=1;1;1;b=H*Xb =1.83331.08330.7833與題中相同 先求出解X，與數(shù)據(jù)1,1,1作比較。H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;b=1.83;1.08;0.783;X=HbX =1.08000.54001.4400與1,1,1!相差較大，矩陣為病態(tài)矩陣用計(jì)算

11、機(jī)繪出上面四個(gè)函數(shù)的圖形。解：MATLAB代碼如下(1) syms x; taylor(si n( x) x=0:0.01*pi:pi plot(x,s in( x)syms x;taylor(x) x=0:0.01*pi:pi/2 plot(x)syms x;taylor(x-xA3/6) fplot(x-xA3/6,0 pi/2) syms x;taylor(x-xA3/6+xA5/120)fplot(x-xA3/6+xA5/120,0 pi/2)ita- i|9 3(叫ftfie卜HA- J1-Al1 b結(jié)果圖形右：x=0:0.1:pi; y=si n( x);plot(x,y,-sk)

12、; hold on x=0:0.1:pi/2;y=x;plot(x,y,-b*) hold on fplot(x-x.A3/6,0,pi/2,0,2,2e-3,-gx) hold onfplot(x-x.A3/6+x.A5/120,0,pi/2,0,2,2e-3,-r.) hold offlege nd(si n(x),x,x-xA3/6,x-A3/6+xA5/120,2) xlabel(x)ylabel(y)title(Taylor approximate n)結(jié)果分析：圖中紅色點(diǎn)線為正弦曲線，藍(lán)色的星線為一階泰勒逼近， 綠色叉線為二階泰勒逼近，黑色正方形線為三階泰勒逼近?？梢娙A泰勒逼近效

13、果最 好，泰勒級(jí)數(shù)越高，逼近效果越好。3. 2繪制飛機(jī)的降落曲線一架飛機(jī)飛臨北京國(guó)際機(jī)場(chǎng)上空時(shí)，其水平速度為 540km/h，飛行高 度為1 000m。飛機(jī)從距機(jī)場(chǎng)指揮塔的橫向距離12 000m處開始降落。 根據(jù)經(jīng)驗(yàn)，一架水平飛行的飛機(jī)其降落曲線是一條三次曲線。 建立直 角坐標(biāo)系，設(shè)飛機(jī)著陸點(diǎn)為原點(diǎn)0,降落的飛機(jī)為動(dòng)點(diǎn)P(x,y)，則x表示飛機(jī)距指揮塔的距離，y表示飛機(jī)的飛行高度，降落曲線為y(x) a0 a-ix a2x23a3X該函數(shù)滿足條件：y(0)0, y(12 000)1000y(0)0, y(12000)0(1)試?yán)脃(x)滿足的條件確定三次多項(xiàng)式中的四個(gè)系數(shù);(2)用所求出的三

14、次多項(xiàng)式函數(shù)繪制出飛機(jī)降落曲線。function S=f(x1,y1,x2,y2,x3,y3,x4,y4) format longa1=1 ,x1,x1A2,x1A3;a2=1以2以2八2以2八3; a3=0,1,2*x3,3*x3A2; a4=0,1,2*x4,3*x4A2; a=a1; a2;a3;a4; b=y1； y2;y3;y4; s=ab;x=-12000:250:0;y=s(3)*x42-s(4)*x.A3 plot(x,y,-d) xlabel(x) ylabel(y)以上為M文件內(nèi)容，在命令窗口鍵入f(0,0,12000,1000,0,0,12000,0)運(yùn)行結(jié)果： ans為

15、多項(xiàng)式系數(shù)ans =1.0e-004 *00WQaoo-TOD-WQ-壬Soo-、42QDun-012000-TDOOO-SOOO-40000.20833333333333-0.000011574074074. 1曾任英特爾公司董事長(zhǎng)的摩爾先生早在 1965年時(shí)，就觀察到 一件很有趣的現(xiàn)象：集成電路上可容納的零件數(shù)量，每隔一年半左右 就會(huì)增長(zhǎng)一倍，性能也提升一倍。因而發(fā)表論文，提出了大大有名的摩爾定律（Moore s LaW,并預(yù)測(cè)未來這種增長(zhǎng)仍會(huì)延續(xù)下去。下面 數(shù)據(jù)中，第二行數(shù)據(jù)為晶片上晶體數(shù)目在不同年代與1959年時(shí)數(shù)目比較的倍數(shù)。這些數(shù)據(jù)是推出摩爾定律的依據(jù)：年代195919621963

16、19641965增加倍數(shù)13456試從表中數(shù)據(jù)出發(fā)，推導(dǎo)線性擬合的函數(shù)表達(dá)式。解：解法一MATLAB 碼： x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p仁polyfit(x,y,1) y1二polyval(p1,x) plot(x,y1,-,x,y,r*) xlabel(x),ylabel(y);運(yùn)行結(jié)果：p1 =1.0e+003 *0.0008 -1.62554.96235.7925y1 =0.81133.3019 4.1321線性擬合的函數(shù)表達(dá)式：Y=0.8302x-1.6255e+003解法二：運(yùn)行結(jié)果：xs =1.0e+003 *-1.6255283

17、018912380.000830188679248x=1959 1962 1963 1964 1965; y=1;3;4;5;6;for i=1:le ngth(x) for j=1:2A(i,j)=x(i)F-1);endend L,U=lu(A*A); xs=U(L(A*y)從而年代y與增加倍數(shù)x之間的關(guān)系為：y=-1625.528301891238+0.830188679248x-235. 1用幾種不同的方法求積分012 dx x的值運(yùn)行結(jié)果：11 =pi12 =313 =3.133314 =3.14165. 5求空間曲線L :z 2 cost sint弧長(zhǎng)公式為L(zhǎng) 0 x2(t)y2(

18、t) z2(t)dt解：運(yùn)行程序?yàn)椋簊yms t;x=diff(cos(t); y=diff(si n( t);z=diff(2-cos(t)-si n( t);y=i nt(xA2+yA2+zA2)A0.5,t,0,2*pi); % digits(14);i=vpa(y)運(yùn)行結(jié)果：i二8.7377525709894符號(hào)函數(shù)積分(1)牛頓-萊布尼茨公式；(2)梯形公式；(3)辛卜生公式；(4)復(fù) 合梯形公式。解：syms xi1=i nt(4心+xT),x,0,1)a=0;b=1;h_b a ；i2=(4/(1+aA2)+4/(1+bA2)/2i3=h/6*(4/(1+aA2)+4*4/(1+

19、(a+b/2)A2)+4/(1+bA2)M=100;h=(b-a)/M;i4=0;for k=1:(M-1)x=a+h*k; i4=i4+4/(1+xA2); end i4=h*(4/(1+aA2)+4/(1+bA2)/2+h*i4結(jié)果分析：牛頓萊布尼茲公式得到精確結(jié)果，采用梯形公式得到的結(jié)果比采用Simpson公式的精確度要低，米用復(fù)化梯形公式在步長(zhǎng)取得越來越小 的狀態(tài)下可以提高精度。x cost y sint6. 1用歐拉公式和四階龍格-庫(kù)塔法分別求解下列初值問題;(i)y0.9y1 2x,y(0)1；x 0,1（2）y芒，）2；X 0,1解：(1)歐拉公式：function t,x=Eu

20、ler(fun,t0,tt,x0,N) h=(tt-t0)/N;t=t0+0:N*h;x(1,:)=x0;for k=1:Nf=feval(fu n, t(k),x(k,:);f=f;x(k+1,:)=x (k, :)+h*f;end以Euler.m保存fun cti on f=Euler_fu n(t,x)f=0.9*x./(1+2*t);_以Euler_fu n.m保存fun cti on main _Euler t,x=Euler(Euler_fu n,0,1,1,20); fh=dsolve(Dx=0.9*x/(1+2*t),x(0)=1,); for k=1:8ft(k)=t(k);

21、fx(k)=subs(fh,ft(k);endt,x/ 亠以main_Euler.m保存輸入：main_Euler四階龍格-庫(kù)塔法function R=rk4(f,a,b,ya,N)%y=f(x,y)%a,b為左右端點(diǎn)%N為迭代步長(zhǎng)%h為步長(zhǎng)%ya為初值h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;Y(1)=ya;for j=1:N運(yùn)行結(jié)果：ans =0 1.00000.05001.04500.10001.08780.15001.12850.20001.16760.25001.20510.30001.24130.35001.27620.4000

22、1.31000.45001.34270.50001.37450.55001.40550.60001.43560.65001.46490.70001.49360.75001.52160.80001.54900.85001.57580.90001.60210.95001.62781.00001.6531運(yùn)行結(jié)果：ans =01.00000.12501.09540.25001.18070.37501.25850.50001.33060.62501.39810.75001.46170.87501.52231.00001.5801k1=h*feval(f,T(j), Y(j); k2=h*feval(f

23、,T(j)+h/2, Y(j)+k1/2); k3=h*feval(f,T(j)+h/2, Y(j)+k2/2); k4=h*feval(f, Y(j)+h, Y(j)+k3);Y (j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endan s=T Y以rk4.m保存fun cti on z=f(x,y)z=0.9*y./(1+2*x);以f.m保存輸入：rk4(f,0,1,1,8)6. 4列出函數(shù)f(x)1 ln (x) (x)在區(qū)間0, e上的函數(shù)值表并作出它的圖象。其中，(X)是初值問題dd-ln (x) 2xdx的解(0) 0解v=dsolve(Dv*log(x)=2*x

24、,v(0)=0,x);f=(1-log(v)*vf =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x) ezplot(f)輸出結(jié)果：f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)7.4 一個(gè)10次項(xiàng)式P(x) x10 a1x9 a2x8 La 9 xa10的系數(shù)為1 a1 a2 a9 a10=1 -55 1320 -81 50 157 773 02 055 341 693 0-840 950 0 127 535 76 -106 286 40 632 880 0試用多項(xiàng)式的求根指令roots求出該10次方程的10個(gè)根

25、，然后修改9次項(xiàng)的系數(shù)-55為-56,得新的10次方程，求解新的方程，觀察根的變化是否很顯著。解：p=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800; roots(p)ans =10.6051 + 1.0127i10.6051 - 1.0127i8.5850 + 2.7898i8.5850 - 2.7898i5.5000 + 3.5058i5.5000 - 3.5058i2.4150 + 2.7898i2.4150 - 2.7898i0.3949 + 1.0127i0.3949 - 1.0

26、127ip=1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800; roots(p)ans =21.73357.3501 + 7.7973i7.3501 - 7.7973i4.3589 + 3.3285i4.3589 - 3.3285i5.18312.4378 + 2.7974i2.4378 - 2.7974i0.3949 + 1.0127i0.3949 - 1.0127i結(jié)果分析：改變系數(shù)之后，根的變化顯著& 2用差分法解常微分方程邊值問題：y y 1 0,x (0,1)y(0) 0,y(1) 0取h=0.1,xj=jh(j=0,1,2/-,10)求y1, y2,y9并與該問題的準(zhǔn)確解y(x) e e x ex 11 e 1 e比較，列出各節(jié)點(diǎn)處的近似解、準(zhǔn)確解和誤差。解析：ode*函數(shù)無能為力Matlab中，提供了 bvp解算器。soli nitbvp in it(x,yin it,params)solbvpsolver(odefu n,bcfu n,soli nit