PAGE1《數值分析NumericalAnalysis》課程設計系專：班級：姓名：學號：指導老師：

實驗一1．1水手、猴子和椰子問題：五個水手帶了一只猴子來到南太平洋的一個荒島上，發現那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個水手醒來后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個水手也陸續起來，和第一個水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問原先共有幾只椰子？試分析椰子數目的變化規律，利用逆向遞推的方法求解這一問題（15621）。問題分析：設最初的椰子數為，即第一個水手所處理之前的椰子數，用分別表示五個水手對椰子動了手腳以后剩余的椰子數目，則根據問題有。再用x表示最后每個水手平分得到的椰子數，于是有。所以，利用逆向遞推的方法，有。由于椰子數為一正整數，用任意的x作為初值遞推出的數據不一定是合適的。這里用for循環語句結合break語句來尋找合適的x和，對任意的x遞推計算出，當計算結果為正整數時，結果正確，否則選取另外的x再次重新遞推計算，直到計算出的結果p0為正整數為止。程序設計：n=input('inputn：');forx=1:np=5*x+1;fork=1:5p=5*p/4+1;endifp==fix(p),break,endenddisp([x,p])運行結果：inputn：1221.0e+003*0.12201.8728inputn：1200102315621inputn：45678901023156211．2設，①從盡可能精確的近似值出發，利用遞推公式：計算機從到的近似值；解：先求出I0的精確值，再用遞推公式symsx;I0=int(1/(5+x),x,0,1);I0=vpa(I0);fori=1:20I0=-5*I0+1/iEnd經過20次迭代，最終結果為I0=0.007997523028232163831452112940177②從較粗糙的估計值出發，用遞推公式：計算從到的近似值；解：先求出的值近似精確值，再用遞推公式symsx;I0=int(0.2*x^30,x,0,1);I0=vpa(I0);fori=30:-1:2I0=-0.2*I0+0.2/iEnd經過30次迭代，最終結果為I0=0.088392216030226868941404253296685③分析所得結果的可靠性以及出現這種現象的原因。解：方法2比較可靠，因為方法1每次都將誤差放大5倍，而方法2則將誤差縮小到原來的0.2倍。

實驗二2．1用高斯消元法的消元過程作矩陣分解。設消元過程可將矩陣A化為上三角矩陣U，試求出消元過程所用的乘數、、并以如下格式構造下三角矩陣L和上三角矩陣U驗證：矩陣A可以分解為L和U的乘積，即A=LU。解：算法分析：由于matlab中含有將A分解成一個下三角矩陣和一個上三角矩陣，所以直接調用。程序及其運行結果：A=[2023;181;2-315][L,U]=lu(A)A=20231812-315L=1.0000000.05001.000000.1000-0.40511.0000U=20.00002.00003.000007.90000.85000015.0443矩陣A可以分解為L和U的乘積2.3驗證希爾伯特矩陣的病態性：對于三階矩陣取右端向量，驗證：（1）向量是方程組的準確解；（2）取右端向量b的三位有效數字得，求方程組的準確解，并與X的數據作比較。說明矩陣的病態性。解：(1)直接驗證即可A=[11/21/3;1/21/31/4;1/31/41/5];B=[1;1;1];b=sym(A*B)結果與上述結果相同(2)先求出解，與數據作比較。A=[11/21/3;1/21/31/4;1/31/41/5];b=[1.83;1.08;0.783];B=A\b結果為B=1.0799999999999990.5400000000000041.439999999999997結果誤差變化較大，說明數值的微小變化導致結果的巨大誤差，矩陣A為病態的

實驗三3．1用泰勒級數的有限項逼近正弦函數用計算機繪出上面四個函數的圖形。解：程序：x=0:pi/100:pi;m=0:pi/100:pi/2;y0=sin(x);y1=m;y2=m-m.^3/6;y3=m-m.^3/6+t.^5/120;plot(x,y0,m,y1,m,y2,m,y3)運行結果：

