matlab數值分析實驗

數值分析實驗董海云數理學院數學實驗教學中心目錄0Matlab介紹入門知識 31緒論 171.1例題解答 171.2Matlab中數值計算精度 202線性方程組的直接解法 222.1例題解答 222.2Matlab解線性方程組常用命令介紹 363線性方程組的迭代解法 383.1例題解答 383.2Matlab迭代解法用到的函數介紹 534方陣特征值和特征向量的計算 554.1例題解答 554.2Matlab關于方陣特征值為特征向量函數介紹 625非線性方程求根 645.1例題解答 645.2Matlab非線性方程求根的命令 856插值法 866.1例題解答 866.2Matlab插值函數介紹 1017數據擬合和最佳平方逼近 1037.1例題解答 1037.2Matlab數據擬合命令介紹 1138數值積分與數值微分 1148.1例題解答 1149常微分方程數值解法 1389.1例題解答 1389.2Matlab常微分方程數值解常用命令介紹 1540Matlab介紹入門知識1.Matlab簡介MATLAB的含義是矩陣實驗室（MATRIXLABORATORY），主要用于方便矩陣的存取，其基本元素是無須定義維數的矩陣.MATLAB自問世以來,就是以數值計算稱.MATLAB進行數值計算的基本單位是復數數組（或稱陣列），這使得MATLAB高度“向量化”.經過十幾年的完善和擴充，現已發(fā)展成為線性代數課程的標準工具.由于它不需定義數組的維數，并給出矩陣函數、特殊矩陣專門的庫函數，使之在求解諸如信號處理、建模、系統(tǒng)識別、控制、優(yōu)化等領域的問題時，顯得大為簡捷、高效、方便，這是其它高級語言所不能比擬的.MATLAB中包括了被稱作工具箱（TOOLBOX）的各類應用問題的求解工具.工具箱實際上是對MATLAB進行擴展應用的一系列MATLAB函數（稱為M文件），它可用來求解各類學科的問題，包括信號處理、圖象處理、控制系統(tǒng)辨識、神經網絡等.隨著MATLAB版本的不斷升級，其所含的工具箱的功能也越來越豐富，因此，應用范圍也越來越廣泛，MATLAB提供的工具箱已覆蓋信號處理、系統(tǒng)控制、統(tǒng)計計算、優(yōu)化計算、神經網絡、小波分析、偏微分方程、模糊邏輯、動態(tài)系統(tǒng)模擬、系統(tǒng)辨識和符號運算等領域.當前它的使用范圍涵蓋了工業(yè)、電子、醫(yī)療、建筑等各行各業(yè).MATLAB中包括了圖形界面編輯GUI，讓使用者也可以象VB、VC、VJ、DELPHI等那樣進行一般的可視化的程序編輯.在命令窗口（matlabcommandwindow）鍵入simulink，就出現(SIMULINK)窗口.以往十分困難的系統(tǒng)仿真問題，用SIMULINK只需拖動鼠標即可輕而易舉地解決問題，這也是近來受到重視的原因所在.MATLAB語言由美國TheMathWorks開發(fā),最早是由C.Moler用Fortran語言編寫的，用來方便地調用LINPACK和EISPACK矩陣代數軟件包的程序.后來他創(chuàng)立了MATHHWORKS公司，對MATLAB作了大量的、堅持不懈的改進.CleveB.Moler是TheMathWork公司的主席和首席科學家.曾任密歇系教授.他在兩個計算機硬件制造商Intel公司的Hypercube組織和ArdenComputers公司工作了五年.他的主要專業(yè)興趣在于數值分析和科學計算.他是MATLAB軟件的創(chuàng)始者，也是著名的矩陣計算軟件包LINPACK和EISPACK的著作這一，已撰寫了三本有相關數值方法的教材.同時，他在SIAM（美國工業(yè)與應用數學學會）歷任期刊編輯、委員會成員和副總裁，并從1996年開始擔任理事會成員.2.Matlab入門知識Matlab變量名是以字母開頭，后接字母、數字或下劃線的字符序列，最多63個字符.在MATLAB中，變量名區(qū)分字母的大小寫.賦值語句:變量=表達式或表達式其中表達式是用運算符將有關運算量連接起來的式子，其結果是一個矩陣.clear命令用于刪除MATLAB工作空間中的變量.who和whos這兩個命令用于顯示在MATLAB工作空間中已經駐留的變量名清單.who命令只顯示出駐留變量的名稱，whos在給出變量名的同時，還給出它們的大小、所占字節(jié)數及數據類型等信息.利用MAT文件可以把當前MATLAB工作空間中的一些有用變量長久地保留下來，擴展名是.mat.MAT文件的生成和裝入由save和load命令來完成.常用格式為：save文件名[變量名表][-append][-ascii]load文件名[變量名表][-ascii]其中，文件名可以帶路徑，但不需帶擴展名.mat，命令隱含一定對.mat文件進行操作.變量名表中的變量個數不限，只要內存或文件中存在即可，變量名之間以空格分隔.當變量名表省略時，保存或裝入全部變量.-ascii選項使文件以ASCII格式處理，省略該選項時文件將以二進制格式處理.save命令中的-append選項控制將變量追加到MAT文件中.(1)向量的創(chuàng)建用步長生成法：數組=初值:步長(增量):終值>>a=1:0.5:3a=1.00001.50002.00002.50003.0000用linspace生成：數組=linspace(初值,終值,等分點數目)>>b=linspace(1,3,5)b=1.00001.50002.00002.50003.0000列向量用分號（;）作為分行標記：>>c=[1;2;3;4;]c=1234若不想輸出結果，在每一條語句后用分號作為結束符，若留空或用逗號結束，則在執(zhí)行該語句后會把結果輸出來.>>a+b;>>a+bans=23456(2)矩陣的創(chuàng)建直接輸入:最簡單的建立矩陣的方法是從鍵盤直接輸入矩陣的元素.具體方法如下：將矩陣的元素用方括號括起來，按矩陣行的順序輸入各元素，同一行的各元素之間用空格或逗號分隔，不同行的元素之間用分號分隔.>>A=[123;456;235]A=123456235利用矩陣函數創(chuàng)建:>>B=magic(3)%魔方陣B=816357492>>C=hilb(3)%3階Hilbert矩陣C=1.00000.50000.33330.50000.33330.25000.33330.25000.2000Matlab中用%引導注釋其它創(chuàng)建矩陣函數還有:eye(m,n):生成m行n列單位矩陣.zeros(m,n):生成m行n列全零矩陣.ones(m,n):生成全1矩陣.rand(m,n):生成隨機矩陣.rand:生成一個隨機數.diag(A):取A的對角線元素.tril(A):取A的下三角元素.triu(A):取A的上三角元素.hilb(n):生成n維Hilbert矩陣.randn(n)：產生均值為0，方差為1的標準正態(tài)分布隨機矩陣.vander(V):生成以向量V為基礎向量的范得蒙矩陣.invhilb(n):求n階的希爾伯特矩陣的逆矩陣.toeplitz(x,y):生成一個以x為第一列，y為第一行的托普利茲矩陣compan(p):生成伴隨矩陣,p是一個多項式的系數向量，高次冪系數排在前，低次冪排在后.pascal(n):生成一個n階帕斯卡矩陣.compan:生成伴隨矩陣(3)矩陣運算MATLAB的基本算術運算有：＋(加)、－(減)、*(乘)、/(右除)、\(左除)、^(乘方).加法：>>A+Bans=939710136127減法：>>B-Aans=7-13-10126-3乘法：>>A*Bans=263826718371456243除法:>>magic(3)/hilb(3)ans=1.0e+003*0.2160-1.17601.14000.0570-0.40800.4500-0.22801.2240-1.1400在MATLAB中，有一種特殊的運算，因為其運算符是在有關算術運算符前面加點，所以叫點運算.點運算符有.*、./、.\和.^.兩矩陣進行點運算是指它們的對應元素進行相關運算，要求兩矩陣的維參數相同.>>A.*Bans=821812254282710MATLAB提供了6種關系運算符：<(小于)、<=(小于或等于)、>(大于)、>=(大于或等于)、==(等于)、～=(不等于).>>A>Bans=010100001MATLAB提供了3種邏輯運算符：&(與)、|(或)和～(非).在邏輯運算中，確認非零元素為真，用1表示，零元素為假，用0表示.設參與邏輯運算的是兩個標量a和b，那么，a&ba,b全為非零時，運算結果為1，否則為0.a|ba,b中只要有一個非零，運算結果為1.～a當a是零時，運算結果為1；當a非零時，運算結果為0.3.矩陣操作和矩陣函數矩陣通過下標引用矩陣的元素，矩陣元素的序號就是相應元素在內存中的排列順序.在MATLAB中，矩陣元素按列存儲，先第一列，再第二列，依次類推.(1)矩陣拆分利用冒號表達式獲得子矩陣.A(:,j)表示取A矩陣的第j列全部元素；A(i,:)表示A矩陣第i行的全部元素；A(i,j)表示取A矩陣第i行、第j列的元素.A(i:i+m,:)表示取A矩陣第i～i+m行的全部元素；A(:,k:k+m)表示取A矩陣第k～k+m列的全部元素，A(i:i+m,k:k+m)表示取A矩陣第i～i+m行內，并在第k～k+m列中的所有元素.此外，還可利用一般向量和end運算符來表示矩陣下標，從而獲得子矩陣.end表示某一維的末尾元素下標.(2)利用空矩陣刪除矩陣的元素在MATLAB中，定義[]為空矩陣.給變量X賦空矩陣的語句為X=[].(3)矩陣的轉置轉置運算符是單撇號(‘).(4)矩陣的旋轉利用函數rot90(A,k)將矩陣A旋轉90o的k倍，當k為1時可省略.(5)矩陣的左右翻轉對矩陣實施左右翻轉是將原矩陣的第一列和最后一列調換，第二列和倒數第二列調換，…，依次類推.MATLAB對矩陣A實施左右翻轉的函數是fliplr(A).(6)矩陣的上下翻轉MATLAB對矩陣A實施上下翻轉的函數是flipud(A).(7)方陣A的逆矩陣inv(A)>>A=magic(3)A=816357492>>B=inv(A)B=0.1472-0.14440.0639-0.06110.02220.1056-0.01940.1889-0.1028>>A*Bans=1.00000-0.0000-0.00001.000000.000001.0000(8)方陣的行列式>>det(A)ans=-360(9)矩陣的跡>>C=trace(A)C=15(10)一些常用的基本初等三角函數三角函數:sin(x),cos(x),tan(x)反三角函數:asin(x),acos(x),atan(x)指數函數:exp(x)自然對數:log(x)常用對數:log10(x)以2為底的對數:log2(x)開平方:sqrt(x)絕對值:abs(x)計算一般函數值:eval(f)求虛部函數:imag(x)求實部函數:real(x)角相位函數:angle(x)共軛復數函數:conj(x)沿零方向取整:fix(x)舍入取整:round(x)沿負無窮大方向取整:floor(x)沿正無窮大方向取整:ceil(x)求除法的余數:rem符號函數:sign(x)最大公約數:gcd()4.圖形可視化(1)二維繪圖指令plotplot函數的基本調用格式為:plot(x,y,)其中x和y為長度相同的向量，分別用于存儲x坐標和y坐標數據.plot(x)plot函數最簡單的調用格式.當x是實向量時，以該向量元素的下標為橫坐標，元素值為縱坐標畫出一條連續(xù)曲線.實際上是繪制折線圖.plot(x1,y1,x2,y2,…,xn,yn)當輸入參數都為向量時，x1和y1，x2和y2，…，xn和yn分別組成一組向量對，每一組向量對的長度可以不同.每一向量對可以繪制出一條曲線，可以在同一坐標內繪制出多條曲線.plotyy(x1,y1,x2,y2)繪制出具有不同縱坐標標度的兩個圖形.holdon/off保持原有圖形還是刷新原有圖形,不帶參數的hold命令在兩種狀態(tài)之間進行切換.plot(x1,y1,選項1,x2,y2,選項2,…,xn,yn,選項n)設置曲線樣式進行繪圖.選項字段見下表:表Matlab常用線形與顏色標記表符號線型符號線型符號線型顏色含義-實線.實點標記^朝上三角g綠色--虛線o圓圈<朝左三角r紅色：點線X叉字符>朝右三角c青色-.點劃線+加號p五角星m洋紅*星號h六角形y黃色s方塊k黑色d菱形w白色v朝下三角b藍色(2)圖形標注:title('圖形名稱'):圖形標題xlabel('x軸說明')ylabel('y軸說明')text(x,y,'圖形說明')legend('圖例1','圖例2',…)gtext('用鼠標確定位置的字符說明')(3)坐標控制axisaxis([xminxmaxyminymaxzminzmax])axis函數功能豐富，常用的格式還有：axisequal：縱、橫坐標軸采用等長刻度.axissquare：產生正方形坐標系(缺省為矩形).axisauto：使用缺省設置.axisoff：取消坐標軸.axison：顯示坐標軸.gridon/off:網格開/關boxon/off:加/不加邊框線上述命令示例如下:>>x=1:length(peaks);>>plot(x,peaks);>>boxon;>>title('繪制混合圖形');>>xlabel('X軸');>>ylabel('Y軸');繪制圖像為:(4)二維數值函數的專用繪圖函數fplotfplot(functionname,[a,b],tol,選項)其中functionname為函數名，以字符串形式出現，[a,b]為繪圖區(qū)間，tol為相對允許誤差，其系統(tǒng)默認值為2e-3.選項定義與plot函數相同.>>fplot(@(x)[tan(x),sin(x),cos(x)],2*pi*[-11-11]);(5)二維符號函數曲線專用命令ezplotf=f(x)時:ezplot(f)：在默認區(qū)間-2π<x<2π繪制f=f(x)的圖形.ezplot(f,[a,b])：在區(qū)間a<x<b繪制f=f(x)的圖形f=f(x,y)時:ezplot(f)：在默認區(qū)間-2π<x<2π和-2π<y<2π繪制f(x,y)=0的圖形.ezplot(f,[xmin,xmax,ymin,ymax])：在區(qū)間xmin<x<xmax和ymin<y<ymax繪制f(x,y)=0的圖形ezplot(f,[a,b])：在區(qū)間a<x<b和a<y<b繪制f(x,y)=0的圖形若x=x(t),y=y(t):ezplot(x,y)：在默認區(qū)間0<t<2π繪制x=x(t)和y=y(t)的圖形.ezplot(x,y,[tmin,tmax])：在區(qū)間tmin<t<tmax繪制x=x(t)和y=y(t)的圖形>>figure;ezplot('cos(tan(pi*x))',[0,1]);(6)圖形窗口的分割subplotsubplot(m,n,p)該函數將當前圖形窗口分成m×n個繪圖區(qū)，即每行n個，共m行，區(qū)號按行優(yōu)先編號，且選定第p個區(qū)為當前活動區(qū).在每一個繪圖區(qū)允許以不同的坐標系單獨繪制圖形.(7)其他坐標系下的二維數據曲線圖對數坐標圖形:semilogx(x1,y1,選項1,x2,y2,選項2,…)semilogy(x1,y1,選項1,x2,y2,選項2,…)loglog(x1,y1,選項1,x2,y2,選項2,…)極坐標圖polar:polar(theta,r,選項)其中theta為極坐標極角，r為極坐標矢徑，選項的內容與plot函數相似.二維統(tǒng)計分析圖:bar(x,y,選項):條形圖stairs(x,y,選項):階梯圖stem(x,y,選項):桿圖fill(x1,y1,選項1,x2,y2,選項2,…):填充圖(8)三維曲線plot3plot3(x1,y1,z1,選項1,x2,y2,z2,選項2,…,xn,yn,zn,選項n)其中每一組x,y,z組成一組曲線的坐標參數，選項的定義和plot函數相同.當x,y,z是同維向量時，則x,y,z對應元素構成一條三維曲線.當x,y,z是同維矩陣時，則以x,y,z對應列元素繪制三維曲線，曲線條數等于矩陣列數.>>t=0:0.1:8*pi;>>plot3(sin(t),cos(t),t);(9)產生三維數據在MATLAB中，利用meshgrid函數產生平面區(qū)域內的網格坐標矩陣.其格式為：[X,Y]=meshgrid(x,y);語句執(zhí)行后，矩陣X的每一行都是向量x，行數等于向量y的元素的個數，矩陣Y的每一列都是向量y，列數等于向量x的元素的個數.(10)繪制三維曲面的函數surf函數和mesh函數的調用格式為:mesh(x,y,z,c)surf(x,y,z,c)一般情況下，x,y,z是維數相同的矩陣.x,y是網格坐標矩陣，z是網格點上的高度矩陣，c用于指定在不同高度下的顏色范圍.(11)標準三維曲面sphere函數的調用格式為：[x,y,z]=sphere(n)cylinder函數的調用格式為：[x,y,z]=cylinder(R,n)MATLAB還有一個peaks函數，稱為多峰函數，常用于三維曲面的演示.(12)其他三維繪圖指令介紹bar3函數繪制三維條形圖，常用格式為bar3(y)bar3(x,y)stem3函數繪制離散序列數據的三維桿圖，常用格式為：stem3(z)stem3(x,y,z)pie3函數繪制三維餅圖，常用格式為：pie3(x)fill3函數等效于三維函數fill，可在三維空間內繪制出填充過的多邊形，常用格式為：fill3(x,y,z,c)5.程序控制結構(1)數據的輸入:A=input(提示信息，選項)其中提示信息為一個字符串，用于提示用戶輸入什么樣的數據.如果在input函數調用時采用's'選項，則允許用戶輸入一個字符串.(2)數據的輸出:disp(輸出項)(3)程序的暫停:pause(延遲秒數)如果省略延遲時間，直接使用pause，則將暫停程序，直到用戶按任一鍵后程序繼續(xù)執(zhí)行.若要強行中止程序的運行可使用Ctrl+C命令.(4)單分支if語句：if條件語句組end當條件成立時，則執(zhí)行語句組，執(zhí)行完之后繼續(xù)執(zhí)行if語句的后繼語句，若條件不成立，則直接執(zhí)行if語句的后繼語句.(5)雙分支if語句：if條件語句組1else語句組2end當條件成立時，執(zhí)行語句組1，否則執(zhí)行語句組2，語句組1或語句組2執(zhí)行后，再執(zhí)行if語句的后繼語句.(6)多分支if語句：if條件1語句組1elseif條件2語句組2……elseif條件m語句組melse語句組nend語句用于實現多分支選擇結構.(7)switch語句:switch表達式case表達式1語句組1case表達式2語句組2……case表達式m語句組motherwise語句組nend(8)try語句語句格式為：try語句組1catch語句組2endtry語句先試探性執(zhí)行語句組1，如果語句組1在執(zhí)行過程中出現錯誤，則將錯誤信息賦給保留的lasterr變量，并轉去執(zhí)行語句組2.(9)for語句for語句的格式為：for循環(huán)變量=表達式1:表達式2:表達式3循環(huán)體語句end其中表達式1的值為循環(huán)變量的初值，表達式2的值為步長，表達式3的值為循環(huán)變量的終值.步長為1時，表達式2可以省略.for語句更一般的格式為：for循環(huán)變量=矩陣表達式循環(huán)體語句end執(zhí)行過程是依次將矩陣的各列元素賦給循環(huán)變量，然后執(zhí)行循環(huán)體語句，直至各列元素處理完畢.(10)while語句while語句的一般格式為：while(條件)循環(huán)體語句end其執(zhí)行過程為：若條件成立，則執(zhí)行循環(huán)體語句，執(zhí)行后再判斷條件是否成立，如果不成立則跳出循環(huán).(11)break語句和continue語句與循環(huán)結構相關的語句還有break語句和continue語句.它們一般與if語句配合使用.break語句用于終止循環(huán)的執(zhí)行.當在循環(huán)體內執(zhí)行到該語句時，程序將跳出循環(huán)，繼續(xù)執(zhí)行循環(huán)語句的下一語句.continue語句控制跳過循環(huán)體中的某些語句.當在循環(huán)體內執(zhí)行到該語句時，程序將跳過循環(huán)體中所有剩下的語句，繼續(xù)下一次循環(huán).(12)循環(huán)的嵌套如果一個循環(huán)結構的循環(huán)體又包括一個循環(huán)結構，就稱為循環(huán)的嵌套，或稱為多重循環(huán)結構.(13)函數文件的基本結構函數文件由function語句引導，其基本結構為function輸出形參表=函數名(輸入形參表)注釋說明部分函數體語句其中以function開頭的一行為引導行，表示該M文件是一個函數文件.函數名的命名規(guī)則與變量名相同.輸入形參為函數的輸入參數，輸出形參為函數的輸出參數.當輸出形參多于一個時，則應該用方括號括起來.(14)函數調用函數調用的一般格式是：[輸出實參表]=函數名(輸入實參表)注意的是，函數調用時各實參出現的順序、個數，應與函數定義時形參的順序、個數一致，否則會出錯.函數調用時，先將實參傳遞給相應的形參，從而實現參數傳遞，然后再執(zhí)行函數的功能.在MATLAB中，函數可以嵌套調用，即一個函數可以調用別的函數，甚至調用它自身.一個函數調用它自身稱為函數的遞歸調用.(15)函數參數的可調性在調用函數時，MATLAB用兩個永久變量nargin和nargout分別記錄調用該函數時的輸入實參和輸出實參的個數.只要在函數文件中包含這兩個變量，就可以準確地知道該函數文件被調用時的輸入輸出參數個數，從而決定函數如何進行處理.(16)全局變量與局部變量全局變量用global命令定義，格式為：global變量名(17)程序調試Debug菜單項:該菜單項用于程序調試，需要與Breakpoints菜單項配合使用.Breakpoints菜單項:該菜單項共有6個菜單命令，前兩個是用于在程序中設置和清除斷點的，后4個是設置停止條件的，用于臨時停止M文件的執(zhí)行，并給用戶一個檢查局部變量的機會，相當于在M文件指定的行號前加入了一個keyboard命令.調試命令:除了采用調試器調試程序外，MATLAB還提供了一些命令用于程序調試.命令的功能和調試器菜單命令類似，具體使用方法請讀者查詢MATLAB幫助文檔.1緒論1.1例題解答例1.1計算,.解:創(chuàng)建符號函數:>>symsx;>>f=sym('sin(x)')f=sin(x)展開至7階泰勒級數：>>h=taylor(f,8,0)h=x-1/6*x^3+1/120*x^5-1/5040*x^7求泰勒級數在處的函數值：>>subs(h,x,0.5)ans=0.479425533234127也可以通過內聯函數來求解:>>H=inline(h)H=Inlinefunction:H(x)=x-1./6.*x.^3+1./120.*x.^5-1./5040.*x.^7>>feval(H,0.5)ans=0.479425533234127例1.2計算積分值.解:解法一:(符號法)：>>I=int('1/(1+x)','x',0,1)I=log(2)解法二:(數值法)：>>x=0:0.2:1;%將[0,1]等分為4等份>>f=1./(1+x);%分別計算每一個等分點的函數值>>I=0;>>fori=1:5I=I+(f(i)+f(i+1))/2*0.2;%將每一小曲邊的梯形累加起來作為積分值End>>vpa(I,9)%取結果的小數精度為9位小數ans=.695634921例1.3略例1.4不用開平方根計算的值.解:解法一(符號法):>>A=sym('a');>>sqrt(A)ans=a^(1/2)解法二(數值法):按以下迭代公式迭代計算近似值:建立函數文件msqrt.mfunctionx=msqrt(x0,a)%用迭代法近似計算平方根%x0為初始迭代值,a為開平方數formatlong;x=zeros(20,1);x(1)=x0;fori=2:20x(i)=1/2*(x(i-1)+a/x(i-1));enddisp(x);用編寫的函數計算,:>>msqrt(2,3);2.0000000000000001.7500000000000001.7320508100147271.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.7320508075688771.732050807568877上述結果為迭代過程計算的中間結果,分析數據可知迭代收斂速度快,只需四次計算即可計算出較為準確的數值.例1.5略.例1.6計算,視已知數為精確數,用4位浮點數計算.解:直接在Matlab中輸入式子:>>1/759-1/760ans=1.7336e-006若先轉化為浮點數再運算可得:>>a=1/759,b=1/760,a-ba=0.0013b=0.0013ans=1.7336e-006可見Matlba在計算時,數據結構都取為雙精度而提高了運算準確度.若以符號運算計算之,有:>>a=sym('1/759'),b=sym('1/760'),c=a-ba=1/759b=1/760c=1/576840可見符號運算準確但耗費運算時間.例1.7略.例1.8解方程.解:符號法解方程:>>x=solve('x^2-18*x+1','x')x=9+4*5^(1/2)9-4*5^(1/2)將結果保留小數點6位:>>vpa(x,6)ans=17.9443.5572e-11.2Matlab中數值計算精度1.Matlab中有三種運算精度,它們分別為數值算法、符號算法和可控精度算法,將它們分別介紹如下:數值算法把每個數取為16位,計算按浮點運算進行,它是運算速度最快的一種算法.符號算法把每個數都變?yōu)榉柫?運算按有理量計算進行,它的優(yōu)點是能夠得到精確結果,缺點是占用空間大,并且運算速度最慢.可控精度算法介于上述兩種算法之間,它能夠使運算在可控的精度下進行計算.2.Matlab的數據顯示格式,列表如下:表Matlab數據顯示格式命令命令意義舉例()formatshort短格式方式，顯示5位定點十進制數3.1416formatlong長格式方式，顯示15位定點十進制數formatshorte最優(yōu)化短格式顯示,5位加指數3.1416e+000formatlonge最優(yōu)格式,15位加指數formatshortg5位定點或浮點格式3.1416formatlongg對雙精度，顯示15位定點或浮點格式，對單精度，顯示7位定點或浮點格式formatshorteng至少5位加3位指數3.1416e+000formatlongeng16位加至少3位指數formathex十六進制格式方式400921fb54442d18formatbank銀行格式.按元、角、分（小數點后具有兩位）的固定格式3.14format++格式，以＋，—和空格分別表示中的正數，負數和零元素+format缺省時為默認短格式方式與formatshort相同3.1416formatrat分數格式形式.用有理數逼近顯示數據355/113formatloose松散格式.數據之間有空行formatcompact緊湊格式.數據之間無空行vpa(date,n)將數據date以n位有效數字顯示vpa(pi,5)=3.1416format并不影響matlab如何計算和存儲變量的值.對浮點型變量的計算，即單精度或雙精度，按合適的浮點精度進行，而不論變量是如何顯示的.對整型變量采用整型數據.整型變量總是根據不同的類（class）以合適的數據位顯示.3.Matlab的特殊變量ans:對最近輸入的反應computer:當前計算機類型eps:浮點精度flops:計算浮點操作次數，現已不再常用i:虛部單位inf:無窮大inputname:輸入參數名j:虛部單位nan:非數值nargin:輸入參數的數目nargout:輸出參數的數目(用戶定義函數)pi:圓周率realmax:最大正浮點數realmin:最小正浮點數vararginvarargout:返回參數數目(matlab函數)cputime:CPU工作時間2線性方程組的直接解法2.1例題解答例2.1用Gauss消元法解方程組:解:直接建立求解該方程組的M文件Gauss.m如下:%求解例題2.1%高斯法求解線性方程組Ax=b%A為輸入矩陣系數,b為方程組右端系數%方程組的解保存在x變量中%先輸入方程系數A=[123;275;149];b=[16-3]';[m,n]=size(A);%檢查系數正確性ifm~=nerror('矩陣A的行數和列數必須相同');return;endifm~=size(b)error('b的大小必須和A的行數或A的列數相同');return;end%再檢查方程是否存在唯一解ifrank(A)~=rank([A,b])error('A矩陣的秩和增廣矩陣的秩不相同,方程不存在唯一解');return;end%這里采用增廣矩陣行變換的方式求解c=n+1;A(:,c)=b;%%消元過程fork=1:n-1A(k+1:n,k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/A(k,k))*A(k,k:c);End%%回代結果x=zeros(length(b),1);x(n)=A(n,c)/A(n,n);fork=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);end%顯示計算結果disp('x=');disp(x);直接運行上面的M文件或在Matlab命令窗口中直接輸入Gauss即可得出結果.在Matlab命令窗口中輸入Gauss得出結果如下:>>Gaussx=2.00001.0000-1.0000擴展:Matlab求解線性方程的幾種命令如下(方程組的一般形式可用矩陣和向量表示成:,但運用下列方法的前提必須保證所求解的方程為恰定方程,即方程組存在唯一的一組解):運用求逆思想:或;左除法,原理上是運用高斯消元法求解,但Matlab在實際執(zhí)行過程中是通過分解法進行的(即先將矩陣A作分解,再回代計算):;用符號矩陣法計算,這種計算方法最接近精確值,但計算速度最慢:通過將矩陣施行初等行變換化成行簡化階梯形的辦法,可以這樣實現之:;.上面四種常用的辦法示例如下:>>A=[123;275;149]%上面示例方程組系數A=123275149>>b=[16-3]'%方程組右端的系數b=16-3>>x1_1=inv(A)*b,x1_2=A^(-1)*b%方法一,求逆思想x1_1=2.00001.0000-1.0000x1_2=2.00001.0000-1.0000>>x2=A\b%方法二,左除思想x2=21-1>>x3=sym(A)\sym(b)%方法三,符號法x3=21-1>>C=[A,b],rref(C)%方法四,行簡化階梯形思想,最后輸出結果的一列為解C=12312756149-3ans=10020101001-1例2.3用消元法思想求的逆陣.解:解法一:直接建立求解的M文件:Gauss_Jordan.m,源程序如下:%Gauss-Jordan法求例2.3clc;A=[1-10;223;-121];A1=A;%先保存原來的方陣A[n,m]=size(A);ifn~=merror('A必須為方陣');return;endA(:,n+1:2*n)=eye(n);%構造增廣矩陣fork=1:n[l,m]=max(abs(A(k:n,k)));%按列選主元ifA(k+m-1,k)==0error('找到列最大的元素為零,錯誤');return;endifm~=1%交換Temp=A(k,:);A(k,:)=A(k+m-1,:);A(k+m-1,:)=Temp;endfori=1:nifi~=kA(i,:)=A(i,:)-A(k,:)*A(i,k)/A(k,k);endendendfori=n:(-1):1A(i,:)=A(i,:)/A(i,i);endA(:,1:n)=[];disp('A=');disp(A1);disp('用Gauss-Jandan算得矩陣A的逆矩陣為:');disp('inv(A)=');disp(A);clearTempiklmn;%清除臨時變量在Matlab命令窗口中輸入Gauss_Jordan回車后得到結果如下:A=1-10223-121用Gauss-Jandan算得矩陣A的逆矩陣為:inv(A)=-41-3-51-36-14解法二:用通過將矩陣施行初等行變換化成行簡化階梯形的辦法,可以借助命令求解,非常方便:輸入矩陣:>>A=[1-10;223;-121]A=1-10223-121>>C=[A,eye(length(A))]C=1-10100223010-121001>>invA=rref(C)invA=100-41-3010-51-30016-14后三列即為其逆陣,檢驗其正確性:>>c=invA(:,4:6)c=-41-3-51-36-14>>d=c*Ad=100010001可見求解正確.例2.4用分解LU的方法求解方程組解:解線性方程組中分解的L,U可以實現矩陣A的三角分解,使得A=L*U.L,U應該是下三角和上三角矩陣的,這樣才利于回代求根.但是MATLAB中的LU分解與解線性方程組中的L,U不一樣.MATLAB的LU分解命令調用格式為:[L,U]=lu(A)MATLAB計算出來的L是"準下三角"(交換L的行后才能成為真正的下三角陣),U為上三角矩陣,但它們還是滿足A=L*U的.先錄入矩陣系數.>>A=[2426;49615;26918;6151840]A=242649615269186151840>>b=[9232247]'b=9232247將A作分解,方法是使用矩陣分解的LU命令即可:>>[L,D]=lu(A)L=0.33331.0000-0.66671.00000.66671.0000000.3333-1.00001.000001.0000000U=6.000015.000018.000040.00000-1.0000-6.0000-11.666700-3.0000-7.0000000-0.3333再檢驗其正確性:>>C=L*UC=242649615269186151840解方程組>>y=L\by=47.0000-8.3333-2.00000.3333解方程組得到方程組的最終解:>>x=U\yx=0.50002.00003.0000-1.0000故方程組的最終解為:例2.5解方程組試用平方根法,改進的平方根法和追趕法分別解之.解:先輸入矩陣:>>A=[610;141;0114]A=6101410114>>b=[624322]'b=624322平方根法:先對A矩陣作分解:>>L=chol(A)L=2.44950.4082001.95790.5108003.7066檢驗其正確性>>L'*Lans=6.00001.000001.00004.00001.000001.000014.0000將L轉化為下三角矩陣:>>L=L'L=2.4495000.40821.9579000.51083.7066解方程組:>>y=L\by=2.449511.747385.2526再解方程組,得到最終解:>>x=L'\yx=1.0000-0.000023.0000故平方根法解上述方程的結果為:改進的平方根法:先對矩陣A作LDL分解:>>[L,D]=ldl(A)L=1.0000000.16671.0000000.26091.0000D=6.00000003.833300013.7391檢驗其分解正確性:>>L*D*L'ans=6101410114解方程組:>>y=L\by=623316解方程組:>>x=(D*L')\yx=1023故改進的平方根法解上述方程的結果亦為:追趕法:編制追趕法求解該方程的程序如下:%pursue.m%三對角線性方程組的追趕法解方程組例2.5%輸入矩陣clc;A=[610;141;0114]f=[624322][n,m]=size(A);%分別取對角元素a=zeros(1,n);a(2:n)=diag(A,-1);c=diag(A,1);%此處用變量d存儲A主對角線上的元素,因已用變量b存儲方程右邊的系數b=diag(A);ifb(1)==0error('主對角元素不能為0');return;end%初始計算,式(2.31)alpha(1)=b(1);beta(1)=c(1)/b(1);%按照公式(2.31)計算fori=2:n-1alpha(i)=b(i)-a(i)*beta(i-1);ifalpha(i)==0error('錯誤:在解方程過程中α為0');return;endbeta(i)=c(i)/alpha(i);end%對最后一行作計算alpha(n)=b(n)-a(n)*beta(n-1);ifalpha(n)==0error('錯誤:在解方程過程中最后一個α為0');return;end%以下按照公式(2.32)計算,解Ly=fy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/alpha(i);end%以下按照公式(2.33)計算,解Ux=yX(n)=y(n);fori=n-1:-1:1X(i)=y(i)-beta(i)*X(i+1);enddisp('X=');disp(X);在Matlab命令窗口輸入pursue計算結果如下:>>pursueA=6101410114f=624322X=1023其中A為系數矩陣,f為矩陣右端的系數,最后計算結果為X.由以上計算可知追趕法解該方程的結果亦為:例2.6,求.解:輸入矩陣:>>x=[10.5-0.3]x=1.00000.5000-0.3000計算其1-范數>>norm_1=norm(x,1)norm_1=1.8000計算其2-范數>>norm_2=norm(x)norm_2=1.1576計算其無窮大范數>>norm_inf=norm(x,inf)norm_inf=1還可以計算其無窮小范數(即各分量絕對值中的最小值)>>norm_minusInf=norm(x,-inf)norm_minusInf=0.3000例2.7求.解:先輸入矩陣:>>A=[-210;1-21;01-2]A=-2101-2101-2求A的1-范數(列和范數):>>norm_1=norm(A,1)norm_1=4求解A的無窮大范數(行和范數):>>norm_inf=norm(A,inf)norm_inf=4求A的2-范數(最大特征值):>>norm_2=norm(A,2)norm_2=3.4142還可以求解A的F-范數:>>norm_F=norm(A,'fro')norm_F=4譜半徑可以通過按其概念進行計算:對其特征值的絕對值取最大值即可:>>R_A=max(abs(eig(A)))R_A=3.4142例2.8n階Hilbert矩陣考查其條件數.解:上述矩陣為抽象矩陣,而計算機只能進行有限次迭代,我們從考查其條件數,為此編制如下程序Hilb_cond10.m:%Hilb_cond10.m%考查從3階至10的矩陣2—條件數fori=3:10h=hilb(i);condA(i)=cond(h,2);enddisp('ncond');fori=3:10s=sprintf('%d%f',i,condA(i));disp(s);end運行后得到如下結果:ncond3524.056778415513.7387395476607.250242614951058.6410057475367356.3703939493153986466.270940結果中左邊為的階數,右邊為對應的條件數,從以上結果也可分析可知:當n的階數稍高時,矩陣即出現嚴重的病態(tài).例2.9求的條件數解:建立2階矩陣:>>H=hilb(2)H=1.00000.50000.50000.3333求其2-條件數>>cond_2=cond(H,2)cond_2=19.2815求其1-條件數>>cond_1=cond(H,1)cond_1=27.0000求其無窮大條件數>>cond_inf=cond(H,inf)cond_inf=27.0000還可求其F-條件數>>cond_f=cond(H,'fro')cond_f=19.33332.2Matlab解線性方程組常用命令介紹1.求秩命令rank在解線性方程組時,為了判斷是否存在解,我們應先判斷矩陣的秩.調用格式為:c=rank(A)2.矩陣零空間向量null當線性方程組的秩時,方程組有無窮多個解,使用x=null(A)可得到滿足的一個解向量,可為符號矩陣或數值矩陣:x=null(A)或x=null(sym(A))3.方陣的LU分解命令lu調用格式為:[L,U]=lu(A)L為準下三角矩陣,U為上三角矩陣,滿足A=LU.[L,U,P]=lu(A)L為下三角矩陣,U為上三角矩陣,P為變換方陣元素位置的換位陣,滿足PA=PL.其它調用格式請用helplu獲得更多信息.4.Cholsky分解cholL=chol(A)其中L為一個下三角矩陣,滿足.必須為正定矩陣.5.條件數condc=cond(A,p)A可為向量或矩陣,P取值為下列之一:1——向量或矩陣的返回1—條件數.2——返回向量或矩陣的2—條件數.inf——返回向量或矩陣的—條件數.'fro'——返回向量或矩陣的F—條件數.6.奇異值分解svd[U,S,V]=svd(A)將A分解為正交矩陣U,對角矩陣S和正交矩陣V的乘積,使得A=USVT.7.線性方程組的符號解linsolveX=linsolve(A,b)它等價于X=sym(A)\sym(b).返回方程組的符號解.3線性方程組的迭代解法3.1例題解答例3.1用Jacobi迭代法解以下方程:解:編制迭代計算的M文件程序如下:%Jacobi迭代法求解例3.1%A為方程組的增廣矩陣clc;A=[10-2-13;-210-115;-1-2510]MAXTIME=50;%最多進行50次迭代eps=1e-5;%迭代誤差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%進入迭代計算disp('迭代過程X的值情況如下:')disp('X=');while1disp(x');fori=1:1:ns=0.0;forj=1:1:nifj~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfori=1:1:nmaxeps=max(0,abs(x(i)-y(i)));%檢查是否滿足迭代精度要求endifmaxeps<=eps%小于迭代精度退出迭代fori=1:1:nx(i)=y(i);%將結果賦給xendreturn;endfori=1:1:n%若不滿足迭代精度要求繼續(xù)進行迭代x(i)=y(i);y(i)=0.0;endk=k+1;ifk>MAXTIME%超過最大迭代次數退出error('超過最大迭代次數,退出');return;endend運行該程序結果如下:A=10-2-13-210-115-1-2510迭代過程X的值情況如下:X=0000.30001.50002.00000.80001.76002.66000.91801.92602.86400.97161.97002.95400.98941.98972.98230.99621.99612.99380.99861.99862.99770.99951.99952.99920.99981.99982.99970.99991.99992.99991.00002.00003.00001.00002.00003.0000容易看出迭代計算最后結果為:例3.2用Gauss-Seidel迭代法計算例3.1并作比較.解:編制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m%Gauss_Seidel迭代法求解例3.2%A為方程組的增廣矩陣clc;formatlong;A=[10-2-13;-210-115;-1-2510][n,m]=size(A);%最多進行50次迭代Maxtime=50;%控制誤差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次數小于最大迭代次數,進入迭代fork=1:Maxtimedisp(x);fori=1:ns=0.0;forj=1:nifi~=js=s+A(i,j)*x(j);%計算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此時迭代的值end%因為方程的精確解為整數,所以這里將迭代結果向整數靠近的誤差作為判斷迭代是否停止的條件ifsum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代結果:');Xformatshort;完成后直接在Matlab命令窗口中輸入Gauss_Seidel回車后可得到如下結果:>>Gauss_SeidelA=10-2-13-210-115-1-2510x=0000.3000000000000001.5600000000000002.6840000000000000.8804000000000001.9444800000000002.9538720000000000.9842832000000001.9922438400000002.9937541760000000.9997021448448001.9998545228697602.9998822381168640.9999591283856381.9999800494888142.9999838454726540.9999943944450281.9999972634362712.9999977842635140.9999992311136061.9999996246490722.9999996960823500.9999998945380491.9999999485158452.9999999583139480.9999999855345641.9999999929383082.9999999942822360.9999999980158851.9999999990314012.9999999992157370.9999999997278541.9999999998671452.9999999998924290.9999999999626721.9999999999817782.9999999999852460.9999999999948801.9999999999975012.9999999999979760.9999999999992981.9999999999996572.9999999999997220.9999999999999041.9999999999999532.9999999999999620.9999999999999871.9999999999999942.9999999999999950.9999999999999981.9999999999999992.9999999999999991.0000000000000002.0000000000000003.000000000000000迭代結果:X=123可見對此題Gauss_Seidel法的收斂速度還是很快的.例3.3取,用Gauss-Seidel迭代法和SOR法計算下列方程組的解:解:Gauss-Seidel迭代法可利用上題中的程序,把輸入矩陣A換掉就可以了,以下編制求解該程序的SOR.m%SOR法求解例3.3%w=1.45%方程組系數矩陣clc;A=[4-2-1;-24-2;-1-23]%方程組右端系數b=[0,-2,3]'w=1.45;%最大迭代次數Maxtime=100;%精度要求Eps=1E-5;%以15位小數顯示formatlong;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代過程:');disp('x=');while1y=x;disp(x');%計算過程fori=1:ns=b(i);forj=1:nifj~=is=s-A(i,j)*x(j);endendifabs(A(i,i))<1E-10|k>=Maxtimeerror('已達最大迭代次數或矩陣系數近似為0,無法進行迭代');return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endifnorm(y-x,inf)<Eps%達到精度要求退出計算break;endk=k+1;enddisp('最后迭代結果:');%最后的結果X=x'%設為默認顯示格式formatshort;為了能有可比性,程序中的誤差控制改用無窮大范數來度量(即將其改為:norm(y-x,inf)<Eps,并在循環(huán)中用y先保存迭代前x的值:y=x),此外,兩個程序中的迭代初始值.以下為兩種方法運行結果:(1)法運行結果如下:>>Gauss_SeidelA=4-2-10-24-2-2-1-233x=1110.7500000000000000.3750000000000001.5000000000000000.9322631785694630.9222812405563841.9256085532274100.9425427585850440.9340756559062271.9368980234658330.9586586469134330.9525664386408791.9545971747317300.9702