計(jì)算機(jī)仿真技術(shù)實(shí)驗(yàn)指導(dǎo)書

1、計(jì)算機(jī)仿真技術(shù)實(shí)驗(yàn)指導(dǎo)書黃劍航 編 莆田學(xué)院電子信息工程學(xué)系2011年8月71 / 75文檔可自由編輯打印 前言計(jì)算機(jī)仿真技術(shù)是利用計(jì)算機(jī)科學(xué)和技術(shù)的成果建立被仿真對象的模型，并在某些實(shí)驗(yàn)條件下對模型進(jìn)行動(dòng)態(tài)實(shí)驗(yàn)的一門綜合性技術(shù)。它具有高效、安全、受環(huán)境條件的約束較少、可改變時(shí)間比例尺等優(yōu)點(diǎn)，已成為分析、設(shè)計(jì)、運(yùn)行、評價(jià)、培訓(xùn)系統(tǒng)（尤其是復(fù)雜系統(tǒng)）的重要工具，是電子信息本科專業(yè)的專業(yè)必修課。計(jì)算機(jī)仿真離不開基本的專業(yè)理論知識和計(jì)算機(jī)編程，因此學(xué)習(xí)本門課程可以使同學(xué)們鞏固相關(guān)專業(yè)領(lǐng)域的基本知識、原理和方法，訓(xùn)練編程能力，掌握仿真的基本概念，并將仿真技術(shù)應(yīng)用在專業(yè)領(lǐng)域的分析，問題解決上面，為今后進(jìn)

2、一步的學(xué)習(xí)、科研打下基礎(chǔ)。本課程注重實(shí)踐和理論結(jié)合，在實(shí)驗(yàn)室上機(jī)完成全部課程。結(jié)合課程和專業(yè)特點(diǎn)，實(shí)驗(yàn)指導(dǎo)書選取了MATLAB語言基礎(chǔ)、基本算法仿真、趣味算法仿真、信號與系統(tǒng)仿真和Simulink仿真這幾部分主要內(nèi)容，作為學(xué)習(xí)仿真技術(shù)的主要內(nèi)容。希望同學(xué)們在實(shí)驗(yàn)課前要做好預(yù)習(xí)工作，掌握基本方法原理，實(shí)驗(yàn)課堂上做好每個(gè)實(shí)驗(yàn)，對實(shí)驗(yàn)結(jié)果要做好必要的記錄和保存，課后認(rèn)真完成實(shí)驗(yàn)報(bào)告。只要認(rèn)真對待每一次實(shí)驗(yàn)練習(xí)，相信同學(xué)們在專業(yè)理論知識和編程仿真實(shí)踐上真正會(huì)有所收獲。目錄第1部分 MATLAB基礎(chǔ). 1實(shí)驗(yàn)1.1 MATLAB入門 . 1實(shí)驗(yàn)1.2 MATLAB數(shù)組及矩陣運(yùn)算 . . . .3實(shí)驗(yàn)1.

3、3 MATLAB編程應(yīng)用.6實(shí)驗(yàn)1.4 MATLAB繪圖. .10實(shí)驗(yàn)1.5 MATLAB匿名函數(shù). 11第2部分 基本算法實(shí)現(xiàn).16實(shí)驗(yàn)2.1 數(shù)值積分算法.16實(shí)驗(yàn)2.2 插值與擬合.18第3部分 趣味算法實(shí)現(xiàn).23實(shí)驗(yàn)3.1 Hill密碼.23實(shí)驗(yàn)3.2 隨機(jī)模擬 .28第4部分 信號與系統(tǒng)仿真.34實(shí)驗(yàn)4.1 連續(xù)時(shí)間信號在MATLAB中的表示. 34實(shí)驗(yàn)4.2 連續(xù)時(shí)間信號在MATLAB中的運(yùn)算.40實(shí)驗(yàn)4.3 傅里葉變換及其性質(zhì).45實(shí)驗(yàn)4.4 連續(xù)時(shí)間LTI系統(tǒng)的時(shí)域分析.51第5部分 Simulink建模仿真.57實(shí)驗(yàn)5.1 Simulink仿真基礎(chǔ).57實(shí)驗(yàn)5.2 連續(xù)系統(tǒng)建

4、模仿真.61第1部分 MATLAB基礎(chǔ)實(shí)驗(yàn)1.1 MATLAB入門1. 實(shí)驗(yàn)?zāi)康氖煜ATLAB工作環(huán)境和基本操作。2. 實(shí)驗(yàn)內(nèi)容2.1熟悉MATLAB工作環(huán)境，掌握clc，clear，cd這幾個(gè)常見命令。 圖1 MATLAB桌面，具體布局可能因機(jī)器和版本的不同而會(huì)有變化2.2 假設(shè)x=3，y=4。利用MTLAB 計(jì)算下列表達(dá)式： ， 要求：（1）在命令窗口直接輸入，得到運(yùn)算結(jié)果；（2）編寫M文檔并在命令窗口執(zhí)行，然后用工作空間平臺查看有哪些變量在當(dāng)前工作區(qū)。2.3 一小球從空中下落的位移公式為： 利用MATLAB 計(jì)算小球在t=5s 時(shí)的位置，已知：. 2.4 通過以下兩種方式得到關(guān)于exp

