南郵MATLAB數學實驗答案全.doc

《南郵MATLAB數學實驗答案全.doc》由會員分享，可在線閱讀，更多相關《南郵MATLAB數學實驗答案全.doc（30頁珍藏版）》請在裝配圖網上搜索。

1、第一次練習教學要求：熟練掌握Matlab軟件的基本命令和操作，會作二維、三維幾何圖形，能夠用Matlab軟件解決微積分、線性代數與解析幾何中的計算問題。補充命令vpa(x,n)顯示x的n位有效數字，教材102頁fplot(f(x),a,b)函數作圖命令，畫出f(x)在區(qū)間a,b上的圖形在下面的題目中為你的學號的后3位（1-9班）或4位（10班以上）1.1 計算與 syms xlimit(902*x-sin(902*x)/x3)ans =366935404/3limit(902*x-sin(902*x)/x3,inf)ans =01.2 ，求 syms xdiff(exp(x)*cos(902*

2、x/1000),2)ans = (46599*cos(451*x)/500)*exp(x)/250000 - (451*sin(451*x)/500)*exp(x)/2501.3 計算dblquad(x,y) exp(x.2+y.2),0,1,0,1)ans =2.13941.4 計算 syms xint(x4/(9022+4*x2) ans = (91733851*atan(x/451)/4 - (203401*x)/4 + x3/121.5 syms xdiff(exp(x)*cos(902*x),10)ans =-356485076957717053044344387763*cos(90

3、2*x)*exp(x)-3952323024277642494822005884*sin(902*x)*exp(x)1.6 給出在的泰勒展式(最高次冪為4). syms xtaylor(sqrt(902/1000+x),5,x) ans = -(9765625*451(1/2)*500(1/2)*x4)/82743933602 +(15625*451(1/2)*500(1/2)*x3)/91733851 -(125*451(1/2)*500(1/2)*x2)/406802 + (451(1/2)*500(1/2)*x)/902 +(451(1/2)*500(1/2)/5001.7 Fibona

4、cci數列的定義是用循環(huán)語句編程給出該數列的前20項（要求將結果用向量的形式給出）。x=1,1;for n=3:20 x(n)=x(n-1)+x(n-2);endxx=Columns 1 through 10 1 1 2 3 5 8 13 21 34 55Columns 11 through 20 89 144 233 377 610 987 1597 2584 4181 67651.8 對矩陣，求該矩陣的逆矩陣，特征值，特征向量，行列式，計算，并求矩陣（是對角矩陣），使得。A=-2,1,1;0,2,0;-4,1,902/1000;inv(A)ans =0.4107 0.0223 -0.455

5、4 0 0.5000 0 1.8215 -0.4554 -0.9107eig(A)ans =-0.5490 + 1.3764i -0.5490 - 1.3764i 2.0000 det(A)ans =4.3920P,D=eig(A)P = %特征向量0.3245 - 0.3078i 0.3245 + 0.3078i 0.2425 0 0 0.9701 0.8944 0.8944 0.0000 D = -0.5490 + 1.3764i 0 0 0 -0.5490 - 1.3764i 0 0 0 2.0000 P*D6*inv(P) %A6的值ans =15.3661 12.1585 + 0.0

6、000i -5.8531 0 64.0000 0 23.4124 -5.8531 + 0.0000i -1.6196 1.9 作出如下函數的圖形（注：先用M文件定義函數，再用fplot進行函數作圖）：m文件: function y=fenduan(x) if x=1/2 y=2*xelse x syms n A=sym(4,2;1,3);x=1;2;P,D=eig(A) %沒有sym下面的矩陣就會顯示為小數P = -1, 2 1, 1 D = 2, 0 0, 5 An=P*Dn*inv(P) An = 2n/3 + (2*5n)/3, (2*5n)/3 - (2*2n)/3 5n/3 - 2n

7、/3, (2*2n)/3 + 5n/3 xn=An*x xn = 2*5n - 2n 2n + 5n3.2 對于練習1中的，求出的通項. syms n A=sym(2/5,1/5;1/10,3/10); x=1;2;P,D=eig(A) P = -1, 2 1, 1 D = 1/5, 0 0, 1/2 An=P*Dn*inv(P) An = (2*(1/2)n)/3 + (1/5)n/3, (2*(1/2)n)/3 - (2*(1/5)n)/3 (1/2)n/3 - (1/5)n/3, (1/2)n/3 + (2*(1/5)n)/3xn = 2*(1/2)n - (1/5)n (1/2)n +

8、 (1/5)n3.3 對隨機給出的，觀察數列.該數列有極限嗎? A=4,2;1,3;a=;x=2*rand(2,1)-1;for i=1:20 a(i,1:2)=x; x=A*x; end for i=1:20 if a(i,1)=0 else t=a(i,2)/a(i,1); fprintf(%g,%gn,i,t); endend 結論：在迭代17次后，發(fā)現數列存在極限為0.53.4 對120頁中的例子，繼續(xù)計算.觀察及的極限是否存在. （120頁練習9） A=2.1,3.4,-1.2,2.3;0.8,-0.3,4.1,2.8;2.3,7.9,-1.5,1.4;3.5,7.2,1.7,-9.

9、0;x0=1;2;3;4;x=A*x0;for i=1:1:100a=max(x);b=min(x);m=a*(abs(a)abs(b)+b*(abs(a) A=2.1,3.4,-1.2,2.3;0.8,-0.3,4.1,2.8;2.3,7.9,-1.5,1.4;3.5,7.2,1.7,-9.0;P,D=eig(A)P = -0.3779 -0.8848 -0.0832 -0.3908 -0.5367 0.3575 -0.2786 0.4777 -0.6473 0.2988 0.1092 -0.7442 -0.3874 -0.0015 0.9505 0.2555D = 7.2300 0 0 0

10、 0 1.1352 0 0 0 0 -11.2213 0 0 0 0 -5.8439結論：A的絕對值最大特征值等于上面的的極限相等，為什么呢？還有，P的第三列也就是-11.2213對應的特征向量和上題求解到的y也有系數關系，兩者都是-11.2213的特征向量。3.6 設，對問題2求出若干天之后的天氣狀態(tài)，并找出其特點（取4位有效數字）. （122頁練習12） A2=3/4,1/2,1/4;1/8,1/4,1/2;1/8,1/4,1/4;P=0.5;0.25;0.25;for i=1:1:20 P(:,i+1)=A2*P(:,i);endPP = Columns 1 through 10 0.5

11、000 0.5625 0.5938 0.6035 0.6069 0.6081 0.6085 0.6086 0.6087 0.6087 0.2500 0.2500 0.2266 0.2207 0.2185 0.2178 0.2175 0.2174 0.2174 0.2174 0.2500 0.1875 0.1797 0.1758 0.1746 0.1741 0.1740 0.1739 0.1739 0.1739 Columns 11 through 200.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.

12、2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 Column 21 0.6087 0.2174 0.1739結論：9天后，天氣狀態(tài)趨于穩(wěn)定P*=（0.6087,0.2174,0.1739）T3.7 對于問題2，求出矩陣的特征值與特征向量，并將特征向量與上一題中的結論作對比. （122頁練習14） A2=3/4,1/2,1/4;1/8,1/4,1/2;1/8,1/4,1/

13、4; P,D=eig(A2)P = -0.9094 -0.8069 0.3437 -0.3248 0.5116 -0.8133 -0.2598 0.2953 0.4695D = 1.0000 0 0 0 0.3415 0 0 0 -0.0915分析：事實上，q=k(-0.9094, -0.3248, -0.2598)T均為特征向量，而上題中P*的3個分量之和為1，可令k(-0.9094, -0.3248, -0.2598)T=1，得k=-0.6696.有q=(0.6087, 0.2174, 0.1739),與P*一致。3.8對問題1，設為的兩個線性無關的特征向量，若，具體求出上述的，將表示成的

14、線性組合，求的具體表達式，并求時的極限，與已知結論作比較. （123頁練習16） A=3/4,7/18;1/4,11/18;P,D=eig(A);syms k pk;a=solve(u*P(1,1)+v*P(1,2)-1/2,u*P(2,1)+v*P(2,2)-1/2,u,v);pk=a.u*D(1,1).k*P(:,1)+a.v*D(2,2).k*P(:,2) pk = -5/46*(13/36)k+14/23 5/46*(13/36)k+9/23或者：p0=1/2;1/2;P,D=eig(sym(A);B=inv(sym(P)*p0 B = 5/46 9/23syms kpk=B(1,1)

15、*D(1,1).k*P(:,1)+B(2,1)*D(2,2).k*P(:,2) pk = -5/46*(13/36)k+14/23 5/46*(13/36)k+9/23 vpa(limit(pk,k,100),10) ans = .6086956522 .3913043478結論：和用練習12中用迭代的方法求得的結果是一樣的。第四次練習教學要求：會利用軟件求勾股數，并且能夠分析勾股數之間的關系。會解簡單的近似計算問題。4.1 求滿足，的所有勾股數，能否類似于（11.8），把它們用一個公式表示出來？解法程序1：for b=1:998 a=sqrt(b+2)2-b2); if(a=floor(a)

16、 fprintf(a=%i,b=%i,c=%in,a,b,b+2) endend運行結果：a=4,b=3,c=5a=6,b=8,c=10a=8,b=15,c=17a=10,b=24,c=26a=12,b=35,c=37a=14,b=48,c=50a=16,b=63,c=65a=18,b=80,c=82a=20,b=99,c=101a=22,b=120,c=122a=24,b=143,c=145a=26,b=168,c=170a=28,b=195,c=197a=30,b=224,c=226a=32,b=255,c=257a=34,b=288,c=290a=36,b=323,c=325a=38,b

17、=360,c=362a=40,b=399,c=401a=42,b=440,c=442a=44,b=483,c=485a=46,b=528,c=530a=48,b=575,c=577a=50,b=624,c=626a=52,b=675,c=677a=54,b=728,c=730a=56,b=783,c=785a=58,b=840,c=842a=60,b=899,c=901a=62,b=960,c=962解法程序2： n=0;m=;for a=1:100 for c=a+1:1000 b=sqrt(c2-a2); if (b=floor(b)&(ba)&(c-b)=2) n=n+1; m(:,l)

18、=a,b,c; end endendm勾股數，的解是： 以下是推導過程：由，有顯然，從而是2的倍數.設，代入上式得到：因為，從而.4.2 將上一題中改為,分別找出所有的勾股數.將它們與時的結果進行比較，然后用公式表達其結果。（1）時通項：a=8,b=6,c=10a=12,b=16,c=20a=16,b=30,c=34a=20,b=48,c=52a=24,b=70,c=74a=28,b=96,c=100a=32,b=126,c=130a=36,b=160,c=164a=40,b=198,c=202a=44,b=240,c=244a=48,b=286,c=290a=52,b=336,c=340a=

19、56,b=390,c=394a=60,b=448,c=452a=64,b=510,c=514a=68,b=576,c=580a=72,b=646,c=650a=76,b=720,c=724a=80,b=798,c=802a=84,b=880,c=884a=88,b=966,c=970（2）5時通項： a=15,b=20,c=25a=25,b=60,c=65a=35,b=120,c=125a=45,b=200,c=205a=55,b=300,c=305a=65,b=420,c=425a=75,b=560,c=565a=85,b=720,c=725a=95,b=900,c=905（3）6時通項a=

20、12,b=9,c=15a=18,b=24,c=30a=24,b=45,c=51a=30,b=72,c=78a=36,b=105,c=111a=42,b=144,c=150a=48,b=189,c=195a=54,b=240,c=246a=60,b=297,c=303a=66,b=360,c=366a=72,b=429,c=435a=78,b=504,c=510a=84,b=585,c=591a=90,b=672,c=678a=96,b=765,c=771a=102,b=864,c=870a=108,b=969,c=975（4）7時通項a=21,b=28,c=35a=35,b=84,c=91a=

21、49,b=168,c=175a=63,b=280,c=287a=77,b=420,c=427a=91,b=588,c=595a=105,b=784,c=791綜上：當c-b=k為奇數時，通項當c-b=k為偶數時，通項4.3 對，（），對哪些存在本原勾股數？（140頁練習12）程序：for k=1:200 for b=1:999 a=sqrt(b+k)2-b2); if(a=floor(a)&gcd(gcd(a,b),(b+k)=1) fprintf(%i,k); break; end endend運行結果：1,2,8,9,18,25,32,49,50,72,81,98,121,128,162,

22、169,200,4.4 設方程(11.15)的解構成數列,觀察數列,.你能得到哪些等式？試根據這些等式推導出關于的遞推關系式. （142頁練習20）解：1000以內解構成的數列 , , , 如下： n 1 2 3 4 5 6 2 7 26 97 362 1351 1 4 15 56 209 780 3 11 41 153 571 2131 4 15 56 209 780 2911 1 3 11 41 153 571我們發(fā)現這些解的關系似乎是：=因為=，所以。有以下結論： （4.1）可以看成一個線性映射,令,（4.1）可寫成：4.5 選取對隨機的，根據的概率求出的近似值。（取自130頁練習7）提

23、示：(1)最大公約數的命令:gcd(a,b)(2)randint(1,1,u,v)產生一個在u,v區(qū)間上的隨機整數程序： m=90200;s=0;for i=1:m a=randint(1,2,1,109); if gcd(a(1),a(2)=1; s=s+1; endendpi=sqrt(6*m/s)pi =3.14314.6 用求定積分的Monte Carlo法近似計算。（102頁練習16）提示：Monte Carlo法近似計算的一個例子。對于第一象限的正方形,內畫出四分之一個圓向該正方形區(qū)域內隨即投點,則點落在扇形區(qū)域內的概率為.投次點,落在扇形內的次數為,則,因此.程序如下n=100000;nc=0;for i=1:n x=rand;y=rand; if(x2+y2=1) nc=nc+1; endendpi=4*nc/n解：程序：a=0;b=1;m=1000;H=1;s=0;for i=1:m xi=rand(); yi=H*rand(); if yisqrt(1-xi2); s=s+1; endendpi=4*H*(b-a)*s/m運行結果：pi = 3.1480四、通過本課程學習，談談你開設對這門課的認識，對教學以及上機實驗提出自己的和建議。