數 學 實 驗 報 告（試 行 稿） 數學試驗作為一門旨在加深學生更深入了解體會所學數學原理內涵并強化應用所學知識能力的新興課程，尚無固定模式與基本內容相對統一的教材。而基于各不同專業不同的實際需要并為適應不斷發展的科學技術，不應也不可能要求此類課出現類似有幾百年歷史而形成的相對穩定的基礎數學教學模式與相對統一的教材。而相應配套的試驗報告更應具有相當的靈活性，包括格式與基本要求，此試行稿僅提供一種模式參考，希望在使用過程中不斷改進，去舊納新。我們熱切期望教師與同學們提出自己的想法與意見，使報告無論從形式、內容與基本要求諸方面不斷完善，而這一過程是無止境的。 在報告中依據需要附加了一些Matlab的使用方法及與試驗相關的部份程序，僅供參考。 學號： 班號： 姓名： 成績： 實驗一 Matlab的使用目的：初步掌握Matlab的常用功能，包括： （1）算術運算及常用函數的數值運算。 （2）矩陣的運算，應特別強調求逆、特征值、特征向量等的運算指令與使用方法。 （3）掌握2維及3維畫圖方法。 （4）初步掌握編寫幾個語句的小程序，并學會編輯、修改、存盤、打開及調用自己編寫的程序的方法。 （5）了解符號數學工具箱。注：教師可依據實際需要，強調或增加實驗內容，重點介紹某些工具箱的使用。內容： 1上機操作第一章書中所給實例，學習掌握各種數據輸入方法。包括Linspace (a,b,n)。這一函數生成從a到b共n個數值的等差數組，公差不再給出。 2完成下列題目，對每一題應附有必須的程序語句（1）稱為陣左除矩陣，其計算結果與大致相同，但二者又有差異。，稱右除，結果基本相同于。用左除與求逆的方法求解方程組，其中 ； （a）方法： 程序語句： 計算結果： （b）方法： 程序語句： 計算結果： （2）設有分塊陣A如下：其中E、R、O、S分別為單位陣、隨機陣、零陣與對角陣，通過數值計算驗證要求：使矩陣生成函數，及分別生成上述4個矩陣，寫出生成語句。 E= ；R= O= ；S=寫出A的生成語句 A=計算，寫出驗證時你所用的指令： 附：diag( )生成一個對角陣，一般調用格式為diag(v, k) 是一個維行向量，是一個整數，時可缺省. 用下列說明其使用方法. 如，diag生成，divg生成. （3）通過Matlab的幫助系統掌握函數，的用法，并用此諸命令解決以下問題。 （a）已知一個多項式零點為求此按降冪排列的多項式要求：生成多項式的語句為： 計算之值的指令與結果為： （4）矩陣的特征方程是關于的多項式，其根為的特征值，現已知多項式，若能確定一個矩陣，使的特征多項式恰為，則稱為的友元矩陣，從而求的根的問題轉化為求的特征值問題。Matlab完成求友元陣的函數為，這里是一個按的降冪排列得到的系數行向量。 利用求友元陣的特征值的方法求 的根，并與求根函數比較所得結果。要求：求的指令與結果為： 求的特征值的語句結果分別為： 的結果為： 結論： （5）作函數與的圖形要求：在一幅圖上用plot畫出這4條曲線，你的指令為： 用subplot畫4個函數圖形的指令為： （6）fplot是專用于繪制一元函數曲線的命令，其自變量取值點步長是通過內部自適應算法而產生的，所以對曲線伏起劇烈的函數，用fplot命令比用一般等距取點的plot命令給出的曲線更光滑準確。fplot的具體使用格式為：要求：創建如下函數： function ； 在上，令誤差， 繪圖 （7）用作圖法求的根的近似值要求：繪圖語句為： 對軸軸加標志并加網格線的程序語句為： （8）作曲面的三維圖形要求：用你所知的方法繪圖 （a） （b）注：在（5）、（6）、（7）、（8）各作圖題目中，所畫圖形若有條件請打印出來，若不具備條件，請用手工繪制略圖，附在你的試驗報告上，以后的繪圖題目亦有類似要求。 （9）建立M函數作以下計算： （a）自然數的階乘；（b）中取的組合。要求： （a）的程序文件為： （b）的程序文件為： （10）input是編程中，人機對話的主要函數之一，用help查詢其用法，編制滿足如下要求的程序。 用input語句輸入一個函數。 求其一階導數與二階導數。 將函數，及其一階與二階導數畫在同一幅圖中。要求：編寫滿足要求的程序： （11）繪制極坐標系下曲線 要求：繪圖語句： 討論參數的影響 （12）空間兩曲面交線 將與的曲面圖形和交線圖形分別畫在兩張圖上（用函數subplot）。其圖附于報告上。 （13）編寫任意函數展開為各階Taylor多項式的程序，并將各階展開畫在同一幅圖中。要求：你的程序 學號： 班號： 姓名： 成績： 實驗二 平面線性映射的迭代目的：明確在線性映射中，特征值與特征向量所起的作用，通過具體實例，感受線性變換A的特征值中模最大者在計算機動化設計中所扮演的角色。內容： （1）已知概率轉移矩陣 計算的特征值，與及屬于它們的特征向量，寫出求它們的Matlab指令。 指令： 結果： （2）已知今天晴、陰、雨的概率為（P11，P21，P31），則第天晴陰雨的概率計算公式為， 其中利用（1）的結果，寫出利用的特征值與特征向量的利于計算的計算公式，你的公式為： 時，你的計算結果為： （3）與的特征向量作比較，你發現了什么？請敘述你的想法。 （4）設有矩陣及4個二維向量如下 ，求的特征值與屬于他們的特征向量的指令為： （5）計算 考察三個點列的軌跡，在Matla上畫出它們的圖形，將圖形附在下面，注意在圖中加標記，最好將三個圖畫在同一個圖中以作比較，結合書中問題5談談你的感想。 （6）對迭代 迭代映射為從某個方便的開始，畫出由生成的點列。下面用表示該點列。 （7） 計算與，迭代足夠多次，觀察點列的圖形。將PLANAR1改寫為Matlab程序，并用該程序完成本練習。PLANAR1.m （8）由試驗結果看出，多次迭代后映射的作用近似于下式 說明的含意 （9）運行本節后所附Planer2的Matlab程序，觀察這17條軌線，結合練習4，體會模最大的特征值在變換中所起的作用。 （10）對矩陣 重復（9）的練習，注意，與的生成圖形作比較，找出它的按絕對值最大的特征值與屬于這一特征值的方向。(1, 1)(1, 1)(1, 1)(1, 1) （11）試編寫一個程序，不妨稱之為動畫演示程序，將所給圖形沿方向拉長，沿方向壓縮，使其變成一個斜置棱形。你的程序為：運行你的程序得到的圖形是： 12對單位圓面，用滿秩矩陣作坐標變換，對得到的數據畫圖，觀察圖形與你所選矩陣的關系。 你所選坐標變換矩陣 C = 你的畫圖程序為： 變換后的圖形為： C的特征值與特征向量分別為 ； ； 依據圖形的變化情形，你的結論是： 按你的意愿將單位圓在某些方向拉伸或壓縮，依據事先的設定，你如何設計你的坐標變換矩陣。 你提的壓縮，拉伸方案是 你設定的矩陣C應如何得到 上機驗證之。 注意：fill作為繪圖函數，繪制并充填二維圖形，請借助help了解它的詳細使用情況。運行下列程序語句 x=0 1 1 0 0 ; y=0 0 1 1 0; % 正方形頂點 fill (x, y, y) % 繪圖并以黃色填充 axis ( 1 2 1 2 ) % 改變作標軸 說明：fill畫出連接點（0，0）（1，0）（1，1）（0，1）（0，0）的封閉折線圖形。 學號： 班號： 姓名： 成績： 實驗3 常微分方程數值解目的： 對實際問題建立其所滿足要求的微分方程并解之，是培養學生應用數學知識能力的重要環節。因一般非線性微分方程，大都不存在解析解，因此采用數值算法，分析所得數據并結合對其相關圖形的研究就成為處理具體課題的重要科研手段。 通過該實驗，應使學生體會在何種情形下可對非線性方程近似轉化為線性方程以便求得解析解，而在何種情形下，這一近似因產生大的誤差被放棄而采用求數值解的方法。 通過試驗，要求同學們熟練掌握Matlab中求解微分方程數值解的函數ode23與ode45的用法。 試驗布置的題目一般為簡化的實際問題，一方面旨在培養學生的建模能力，另一方面亦要求其有一定的編程，求解方程及對相關數據包括圖形的分析能力，這正是按排此一實驗的初忠。內容： （1）通過課堂學習，明了ode23與ode45兩函數其程序內核的數學方法分別近似為二階與四階龍格庫塔法。需要強調的是盡管ode45的誤差精度高于ode23，但在一些情形下，ode23處理實際問題的效果要優于ode45。 （2）對方程（單擺問題） 令，方程化為與其等價的一階方程組 求解程序為（取）（弧度） 1. function xdot=danbai (, ) 2. g=9.8; 3. l=25； 4. xdot(1)=x(2); 5. xdot(2)=g/l * sin (x (1) ; 6. xdot=xdot(1); xdot(2); 7. t, x=ode23(danbai, 0, 10, 0.1745, 0針對上敘程序，回答如下問題： a) t的維數由什么確定？運行程序試一試。 b) 寫出的形式，說明矩陣兩列的含意。 c) 去掉語句6中分號，觀察出現的情形，并解釋說明之。 d) t的值在第幾條語句中賦值？ 仔細體會程序的運作原理。 （3）對改進的歐拉公式 a) 編寫求微分方程數值解的程序 你的程序為： b) 用上述程序與ode23解下列微分方程 （1） 或 （2） ， （4）一只小船渡過寬為的河流，目標是起點正對著的另一岸上點，已知河水流速與船在靜水中的速度之比為 a) 建立小船航線的方程，求其解析解。 b) 設，=，用數值解法求渡河所需時間，任意時刻小船的位置及航行曲線，作圖并與解析解比較。 （5）研究種群競爭模型，當甲、乙兩個種群各獨自生存時，數量演變服從下面規律。其中，分別為時刻甲、乙兩個種群的數量，為其固有增長率，為它們的最大容量，而當這兩個種群在同一環境中生存時，由于乙消耗有限資源對甲的增長產生影響，將甲的方程修改為 （1）這里的含義是，對于供養甲的資源而言，單位數量乙（相對）的消耗率為單位數量甲（相對）消耗率的倍，類似，甲的存在也影響了乙的增長，乙的方程應改為 （2）給定種群初始值為及參數后方程（1）與（2）確定了兩種群的變化規律，因其解析解不存在，試用數值解法研究以下問題： a) 設，計算，畫出它們的圖形及相圖，說明時間充分大以后的變化趨勢。 b) 改變，但不變（或保持）并分析所得結果，若，再分析結果，由此你得到什么結論，請用各參數生態學上的含意作出 c) 試驗當，時會有什么結果；當，時又會出現什么結果，能解釋這些結果嗎？ （6）前言中，曾就導彈B追蹤飛機A的軌跡問題展開討論，那里，為簡化問題，假設B與A的距離不變，其中，為的飛行軌跡的參數方程，待求函數。，是的飛行軌跡，是已知的。下面，我們繼續討論，并去掉不變的條件，而令。 （a）證明，B的飛行路線滿足微分方程其中為B的恒定速率。 （b）下面的4個程序動態顯示了B追蹤A并最終追上A的過程。讀懂并運行程序。A的逃逸路線： function S=A(t) If t6 S=8*t; 0; else S=8*(12-t); 1; end說明：A先沿軸正向飛行，后當時，調頭回飛。B的追蹤路線（待求解方程）： function zs, isterminal, direction=B(t, z, flag) global w % speed of the B X=A(t); h=X-z; nh=norm(h); % nh= if nargin3 | isempty (flag) zs=(w/nh)*h; else switch (flag) case events zs=nh-1e-3; isterminal=1; direction=0; otherwise error (unknow flag:, flag); end end說明：由(a)，待解方程中存在；為在適當時間（變的很小時）中止運算，求解方程的積分器ode23及ode45借助函數odeset( )設定events選項為on時，當isterminal=1，且direction=0時，將終止迭代運算。詳見help odeset. 求解方程及圖示追蹤過程的主程序如下： % 以 main1-m存盤 global w yo=60; 70; % initial starting point of B w=10; % speed of the B options=odeset (RelTol, 1e-5, Events, on); t, y=ode23 (B, 0, 20, yo, options); j= ; for h=1: lengthe (t) w=A(t(h); j=j; w;end xmin=min(min(y(:, 1), min(j(:, 1); xmax=max(max(y(:, 1), max(j(:, 1); ymin=min(min(y(:, 2), min(j(:, 2); ymax=max(max(y(:, 2), max(j(:, 2); clf; hold on ; axis=(xmin, xmax, ymin, ymax); for h=1: length(t)-1 plot (y(h, 1), y (h+1, 1), y(h, 2), y(h+1, 2), -, color, red, EraseMode, none); plot (j(h, 1), j(h+1, 1), j(h, 2), j(h+1, 2), :, color, black, EraseMode, none); drawnow; pause (0.3); end p=max(size(y); cross_1(y(p, 1), y(p, 2);hold off說明：cross_1是在B追上A的那點處作一個標記。下面是這一函數的代碼。 function cross_1 (cx, cy) v=2; kx=cx, cx, cx, cx-v, cx+v; ky=cy, cy+2.5*v, cy+1.5*v, cy+1.5*v, cy+1.5*v; plot (kx, ky) plot (cx, cy, 0) （c）如果建立微分方程很困難，則可用仿真的方法在計算機上一步步模擬B追蹤A的實際過程。請給出算法，并編寫程序動態顯示B追蹤A的過程。 你的理論依據是： 你的程序。 學號： 班號： 姓名： 成績： 實驗四 Matlab程序編寫目的： 程序的編寫調試，是建模的重要組成部分。為進一步提高同學們的編程能力，本節試驗給出了一系列具有一定難度的實用程序，要求同學們讀懂并上機運行之。 關于程序的功能說明，除必要的文字簡述，大部份附加于程序的幫助系統內。而更深入了解程序的運作機制，則需同學們進一步閱讀相關書籍。 程序的選擇力求其實用性，并重視動畫的演示效果，希望這些程序的給出，對激發同學們編寫高質量應用程序的熱情起到拋磚引玉的效果。 一、options=odeset( )的設定形式 以求解高階微分方程為例，重點說明結構變量options中事件變量events的使用。 （1）求解高階方程時，還須將其等價的轉化為一階方程組，已知待求解方程為 與其等價的一價微分方程組為： （2）求解范德堡（Vandor Pol）方程 與其等價的一階微分方程組為：在前項給出的結果下，編寫待求解程序為：%以vdpol.m存盤 if navginodeset % 顯示odeset的設定內容Abstol:(絕對誤差),缺省值le-6Reltol:(相對誤差),缺省值le-3Events: function在本問題中，結構options除Events的選項設定為函數Vdpolevents的地址(指針)外,其佘全部用缺省值, 在調用ode45時，積分器指向vdpol, 由上一步值，遞歸的計算下一步值，同時調用events所指向的函數vdpolevents,并將新計算的函數值傳給vdpolevents, 之后計算value是否為0, 若不為0, 返回ode45繼續計算, 否則, 檢查is terminal的值,若此值為0, 則ode45繼續計算, 并將使Value為0的值返給te與ye, 但若isterminal值為1, 則終止ode45的調用,并將使Value=0的t與y值返給te與ye, 對direction的賦值為0, 1, -1, 值為0時，忽略Value是大于0還是小于0趨于0，賦值為1時，則當Value大于0趨于0時，終止運算，。在實驗三的問題（6）中，給出了導彈追綜并追上擊中飛機的追綜路線求解程序，在那里終止ode45運算的事件程序混編在待求解方程中。重新改寫程序，將終止運算程序與待求解方程分離。你的程序為： （a）待求解的B的追綜路線程序： （b）終止ode45調用的事件函數： （c）求解方程的主程序為：二、應用程序1：ex_draw(arg)這一程序在當前坐標軸下用鼠標拖動畫線，在畫線起點，按鼠標左鍵，出現一個“+”字，之后按住左鍵（不放鍵），拖動鼠標，則可畫出一條紅線，指針到哪兒，紅線跟到哪兒，放開鼠標，一條線畫線結束，之后可重復畫第二，第三條線，。按任意鍵結束畫線。函數以es_draw.m存盤。以下為程序代碼。%遞歸調用的畫線程序function ex_draw(arg)global DRAW_HL % 全局變量一般用大寫字母。if nargin=0arg=ex_line;endif isstr(arg) % initial call,clear all % set things up.Hf=gcf;%若當前有圖形窗口，將窗口句柄賦于 Hf，若無圖形窗口，則建立%新窗口對象，并將其句柄值存于Hf中。Set (Hf , Pointer,crossh,;%set up callback for line startBackingStore,off WindowButtonDownFcn,ex_draw(1)figure(Hf)elseif arg=1 % callback is line start pointfp=get(gca,CurrentPoint); % start of line pointset(gca,Userdata,fp(1,1:2) % store in axes userdataser(gcf,WindowButtonMotionFcn,ex_draw(2),WindowButtonUpFcn,ex_draw(3)elseif arg=2 % callback is mouse motioncp=get(gca,CurrentPoint);cp=cp(1,1:2);fp=get(gca,Userdata);Hl=line(Xdata,fp(1);cp(1),Ydata,fp(2);cp(2),EraseMode,xor,Color,r,LineStyle,-,Clipping,off);if isempty(DRAW_HL) % delete prior line if it existsdelete(DRAW_HL)endDRAW_HL=H1; % store current line handleelseif arg=3 % callback is line end point, finish upset(gcf,Pointer,arrow,BackingStore,onWindowButtonDownFcn, ,WindowButtonMotionFcn, ,WindowButtonUpFcn, )set(gca,Userdata,);set(DRAW_HL,EraseMode,normal) % render line betterkey=waitforbuttonpress;if keygtext(finish);returnelseex_draw;endend說明：nargin 與nargout 為Matlab提供的檢測函數輸入輸出參數個數的函數。arg則表明函數有輸入參數，但參數形式不定。亦可無參數調用。注意程序的遞歸調用。說明：arg=1時的調用情形：arg=2 ：arg=3 ：說明Key=waitforbuttonpress 的功能：三、應用程序2：eigshow(arg)eigshow( ) 動態顯示線性變換中特征值與特征向量的作用。注意各子程序的功能。讀懂并上機運行。簡述各程序模塊的功能。initialize:initv:action:movev:setmode:function eigshow(arg) %顯示二維矩陣的特征值的本質.if nargin= =0; initializeelseif arg= =0 actionelseif arg0 setmode(arg)else initialize(arg);endfunction initialize(arg)if nargin= =0 arg=6;endif isequal(get(gcf,tag),eigshow); h=get(gcf,userdata); mats=h.mats;else set(gcf,numbertitle,off,menubar,none) h.svd=0; mats= %二階方陣數組.5/4 0; 0 3/45/4 0; 0 3/41 0; 0 10 1; 1 00 1; -1 01 3; 4 2/41 3; 2 4/43 1; 4 2/43 1; -2 4/42 4; 2 4/42 4; -1 -2/46 4; -1 2/4randn(2,2);end%根據參數的不同確定變換矩陣.if all(size(arg)= =1) if (arglength(mats)mindex=arg;A=eval(matsmindex); %變換矩陣由二階方陣數組中的某一個. else A=randn(2,2); S=sprintf(%4.2f %4.2f; %4.2f %4.2f,A.); %一個二維隨機陣. mindex=length(mats); mats=mats(1:mindex-1); S; mats(mindex); endelse A=arg; %變換矩陣直接在參數中指定. if isstr(A);S=A;A=eval(A); else S= sprintf( %4.2f %4.2f; %4.2f %4.2f,A.); end if any(size(A) =2)error(Matrix must be 2-by-2) end mats=S; mats; mindex=1;endclfif h.svd, t=svd/ (eig); else, t=eig/ (svd);enduicontrol( %界面控制按鈕style,pushbutton,units,normalized,position,.86 .60 .12 .06,string,t,value,h.svd,callback,eigshow(-1);uicontrol( style,pushbutton,units,normalized,position, .86 .50 .12 .06,string,help,callback,helpwin eigshow)uicontrol( style,pushbutton,units,normalized,position,.86 .40 .12 .06,string,close,callback,close(gcf)uicontrol( style,popup, %變換矩陣選擇下拉菜單.units,normalized,position,.28 .92 .48 .08,string,mats,tag,mats, fontname,courier fontweight,bold,fontsize,14, value,mindex, callback,eigshow(get(gco,value);s=1.1*max(1,norm(A); %圖形繪制.axis(-s s s s)axis squarexcolor=0 .6 0Axcolor=0 0 .8h.A=A;h.mats=mats;h.x=initv(1 0,x,xcolor);h.Ax=initv(A(:, 1),Ax,Axcolor);if h.svd h.y=initv(0 1,y,xcolor); h.Ay=initv(A(:, 2),Ay,Axcolor); xlabel(Make A*x pependicular to A*y,fontweight,bold) set(gcf,name,svdshow)else xlabel(Make A*x parallel to x,fontweight,bold) set(gcf,name,eigshow)endset(gcf,tag,eigshow, userdata,h, windowbuttondownfcn, eigshow(0);set(gcf, windowbuttonmotionfcn , eigshow(0) ), windowbuttonupfcn, set(gcf, windowbuttonmotionfcn , )%-function h=initv(v,t,color) %動畫演示程序體.h.mark=line(v(1),v(2),marker,.,erase,none,color,color);h.line=line(0 v(1),0 v(2),erase,xor,color,color);h.text=text(v(1)/2,v(2)/2,t,fontsize,12,erase,xor,color,color);%-function actionh=get(gcf,userdata); %與用戶的交互作用程序.pt=get(gca,currentpoint);x=pt(1,1:2);x=x/norm(x);movev(h.x,x);A=h.A;movev(h.Ax,A*x);if h.svd y=-x(2);x(1); movev(h.y,y); movev(h.Ay,A*y);end%-function movev(h,v)set(h.mark,xdata,v(1),ydata,v(2);set(h.line, xdata,0 v(1),ydata,0 v(2);set(h.text,pos,v/2);%-function setmode(arg)h=get(gcf,userdata);h.svd=h.svd;set(gcf,userdata,h)initialize(get(findobj(gcf,tag,mats),value)四、閱讀與理解1：給出三段程序。mmgetxy( )調用mmcxy( ),mmcxy( )調用mmgcf，程序代碼如下：-function Hf=mmgcfHf=get(0,Children);if isempty(Hf)returnelse Hf=get(0,CurrentFigure);end-function out=mmcxy(arg)global MMCXY_OUTifnarginHf=mmgcf;if isempty(Hf);error(No Figure Available.),endHa=findobj(Hf,Type,axes);if isempty(Ha),error(No Axes in Current Figure.),endHu=uicontrol(Hf,Style,text, units,pixels, Position,1 1 140 15, HorizontalAlignment,left);set(Hf,Pointer,crossh, WindowButtonMotionFcn,mmcxy( move ), WindowButtonDownFcn,mmcxy( end ), Userdata,Hu);figure(Hf) %bring figure forwardif nargout %must return x-y data key=waitforbuttonpress; %pause until mouse is pressed if dey, out= ; %return empty if aborted mmcxy(end); %clean things up else out=MMCXY_OUT; %now that move is complete return point end returnend elseif strcmp(arg,move) %mouse is moving in figure windowcp=get(gca,CurrentPoint); %get current mouse posstionMMCXY_OUT=cp(1,1:2);xystr=sprintf(%.3g,%.3g,MMCXY_OUT);Hu=get(gcf,Userdata);set(Hu,String,xystr); %put x-ycoordinates in text boxelseif strcmp(arg,end) %mouse click occurred,clean things upHu=get(gcf,Userdata);delete(Hu)set(gcf,Pointer,arrow, WindowButtonMotionFcn, , WindowButtonDownFcn, , Userdata, ); endfunction xy=mmgctxy(n)if nargin= =0,n=5;end % nobody wants more points!xy= ;s=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,;axis(axis);hold onfor i=1:ntmp=mmcxy;if isempty(tmp) returnelse xy=xy;tmp; plot(xy(i,1),xy(i,2),r*) text(xy(i,1),xy(i,2)+0.05,si)endendfor j=1:2for i=1:2waitforbuttonpress;dot=get(gca,CurrentPoint);xdata(i)=dot(1,1);ydata(i)=dot(1,2);end line(xdata,ydata);endhold off