5、 函數(shù)的幫助 ：（1）在命令窗口中輸入help exp 命令。（2）運(yùn)用幫助空間窗口。2.5 計(jì)算函數(shù)在為如下取值時(shí)候的函數(shù)值，取4.9 3.2, 100, 1.5, 9.75, 2.56，16， 4.9， 10。3. 思考題3.1 MATLAB 的命令窗口的作用是什么？編輯/調(diào)試窗口？圖像窗口？ 3.2 列出幾種不同的得到MATLAB 幫助的方法。 3.3什么是工作區(qū)？在同一工作區(qū)內(nèi)，你如何決定它里面存儲了什么？ 3.4 你怎樣清空MATLAB 工作區(qū)內(nèi)的內(nèi)容？ 實(shí)驗(yàn)1.2 MATLAB數(shù)組及矩陣運(yùn)算1. 實(shí)驗(yàn)?zāi)康氖煜ATLAB數(shù)組及矩陣運(yùn)算。 2. 實(shí)驗(yàn)內(nèi)容2.1 確定下列數(shù)組

6、的大小，通過whos 或工作空間窗口（The workspace browser）檢查你 的答案。注意在本練習(xí)中后面的數(shù)組可能要用到前面數(shù)組的定義。 (1) u=10 20*i 10+20(2) v=-1;20;3(3) w=1 0 -9;2 -2 0;1 2 3(4) x=u' v(5) y(3,3)=-7(6) z=zeros(4,1) ones(4,1) zeros(1,4)'(7) v(4)=x(2,1) 2.2 w(2,1)的值是多少？  2.3 x(2,1)的值是多少？  2.4 y(2,1)的值是多少? 2.5 當(dāng)語句3(7

7、)執(zhí)行后，v(3)的值是多少？  2.6 c 數(shù)組的定義如下，寫出下面子數(shù)組的內(nèi)容。 c =1.1000   -3.2000   3.4000   0.60000.6000   1.1000    -0.6000   3.10001.3000    0.6000    5.5000     0 (1) c(2,:)  (2) c(:,en

8、d)  (3)  c(1:2,2:end)  (4) c(6)(5) c(4:end) (6) c(1:2,2:4) (7) c(1 4,2) (8) c(2 2,3 3) 2.7 當(dāng)賦值語句執(zhí)行后，下列數(shù)組的內(nèi)容是多少？ (1) a=1 2 3; 4 5 6; 7 8 9;a(3 1,:)=a(1 3,:); (2) a=1 2 3; 4 5 6; 7 8 9;a(1 3,:)=a(2 2,:); (3) a=1 2 3; 4 5 6; 7 8 9;a=a(2 2,:); 2.8 當(dāng)數(shù)組執(zhí)行后，下列數(shù)組a 的內(nèi)容是多少？

9、(1) a=eye(3,3);b=1 2 3;a(2,:)=b; (2) a=eye(3,3);b=4 5 6 ;a(:,3)=b' (3) a=eye(3,3);b=7 8 9;a(3,:)=b(3 1 2); (4) a=eye(3,3);b=7 8 9;a(3,:)=b(3 1 2); 2.9 假設(shè)a,b,c 和d 的定義如下： 分別運(yùn)行出下列表達(dá)式的運(yùn)算結(jié)果，并思考。 (1) a + b  (2) a .* c  (3) a * b  (4) a * c(5) a + c  (6) a + d  (

10、7) a .* d  (8) a * d3. 思考題3.1 數(shù)組,矩陣,向量有什么區(qū)別?3.2 回答關(guān)于下列矩陣的有關(guān)問題 (1) C 的大小是多少？ (2) C(2,3)的值是多少？ (3) 列出值為0.6 的元素的下標(biāo)實(shí)驗(yàn)1.3 MATLAB編程應(yīng)用1. 實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)函數(shù)的編制，掌握MATLAB的實(shí)際應(yīng)用。通過不同的程序結(jié)構(gòu)和不同的實(shí)際應(yīng)用問題，掌握MATLAB的編程方法。2. 實(shí)驗(yàn)內(nèi)容2.1  畫出函數(shù) 在區(qū)間0,10的波形。  用MATLAB 編輯器創(chuàng)建一個(gè)新的M 文件，把上面的語句寫入這個(gè)文件并命名為 ： test1.m。然后在命令窗口中輸入test1

11、執(zhí)行這個(gè)文件。看得到什么結(jié)果？2.2 編寫一個(gè)MATLAB 程序，要求輸入圓柱體的半徑和高，然后輸出體積。2.3 一個(gè)程序?qū)嵗龑W(xué)習(xí)： （溫度轉(zhuǎn)換）設(shè)計(jì)一個(gè)MATLAB 程序，讀取一個(gè)華氏溫度的輸入，輸出開爾文溫度。 華氏溫度和開爾文溫度的轉(zhuǎn)換關(guān)系式可在物理學(xué)課本中找到。其關(guān)系式為： 在物理學(xué)參考書中舉了一些例子，我們可以用來檢驗(yàn)我們程序是否正確。例如 我們設(shè)計(jì)程序的步驟如下 ：（1）提示用戶鍵入華氏溫度值 （2）讀取輸入值 （3）通過關(guān)系式轉(zhuǎn)換為開氏溫度 （4）輸出結(jié)果，結(jié)束 我們將會(huì)用input 函數(shù)輸入華氏溫度，用fprintf 函數(shù)輸出結(jié)果。 % Script % Purpose:% T

12、o convert an input temperature from degrees Fahrenheit to% an output temperature in kelvins.% Record of revisions:% Date Programmer Description of change% = = =% 12/01/97 S.J.Chapman Original code%Define variables:% temp_f -Temperature in degrees Fahrenheit% temp_k -Temperature in kelvins%Prompt the

13、 user for the input temperature.temp_f=input('Enter the temperature in degrees Fahrenheit:');%Converttokelvins.temp_k=(5/9)*(temp_f-32)+273.15;%Writeouttheresult.fprintf('%6.2f degrees Fahrenheit = %6.2f kelvins.n',.temp_f,temp_k);我們輸入上面的例子中的華氏溫度值，以檢測程序的正確性。注意用戶的輸入值已用黑 體字標(biāo)出。 >>