3．2繪制飛機的降落曲線一架飛機飛臨北京國際機場上空時，其水平速度為540km/h，飛行高度為1000m。飛機從距機場指揮塔的橫向距離12000m處開始降落。根據經驗，一架水平飛行的飛機其降落曲線是一條三次曲線。建立直角坐標系，設飛機著陸點為原點O，降落的飛機為動點，則表示飛機距指揮塔的距離，表示飛機的飛行高度，降落曲線為該函數滿足條件：（1）試利用滿足的條件確定三次多項式中的四個系數；（2）用所求出的三次多項式函數繪制出飛機降落曲線。解：（1）程序：A=[1000;11200012000^212000^3;0100;01240003*(12000^2)];b=[0,1000,0,0]';inv(A)*b運行結果：ans=1.0e-004*000.2083-0.0000（2）程序：x=0:50:12000;y=0.207e-004*(x.^2)+-0.00001158e-004*(x.^3);plot(x,y)運行結果：實驗四4．1曾任英特爾公司董事長的摩爾先生早在1965年時，就觀察到一件很有趣的現象：集成電路上可容納的零件數量，每隔一年半左右就會增長一倍，性能也提升一倍。因而發表論文，提出了大大有名的摩爾定律（Moore’sLaw），并預測未來這種增長仍會延續下去。下面數據中，第二行數據為晶片上晶體數目在不同年代與1959年時數目比較的倍數。這些數據是推出摩爾定律的依據：年代19591962196319641965增加倍數13456試從表中數據出發，推導線性擬合的函數表達式。程序設計：y=[13456];x2=[19591962196319641965];P=polyfit(x2,y,1);disp(P);運行結果為：1.0e+003*0.0008-1.6255結果分析：由運行結果可知線性擬合的函數表達式為y=0.0008x-1.6255，由于是近似值，所以不免有誤差。實驗五5．1用幾種不同的方法求積分的值。（1）牛頓-萊布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）復合梯形公式。解：（1）牛頓-萊布尼茨公式輸入：symsx;y=4/(1+x^2);int_y=int(y,x,0,1)輸出：int_y=pi（2）梯形公式functiony=dada(x)y=4./(1+x.^2);以'dada.m'保存輸入：T=(1-0)/2*(dada(1)+dada(0))結果：T=3（3）辛卜生公式輸入：S=quad('dada',0,1)結果：S=3.141592682924567（4）復合梯形公式輸入d=1/100;x=0:d:1;f=dada(x);FT=trapz(f)*d結果：FT=3.1415759869231295．2設X為標準正態隨機變量，即X~N（0，1）。現分別取，試設計算法計算30個不同的概率值；，并將計算結果與概率論教科書中的標準正態分布函數表作比較。（提示：解：m=0.1:0.1:3;P=zeros(30,1);fori=1:30P(i)=1/sqrt(2*pi)*int(sym('exp(-x^2/2)'),m(i),4);endPm=0.1:0.1:3;P=zeros(30,1);fori=1:30P(i)=1/sqrt(2*pi)*int(sym('exp(-x^2/2)'),-4,m(i));endPP=0.46010.42070.38210.34450.30850.27420.24190.21180.18400.15860.13560.11500.09680.08070.06680.05480.04450.03590.02870.02270.01780.01390.01070.00820.00620.00460.00340.00250.00180.0013P=0.53980.57920.61790.65540.69140.72570.75800.78810.81590.84130.86430.88490.90320.91920.93320.94520.95540.96400.97130.97720.98210.98610.98920.99180.99380.99530.99650.99740.99810.99865．5求空間曲線L：弧長公式為解：運行程序為：symst;a=diff(cos(t));b=diff(sin(t));c=diff(2-cos(t)-sin(t));y=int((a^2+b^2+c^2)^0.5,t,0,2*pi);vpa(y)運行結果：ans=8.737752570984804741619351957708

