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

數(shù)值 課分 程析 設(shè)計(jì)專業(yè)班級(jí)：信息與計(jì)算科學(xué)09-1班姓 名：陳育偉 學(xué) 號(hào)：20096361 實(shí)驗(yàn)一11 水手、猴子和椰子問題：五個(gè)水手帶了一只猴子來到南太平洋的一個(gè)荒島上，發(fā)現(xiàn)那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個(gè)水手醒來后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個(gè)水手也陸續(xù)起來，和第一個(gè)水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問原先共有幾只椰子？試分析椰子數(shù)目的變化規(guī)律，利用逆向遞推的方法求解這一問題（15621）?！締栴}分析】首先分析問題的解答方法，應(yīng)該采用逆推的方式解決問題。每個(gè)水手起來的椰子數(shù)量等于前一個(gè)水手醒來時(shí)椰子數(shù)量少一的五分之四。最后每個(gè)水手得到的椰子數(shù)量會(huì)等于最后一堆椰子少一的五分之一。根據(jù)這個(gè)逆推，由于椰子的數(shù)量是整數(shù)，所以利用循環(huán)語句知道整數(shù)的解答為止?！境绦?qū)崿F(xiàn)】n=input(input n：)；for x=1：n p=5*x+1； for k=1：5 p=5*p/4+1; end if p= =fix(p)，break，end end disp(x，p)運(yùn)行這段程序后，屏幕出現(xiàn)要求從鍵盤輸入 x 數(shù)據(jù)的信息input n，輸入1200后，MATLAB計(jì)算出合適的 x 和 p0 的值為 1023 15621 【結(jié)果分析】要使得最初的椰子數(shù)p0為整數(shù)，必須取 (x +1) 為 4 5( =1024)的倍數(shù)，一種簡單的處理可取 x = 1023。12設(shè)，（1）從盡可能精確的近似值出發(fā)，利用遞推公式：計(jì)算機(jī)從到的近似值；（2）從較粗糙的估計(jì)值出發(fā)，用遞推公式：計(jì)算從到的近似值；（3）分析所得結(jié)果的可靠性以及出現(xiàn)這種現(xiàn)象的原因【算法分析】 syms x n; int(1/(x+5),0,1) ans log(2)+log(3)-log(5)eval(ans)ans = 0.1823Now we use to compute In where n rang from 1 to 20.s=0.1823for n=1:20s=-5*s+1/n;disp(s)ends = 0.1823 0.0885 0.0575 0.0458 0.0208 0.0958 -0.3125 1.7054 -8.4018 42.1200 -210.5002 1.0526e+003 -5.2629e+003 2.6314e+004 -1.3157e+005 6.5786e+005 -3.2893e+006 1.6447e+007 -8.2233e+007 4.1116e+008 -2.0558e+009（2）首先計(jì)算I（30）的近似值 syms x n; int(x30/(x+5),0,1)eval (ans)0 /即I（30）0s=0;for n=30:-1:2s=-s/5+1/5*n;disp(s)ends = 0 6 4.6000 4.6800 4.4640 4.3072 4.1386 3.9723 3.8055 3.6389 3.4722 3.3056 3.1389 2.9722 2.8056 2.6389 2.4722 2.3056 2.1389 1.9722 1.8056 1.6389 1.4722 1.3056 1.1389 0.9722 0.8056 0.6389 0.4722 0.3056【結(jié)果分析】第一小題當(dāng)中首先算出較為精確的I0，之后經(jīng)過代入的方法計(jì)算出 I（20），結(jié)果是比較準(zhǔn)確的。但是第二小題當(dāng)中先是算出必是很精確的I（30），本來所得誤差比第一步就稍微大一點(diǎn)了，然后再用這個(gè)誤差稍微大的I（30）回代計(jì)算I（20）I（1），而從I（30）算到I（20）的時(shí)候誤差這時(shí)已經(jīng)夠大了，采用這個(gè)誤差更大的結(jié)果去計(jì)算我們想要的值，當(dāng)然誤差大得驚人啦，顯然比第一種誤差大得多。所以我們使用第一種方法得到的結(jié)果較接近準(zhǔn)確值。實(shí)驗(yàn)二21用高斯消元法的消元過程作矩陣分解。設(shè)消元過程可將矩陣A化為上三角矩陣U，試求出消元過程所用的乘數(shù)、并以如下格式構(gòu)造下三角矩陣L和上三角矩陣U驗(yàn)證：矩陣A可以分解為L和U的乘積，即A=LU。【算法分析】(1) The process of Gaussian elimination A= 20, 2,3;1,8,1;2,-3,15;B(1)=0,0,0;-1,-2/20,-3/20;0,0,0;C(1)=A+B(1)得到結(jié)果為：C(1) = 20.0000 2.0000 3.0000 0 7.9000 0.8500 2.0000 -3.0000 15.0000B(2)=0,0,0; 0,0,0; -20/10,-2/10,-3/10;C(2)=C(1)+B(2)得到結(jié)果為：C(2) = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 -3.2000 14.7000B(3)=0,0,0; 0,0,0; 0,0.32,0.85*3.2/7.9;C(3)=B(3)+C(2)得到結(jié)果為：C(3)= 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443So U= 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443And from the Gaussion elimination we have thatL= 1 0 0 1/20 1 0 1/10 -3.2/7.9 1 驗(yàn)證：L=1，0，0；0.05，1，0；0.1，-0.405063，1; U=20，2，3；0，7.9，0.85；0，0，15.0443; A*Bans = 20.0000 2.0000 3.0000 1.0000 8.0000 1.0000 2.0000 -2.9995 14.9998 所以LU = 1 0 0 20.0000 2.0000 3.0000 1/20 1 0 0 7.9000 0.8500 1/10 -3.2/7.9 1 0 0 15.0443 L*U 20 2 3 1 8 1 2 -3 15【結(jié)果分析】因?yàn)樵谶\(yùn)算過程中遇到除不盡的分?jǐn)?shù)采用的是四舍五入的方法來記錄數(shù)據(jù)，所以所得結(jié)果是一個(gè)近似值，但是結(jié)果非常的接近準(zhǔn)確值，所以說這種方法是正確的。22 用矩陣分解方法求上題中A的逆矩陣。記分別求解方程組由于三個(gè)方程組系數(shù)矩陣相同，可以將分解后的矩陣重復(fù)使用。對第一個(gè)方程組，由于A=LU，所以先求解下三角方程組,再求解上三角方程組，則可得逆矩陣的第一列列向量；類似可解第二、第三方程組，得逆矩陣的第二列列向量的第三列列向量。由三個(gè)列向量拼裝可得逆矩陣?！境绦?qū)崿F(xiàn)】A = 20 2 3;1 8 1;2 -3 15;b1 = 1;0;0;b2 = 0;1;0;b3 = 0;0;1;L,U = lu(A);Y = inv(L)*b1;a1 = inv(U)*Y;Y = inv(L)*b2;a2 = inv(U)*Y;Y = inv(L)*b3;a3 = inv(U)*Y;a1 = 0.0517 -0.0055 -0.0080a2 = -0.0164 0.12370.0269a3 = -0.0093 -0.00720.0665【結(jié)果分析】所以A的逆矩陣是0.0517 -0.0164 -0.0093-0.0055 0.1237 -0.0072-0.0080 0.0269 0.0665同時(shí)，a1、a2、a3分別是三個(gè)方程組的解。實(shí)驗(yàn)三31 用泰勒級(jí)數(shù)的有限項(xiàng)逼近正弦函數(shù)用計(jì)算機(jī)繪出上面四個(gè)函數(shù)的圖形?！舅惴ǚ治觥?1) y1=x,x 0,pi/2則算法如下： x=0:.1:pi/2; y1=sin(x); y2=x; plot(x,y1,x,y2)得到的圖形如右圖：（2）y(x)=sin x ,x0,；則y(x)x- x/3！具體算法如下： x=0:.1:pi;y1=sin(x);y2=x-x.3/6;plot(x,y1,x,y2)得到的圖形如右圖：（3）y3=x-x3/6+x5/120,x0，pi/2具體算法如下：x=0:.1:pi/2;y1=sin(x);y2=x-x.3/6+x.5/120;plot(x,y1,x,y2)得到的圖形如右圖：【結(jié)果分析】函數(shù)采用taloy展開項(xiàng)來逼近原函數(shù)，依次采用一階、二階、三階的泰勒公式，顯然有限項(xiàng)越多，逼近效果越好，從圖形就可以很明顯的看出來了。33追趕曲線的計(jì)算機(jī)模擬問題描述：歐洲文藝復(fù)興時(shí)期的著名人物達(dá)芬奇曾經(jīng)提出一個(gè)有趣的“狼追兔子”問題，當(dāng)一只兔子正在它的洞穴南面60碼處覓食時(shí)，一只餓狼出現(xiàn)在兔子正東的100碼處。兔子急忙奔向自己的洞穴，狼立即以快于兔子一倍的速度緊追兔子不放。兔子一旦回到洞穴便逃脫厄，問狼是否會(huì)追趕上兔子？這一問題的研究方法可以推廣到如魚雷追擊潛艇、地對空導(dǎo)彈擊飛機(jī)等問題上去。在對真實(shí)系統(tǒng)做實(shí)驗(yàn)時(shí)，可能時(shí)間太長、費(fèi)用太高、危險(xiǎn)太大、甚至很難進(jìn)行。計(jì)算機(jī)模擬是用計(jì)算機(jī)模仿實(shí)物系統(tǒng)，對系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行動(dòng)態(tài)演示，以評(píng)價(jià)或預(yù)測系統(tǒng)的行為效果。根據(jù)模擬對象的不同特點(diǎn)，分為確定性模擬和隨機(jī)性模擬兩大類。模擬通常所用的是時(shí)間步長法，即按照時(shí)間流逝的順序一步一步對所研究的系統(tǒng)進(jìn)行動(dòng)態(tài)演示，以提取所需要的數(shù)據(jù)。問題分析：首先計(jì)算狼的初始位置到兔子洞穴的直線距離：由于狼奔跑的速度是兔子速度的兩倍，兔子跑60碼的時(shí)間狼可以跑120碼。如果狼沿直線奔向兔窩，應(yīng)該是可以追上兔子的。但是，有人推導(dǎo)出狼在追趕兔子過程中的運(yùn)動(dòng)曲線為根據(jù)曲線方程，當(dāng)時(shí)，。也就是說，在沒有兔窩的情況下兔子一直往北跑，在跑到大約66碼處將被狼追上。由此可知，在有兔窩時(shí)狼是追趕不上兔子的。用計(jì)算機(jī)模擬的方法也可以得到同樣的結(jié)論。取時(shí)間步長為1s，隨時(shí)間步長的增加，考慮這一系統(tǒng)中的各個(gè)元素（狼和兔子）所處的位置變化規(guī)律，用計(jì)算機(jī)作出模擬。最后，根據(jù)第60s時(shí)狼所在的位置的坐標(biāo)，判斷狼是否能追上兔子。問題思考與實(shí)驗(yàn)：（1）設(shè)兔子奔跑的速度為，則狼運(yùn)動(dòng)的速度為。建立平面直角坐標(biāo)系，若當(dāng)時(shí)刻，兔子位于點(diǎn)處，狼位于點(diǎn)處。試根據(jù)，的坐標(biāo)確定一個(gè)單位向量描述狼在時(shí)段內(nèi)的運(yùn)動(dòng)方向。（2）根據(jù)狼的運(yùn)動(dòng)方向和速度推導(dǎo)到的坐標(biāo)的具體表達(dá)式；（3）用計(jì)算機(jī)繪制追趕曲線的圖形（包括靜態(tài)和動(dòng)態(tài)的圖形）?！締栴}求解】 由于P點(diǎn)的運(yùn)動(dòng)方向始終指向Q,設(shè)在t=tk時(shí)刻的位置為Pk(xk,yk) ,則P點(diǎn)的運(yùn)動(dòng)方向可以由下面的單位向量表示Ek=PkQkPkQk = 0-xk(0-xk)2+(tk-yk)2，tk-yk(0-xk)2+(tk-yk)2P點(diǎn)的速度是2m/s，取時(shí)間步長為t ，設(shè)在t=tk+1時(shí)刻，P點(diǎn)的位置為Pk+1(xk+1,yk+1)，于是P點(diǎn)的位置變化規(guī)律為xk+1=xk+2t（0-xk）(0-xk)2+(tk-yk)2yk+1=yk+2t(tk-yk)(0-xk)2+(tk-yk)2【編程實(shí)現(xiàn)】建立M腳本文件 m3_3.mx2(1)=100;y2(1)=0;t(1)=0;tstep=1;for k=1:60 t(k+1)=t(k)+tstep; x2(k+1)=x2(k)+2*(0-x2(k)/sqrt(0-x2(k)2+(t(k)-y2(k)2); y2(k+1)=y2(k)+2*(t(k)-y2(k)/sqrt(0-x2(k)2+(t(k)-y2(k)2); x1(k)=0; y1(k)=t(k);endaxis(0,100,0,60)hold onfor k=1:61 plot(x2(k),y2(k),*,x1(k),y1(k),o) pause(0.5)endhold offd=sqrt(0-x(61)2+(t(61)-y(61)2);fprintf(狼與兔子的距離為%5.4n,d)運(yùn)行程序可以出現(xiàn)動(dòng)態(tài)的圖形去掉語句pause(0.5)后，可以顯示靜態(tài)的圖形如下圖：解決本題時(shí)，使用了計(jì)算機(jī)模擬的方法，動(dòng)態(tài)演示了系統(tǒng)的過程，使得問題的解決大大地簡化了。實(shí)驗(yàn)四43 神經(jīng)元模型用于蠓的分類識(shí)別（MCM1989A題）問題描述：生物學(xué)字試圖對兩類蠓蟲（Af與Apf）進(jìn)行鑒別，依據(jù)的資料是蠓蟲的觸角和翅膀的長度，已經(jīng)測得9只Af和6只Apf的數(shù)據(jù)（觸角長度用x表示，翅膀長度用y表示）Af數(shù)據(jù)x124136138138138140148154156Y127174164182190170182182208Apf數(shù)據(jù)x1.141.181.201.261.281.30Y1.781.961.862.002.001.96現(xiàn)需要解決三個(gè)問題：（1）如何憑借原始資料（15對數(shù)據(jù)，被稱之為學(xué)習(xí)樣本）制定一種方法，正確區(qū)分兩類蠓蟲；（2）依據(jù)確立的方法，對題目提供的三個(gè)樣本：（1.24,1.80）,(1.28,1.84),(1.40,2.04)加以識(shí)別；（3）設(shè)Af是寶貴的傳粉益蟲，Apf是某種疾病的載體，是否應(yīng)該修改分類方法?！締栴}分析】首先畫出15對數(shù)據(jù)的散點(diǎn)圖，其中，Af 用 * 標(biāo)記，Apf 用 標(biāo)記。觀察圖4-5 可以發(fā)現(xiàn)，代表 Af 的點(diǎn)位于圖中偏右，而代表 Apf 的點(diǎn)位于圖中偏左。應(yīng)該存在一條位于兩類點(diǎn)之間的直線 L， 作為 Af 和 Apf 分 界線，這條直線 L 應(yīng)依據(jù)問題所給的數(shù)據(jù)，即學(xué)習(xí)樣本來確定。設(shè)這條直線的方程為。對于平面上任意一點(diǎn) P(x，y)，如果該點(diǎn)在直線上，將其坐標(biāo)代入直線方程則使方程成為恒等式，即使方程左端恒為零；如果點(diǎn) P(x，y) 不在直線上，將其坐標(biāo)代入直線方程，則方程左端不為零。由于 Af 和 Apf 的散點(diǎn)都不在所求的直線上，故將問題所提供的數(shù)據(jù)代入直線方程左端所得到的表達(dá)式的值應(yīng)有大于零或者小于零兩種不同的結(jié)果。【算法分析】建立一個(gè)判別系統(tǒng)，引入判別函數(shù) g(x，y)，當(dāng) (x，y) 代表 Af 類時(shí)，令g(x，y)0，否則 g(x，y)0。 為了對判別系統(tǒng)引入學(xué)習(xí)機(jī)制，在學(xué)習(xí)過程中將兩種不同的狀態(tài)，以“1”和“1”表示。取 則由所給數(shù)據(jù)形成約束條件，這是關(guān)于判別函數(shù)中的三個(gè)待定系數(shù)值 的線性方程組： 此為包括三個(gè)未知數(shù)共15個(gè)方程的超定方程組，可以求方程組的最小二乘解。 下列程序用于求超定方程組的最小二乘解，并繪制出分類邊界曲線的圖解。xy=1.24 1.27;1.36 1.74;1.38 1.64;1.38 1.82;1.38 1.90; 1.40 1.70;1.48 1.82;1.54 1.82;1.56 2.08;1.14 1.78; 1.18 1.96;1.20 1.86;1.26 2.00;1.28 2.00;1.30 1.96; %學(xué)習(xí)樣本數(shù)據(jù) z=1;1;1;1;1;1;1;1;1;-1;-1;-1;-1;-1;-1; x=xy(:,1);y=xy(:,2);x1=x(1:9);y1=y(1:9);x2=x(10:15);y2=y(10:15); plot(x1,y1,*,x2,y2,x),pause %繪制原始數(shù)據(jù)散點(diǎn)圖 X=x y ones(x); A=X*X;B=X*z;w=AB; %求解正規(guī)方程組 a=w(1)/w(2);b=w(3)/w(2); t=1.10:0.02:1.60;u=a*x+b ; %確定分類直線數(shù)據(jù) plot(x1,y1,*,x2,y2,x,t,u) %在散點(diǎn)圖中畫分類直線 運(yùn)行上面程序可求出超定方程組的最小二乘解并畫出分類邊界曲線。為了由所給數(shù)據(jù)用判別函數(shù)判別三個(gè)新蠓蟲的類屬，即當(dāng) 時(shí)，判為 Af 類；當(dāng) 時(shí)，判為 Apf 類。運(yùn)行上面程序后,鍵入下面命令 xx=1.24 1.80 1;1.28 1.84 1;1.40 2.04 1； xx*w plot(t，u，xx(:,1)，xx(:,2)，o) 求得 ans = -0.3877 -0.2384 -0.0235 這說明，所給數(shù)據(jù)反映出三個(gè)蠓蟲均屬于Apf類?！窘Y(jié)果分析】上面的解決問題方案是在學(xué)習(xí)過程中用 +1 和 1 分別代表正數(shù)和負(fù)數(shù)來完成的。這只是一種人為的規(guī)定，并不是一成不變的。當(dāng)Apf是害蟲時(shí)，可以修改超定方程組的右端項(xiàng)中“-1”為“-0.6”，或者將 1 改為其它的負(fù)數(shù)以重新求超定方程組的最小二乘解獲得分類邊界直線的方程。這樣將與前面所求分類邊界直線的方程不一樣，當(dāng)然對新給定的蠓蟲的翅膀和觸角長度數(shù)據(jù)來做判斷其結(jié)果也是不同的。實(shí)驗(yàn)五51 用幾種不同的方法求積分的值。（1）牛頓-萊布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）復(fù)合梯形公式?！境绦?qū)崿F(xiàn)】function t=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f) 由定步長梯形法求得值為 3function q,srm,err = naadapt( fname,a,b,e)srm=zeros(30,6);it=0;done=1;m=1;sr=simpson(fname,a,b,e);srm(1,1:6)=sr;state=it;while (state=it) n=m; for j=n:-1:1 p=j; sr0=srm(p,:); err=sr0(5);e=sr0(6); if e=err state=done; a=sr0(1);b=sr0(2);err=sr0(5);e=sr0(6); c=0.5*(a+b);e2=0.5*e; sr1=simpson(fname,a,c,e2); sr2=simpson(fname,c,b,e2); err=abs(sr0(3)-sr1(3)-sr2(3)/10;if err x,y=ode23(lore,4,5,0)x = 4.00000000000000 4.00000093280783 4.00000559684697 4.00002891704269 4.00014551802130 4.00072852291430 4.00364354737933 4.01821866970447 4.09109428133017 4.19109428133017 4.29109428133017 4.39109428133017 4.49109428133017 4.59109428133017 4.69109428133017 4.79109428133017 4.89109428133017 5.00000000000000y = 1.0e+002 * 0 0.00000080000028 0.00000480001007 0.00002480026893 0.00012480681052 0.00062497072747 0.00312907383321 0.01573209090479 0.08086238316915 0.17628097287779 0.27924859887657 0.39042430504800 0.51052612335396 0.64033652423016 0.78070837535945 0.93257145727939 1.09693958889058 1.29146855032127【結(jié)果分析】第二個(gè)積分的求解問題是通過轉(zhuǎn)化為常微分初值問題來解決的，所得結(jié)果不是一個(gè)近似值，而是一個(gè)精確值，由此可知，積分和微分的求解問題是可以相互轉(zhuǎn)化的，數(shù)學(xué)個(gè)問題之間是相互貫通的。64 列出函數(shù)在區(qū)間0，e上的函數(shù)值表并作出它的圖象。其中，是初值問題的解?！締栴}分析】從題目中看出，In(v)在零點(diǎn)無定義，所以計(jì)算f：dfdx=dfdvdvdx=-Invdvdx=-2xf0=limv0(1-In(v)v=0【程序?qū)崿F(xiàn)】 clear;simplify(dsolve(Dy=2*x/log(y),x) ans = (x2+2*C1)/lambertw(x2+2*C1)*exp(-1) fun=inline(2*x/log(y),x,y) ode23(fun,1,2.7,1)fun = Inline function: fun(x,y) = 2*x/log(y)Warning: Divide by zero.作圖得：實(shí)驗(yàn)七7.4一個(gè)10次項(xiàng)式的系數(shù)為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試用多項(xiàng)式的求根指令roots求出該10次方程的10個(gè)根，然后修改9次項(xiàng)的系數(shù)-55為-56，得新的10次方程，求解新的方程，觀察根的變化是否很顯著。【問題解答】此題應(yīng)用roots(p)函數(shù)。p=1 -551320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;所以roots(p)= 551320.032920469 + 0.00000000000000i-2.24160771698781 + 0.00000000000000i-0.626975644192937 + 1.55805673614711i-0.626975644192937 - 1.55805673614711i1.21700668114928 + 0.00000000000000i0.820510882955133 + 0.840751943017165i0.820510882955133 - 0.840751943017165i0.302305044859346 + 0.994768021247990i0.302305044859346 - 0.994768021247990ip=1 -561320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;所以roots(p)= 561320.032333995 + 0.00