14、 temp_conversionEnter the temperature in degrees Fahrenheit:212212.00 degrees Fahrenheit = 373.15 kelvins.>> temp_conversionEnter the temperature in degrees Fahrenheit:-110-110.00 degrees Fahrenheit = 194.26 kelvins.這個(gè)結(jié)果和物理教科書的結(jié)果相同。  2.4 編寫一個(gè)函數(shù)文件，計(jì)算出坐標(biāo)系中用戶指定兩點(diǎn)(X1,Y1)和(X2,Y2)之間的距離。要求有輸入、輸出

15、及其相關(guān)提示。2.5 雙曲余弦的定義如下： ，編寫一個(gè)程序，計(jì)算出用戶指定的x 的值對應(yīng)的雙曲余弦值。用這個(gè)程序計(jì)算雙曲余弦值的若干值，并和MATLAB 中的內(nèi)建函數(shù)cosh(x)得到的值比較看看是否完全相同。并用MATLAB 打印出這個(gè)函數(shù)的圖象。2.6 電子工程：負(fù)載的最大輸出功率一個(gè)內(nèi)阻Rs50，電動(dòng)勢V120V 的電源驅(qū)動(dòng)一個(gè)負(fù)載RL。當(dāng)RL 為多少時(shí)，RL 的功率最大？在這種情況下，功率為多少？畫以RL 為自變量的RL 功率圖。2.7 利用公式求的近似值，直到最后一項(xiàng)的絕對值小于為止。2.8 Fibonacci(斐波納契）序列的元素滿足Fibonacci 規(guī)則： 且；現(xiàn)要求該序列中第

16、一個(gè)大于20000 的元素，并指明該元素是序列的第幾項(xiàng)。 2.9 在田徑比賽中，一個(gè)身高為一米八零的鉛球運(yùn)動(dòng)員，大概以多大的角度（和水平方向夾角）推鉛球，才能使鉛球推得最遠(yuǎn)，并求出最遠(yuǎn)距離。不計(jì)空氣阻力，假設(shè)鉛球出手點(diǎn)和運(yùn)動(dòng)員高度相等，且鉛球出手瞬間初始速度大小為14m/s，重力加速度取g=10m/s2。2.10 打印出所有的“水仙花數(shù)”，所謂“水仙花數(shù)”是指一個(gè)三位數(shù)，其各位數(shù)字立方和等于該數(shù)本身。例如153是一個(gè)水仙花數(shù)，滿足15313+53+33。2.11 有一個(gè)函數(shù)：，寫一程序，輸入x，輸出y值。2.12 輸入3個(gè)整數(shù)，要求按大小輸出。3. 思考題3.1 MATLAB 啟動(dòng)時(shí)

17、，當(dāng)前工作路徑是什么？在任意路徑下新建一個(gè)名字為mynewdir的新文件夾， 把當(dāng)前工作路徑改為mynewdir，然后打開一個(gè)編輯窗口，增加以下語句：  % create an input array from -2*pi to 2*pit = -2*pi:pi/10:2*pi;% calculate|sin(t)|x=abs(sin(t);%plot resultplot(t,x); 把此文件以test2.m 為文件名保存，然后在命令窗口中輸入test2 執(zhí)行此文件，觀察結(jié)果。關(guān)閉畫圖窗口，然后返回到原始工作路徑，再在命令窗口中輸入test2?？从泻吻闆r發(fā)生，為什么？3.

18、2 MATLAB中的函數(shù)文件和腳本文件有什么區(qū)別？實(shí)驗(yàn)1.4 MATLAB繪圖1. 實(shí)驗(yàn)?zāi)康恼莆誐ATLAB的幾種繪圖方法和相關(guān)命令，包括plot，meshgird等命令。掌握結(jié)果圖的一些參數(shù)設(shè)置（如坐標(biāo)范圍，坐標(biāo)軸標(biāo)示，曲線類型等等）。2. 實(shí)驗(yàn)內(nèi)容2.1 繪制正切函數(shù) 在0x1 上的圖象，并為x 軸和y 軸添加標(biāo)簽。增量設(shè)為0.1。2.2 顯示同樣的圖象，把sinx添加到第二個(gè)圖形中作為第二條曲線。2.3 使用plot3 函數(shù)繪制曲線、和的圖象，不要添加坐標(biāo)軸標(biāo)簽，但要打開網(wǎng)格。2.4 繪制以參數(shù)方程表示的空間曲線。2.5 繪制以極坐標(biāo)方程表示的心臟線。2.6 分別在，的條件下使用命令fp

19、lot做出函數(shù)的圖形。2.7 試在矩形區(qū)域，上分別繪制函數(shù)與對應(yīng)的三維網(wǎng)格表面圖和三維曲面圖。3. 思考題3.1 執(zhí)行完 x,y=meshgrid(-10:0.5:10)之后的x、y值是多少？3.2 說明在計(jì)算諸如 z=x2+y2 這樣的f(x,y)函數(shù)時(shí)，必須先使用meshgrid函數(shù)的原因。實(shí)驗(yàn)1.5 MATLAB匿名函數(shù)1. 實(shí)驗(yàn)?zāi)康腗ATLAB7.0以后引入了一種新的函數(shù)類型-匿名函數(shù)(anonymous function)。匿名函數(shù)不僅完全實(shí)現(xiàn)了內(nèi)聯(lián)函數(shù)(inline function)的所有功能，而且一些用內(nèi)聯(lián)函數(shù)實(shí)現(xiàn)起來不大方便或者不簡潔的功能，用匿名函數(shù)也非常容易實(shí)現(xiàn)。除此之外