實驗六6．1用歐拉公式和四階龍格-庫塔法分別求解下列初值問題；解：（1）程序：歐拉公式：function[t,x]=Euler(fun,t0,tt,x0,N)h=(tt-t0)/N;t=t0+[0:N]'*h;x(1,:)=x0';fork=1:Nf=feval(fun,t(k),x(k,:));f=f';x(k+1,:)=x(k,:)+h*f;end以'Euler.m'保存functionf=Euler_fun(t,x)f=0.9*x./(1+2*t);以'Euler_fun.m'保存functionmain_Euler[t,x]=Euler('Euler_fun',0,1,1,20);fh=dsolve('Dx=0.9*x/(1+2*t)','x(0)=1');fork=1:21ft(k)=t(k);fx(k)=subs(fh,ft(k));end[t,x]以'main_Euler.m'保存輸入：main_Eulerans=01.00000.05001.04500.10001.08780.15001.12850.20001.16760.25001.20510.30001.24130.35001.27620.40001.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四階龍格-庫塔法functionR=rk4(f,a,b,ya,N)h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;Y(1)=ya;forj=1:Nk1=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'保存functionz=f(x,y)z=0.9*y./(1+2*x);以'f.m'保存輸入：rk4('f',0,1,1,20)結果：ans=01.00000.05001.03920.10001.07650.15001.11210.20001.14630.25001.17910.30001.21080.35001.24150.40001.27120.45001.30000.50001.32810.55001.35540.60001.38210.65001.40820.70001.43370.75001.45860.80001.48310.85001.50710.90001.53060.95001.55381.00001.5765（2）程序：歐拉公式function[t,x]=Eulers(fun,t0,tt,x0,N)h=(tt-t0)/N;t=t0+[0:N]'*h;x(1,:)=x0';fork=1:Nf=feval(fun,t(k),x(k,:));f=f';x(k+1,:)=x(k,:)+h*f;end以'Eulers.m'保存functionf=Euler_funs(t,x)f=-t*x/(1+t^2);以'Euler_funs.m'保存functionmains_Euler[t,x]=Euler('Euler_funs',0,1,2,20);fh=dsolve('Dx=0.9*x/(1+2*t)','x(0)=1');fork=1:21ft(k)=t(k);fx(k)=subs(fh,ft(k));end[t,x]以'mains_Euler.m'保存輸入：mains_Euler結果：ans=02.00000.05002.00000.10001.99500.15001.98510.20001.97060.25001.95160.30001.92870.35001.90210.40001.87250.45001.84020.50001.80580.55001.76960.60001.73230.65001.69410.70001.65540.75001.61650.80001.57770.85001.53920.90001.50120.95001.46391.00001.4274四階龍格-庫塔法functionR=rk4s(f,a,b,ya,N)h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;Y(1)=ya;forj=1:Nk1=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']以'rk4s.m'保存functionz=fs(x,y)z=-x*y/(1+x^2);以'fs.m'保存輸入：rk4s('fs',0,1,2,20)結果：ans=02.00000.05001.99180.10001.97950.15001.96320.20001.94330.25001.92000.30001.89360.35001.86450.40001.83310.45001.79980.50001.76500.55001.72910.60001.69230.65001.65510.70001.61760.75001.58010.80001.54290.85001.50600.90001.46970.95001.43401.00001.39906.2用求解常微分方程初值問題的方法計算積分上限函數的值，實際上是將上面表達式兩端求導化為常微分方程形式，并用初值條件。試用MATLAB中的指令ode23解決定積分的計算問題。解：odefun=@(x,y)2*pi*exp(x)/x;%定義函數tspan=[45];%求解區間y0=0;%初值[x,y]=ode23(odefun,tspan,y0)x=4.00004.00004.00004.00004.00014.00074.00364.01824.09114.19114.29114.39114.49114.59114.69114.79114.89115.0000y=00.00010.00050.00250.01250.06250.31291.57328.086217.628127.924939.042451.052664.033778.070893.2571109.6940129.1469當x=5時，y=129.1469

