1、課 程 設 計 設計題目 數值分析 學生姓名 李飛吾 學 號 xxxxxxxx 專業班級 信息計xxxxx班 指導教師 設 計題 目共15題如下成績課 程 設 計 主 要 內 容設計目的：通過不同題目的理解，進行算法分析。通過MATLAB軟件進行編程對題目進行解決。個別題目設計驗證，加深對數值分析的理解。函數的圖像繪制的運用設計題目：題11 利用逆向遞推的方法求解問題，通過條件終止地推題12 從某個初始值開始，利用遞推公式進行積分估值題13 繪制Koch分形曲線，節點之間的關系與坐標變換題21 用高斯消元法的消元過程作矩陣分解，LU分解題22 矩陣分解方法求上題中A的逆矩陣，針對不同的b,而重

2、復利用已知的LU題23 驗證希爾伯特矩陣的病態性,矩陣基本運算題31 用泰勒級數的有限項逼近正弦函數，由圖像觀察逼近效果題32 繪制飛機的降落曲線，線性方程組求解，與繪圖題41 線性擬合的函數表達式的推導，使用了兩種代碼方法題51 用幾種不同的方法求積分，觀察數值積分的逼近效果題55 求空間曲線L弧長。求導后使用符號函數積分計算題61 用歐拉公式和四階龍格-庫塔法分別求解下列初值問題，代碼搜索內容。題64 常微分方程的解，dsolve()函數使用題82 差分法解常微分方程邊值問題，ode函數無能為力，Matlab中提供bvp解算器。 solinit=bvpinit(x,yinit,params

3、)sol= bvpsolver(odefun,bcfun,solinit,options) 題83 求解圓的半徑，圓心。線性方程組解參數設計總結：（1） 算法是題目的解題核心，好的算法可以使計算更加精確 （題5.1）（2） 圖形繪制在今后的課程設計，或者是論文中可以用到。（3） 無法解決的問題，需要請教室友，或者上網查閱。（4） MATLAB是一個很強大的軟件，提供了很多內置的數學函數，直接進行解題。查閱資料時了解到一些MATLAB論壇。通過帖子閱讀，了解到了MATLAB在科學計算方面，模擬，圖形，視頻等起到的作用。增加了對“計算科學“的理解。指 導 老 師 評 語建議：從學生的工作態度、工作

4、量、設計(論文的)創造性、學術性、使用性及書面表達能力等方面給出評價。 簽名： 20 年 月 日 數值分析課程設計11 水手、猴子和椰子問題：五個水手帶了一只猴子來到南太平洋的一個荒島上，發現那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個水手醒來后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個水手也陸續起來，和第一個水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問原先共有幾只椰子？（15621）試分析椰子數目的變化規律，利用逆向遞推的方法

5、求解這一問題解：算法分析：解該問題主要使用遞推算法，關于椰子數目的變化規律可以設起初的椰子數為，第一至五次猴子在夜里藏椰子后，椰子的數目分別為再設最后每個人分得x個椰子,由題: （k=0,1,2,3,4）所以,利用逆向遞推方法求解 （k=0,1,2,3,4）輸出結果 ： 1023 15621MATLAB代碼：n=input('n= ');n= 15621for x=1:n p=5*x+1; for k=1:5p=5*p/4+1; end結果分析：此題的解題思想是在迭代法中，判斷p為整數時，輸出與p if p=fix(p),breakendenddisp(x,p)12 設，（1）

6、從盡可能精確的近似值出發，利用遞推公式：計算機從到的近似值；（2）從較粗糙的估計值出發，用遞推公式：計算從到的近似值； 解：首先第一步，估計和的值：syms x n; int (x0/(5+x),0,1)ans=log(2)+log(3)-log(5)eval(ans)ans=0.1823則取為0.18syms x n;int(x30/(5+x),0,1)ans =931322574615478515625*log(2)+931322574615478515625*log(3)-931322574615478515625*log(5)-79095966183067699902965545527

7、073/465817912560eval(ans)ans = 0MATLAB中小數點后保留四位，由上面計算知道積分的值不為了零。 所以的取值為0.00001-0.0001MATLAB代碼：i=input('i=');i=0.18;>> if i>=0.1&&i<=0.2 for n=1:1:20 i=-5*i+1/n endelseif i>0&&i<=0.0001 for n=30:-1:2 i=(-1/5)*i+1/(5*n) endend i = 1.1336e+005i = -5.6679e+005i

8、= 2.8339e+006i = -1.4170e+007i = 7.0848e+007i = -3.5424e+008i = 1.7712e+009i = -8.8560e+009i = 4.4280e+010i = -2.2140e+011同理輸入積分初始值i=0時可以得 i=0.0884結果分析：第二種方法所得的結果相對來說比較精確一些，也比較可靠因為第一種方法每一迭代都將最初的誤差放大了五倍，使得最終的誤差越來越大；而第一種方法經過每一次迭代都將誤差縮小為初始誤差的五分之一，使得最終的誤差越來越小，因此相對來說比較可靠，性能較好。13 繪制Koch分形曲線問題描述：從一條直線段開始，將