20、，匿名函數(shù)的調(diào)用效率要遠(yuǎn)高于內(nèi)聯(lián)函數(shù)。本實(shí)驗(yàn)的主要目的就是熟悉匿名函數(shù)的用法。2. 匿名函數(shù)的定義2.1 匿名函數(shù)的基本定義在程序編寫過程中，經(jīng)常需要計(jì)算諸如，等這樣一些數(shù)學(xué)表達(dá)式的值，如何方便地構(gòu)造這些函數(shù)？寫成M文件再調(diào)用肯定可以，但是從簡潔以及管理的方便性考慮這不是好辦法。這里引入了匿名函數(shù)，基本定義如下：Fhandle=(arglist) expr其中，expr是具體的函數(shù)表達(dá)式，arglist是指定的函數(shù)的自變量。具體見下面示例：f=(x) x.2;fx=f(1:10) fx = 1 4 9 16 25 36 49 64 81 100g=(x,y) x.2+y.2;gxy=g(1:1

21、0,2:11) gxy = 5 13 25 41 61 85 113 145 181 2212.2 匿名函數(shù)的種類這里按照自變量的個(gè)數(shù)以及層數(shù)可以分為以下幾種：單變量匿名函數(shù)、多變量匿名、單重匿名函數(shù)以及多重匿名函數(shù)。1.單變量匿名函數(shù)單變量匿名函數(shù)是最簡單的匿名函數(shù)，只含有一個(gè)變量，如：f=(x) x.2就是單變量匿名函數(shù)，除此之外，含有參數(shù)，參數(shù)值已知的單個(gè)自變量的匿名函數(shù)也是單變量匿名函數(shù)，比如：a=10;b=15;f=(x) a*x+b;f(1:5) ans = 25 35 45 55 652.多變量匿名函數(shù)多變量匿名函數(shù)含有兩個(gè)或者兩個(gè)以上的自變量，如g=(x,y) x.2+y.2。

22、同單變量匿名函數(shù)一樣，多變量匿名函數(shù)也可以有自己的參數(shù)值。比如：a=1;b=2;g=(x,y) a*x+y.b;g(1:5,1:5) ans = 2 6 12 20 303.單重匿名函數(shù)到目前為止，上面列的匿名函數(shù)，無論單變量還是多變量匿名函數(shù)，都屬于單重匿名函數(shù)。這類匿名函數(shù)的特點(diǎn)是：只有一個(gè)符號引導(dǎo)，符號之后就是具體的函數(shù)表達(dá)式。自變量輸入單重匿名函數(shù)后，得到的就是具體的數(shù)值結(jié)果。除了單重匿名函數(shù)外，還有二重、多重匿名函數(shù)。這些多重匿名函數(shù)在參數(shù)傳遞方面非常方便。4.多重匿名函數(shù)下面以二重匿名函數(shù)為例來介紹多重匿名函數(shù)。例子如下：f=(a,b) (x) a*x+b f = (a,b)(x)

23、a*x+b其中a，b是外層變量，x是內(nèi)層變量?？梢赃@樣理解這個(gè)表達(dá)式：每個(gè)符號后面括號里的變量作用域一直到表達(dá)式結(jié)尾。這樣a，b的作用域就是(x) a*x+b，而x的作用域就是 a*x+b。因此，若有個(gè)表達(dá)式為fab=f(a,b)，意思是這是一個(gè)在給定的參數(shù)a，b情況下，以x為變量的單層的單變量匿名函數(shù)?？梢岳胒unctions函數(shù)觀察建立的匿名函數(shù)的信息，如下：f=(a,b) (x) a*x+b f = (a,b)(x)a*x+b f23=f(2,3) f23 = (x)a*x+b3. 應(yīng)用實(shí)例3.1 匿名函數(shù)可以非常方便地表示所求方程，并供fzero等求解函數(shù)調(diào)用。求下列方程的根：求解代