實驗七7.1用二分法求下列方程在指定區間[a，b]上的實根近似值：（要求誤差不超過0.01）(1)x-sinx-1=0，[a,b]=[1,3]；(2)xsinx=1，[a,b]=[0,1.5].解；（1）運行程序：a=1;b=3;k=0; y1=a-sin(a)-1; while(b-a)>0.01 k=k+1; x0=(a+b)/2;x(k)=x0; y0=x0-sin(x0)-1;y(k)=y0; ify0*y1<0 b=x0; else a=x0; y1=y0; end end disp([x'y'])運行結果：[2,408488074749405/4503599627370496][3/2,-4481036472577419/9007199254740992][7/4,-1053779023151395/4503599627370496][15/8,-712341393175443/9007199254740992][31/16,8975041611503/2251799813685248][61/32,-85593294866561/2251799813685248][123/64,-154268929443443/9007199254740992][247/128,-7430218095237/1125899906842624][495/256,-11835035444387/9007199254740992][991/512,93879030099/70368744177664][1981/1024,10840787111/1125899906842624][3961/2048,-5875158237047/9007199254740992][7923/4096,-723616715653/2251799813685248][15847/8192,-701966501541/4503599627370496][31695/16384,-164654758197/2251799813685248]運行程序：a=0.1;b=1.5;k=0;y1=a*sin(a)-1;while(b-a)>0.01k=k+1;x0=(a+b)/2;x(k)=x0;y0=x0*sin(x0)-1;y(k)=y0;ify0*y1<0b=x0;elsea=x0;y1=y0;endenddisp([x'y'])運行結果：[4/5,-3838103856873717/9007199254740992][23/20,55933053762297/1125899906842624][39/40,-1738305320249353/9007199254740992][17/16,-323478390340705/4503599627370496][177/160,-98943562990273/9007199254740992][361/320,21826383383843/1125899906842624][143/128,18951300402293/4503599627370496][1423/1280,-7626403754673/2251799813685248][495/256,-11835035444387/9007199254740992][991/512,93879030099/70368744177664][1981/1024,10840787111/1125899906842624][3961/2048,-5875158237047/9007199254740992][7923/4096,-723616715653/2251799813685248][15847/8192,-701966501541/4503599627370496][31695/16384,-164654758197/2251799813685248]7.2方程xsinx=1的根實際上是兩個函數y1(x)=sinx,y2(x)=的交點，用計算機繪出兩個函數在區間[0.2，6]的圖形如下，觀察圖形，分析它們的交點分布規律及特點，試寫出方程的全部實根所在的隔根區間，并給出每一根的近似值。用計算機繪出區間[-6，0.2]上兩個函數的圖形，觀察交點所在位置驗證你的分析結果。解：1，運行程序：function[x,k]=ntf(fun,dfun,x0,eps,N)ifnargin<5N=500;endifnargin<4eps=1e-5;endk=1;whilek<Nx=x0-feval(fun,x0)/feval(dfun,x0);ifabs(x-x0)<epsbreak;endx0=x;k=k+1;endifk==Nwarning('已迭代次數上限');endformatlongfun=inline('x*sin(x)-1');dfun=inline('sin(x)+x*cos(x)');[x,k]=ntf(fun,dfun,1,1e-8)formatlongfun=inline('x*sin(x)-1');dfun=inline('sin(x)+x*cos(x)');[x,k]=ntf(fun,dfun,1,1e-8)運行結果：（1）x=1.114157140871930k=（2）x=2.772604708265991k=52，畫圖驗證運行程序：x=[-6:0.01*pi:0.2];plot(x,sin(x))holdonsymstfplot('1/t',[-6,0.2-51])實驗八8.1分別用直接法、雅可比迭代法、賽德爾迭代法求解線性方程組AX=b。（1），（2），解：（1）a=ones(9);c=4*ones(10);A=diag(diag(c))+diag(diag(a),1)+diag(diag(a),-1);b=[2;1;1;1;1;1;1;1;1;2];B=[A,b]；直接法A\bans=0.4792719919110210.0829120323559150.1890798786653190.1607684529828110.1678463094034380.1678463094034380.1607684529828110.1890798786653190.0829120323559150.479271991911021Jacobi法Jocobi.m文件functionJacobi(A,b,X0,P,error,max1)[nn]=size(A);X=zeros(n,1);fork=1:max1forj=1:nX(j)=(b(j)-A(j,[1:j-1,j+1:n])*X0([1:j-1,j+1:n]))/A(j,j);endXerrX=norm(X-X0,P);X0=X;X1=A\b;if(errX<error)disp('gvgfvgvg')kX1Xreturnendendif(errX>=error)disp('fgjsdj')endx0=[0;0;0;0;0;0;0;0;0;0];Jacobi(A,b,x0,inf,0.001,100);X=0.4793395996093750.0831260681152340.1892738342285160.1610832214355470.1681365966796880.1681365966796880.1610832214355470.1892738342285160.0831260681152340.479339599609375Gauss-Seidel法GS.m文件functionGS(A,b,X0,p,error,max1)[nn]=size(A);X=zeros(n,1);fork=1:max1forj=1:nXX=0;fori=1:nifi<jXX=XX+A(j,i)*X(i);endifi>jXX=XX+A(j,i)*X0(i);endendX(j)=(b(j)-XX)/A(j,j);endXerrX=norm(X-X0,p);X0=X;X1=A\b;if(errX<error)disp('fgdfgh')kX1Xreturnendendif(errX>=error)disp('fghdgh')endx0=[0;0;0;0;0;0;0;0;0;0];GS(A,b,x0,inf,0.001,100）；X=0.4792585372924800.0828957557678220.1890988051891330.1606427729129790.1680229734629390.1677014031447470.1608548912918200.1890394594811370.0829268395536930.479268290111577（2）b=[1;zeros(19,1)