9、線段中間的三分之一部分用一個等邊三角形的另兩條邊代替，形成具有5個結點的新的圖形（圖1-4）；在新的圖形中，又將圖中每一直線段中間的三分之一部分都用一個等邊三角形的另兩條邊代替，再次形成新的圖形（圖1-5），這時，圖形中共有17個結點。這種迭代繼續進行下去可以形成Koch分形曲線。問題分析：考慮由直線段（2個點）產生第一個圖形（5個點）的過程，設和分別為原始直線段的兩個端點。現在需要在直線段的中間依次插入三個點產生第一次迭代的圖形（圖1-4）。顯然，位于點右端直線段的三分之一處，點繞旋轉60度（逆時針方向）而得到的，故可以處理為向量經正交變換而得到向量，形成算法如下：（1）；（2）；（3）；在

10、算法的第三步中，A為正交矩陣。這一算法將根據初始數據（和點的坐標），產生圖1-4中5個結點的坐標。這5個結點的坐標數組，組成一個5×2矩陣。這一矩陣的第一行為為的坐標，第二行為的坐標，第二行為的坐標第五行為的坐標。矩陣的第一列元素分別為5個結點的x坐標 ，第二列元素分別為5個結點的y坐標。問題思考與實驗：（1）考慮在Koch分形曲線的形成過程中結點數目的變化規律。設第k次迭代產生結點數為，第迭代產生結點數為，試寫出和之間的遞推關系式；（2）參考問題分析中的算法，考慮圖1-4到圖1-5的過程，即由第一次迭代的5個結點的結點坐標數組，產生第二次迭代的17個結點的結點坐標數組的算法；（3）

11、考慮由第k次迭代的個結點的結點坐標數組，產生第次迭代的個結點的結點坐標數組的算法；（4）設計算法用計算機繪制出如下的Koch分形曲線（圖1-6）解：(1) (2)(3)算法及(4)代碼分析：p=0 0;10 0; %P為初始兩個點的坐標,第一列為x坐標,第二列為y坐標n=2; %n為結點數A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %旋轉矩陣for k=1:4 d=diff(p)/3; %diff計算相鄰兩個點的坐標之差,得到相鄰兩點確定的向量 %則d就計算出每個向量長度的三分之一,與題中將線段三等分對應 m=4*n-3; %迭代公式 q=p(1:n

12、-1,:); %以原點為起點,前n-1個點的坐標為終點形成向量 p(5:4:m,:)=p(2:n,:); %迭代后處于4k+1位置上的點的坐標為迭代前的相應坐標 p(2:4:m,:)=q+d; %用向量方法計算迭代后處于4k+2位置上的點的坐標 p(3:4:m,:)=q+d+d*A' %用向量方法計算迭代后處于4k+3位置上的點的坐標 p(4:4:m,:)=q+2*d; %用向量方法計算迭代后處于4k位置上的點的坐標 n=m; %迭代后新的結點數目end plot(p(:,1),p(:,2) %繪出每相鄰兩個點的連線axis(0 10 0 3.5)21 用高斯消元法的消元過程作矩陣分解