24、碼如下：>> f=(x) exp(x)+x2+x(sqrt(x)-100f = (x)exp(x)+x2+x(sqrt(x)-100 >> x0=fzero(f,3)x0 = 4.1635>>f(x0) ans = 2.8422e-0143.2 顯式表示下列y關(guān)于x的隱函數(shù)：，利用匿名函數(shù)，可以在MATLAB中顯式寫出y和x的關(guān)系式如下：y=(x) fzero(y) (exp(y)+xy)(1/y)-x2*y,1)這樣，對于任意的x，只需調(diào)用，就能求得相應(yīng)的y值，如：>>y1=y(1)y1 = 2.7779>>y2=y(2)y2 =

25、1.10553.3 顯式表示下列z關(guān)于x，y的隱函數(shù)：借助匿名函數(shù)，可以在MATLAB中顯式寫出z關(guān)于x，y的關(guān)系式如下：>> z=(x,y) fzero(z) z-sin(z*x-0.5)2+x*2*y2-z/10)*.exp(-(x-0.5-exp(-y+z)2+y2-z/5+3),rand);其中，fzero函數(shù)求解該隱函數(shù)的初值為隨機(jī)值rand。利用上述匿名函數(shù)，可以畫出z關(guān)于x，y的圖形，如圖1所示：X,Y=meshgrid(-1:0.1:7,-2:0.1:2);Z=arrayfun(x,y) z(x,y),X,Y);surf(X,Y,Z)xlabel('font

26、size15fontnametimes new romanx','color','b')ylabel('fontsize15fontnametimes new romany','color','b')zlabel('fontsize15fontnametimes new romanz','color','b')title('fontsize15fontnametimes new romanz(x,y)','color','

27、;b') 圖1 實(shí)例3.3結(jié)果圖3.4 求，k取值0,5區(qū)間不同值時(shí)的函數(shù)值，并畫出函數(shù)圖像。分別用inline內(nèi)聯(lián)函數(shù)類型和匿名函數(shù)類型實(shí)現(xiàn)。該例對于一些復(fù)雜的、無顯式積分表達(dá)式的帶參數(shù)積分問題具有通用性。解決代碼如下：tic;k=linspace(0,5);y1=zeros(size(k);for i=1:length(k) kk=k(i); fun=inline('sin(',num2str(kk),'*x).*x.2'); y1(i)=quadl(fun,0,5);endtime=toc;disp('Inline method costs

28、 time:',num2str(time),'seconds!'); tic;f=(k) quadl(x) sin(k.*x).*x.2,0,5);kk=linspace(0,5);y2=zeros(size(kk);for ii=1:length(kk) y2(ii)=f(kk(ii);endtime=toc;disp('Anonymous function method costs time:',num2str(time),'seconds!'); plot(kk,y2);xlabel('k');ylabel('

29、;f(k)');grid on;得到的運(yùn)行結(jié)果圖如下圖2所示：圖2 實(shí)例3.4結(jié)果圖第2部分 基本算法實(shí)現(xiàn)實(shí)驗(yàn)2.1 數(shù)值積分算法 1. 實(shí)驗(yàn)?zāi)康臑榱搜芯窟B續(xù)系統(tǒng)的性質(zhì)，需要建立連續(xù)系統(tǒng)的數(shù)學(xué)模型，然后進(jìn)行仿真。連續(xù)系統(tǒng)通?？梢杂贸Ｎ⒎址匠虂砻枋?。對于連續(xù)系統(tǒng)的仿真，就是用計(jì)算機(jī)求解這些數(shù)學(xué)模型，也就是對微分方程求數(shù)值解。各種數(shù)值積分算法在常微分方程的求解中起到了關(guān)鍵作用。通過算法仿真實(shí)驗(yàn)，達(dá)到對數(shù)值積分算法本質(zhì)更好的理解和應(yīng)用的目的。2. 算法原理歐拉法，改進(jìn)的歐拉法，龍格-庫塔法，線性多步法等，詳見課本。3. 實(shí)驗(yàn)內(nèi)容3.1 設(shè)一微分方程為，初始條件。試編寫一個(gè)程序用歐拉法求其數(shù)值

30、解。3.2 已知方程的解析解為，且時(shí)，。取計(jì)算步長為，試用歐拉法，改進(jìn)歐拉法和四階龍格-庫塔法求取值是05之間的數(shù)值結(jié)果，并和解析解做一比較。3.3 對于線性微分方程，分別用改進(jìn)歐拉法和龍格-庫塔法進(jìn)行求解，然后將結(jié)果與用MATLAB中現(xiàn)成函數(shù)求得的結(jié)果進(jìn)行比較。3.4 設(shè)連續(xù)系統(tǒng)的數(shù)學(xué)模型用下列微分方程描述：試計(jì)算當(dāng)時(shí)，的時(shí)間響應(yīng)。4. 思考題4.1 對于高階微分方程的數(shù)值積分求解，如何進(jìn)行預(yù)處理，轉(zhuǎn)換成為適合歐拉法或者龍格-庫塔法求解的一階形式？4.2 已知一階系統(tǒng)的微分方程為，初始條件，取仿真步長，分別用歐拉法，改進(jìn)歐拉法和龍格-庫塔法計(jì)算該系統(tǒng)仿真第一步的值。實(shí)驗(yàn)2.2 插值與擬合 1

31、. 實(shí)驗(yàn)?zāi)康睦斫獠逯蹬c擬合各自的基本定義，算法原理。實(shí)現(xiàn)算法在實(shí)際中的具體應(yīng)用。掌握相關(guān)的MATLAB函數(shù)用法2. 算法原理2.1 插值方法原理在一些實(shí)際應(yīng)用中，會(huì)得到一些數(shù)據(jù)，通常將這些數(shù)據(jù)列成數(shù)據(jù)表格的形式。這種用數(shù)據(jù)表格形式給出的函數(shù)稱為列表函數(shù)，即，其中點(diǎn)成為結(jié)點(diǎn)。根據(jù)函數(shù)已有的數(shù)據(jù)表格來計(jì)算函數(shù)在一些新的點(diǎn)處的函數(shù)值，這就是插值所要解決的問題。一種最簡單的插值方法就是拉格朗日插值，利用拉格朗日多項(xiàng)式即可計(jì)算相應(yīng)區(qū)間任意點(diǎn)的插值結(jié)果。候拉格朗日插值多項(xiàng)式也被稱為拉格朗日公式，如下： 拉格朗日公式常用于以下兩種簡單的情形。（1）線性插值設(shè)函數(shù)在處的函數(shù)值分別為，則線性插值公式為利用線性插

32、值公式就可以計(jì)算在之間任意一個(gè)新點(diǎn)x處的近似值。（2）二次（拋物線）插值設(shè)函數(shù)在處的函數(shù)值分別為，則二次插值公式為如果這三點(diǎn)不在一條直線上，則上述插值公式得到的是一個(gè)二次函數(shù)，通過這三點(diǎn)的曲線是拋物線。利用三點(diǎn)的二次插值公式就可以計(jì)算在所在區(qū)間的任意一個(gè)新點(diǎn)x處的近似值。2.2 擬合方法原理給定平面上的一組數(shù)據(jù)，要求確定一個(gè)初等函數(shù)，使給定的數(shù)據(jù)點(diǎn)（在某種距離下）與函數(shù)表示的曲線整體上最為接近，這類問題稱為數(shù)據(jù)擬合。設(shè)有經(jīng)驗(yàn)公式（這里x和c均可為向量），要求根據(jù)數(shù)據(jù)，確定參數(shù)c，這樣的問題稱為曲線擬合，其基本原理是最小二乘法，即求使誤差平方和達(dá)到最小的c。MATLAB中有關(guān)擬合的命令如下：p=

33、polyfit(x,y,k) 用k次多項(xiàng)式擬合向量數(shù)據(jù)(x,y)，返回多項(xiàng)式的降冪系數(shù)向量。yi=polyval(p,xi) 根據(jù)多項(xiàng)式p計(jì)算出在xi處插值函數(shù)值yi.c=lsqnonlin(fun,c0) 使用迭代法搜索最優(yōu)參數(shù)c，其中fun是以參數(shù)c（可以是向量）為自變量的表示誤差 (x,y是數(shù)據(jù)向量)的函數(shù)，c0是c的近似值，作為迭代初始值。c=lsqcurvefit(fun,c0,x,y) 根據(jù)數(shù)據(jù)向量x,y 確定經(jīng)驗(yàn)公式（c可以是向量）中的參數(shù)c。這里fun為函數(shù)，c0是c的初始值。與lsqnonlin不同，使用lsqcurvefit時(shí)數(shù)據(jù)x,y從外部輸入。3. 實(shí)例分析3.1 以下

34、是美國人口兩個(gè)世紀(jì)以來的統(tǒng)計(jì)數(shù)據(jù)，試以此預(yù)測2010年，2020年的美國人口。年份180018101820183018401850186018701880人口/百萬5.37.29.612.917.123.231.438.650.2年份189019001910192019301940195019601970人口/百萬62.976.092.0106.5123.2131.7150.7179.3204.0年份198019902000人口/百萬226.5251.4275.0解：先做散點(diǎn)圖，為了方便，也為了避免數(shù)據(jù)相差太大造成誤差，年份均減去1800.t=0:10:200;y=5.3 7.2 9.6 12

35、.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0. 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 275.0;subplot(2,2,1);plot(t,y,'bo');axis(0,300,0,400);text(20,350,'Population');從數(shù)值上看給定數(shù)據(jù)接近拋物線或者指數(shù)曲線。圖1 美國人口散點(diǎn)圖及預(yù)測結(jié)果圖（1）考慮用二次多項(xiàng)式擬合，代碼如下：p=polyfit(t,y,2)Q1=sqrt(sum(y-polyval(p,t).2)ti=0:2:220;

36、yi=polyval(p,ti);yj=polyval(p,210,220)subplot(2,2,2)plot(t,y,'bo',ti,yi,'r',210,220,yj,'ro')text(20,350,'Bipoly');結(jié)果：擬合曲線見圖1右上角，均方誤差=12.5425，預(yù)測人口=303.0567 331.8160.（2）采用馬爾薩斯模型，人口按指數(shù)增長，因此可設(shè)人口曲線為，其中N0, r是待擬合的參數(shù)，相關(guān)代碼如下：fun1=(c,t) c(1)*exp(c(2)*t);c= lsqcurvefit(fun1,6,0.

37、1,t,y)Q2=sqrt(sum(y-c(1)*exp(c(2)*t).2)ti=0:2:220;yi=c(1)*exp(c(2)*ti);yj=c(1)*exp(c(2)*210,220)subplot(2,2,3);plot(t,y,'bo',ti,yi,'r',210,220,yj,'ro')text(20,350,'Exp curve');結(jié)果：擬合曲線見圖1左下角，均方誤差=48.5696，預(yù)測人口=338.3433 389.2984.（3）采用Logistic模型，人口增長規(guī)律為 ，其中Nm,N0,r是待擬合的參數(shù)，

38、相關(guān)代碼如下：fun2=(c,t) c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*t);c=lsqcurvefit(fun2,500,6,0.2,t,y)Q3=sqrt(sum(y-c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*t).2)ti=0:2:220; yi=c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*ti);yj=c(1)./(1+(c(1)/c(2)-1)*exp(-c(3)*210,220)subplot(2,2,4)plot(t,y,'bo',ti,yi,'r',210,220,yj,

39、'ro')text(20,350,'Logistic');結(jié)果：擬合曲線見圖1右下角，均方誤差=19.8133，預(yù)測人口=292.9388 311.8921從以上結(jié)果看，二次多項(xiàng)式擬合均方誤差最小，Logistic模型次之，馬爾薩斯模型最差，但從機(jī)理上考慮，二次多項(xiàng)式和馬爾薩斯模型確定的人口最終會(huì)趨于無窮大，不合理，它們都只能用作近期預(yù)測。相對來說，Logistic模型較為合理。.4. 實(shí)驗(yàn)內(nèi)容4.1 設(shè)函數(shù)在直角坐標(biāo)平面內(nèi)四個(gè)點(diǎn)(0,0), (0,1), (1,0), (1,1)的所對應(yīng)的函數(shù)值分別為1，3，2，5，試用兩次線性插值（雙線性插值）求坐標(biāo)為(0.

40、2,0.4)該點(diǎn)處的函數(shù)值。4.2 設(shè)x0=0.4,x2=0.8,而x1可取0.5,0.6,0.7中任何一點(diǎn)，根據(jù)下表使用二次拉格朗日插值求sin0.6389的近似值時(shí)，問如何選節(jié)點(diǎn)x1，可使其近似值的誤差較?。吭囉?jì)算說明。x0.40.50.60.70.8sinx0.38940.47940.56460.64420.71744.3 用電壓U=10v的電池給電容充電，電容上時(shí)刻t的電壓為：，其中V0是電容的初始電壓，是充電常數(shù)，試由下面一組t，U數(shù)據(jù)確定U0和。t/s0.51234579U/V6.366.487.268.228.668.999.439.635. 思考題5.1 說明雙線性插值的原理。