13、。設消元過程可將矩陣A化為上三角矩陣U，試求出消元過程所用的乘數、并以如下格式構造下三角矩陣L和上三角矩陣U驗證：矩陣A可以分解為L和U的乘積，即A=LU。矩陣LU分解MATLAB代碼:function hl=zhjLU(A)n,n=size(A);RA=rank(A);if RA=n disp('因為A的n階行列式hl等于零，所以A不能進行LU分解.A的秩RA如下：'); RA,hl=det(A); returnendif RA=n for p=1:n h(p)=det(A(1:p,1:p); end hl=h(1:n); for i=1:n if h(1,i)=0 disp

14、('因為A的各階主子式等不等于零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：'); RA,hl return end end if h(1,i)=0 disp('因為A的各階主子式都不等于零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl如下：'); for j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>j L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1); L(i

15、,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end RA,hl,U,L· endend以上上代碼保存為M文件，并在命令窗口輸入A=20 2 3;1 8 1; 2 -3 15;b=0 0 0'h=zhjLU(A)程序運行結果：L = U= 1.0000 0 0 20.0000 2.0000 3.0000 0.0500 1.0000 0 0 7.9000 0.8500 0.1000 -0.4051 1.0000 0 0 15.0443

16、實驗驗證： 可以直接使用MATLAB內置LU分解A=20 2 3;1 8 1; 2 -3 15;L U=lu(A) 輸出結果與上程序輸出結果一致。22 用矩陣分解方法求上題中A的逆矩陣。記 分別求解方程組 由于三個方程組系數矩陣相同，可以將分解后的矩陣重復使用。對第一個方程組，由于A=LU，所以先求解下三角方程組,再求解上三角方程組，則可得逆矩陣的第一列列向量；類似可解第二、第三方程組，得逆矩陣的第二列列向量的第三列列向量。由三個列向量拼裝可得逆矩陣。解: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

17、;0.05,1,0;0.1,-0.4051,1;U=20 2 3;0 7.9 0.85;0 0 15.0443;Y1=Lb1 X1=UY1 Y2=Lb2X2=UY2Y3=Lb3X3=UY3Y1 = 1.0000 -0.0500 -0.1203X1 = 0.0517 -0.0055 -0.0080Y2 =實驗驗證： X1 X2 X3ans = 0.0517 -0.0164 -0.0257 -0.0055 0.1237 0.1165 -0.0080 0.0269 0.0934而：inv（A）=X1 X2 X3 得證 0 1.0000 0.4051X2 = -0.0164 0.1237 0.0269

18、Y3 = 0 1.0000 1.4051X3 = -0.0257 0.11650.0934- 19 -23 驗證希爾伯特矩陣的病態性：對于三階矩陣取右端向量，驗證：（1）向量是方程組的準確解；（2）取右端向量b的三位有效數字得，求方程組的準確解，并與X的數據作比較 。說明矩陣的病態性。解：(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.8333 1.08330.7833 與題中相同(2) 先求出解X，與數據作比較 。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=

19、HbX = 1.0800 0.54001.4400與相差較大，矩陣為病態矩陣31 用泰勒級數的有限項逼近正弦函數用計算機繪出上面四個函數的圖形。解：MATLAB代碼如下(1) syms x;taylor(sin(x)x=0:0.01*pi:piplot(x,sin(x) (2) syms x;taylor(x)x=0:0.01*pi:pi/2plot(x) (3) syms x;taylor(x-x3/6)fplot('x-x3/6',0 pi/2) (4) syms x;taylor(x-x3/6+x5/120)fplot('x-x3/6+x5/120',0

20、pi/2) 結果圖形右： x=0:0.1:pi;y=sin(x);plot(x,y,'-sk');hold onx=0:0.1:pi/2;y=x;plot(x,y,'-b*')hold onfplot('x-x.3/6',0,pi/2,0,2,2e-3,'-gx')hold onfplot('x-x.3/6+x.5/120',0,pi/2,0,2,2e-3,'-r.')hold offlegend('sin(x)','x','x-x3/6','

21、x-3/6+x5/120',2)xlabel('x')ylabel('y')title('Taylor approximation')結果分析：圖中紅色點線為正弦曲線，藍色的星線為一階泰勒逼近，綠色叉線為二階泰勒逼近，黑色正方形線為三階泰勒逼近。可見三階泰勒逼近效果最好，泰勒級數越高，逼近效果越好。32 繪制飛機的降落曲線一架飛機飛臨北京國際機場上空時，其水平速度為540km/h，飛行高度為1 000m。飛機從距機場指揮塔的橫向距離12 000m處開始降落。根據經驗，一架水平飛行的飛機其降落曲線是一條三次曲線。建立直角坐標系，設飛機著陸點

22、為原點O，降落的飛機為動點，則表示飛機距指揮塔的距離，表示飛機的飛行高度，降落曲線為該函數滿足條件：（1）試利用滿足的條件確定三次多項式中的四個系數；（2）用所求出的三次多項式函數繪制出飛機降落曲線。function s=f(x1,y1,x2,y2,x3,y3,x4,y4)format longa1=1 ,x1,x12,x13;a2=1 ,x2,x22,x23;a3=0,1,2*x3,3*x32;a4=0,1,2*x4,3*x42;a=a1; a2;a3;a4;b=y1; y2;y3;y4;s=ab;x=-12000:250:0;y=s(3)*x.2-s(4)*x.3plot(x,y,'

23、;-d')xlabel('x')ylabel('y')以上為M文件內容，在命令窗口鍵入f(0,0,12000,1000,0,0,12000,0)運行結果： ans為多項式系數ans = 1.0e-004 * 0 0 0.20833333333333 -0.0000115740740741 曾任英特爾公司董事長的摩爾先生早在1965年時，就觀察到一件很有趣的現象：集成電路上可容納的零件數量，每隔一年半左右就會增長一倍，性能也提升一倍。因而發表論文，提出了大大有名的摩爾定律（Moores Law），并預測未來這種增長仍會延續下去。下面數據中，第二行數據為晶片

24、上晶體數目在不同年代與1959年時數目比較的倍數。這些數據是推出摩爾定律的依據：年代19591962196319641965增加倍數13456試從表中數據出發，推導線性擬合的函數表達式。解：解法一MATLAB代碼：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=polyfit(x,y,1)y1=polyval(p1,x)plot(x,y1,'-',x,y,'r*')xlabel('x'),ylabel('y');運行結果：p1 = 1.0e+003 * 0.0008 -1.6255y1 =0.8

25、113 3.3019 4.1321 4.9623 5.7925線性擬合的函數表達式： Y=0.8302x-1.6255e+003解法二：運行結果：xs = 1.0e+003 * -1.625528301891238 0.000830188679248x=1959 1962 1963 1964 1965;y=1;3;4;5;6;for i=1:length(x) for j=1:2 A(i,j)=x(i)(j-1); endendL,U=lu(A'*A);xs=U(L(A'*y)從而年代y與增加倍數x之間的關系為：y=-1625.528301891238+0.8301886792

26、48x51 用幾種不同的方法求積分的值。（1）牛頓-萊布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）復合梯形公式。解：syms xi1=int(4/(1+x2),x,0,1)運行結果：i1 =pii2 = 3i3 = 3.1333i4 = 3.1416a=0;b=1;h=b-a;i2=(4/(1+a2)+4/(1+b2)/2i3=h/6*(4/(1+a2)+4*4/(1+(a+b/2)2)+4/(1+b2)M=100;h=(b-a)/M;i4=0;for k=1:(M-1) x=a+h*k; i4=i4+4/(1+x2);endi4=h*(4/(1+a2)+4/(1+b2)/2+h*i4

27、結果分析：牛頓萊布尼茲公式得到精確結果，采用梯形公式得到的結果比采用Simpson公式的精確度要低，采用復化梯形公式在步長取得越來越小的狀態下可以提高精度。55 求空間曲線L：弧長公式為解：運行程序為：syms t;x=diff(cos(t);y=diff(sin(t);z=diff(2-cos(t)-sin(t);y=int(x2+y2+z2)0.5,t,0,2*pi); %符號函數積分digits(14);i=vpa(y)運行結果：i = 8.737752570989461 用歐拉公式和四階龍格-庫塔法分別求解下列初值問題；解：（1）<1>歐拉公式：function t,x=E

28、uler(fun,t0,tt,x0,N)h=(tt-t0)/N;t=t0+0:N'*h;x(1,:)=x0'運行結果：ans = 0 1.0000 0.0500 1.0450 0.1000 1.0878 0.1500 1.1285 0.2000 1.1676 0.2500 1.2051 0.3000 1.2413 0.3500 1.2762 0.4000 1.3100 0.4500 1.3427 0.5000 1.3745 0.5500 1.4055 0.6000 1.4356 0.6500 1.4649 0.7000 1.4936 0.7500 1.5216 0.8000 1

29、.5490 0.8500 1.5758 0.9000 1.6021 0.9500 1.6278 1.0000 1.6531for k=1:N f=feval(fun,t(k),x(k,:); f=f' x(k+1,:)=x(k,:)+h*f;end以'Euler.m'保存function f=Euler_fun(t,x)f=0.9*x./(1+2*t);以'Euler_fun.m'保存function main_Eulert,x=Euler('Euler_fun',0,1,1,20);fh=dsolve('Dx=0.9*x/(1+

30、2*t)','x(0)=1');for k=1:8 ft(k)=t(k); fx(k)=subs(fh,ft(k);endt,x以'main_Euler.m'保存輸入：main_Euler運行結果：ans = 0 1.0000 0.1250 1.0954 0.2500 1.1807 0.3750 1.2585 0.5000 1.3306 0.6250 1.3981 0.7500 1.4617 0.8750 1.5223 1.0000 1.5801<2>四階龍格-庫塔法function R=rk4(f,a,b,ya,N)%y'=f(x,

31、y)%a,b為左右端點%N為迭代步長%h為步長%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 k1=h*feval(f,T(j),Y(j); k2=h*feval(f,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;endans=T' Y' 以'rk4.m'保存function z=f

32、(x,y)z=0.9*y./(1+2*x);以'f.m'保存輸入：rk4('f',0,1,1,8)64 列出函數在區間0，e上的函數值表并作出它的圖象。其中，是初值問題 的解。解v=dsolve('Dv*log(x)=2*x','v(0)=0','x');f=(1-log(v)*v f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x) ezplot('f')輸出結果：f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)7.4

33、一個10次項式的系數為1 a1 a2 a9 a10=1 55 1320 181 50 157 773 902 055 341 693 0 -840 950 0 127 535 76 -106 286 40 632 880 0試用多項式的求根指令roots求出該10次方程的10個根，然后修改9次項的系數-55為-56，得新的10次方程，求解新的方程，觀察根的變化是否很顯著。解： p=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 10.6051

34、+ 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127ip=1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2