41、5.2 結(jié)合例子說明參數(shù)擬合時(shí)候如何選擇參數(shù)的初始值。第3部分 趣味算法實(shí)現(xiàn)實(shí)驗(yàn)3.1 Hill密碼 1. 實(shí)驗(yàn)?zāi)康牧私饧用芎徒饷茉?，掌握模運(yùn)算下的矩陣運(yùn)算以及數(shù)字與字符串之間的轉(zhuǎn)換。2. 算法原理保密通信具有悠久的歷史，現(xiàn)在更是被廣泛應(yīng)用于軍事、經(jīng)濟(jì)、商業(yè)等各行各業(yè)。在保密通信中獎(jiǎng)原信息稱為明碼，加密后的信息稱為密碼。如果不知道加密方法，一般人無法知道明文內(nèi)容，這樣就起到了保密的作用。加密和解密的過程可以抽象為一個(gè)數(shù)學(xué)模型。移位加密法 移位加密法是一種簡單的加密方法。它通過將明文的字母按字母表中的次序平移若干位實(shí)現(xiàn)加密。如加密方法是字母平移5位，則明文字母和密文字母對應(yīng)關(guān)系如下：明文字母

42、ABCDEFGHIJKLMNOPQRSTUVWXYZ密文字母 FGHIJKLMNOPQRSTUVWXYZABCDE明文THEOLYMPICGAMES 加密后就成了 YMJTQDRUNHLFRJX. 一般人不知道它是什么意思，這就起到了加密的作用。如果知道它是由字母表平移5位得來的，那就很容易獲得原文。數(shù)字5是解開密碼的一把鑰匙，稱為密鑰。在移位加密法中，明文字母和密文字母之間的對應(yīng)關(guān)系是固定的，這種加密方法可以通過分析詞頻，利用統(tǒng)計(jì)方法破譯。下面要介紹的Hill密碼是用矩陣運(yùn)算實(shí)現(xiàn)加密的，它不保持明文字母和密文字母之間固定的對應(yīng)關(guān)系，破譯起來比移位法困難。模n運(yùn)算 在模n運(yùn)算下，參與運(yùn)算的只有

43、0，1，2，n-1這n個(gè)元素。規(guī)定兩個(gè)元素的和（積）是按普通加法（乘法）運(yùn)算的結(jié)果減去或加上n的某個(gè)倍數(shù)后得到的在0和n-1之間的數(shù)。例如在模15運(yùn)算下4+13=2，4*13=7。通常記做4+13=2(mod 15),4*13=7(mod 15)。在模n運(yùn)算下，元素m的加法負(fù)元是n-m。如果m*k=1(mod n)，則稱k是m的乘法逆元。如在模15運(yùn)算下，13的逆元是7，4的逆元還是4。值得注意的是在模15運(yùn)算下，除了0沒有逆元外，3，5及3，5的倍數(shù)6，9，10，12也沒用逆元?？梢宰C明：若m和n的最大公因子等于1，即gcd(m,n)=1，則在模n運(yùn)算下元素m存在唯一的逆元。特別，當(dāng)n是素?cái)?shù)

44、時(shí)，元素1,2,n-1都有逆元，此時(shí)數(shù)0,1,2,n-1在模n下可以進(jìn)行加減乘除(0除外)四則運(yùn)算，它們構(gòu)成一個(gè)數(shù)域，稱為有限域。模n下矩陣加法、乘法以及數(shù)乘矩陣運(yùn)算與模n下數(shù)字的加法、乘法運(yùn)算相似，先對矩陣做普通加法、乘法和數(shù)乘矩陣運(yùn)算，然后所有數(shù)字以n取模即可。定義 兩個(gè)方陣A、B如滿足AB=BA=E(mod n)，則稱矩陣A模n可逆，B稱為A的逆矩陣，記做(mod n)。矩陣A的逆矩陣可以利用伴隨矩陣，按求得。Hill加密法（1） 先將英文字母變換成數(shù)字；A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 1 2 3 4 5 6 7 8

45、 9 10 24 25 0稱其為字母表的表值。（2） 將明文按n個(gè)字母分組，并用對應(yīng)的數(shù)字取代字母，構(gòu)成一個(gè)個(gè)n維向量。（3） 取一個(gè)在模n運(yùn)算下可逆的n階矩陣A，用A左乘（2）中的向量得到新向量；把新向量中的數(shù)字按（1）中的對應(yīng)關(guān)系換回成字母即得到密文。（4） 按例中的方法求出可逆矩陣A在模n下的逆矩陣B，用B左乘（3）中得到的新向量則重新得到（2）的向量，從而得到原始明文。矩陣A是解密的關(guān)鍵，稱為密鑰。3. 實(shí)例分析3.1 在模26運(yùn)算下，求矩陣的逆矩陣。解：由于det(A)=3，且(mod 26)，故3.2 采用Hill2加密法（即明文2個(gè)2個(gè)分組），加密矩陣取，對 “THE OLYMP

46、IC GAMES” 進(jìn)行加密，再將結(jié)果解密。解：（1）為了保留英文單詞之間的空格，可以增加逗號“，”，句號“.”和空格“ ”當(dāng)做字母，連同26個(gè)英文字母一起與數(shù)字1,2,28,0之間建立對應(yīng)關(guān)系，其中126分別對應(yīng)AZ，“，”對應(yīng)27，“.”對應(yīng)28，“ ”對應(yīng)0。（2）將“THE OLYMPIC GAMES”兩兩分組，最后不足兩個(gè)字母時(shí)用“空格”補(bǔ)足，并用數(shù)字代替，得（3）取矩陣，輸入A=2 7;4 5;x=20 8; 5 0;15 12;25 13; 16 9; 3 0; 7 1;13 5; 19 0'y=mod(A*x,29)結(jié)果為：y = 9 10 27 25 8 6 21 3

47、 9 4 20 4 20 22 12 4 19 18這個(gè)換成字母就是“IDJT,DYTHVFLUDCSIR”,這就是密文。注意，明文中的兩個(gè)“E”，在密文中分別變成了“J”和“S”；而密文中的三個(gè)“D”對應(yīng)明文中的“H”，“L”和“A”，它們并不一一對應(yīng)。（4）如果知道密鑰A，解密是一件很容易的事情。A是模29下的可逆矩陣，其逆矩陣，用B左乘（3）中的向量重新得到20 5 15 25 16 3 7 13 198 0 12 13 9 0 1 5 0換成字母就是“THE OLYMPIC GAMES”，即原始明文。以下是本題加密和解密的一個(gè)完整的MATLAB程序，程序名為Hill2.mclear;c

48、lc;close all;A=2 7;4 5;x='THE OLYMPIC GAMES' % encryption x2=double(x);x3=(x2-64).*(x2>64)+(x2-17).*(x2=44)+(x2-18).*(x2=46)+(x2-32).*(x2=32);if rem(length(x2),2) =1 x3=x3,0;endx4=reshape(x3,2,length(x3)/2);x5=mod(A*x4,29);x6=x5(:)'x7=(x6+64).*(x6<27 & x6)+(x6+17).*(x6=27)+(x6+

49、18).*(x6=46)+(x6+32).*(x6=0);y=char(x7) % decryption a=mod(det(A),29);for k=1:28 if rem(k*a,29)=1 b=k; break; endend B=mod(b*A(2,2),-A(1,2);-A(2,1),A(1,1),29);y2=double(y);y3=(y2-64).*(y2>64)+(y2-17).*(y2=44)+(y2-18).*(y2=46)+(y2-32).*(y2=32);y4=reshape(y3,2,length(y3)/2);y5=mod(B*y4,29);y6=y5(:)

50、'if y6(end)=0 y6(end)=;endy7=(y6+64).*(y6<27 & y6>0)+(y6+17).*(y6=27)+(y6+18).*(y6=28)+(y6+32).*(y6=0);z=char(y7)4. 實(shí)驗(yàn)內(nèi)容4.1 根據(jù)移位加密法編寫一段加密和解密的程序，并將所編程序應(yīng)用于明文 “HELLOWORLD” ，密鑰取3.4.2 甲方收到與之有秘密通信往來的另一方的一個(gè)密文信息，密文內(nèi)容為：JCFOAMOOLABUOHFAXOIKVLYWQSTPKPYABHROARXAMO按照甲方與乙方的約定，他們之間的密文通信采用Hill2密碼，密鑰為矩

51、陣， 26個(gè)字母的表值如正文所示，問這段密文的原始明文是什么？實(shí)驗(yàn)3.2 隨機(jī)模擬 1. 實(shí)驗(yàn)?zāi)康牧私怆S機(jī)模擬的基本概念和原理，通過實(shí)驗(yàn)內(nèi)容，掌握利用計(jì)算機(jī)進(jìn)行隨機(jī)模擬的方法。2. 背景介紹隨機(jī)模擬是一種隨機(jī)試驗(yàn)的方法，也稱蒙特卡洛(Monte Carlo)方法。這種方法源于美國第二次世界大戰(zhàn)期間研制原子彈的“曼哈頓計(jì)劃”，該計(jì)劃的主持人之一，馮·諾依曼用馳名世界的賭城-摩納哥的蒙特卡洛來命名這種方法，使它蒙上了一層神秘的色彩。設(shè)計(jì)一個(gè)隨機(jī)試驗(yàn)，只要使一個(gè)事件的概率與某個(gè)未知數(shù)有關(guān)，然后通過重復(fù)試驗(yàn)，以頻率近似表示概率，即可求出該未知數(shù)的近似解?，F(xiàn)在，隨著計(jì)算機(jī)的發(fā)展，已按照上述思路

52、建立起一類新的方法-隨機(jī)模擬方法。計(jì)算機(jī)產(chǎn)生的隨機(jī)數(shù)都是按照某種確定的算法產(chǎn)生的，它遵循一定的規(guī)律，一旦初值確定，所有隨機(jī)數(shù)也就隨之確定，這顯然不滿足真正隨機(jī)數(shù)的要求，因此我們稱這種隨機(jī)數(shù)為“偽隨機(jī)數(shù)”。在仿真中，通常，它們能當(dāng)做真正的隨機(jī)數(shù)那樣使用。符合各種統(tǒng)計(jì)規(guī)律的隨機(jī)數(shù)可以通過MATLAB的隨機(jī)數(shù)生成函數(shù)很方便地調(diào)用。3. 實(shí)例分析3.1 的模擬計(jì)算建立一個(gè)概率模型，它和有關(guān)，然后設(shè)計(jì)適當(dāng)?shù)碾S機(jī)試驗(yàn)，并通過這個(gè)試驗(yàn)的結(jié)果來確定。通過這個(gè)實(shí)驗(yàn)理解模擬計(jì)算的基本思想。大家知道圓周率的值本身沒有解析解。現(xiàn)在用隨機(jī)模擬的方法設(shè)計(jì)一種求的近似值的方法，并計(jì)算它的近似值。在平面上畫一組距離為a的平行線，把長度為l(l<a)的針投到平面上，則針與平行線相交的概率為。這是著名的蒲豐投針問題，如圖