MATLAB8.X程序設(shè)計及典型應(yīng)用第四章

第四章數(shù)值計算數(shù)值計算在現(xiàn)代科學研究和工程技術(shù)領(lǐng)域中應(yīng)用最為廣泛。MATLAB正是憑借其卓越的數(shù)值計算功能而稱雄一方?，F(xiàn)在MATLAB擁有的數(shù)值計算功能已經(jīng)非常強大。本章重點講述在MATLAB工作環(huán)境下解決具體的數(shù)值計算問題，包含高等數(shù)學、線性代數(shù)、數(shù)據(jù)分析中的很多重要內(nèi)容：矩陣的計算和分解、多項式運算、函數(shù)的零極點、數(shù)值微積分等。矩陣的多種運算方法MATLAB實現(xiàn)計算矩陣的秩、特征值及其對應(yīng)的特征向量利用矩陣操作求解線性方程組利用多項式實現(xiàn)數(shù)據(jù)點的插值和擬合利用ode()指令求初值問題的常微分方程數(shù)值解本章主要內(nèi)容包括：4.1矩陣的計算

4.1.1矩陣的結(jié)構(gòu)變換

1.矩陣轉(zhuǎn)置的MATLAB運算符為單引號“’”。?A'?求A的轉(zhuǎn)置，其中A可以是行向量、列向量和矩陣?！纠?-1】矩陣轉(zhuǎn)置的MATLAB實例。>>clear,x1=[1;3;5];%3×1的列向量>>y1=x1'%轉(zhuǎn)置結(jié)果為1×3的行向量y1=135>>x2=[123;456];%2×3的矩陣>>y2=x2'%轉(zhuǎn)置結(jié)果為3×2的矩陣y2=142536>>x3=[1-2*i3+4*i;1-4*i2+6*i;3-i-5+4*i];%3×2的復(fù)數(shù)矩陣>>y3=x3'%轉(zhuǎn)置結(jié)果為2×3的共軛矩陣y3=1.0000+2.0000i1.0000+4.0000i3.0000+1.0000i3.0000-4.0000i2.0000-6.0000i-5.0000-4.0000i>>x3=[1-2*i3+4*i;1-4*i2+6*i;3-i-5+4*i];%3×2的復(fù)數(shù)矩陣>>y3=x3.‘%的點轉(zhuǎn)置結(jié)果為2×3的非共軛矩陣y3=1.0000-2.0000i1.0000-4.0000i3.0000-1.0000i3.0000+4.0000i2.0000+6.0000i-5.0000+4.0000i4.1.1矩陣的結(jié)構(gòu)變換對于復(fù)數(shù)矩陣，MATLAB完成的是復(fù)數(shù)矩陣的共軛轉(zhuǎn)置。復(fù)數(shù)矩陣的非共軛轉(zhuǎn)置的運算符為“.‘”。注意：矩陣轉(zhuǎn)置的操作符必須在英文狀態(tài)下輸入。4.1.1矩陣的結(jié)構(gòu)變換指令flipud()和fliplr()完成矩陣的對稱變換，調(diào)用格式如下：?B=flipud(A)?上下方向翻轉(zhuǎn)的矩陣。如果是列向量，返回相反順序的向量；如果是行向量，返回原向量。?B=fliplr(A)?水平方向翻轉(zhuǎn)的矩陣。如果是行向量，返回相反順序的向量；如果是列向量，返回原向量。另外，MATLAB還提供了一個以指定維翻轉(zhuǎn)的函數(shù)flipudim().格式為：?B=flipdim(A,dim)?返回以指定的維翻轉(zhuǎn)矩陣。dim=1，以行方向翻轉(zhuǎn)；dim=2，以列方向翻轉(zhuǎn)?？梢奻lipdim(A,1)操作效果等同于flipud(A),flipdim(A,2)操作效果等同于fliplr(A)。2.對稱變換>>clear,B=[1234;5678;9101112;13141516]B=12345678910111213141516>>flipud(B)%將矩陣B上下方向翻轉(zhuǎn)矩陣ans=13141516910111256781234>>fliplr(B)%將矩陣B水平方向翻轉(zhuǎn)ans=43218765121110916151413【例4-2】矩陣的結(jié)構(gòu)變換實例。指令rot90()可以完成矩陣逆時針旋轉(zhuǎn)90度。

?rot90(A,k)

?矩陣A逆時針翻轉(zhuǎn)k×90度?！纠?-3】矩陣旋轉(zhuǎn)變換實例。>>clear,x=[12;34;56];>>A1=rot90(x)%逆時針旋轉(zhuǎn)90度A1=246135>>A2=rot90(x,2)%逆時針旋轉(zhuǎn)2×90度A2=654321>>A3=rot90(x,3)%逆時針旋轉(zhuǎn)3×90度A3=5316423.旋轉(zhuǎn)MATLAB中提取上三角矩陣的函數(shù)為triu()，提取下三角矩陣的函數(shù)為tril()，格式如下：?triu(X,K)?K=0時，提取X的主對角線及以上的元素。K>0時，提取矩陣X主對角線上方的第K條對角線以上的元素。K<0時，提取矩陣X的主對角線下方第-K條主對角線以上的元素。K缺省時默認為0。?tril(X,K)?K=0時，提取X的主對角線及以下的元素。其中K>0時，提取矩陣X主對角線上方的第K條對角線以下的元素。K<0時，提取矩陣X的主對角線下方第-K條主對角線以下的元素。K缺省時默認為0。4.提取三角陣>>clear,x=magic(3);%創(chuàng)建3階魔方矩陣>>A=triu(x)%提取主對角線上三角矩陣A=816057002>>B=tril(x)%提取主對角線下三角矩陣B=800350492>>B1=tril(x,1)%提取主對角線上方第一條對角線下三角矩陣B1=810357492【例4-6】提取三角陣指令實例。4.1.2矩陣分析矩陣分析是線性代數(shù)中關(guān)于矩陣運算的重要環(huán)節(jié)。矩陣分析包括矩陣的秩、矩陣對應(yīng)的行列式和逆矩陣等運算，MATLAB提供了相應(yīng)的指令。用戶在計算中只需要正確調(diào)用這些指令，就可以快速地得到結(jié)果。矩陣中線性無關(guān)的行數(shù)和列數(shù)成為矩陣的秩。在MATLAB中求秩函數(shù)為rank（）。rank的格式為：

?rank(A)

?計算矩陣A的秩。

【例4-7】求下列矩陣的秩（1）（2）（1）>>clear,a=[142;672;5-4-8];>>R1=rank(a)R1=3（2）>>b=[54-2;452;-228];>>R2=rank(b)R2=21秩MATLAB還提供了一個化簡函數(shù)rref()，其格式如下：

?rref(A)

?將矩陣A化為行階梯陣。

根據(jù)MATLAB的輸出結(jié)果，用戶可以得到矩陣的秩，即行階梯陣中非零行的行數(shù)【例4-8】（續(xù)例4-7）行階梯陣化簡函數(shù)rref（）應(yīng)用實例。>>rref(a)%將矩陣a化簡為行階梯陣ans=100010001>>rref(b)%將矩陣b化簡為行階梯陣ans=10-2012000由化簡結(jié)果可知，矩陣a的秩為3，矩陣b的秩為2.。1秩把方陣A看做行列式|A|，對其按照行列式的規(guī)則求值，稱該值為行列式的值。MATLAB的實現(xiàn)指令為det()，其格式為：?det(A)?計算方陣A對應(yīng)的行列式的值行列式的值為0時，相應(yīng)的矩陣成為奇異矩陣，否則稱為非奇異矩陣（或滿秩矩陣）。對于非奇異矩陣A，有與其同型的非奇異矩陣B，使得A×B=B×A=I，其中I為單位矩陣，則A和B互為逆矩陣。求方陣A逆矩陣的MATLAB函數(shù)為inv()，格式為：?B=inv(A)?求非奇異矩陣A的逆矩陣B2．行列式和逆矩陣>>det(a)ans=66>>det(b)%計算矩陣對應(yīng)行列式的值ans=0>>inv(a)%矩陣a為非奇異矩陣，存在逆矩陣ans=-0.72730.3636-0.09090.8788-0.27270.1515-0.89390.3636-0.2576>>pinv(b)%矩陣b為奇異矩陣，存在偽逆矩陣ans=0.06170.0494-0.02470.04940.06170.0247-0.02470.02470.0988【例4-9】（續(xù)例4-7）矩陣對應(yīng)行列式的值、逆矩陣和廣義逆矩陣運算實例。對于奇異矩陣，調(diào)用MATLAB指令inv()也可以獲得計算結(jié)果.此時MATLAB會給出警告：Warning:Matrixisclosetosingularorbadlyscaled.用戶不要將輸出矩陣認為是該奇異矩陣的逆矩陣，必須通過逆矩陣定義加以驗證。【例4-10】逆矩陣的應(yīng)用：設(shè)A=,B=,C=,求矩陣X，使?jié)M足:?！菊f明】由線性代數(shù)可知：，或者。本題給出兩種計算程序，并比較計算結(jié)果和耗時。編寫文件名為exm4_10的腳本文件：clear,clc,A=[123;234;343];B=[2,5;43];C=[13;40;21];%顯示利用逆矩陣求解時的耗時tic%開啟計時器X1=inv(A)*C*inv(B)toc%結(jié)束計時，并輸出時間%顯示利用矩陣除法求解時的耗時ticX2=A\C/Btoc【例4-10】逆矩陣的應(yīng)用在指令窗中執(zhí)行exm4_11：X1=-4.75004.25004.3571-3.9286-1.10711.1786Elapsedtimeis0.000201seconds.X2=-4.75004.25004.3571-3.9286-1.10711.1786Elapsedtimeis0.000137seconds.【例4-10】逆矩陣的應(yīng)用利用逆矩陣求解方程組時需要兩步計算，即求逆矩陣和矩陣乘法，因此較左除法求解方程組速度要慢。程序在不同計算機上執(zhí)行時耗時不唯一，但時長順序不會改變。如果矩陣A不是方陣或是奇異矩陣時，該矩陣存在與A的轉(zhuǎn)置矩陣同型的矩陣B，使A×B×A=A,B×A×B=B成立，則B稱為A的逆，即偽逆矩陣。求偽逆矩陣的MATLAB函數(shù)為pinv(),其格式為：?B=pinv(A)?求矩陣A的偽逆矩陣B【說明】分析：該方程組為欠定方程組，系數(shù)矩陣不是方陣，因此只能求其偽逆矩陣。編寫文件名為exm4_11的腳本文件：clear,clcA=[2171;8-326];b=[35]';X=pinv(A)*b在指令窗中執(zhí)行exm4_11：X=0.3426-0.06910.30630.2400【例4-11】求方程組的最小范數(shù)解。矩陣的跡：矩陣的對角線元素之和，即矩陣的特征值之和。MATLAB中求跡函數(shù)是trace()，格式如下：?trace(A)?求矩陣A的跡。向量種度量形式。使用范數(shù)可以測量兩個向量或矩陣之間的距離。范數(shù)有多種定義形式。MATLAB計算向量或矩陣的范數(shù)指令為norm()，格式如下：?norm(A,k)?計算向量或矩陣A的k—范數(shù):k可以取1、2、Inf、'fro',分別為計算矩陣的1-范數(shù)、2-范數(shù)、無窮大-范數(shù)、Frobenius-范數(shù),k缺省時為計算2-范數(shù)。3.矩陣的跡和范數(shù)>>clear,A=[86-11;43-12;53-34;34-7-2];>>trace(A)%計算矩陣的跡ans=6>>norm(A,1)%計算矩陣的1-范數(shù)ans=20>>norm(A)%計算矩陣的2-范數(shù)ans=14.6918【例4-12】計算矩陣的跡和范數(shù)。>>norm(A,inf)%計算矩陣的無窮大-范數(shù)ans=16>>norm(A,'fro')%計算矩陣的Frobenius-范數(shù)ans=16.4012>>norm(A(:,2),1)%計算向量的1-范數(shù)ans=16【例4-12】計算矩陣的跡和范數(shù)。矩陣的條件數(shù)是解線性方程組是判斷系數(shù)矩陣的變化對解的影響程度的一個參數(shù)。矩陣的條件數(shù)的定義依賴于范數(shù)。矩陣的條件數(shù)接近1時表示方程為良性的，否則為病態(tài)。MATLAB計算矩陣條件數(shù)的指令為cond(),格式如下：?cond(A,p)?計算矩陣A的p—范數(shù):p可以取1、2、Inf、'fro',分別為計算矩陣的1-范數(shù)條件數(shù)、2-范數(shù)條件數(shù)、無窮大-范數(shù)條件數(shù)、Frobenius-范數(shù)條件數(shù),p缺省時為計算2-范數(shù)條件數(shù)。4.條件數(shù)>>cond(A,1)%計算矩陣的1-范數(shù)條件數(shù)ans=108.1720>>cond(A)%計算矩陣的2-范數(shù)條件數(shù)ans=51.7487>>cond(A,inf)%計算矩陣的無窮大-范數(shù)條件數(shù)ans=67.0968>>cond(A,'fro')%計算矩陣的Frobenius-范數(shù)條件數(shù)ans=58.0285【例4-13】（續(xù)例4-12）計算矩陣A的條件數(shù)。4.1.3矩陣的特征值分析設(shè)n階方陣A，如果數(shù)和n維非零列向量，使關(guān)系式成立，則數(shù)稱為方陣A的特征值，非零向量稱為矩陣A的對應(yīng)于特征值的特征向量，MATLAB計算矩陣特征值和特征向量的函數(shù)為eig()，格式為：?[V,D]=eig(A)?計算矩陣A的特征值D和對應(yīng)的特征向量V，使得AV=VD。如果函數(shù)只有一個輸出宗量，則只給出特征值。>>clear,A=magic(3);%創(chuàng)建3階魔方矩陣>>d=eig(A)d=15.00004.8990-4.8990>>[V,D]=eig(A)V=-0.5774-0.8131-0.3416-0.57740.4714-0.4714-0.57740.34160.8131D=15.00000004.8990000-4.8990【例4-14】計算3階魔方矩陣的特征值和特征向量。4.1.4矩陣的分解矩陣分解是線性代數(shù)中比較重要的內(nèi)容，針對不同的分解方法，MATLAB提供了多個關(guān)于矩陣分解的函數(shù)，具體見表4-1所示。分解類別調(diào)用格式功能LU分解[L,U]=lu(A)計算上三角矩陣U和交換下三角矩陣L，使LU=A[L,U,P]=lu(A)計算上三角矩陣U、有單位對角線的下三角矩陣L和交換矩陣P，滿足LU=PACholesky因式分解R=chol(A)計算正定矩陣A的Cholesky因子，是一個上三角矩陣，滿足R'*R=A。如果A不是一個正定矩陣，則給出一個錯誤信息。[G,p]=chol(A)計算矩陣A的Cholesky因子G。如果A不是一個正定矩陣，則不給出錯誤信息，而是將p設(shè)為正整數(shù)。qr因式分解[Q,R]=qr(A)計算m×n矩陣A分解得到的m×n酉矩陣Q和m×n的上三角矩陣R。Q的列形成了一個正交基。Q和R滿足A=Q*R[Q,R,P]=qr(A)計算矩陣Q、上三角矩陣R和交換矩陣P。Q的列形成一個正交基，R的對角線元素按大小降序排列，滿足A*P=Q*RSVD分解s=svd(A)計算一個包含矩陣A的奇異值的向量s[U,S,V]=svd(A)計算奇異值為對角線且與A同型的矩陣S，以及m×m和n×n的酉矩陣U和V，使A=U*S*V'表4-1矩陣分解函數(shù)及其功能表分析：方程組形式為AX=b，A=LU，即LUX=b，由矩陣除法可知X=U\(L\b)。>>clear,clc,>>A=[32-1;530];b=[2;7];[L,U]=lu(A);>>X=U\(L\b)X=1.400002.2000【例4-15】利用LU分解法求欠定方程組的一個特解。4.1.5線性方程組的求解線性方程組分為齊次方程組和非齊次方程組兩類。對于系數(shù)矩陣為滿秩矩陣的方程組，可以將方程組表述為Ax=b,其解為x=A-1b.MATLAB語句為x=inv(A)*b或者x=A\b.【例4-16】解方程組>>clear,A=[21-11;53-12;23-24;33-32];>>det(A)%由行列式值是否為零判斷A是否存在逆矩陣ans=16>>b=[2353]';>>x=inv(A)*bx=0.5000-1.0000-0.50001.5000由執(zhí)行結(jié)果可知，方程組的解為：x1=0.5;x2=-1;x3=-0.5;x4=1.5.該題也可以通過執(zhí)行語句x=A\b求得?！纠?-16】解方程組齊次線性方程組求解對于系數(shù)陣A為奇異陣的齊次線性方程組，可以利用函數(shù)rref(A)將系數(shù)陣化為行階梯陣，對照化簡結(jié)果可以直接寫出方程組的解?！纠?-17】解方程組對于系數(shù)矩陣為非滿秩矩陣，即奇異矩陣，MATLAB提供了針對齊次線性方程組和非齊次線性方程組的不同求解函數(shù)。MATLAB求解過程如下：>>clear,A=[211-1;422-2;132-4;1222];>>formatrat%有理格式輸出>>x=rref(A)x=100-4/3010-600123/30000由執(zhí)行結(jié)果可知，該方程組有無窮多解，其通解為：【例4-17】解方程組MATLAB還提供了函數(shù)null（）用來求解齊次線性方程組解空間的一組基（即基礎(chǔ)解系），格式如下：?X=null（A，‘r’）?求系數(shù)矩陣為A的齊次方程組一組基礎(chǔ)解系（即一組基）。r是可選參數(shù)：當有r時，輸出有理基；當無r時，輸出正交規(guī)范基?！纠?-18】（續(xù)例4-17）利用函數(shù)null（）求解齊次線性方程組示例。>>x=null(A,'r')x=4/36-23/31則方程組的解為。結(jié)果與【例4-14】一致?！纠?-17】解方程組【例4-19】利用函數(shù)null（）求解齊次線性方程組>>clear,A=[1221;21-2-2;1-1-4-3];>>formatrat%有理式格式輸出>>x=null(A,'r')%求解空間的有理基x=25/3-2-4/31001則方程組的解為。非齊次線性方程組可以表述為AX=b，其中A為系數(shù)矩陣，其解存在3種可能：（1）無解：系數(shù)矩陣A的秩n小于增廣陣(A,b)的秩，則方程組無解；（2）唯一解：系數(shù)陣A的秩n等于增廣陣(A,b)的秩且等于未知變量的個數(shù)m（3）無窮多解：系數(shù)陣A的秩n等于增廣陣(A,b)的秩且小于未知變量的個數(shù)m，則方程組有無窮多解。對于情況（3），MATLAB必須要給出該方程的一個特解和對應(yīng)齊次方程的(n-m)個通解

求解通解用函數(shù)null()，特解可以利用矩陣LU分解法求得。2.非齊次線性方程組求解functionX=equs(A,b,n)%EQUSSolvingthehomogeneousorinhomogeneouslinearequations%A方程組的系數(shù)矩陣%b常數(shù)項對應(yīng)的列向量%n方程組中變量個數(shù)%X線性方程組的解ifnargin>3error('輸入宗量太多，請重新輸入')endB=[A,b];%增廣矩陣R_A=rank(A);%求系數(shù)矩陣的秩R_B=rank(B);%求增廣矩陣的秩ifR_A>n|R_B>nerror('輸入變量個數(shù)有誤，請重新輸入')end編寫一個文件名為equs()函數(shù)文件:formatrat%有理格式顯示ifR_A~=R_B%是否無解

X='Equationshavenosolution!';elseifR_A==R_B&R_A==n%是否是唯一解

X=A\b;else[L,U]=lu(A);X=U\(L\b);%求特解

C=null(A,'r')%求通解endend用戶利用該函數(shù)文件，可以用來求解任意類型的齊次或者非齊次線性方程組。編寫一個文件名為equs()函數(shù)文件:在指令窗中執(zhí)行：>>clear,>>A=[11-3-1;3-1-34;15-9-8];b=[140]';n=4;>>X=equs(A,b,n)Warning:Rankdeficient,rank=2,tol=9.0189e-015.Inequsat23C=3/2-3/43/27/41001X=00-8/153/5則該方程組的解為：【例4-20】求解方程組4.2多項式多項式運算是數(shù)學運算中最基本運算方式之一，也是線性代數(shù)分析和設(shè)計中的重要內(nèi)容。MATLAB提供了多個函數(shù)可以幫助用戶完成多項式的運算。4.2.1多項式的表達和創(chuàng)建MATLAB中將n階多項式p(x)存儲在長度為n+1的行向量p中，行向量p的元素為多項式的系數(shù)，并按x的降冪排列。行向量代表的多項式為：注意：多項式中缺少某個冪次項，則在創(chuàng)建行向量時將該冪次項的系數(shù)寫為0.【例4-21】創(chuàng)建多項式>>p=[4,-3,0,0,1,0]%三次項、二次項、常數(shù)項系數(shù)均為0p=4-300104.2.2多項式的運算

1.多項式的加法和減法運算多項式加減運算實質(zhì)上為兩個向量的加減運算。如果兩個多項式向量長度相同，標準的向量加法和減法是有效的。當兩個多項式階次不同時，必須將較短的（即低階多項式）向量前面用若干個零填補，使得兩個向量長度相同，然后再進行運算。【例4-22】將多項式與相加減。>>a=[1-206];>>b=[0123];%將低階多項式對于行向量補零>>c=a+b%進行加法運算c=1-129>>sum=poly2str(c,'x')%標準形式輸出sum=x^3-1x^2+2x+9>>d=a-b%進行減法運算d=1-3-23>>dif=poly2str(d,'x')%標準形式輸出dif=x^3-3x^2-2x+3【例4-22】將多項式與相加減。兩個多項式乘法運算實質(zhì)是兩個多項式系數(shù)的卷積運算。MATLAB中計算多項式卷積的指令為conv(),其調(diào)用格式為：?conv(a,b)?計算兩個多項式的乘積，其中a,b分別為兩個多項式的系數(shù)向量。多項式除法為乘法的逆運算，MATLAB中的除法運算指令為deconv()，其調(diào)用格式為：?[q,r]=deconv(a,b),?計算兩個多項式的除法，其中a,b為被除數(shù)多項式和除數(shù)多項式的系數(shù)向量，q、r分別為商多項式和余數(shù)多項式的系數(shù)向量。2.多項式的乘法和除法運算>>a=[1-206];>>b=[123];>>m=conv(a,b)%進行乘法運算m=10-101218>>pro=poly2str(m,'x')%標準形式輸出pro=x^5-1x^3+12x+18【例4-23】（續(xù)例4-22）運用函數(shù)conv()和deconv()進行多項式運算實例。>>[q,r]=deconv(a,b)%進行除法運算q=1-4r=00518>>q=poly2str(q,'x')%輸出商多項式的標準形式q=x-4>>r=poly2str(r,'x')%輸出余多項式的標準形式r=5x+18【例4-23】（續(xù)例4-22）運用函數(shù)conv()和deconv()進行多項式運算實例。MATLAB中求多項式全部根的函數(shù)為roots()，其格式為：?x=roots(p)?計算多項式的全部根，并以列向量形式賦給變量x。其中p為多項式的系數(shù)向量?！菊f明】在MATLAB中，多項式和根都表示為向量，其中多項式是行向量，根為列向量。如果已知多項式的根x，可以用函數(shù)poly()創(chuàng)建多項式：?p=poly(x)?由多項式的根向量x創(chuàng)建多項式系數(shù)向量p，其中p的第一個元素為1，對應(yīng)多項式最高階系數(shù)為1。3.多項式求根（1）計算f（x）=0的全部根。（2）由方程f（x）=0的根構(gòu)造一個多項式g（x），并與f（x）進行對比。>>p=[2032-5];>>x=roots(p)%求方程f(x)=0的根x=0.1414+1.5930i0.1414-1.5930i-1.14020.8573由結(jié)果可知，方程有4個根，其中有兩個實根，還有一對共軛復(fù)根。>>g=poly(x)%求多項式g(x)g=1.00000.00001.50001.0000-2.5000>>g=poly2str(g,'x')%標準形式輸出g=x^4+7.7716e-016x^3+1.5x^2+1x-2.5【例4-24】已知。計算多項式的值函數(shù)為polyval()和polyvalm()，調(diào)用格式如下：?y=polyval(P,x)?計算多項式的值。如果x為常數(shù)，則求多項式P在該點的值。計算表達式為；若x為向量或矩陣，則對向量或矩陣中的每個元素求多項式P的值，返回值為與自變量同型的向量或矩陣。?y=polyvalm(P,A)?以方陣A為自變量求多項式P的值4.多項式求值>>clear,clc>>p=[2-300-5];x1=1.2;x2=[1.215;3-29;-12/34];y1=polyval(p,x1)%計算標量對應(yīng)的函數(shù)值y1=-6.0368>>y2=polyval(p,x2)%計算矩陣中各元素對應(yīng)的函數(shù)值y2=1.0e+004*-0.0006-0.00060.08700.00760.00511.09300-0.00050.0315>>y3=polyvalm(p,x2)%計算矩陣對應(yīng)的函數(shù)值y3=-201.916833.2427838.0000-471.7920343.4400606.0000-53.2960-18.6133206.0000【例4-25】polyval()和polyvalm()使用實例：已知

,分別計算x=1.3和時的值。執(zhí)行部分分式展開的函數(shù)是residue()，格式為：?[r,p,k]=residue（b,a）?[b,a]=residue（r,p,k）?對多項式進行部分分式的展開和其逆運算。a為分式的分母系數(shù)向量，b為分式的分子系數(shù)向量。p為極點(Pole)，r為留數(shù)(Residue)，k(x)為直項(Directterm)。對于多項式b(x)和不含重根的n階多項式a(x)之比，滿足：如果a(x)有m階重根，則對應(yīng)的部分分式可以寫成：5.多項式部分分式展開>>clear,clc>>b=[23-26];>>a=[-4053];>>[r,p,k]=residue(b,a)%進行部分分式展開r=-0.81440.0322-1.5573i0.0322+1.5573ip=1.3445-0.6723+0.3254i-0.6723-0.3254ik=-0.5000【例4-26】對有理多項式進行部分分式展開>>[b1,a1]=residue(r,p,k)%進行部分分式逆運算，返回有理多項式形式b1=-0.5000-0.75000.5000-1.5000a1=1.00000.0000-1.2500-0.7500【說明】運用逆運算返回的多項式分母系數(shù)進行了歸一化處理，即將原有理分式分子分母都除以-4即可。【例4-26】對有理多項式進行部分分式展開4.3多項式插值和擬合

在大量的應(yīng)用領(lǐng)域中，人們往往需要依據(jù)離散的數(shù)據(jù)點（即測量值）得出一個解析函數(shù)的表達式：如果假定這些數(shù)據(jù)點是完全正確的，要求以某種方式描述數(shù)據(jù)點之間所分布點的函數(shù)表達式，為插值問題；倘若根據(jù)這些數(shù)據(jù)點，用戶得到某條光滑曲線，它不需要完全經(jīng)過這些點，但能最佳擬合數(shù)據(jù)點，這就是曲線擬合或者回歸。4.3.1.多項式擬合

基于離散的數(shù)據(jù)點，擬合方法不同給出的曲線表達式也不一樣。多項式擬合是用一個多項式來逼近這些給定的數(shù)據(jù)點。擬合的準則是函數(shù)在數(shù)據(jù)點處的誤差平方和為最小,即最小二乘多項式擬合.

MATLAB實現(xiàn)最小二乘多項式擬合的的函數(shù)為polyfit()，調(diào)用格式：?p=polyfit(x,y,m)?根據(jù)數(shù)據(jù)點（x,y），產(chǎn)生一個m階多項式p。其中x,y是兩個等長的向量（x要求按遞增或遞減次序排列），p是一個長度為m+1的多項式系數(shù)行向量。%輸入測量數(shù)據(jù)點clear,clcx=0:.1:1;y=[-.4601.683.286.187.187.357.669.569.409.3011.2];p1=polyfit(x,y,1);%線性擬合y1=poly2str(p1,'x')%化為標準多項式形式>>p2=polyfit(x,y,2);%二階多項式擬合>>y2=poly2str(p2,'x')在指令窗中執(zhí)行文件exm4_27.m，結(jié)果為：y1=10.3982x+1.3764y2=-10.1632x^2+20.5614x-0.14811【例4-27】多項式擬合實例。編寫文件名為exm4_27的腳本文件：圖4.1線性擬合直線和10次擬合曲線圖【例4-27】多項式擬合實例。擬合曲線能最佳擬合數(shù)據(jù)，反應(yīng)數(shù)據(jù)點中隱藏的內(nèi)在規(guī)律，擬合曲線上的函數(shù)值代表實驗測量的真值，因此，測量值與擬合線上對應(yīng)真值的差即為殘差。最小二乘原理就是要保證通過擬合直線得到的殘差平方和為最小。給定n點數(shù)據(jù)，其最大擬合階次為n-1。圖4.1給出了最高階擬合曲線，此時擬合曲線將全部通過數(shù)據(jù)點，但在左右兩邊的極值處，曲線出現(xiàn)了大的波紋，此時多項式表達的性質(zhì)不明顯。因此，計算擬合曲線不是階數(shù)越高越好。常見的多項式擬合選擇線性擬合，即選擇1階?！菊f明】4.3.2.多項式插值插值為對數(shù)據(jù)點之間函數(shù)的估值方法，即用戶要找到一個解析函數(shù)能描述自變量相鄰的兩個點之間任一位置x處y的值。MATLAB中插值函數(shù)有一維、二維和多維。1.一維數(shù)值插值一維數(shù)值插值函數(shù)為interp1()。調(diào)用格式如下：?y0=interp1(x,y,x0,'method')?根據(jù)數(shù)據(jù)點(x,y)的值，計算函數(shù)在x0處的值y0。x0是一個向量或標量，為即將插值的點；y0是一個與x0等長的插值結(jié)果。method是插值方法，可用的方法有：(1)′linear′線性插值法，如果沒有指定插值方法，MATLAB皆以這種方法進行插值。(2)′nearest′最鄰近插值法。(3)′cubic′立方插值法，要求x的值等距離。(4)“spline”分段立方樣條插值(SPLINE)注意：所有插值方法均要求x是單調(diào)的（單調(diào)遞增或單調(diào)遞減），且待插值的數(shù)據(jù)點橫坐標不能超出x給定的范圍，否則，插值結(jié)果會出現(xiàn)非數(shù)(NaN).4.3.2.多項式插值【例4-28】在區(qū)間[05]內(nèi)，取正弦函數(shù)曲線上均勻分布的6個數(shù)據(jù)點作為原始數(shù)據(jù)。再在該區(qū)間上選取均勻分布的21個點作為自變量，利用插值方法分別計算用線性插值、三次樣條插值和三次插得到的函數(shù)應(yīng)變量的值，并進行誤差比較。>>x=linspace(0,5,6);y=sin(x);%創(chuàng)建原始數(shù)據(jù)點>>x0=linspace(0,5,21);%待插值點>>y0=sin(x0);%精確解>>y1=interp1(x,y,x0);%線性插值法>>y2=interp1(x,y,x0,'spline');%三次樣條插值法>>y3=interp1(x,y,x0,'cubic');%立方插值法

4.3.2.多項式插值%插值結(jié)果與精確解之間的殘差>>err=[y1-y0;y2-y0;y3-y0];

%不同插值方法帶來的殘差標準方差>>s=[std(err(1,:)),std(err(2,:)),std(err(3,:))]s=0.06410.01020.0522【說明】函數(shù)std()計算序列的標準方差，插值殘差的標準方差可以衡量數(shù)值插值效果：殘差的標準方差越大，表明插值誤差越大，插值效果越不理想。在四種插值方法中，樣條插值和立方插值的效果比較好，而線性插值的效果比較差。4.3.2.多項式插值二維插值函數(shù)是interp2()，調(diào)用格式如下：?Z1=interp2(X,Y,Z,X1,Y1,'method')?根據(jù)數(shù)據(jù)點(X,Y,Z)的值，計算函數(shù)在(X1,Y1)處由相應(yīng)的插值方法得到的插值結(jié)果Z1。其中X、Y是兩個向量，分別描述兩個參數(shù)的采樣點，X1、Y1是兩個向量或標量，描述欲插值的點。method的取值與一維插值函數(shù)相同。X、Y、Z也可以是矩陣形式。X1、Y1的取值范圍不能超出X、Y的給定范圍，否則，插值結(jié)果會出現(xiàn)非數(shù)(NaN)。MATLAB正是利用插值方法實現(xiàn)二維繪圖和三維繪圖。2.二維數(shù)值插值4.4函數(shù)的零點和極值點

4.4.1函數(shù)的零點

對于任意函數(shù)f(x)，它可能存在零點，也可能沒有零點；即便存在零點，函數(shù)也可能只有一個零點，或者有多個甚至無窮多個零點。對于應(yīng)用解析方法尋找零點較困難，或者解析方法尋找零點無能為力的情況下，MATLAB提供了找尋函數(shù)零點的數(shù)值解法。MATLAB中找尋函數(shù)零點的數(shù)值計算方法為：先猜測一個初始零點或者零點所在的區(qū)間，然后通過不斷迭代，使得猜測值不斷精確化，或區(qū)間不斷收縮，直至達到預(yù)先指定的精度，終止計算輸出結(jié)果。MATLAB中找尋的函數(shù)零點的指令為fzero(),調(diào)用格式如下：?[xx,yy]=fzero(fun,x0)?計算函數(shù)fun在x0附近的零點；x0為初始猜測零點：若x0取標量，該指令將在它兩側(cè)尋找一個與之最靠近的零點；x0為區(qū)間時，指令將在區(qū)間內(nèi)尋找零點。輸出宗量xx、yy分別為零點的橫坐標和縱坐標。如果函數(shù)只有一個輸出宗量，則僅給出零點橫坐標。【說明】該函數(shù)有不止一個零點，但函數(shù)fzero()只能找出一個零點。為了能夠有目的的找尋多個零點中的一個，用戶可以先借助圖形找到初始猜測零點的值。1.一元函數(shù)的零點（1）使用內(nèi)聯(lián)函數(shù)，并計算0.5附近的零點.>>clear,f=inline('(sin(t))^2*exp(-0.1*t)-0.5*abs(t)')>>[t,y]=fzero(f,0.5)%計算0.5附近的零點t=0.5993y=-5.5511e-017為了獲得更多的零點信息，建議用戶用作圖法觀察函數(shù)零點的分布情況，并借助相關(guān)指令來獲取多個零點的近似值?！纠?-29】求的零點。（2）作圖法觀察零點分布。編寫文件名為exm4_29的腳本文件：t=linspace(-8,8,200);y=(sin(t)).^2.*exp(-0.1*t)-0.5*abs(t);plot(t,y,t,zeros(size(t)))%繪圖指令xlabel('t'),ylabel('y(t)')%添加坐標軸名稱在指令窗中輸入exm4_29并執(zhí)行，結(jié)果如圖4-2所示。

圖4.2函數(shù)零點分布觀察圖>>zoomon%執(zhí)行該指令后，可獲得圖形的局部放大圖>>[tt,yy]=ginput(5);zoomoff%利用鼠標獲取5個零點猜測值>>tt%送出零點猜測初始值tt=-2.0046-0.5392-0.02300.58531.6820圖4.3局部放大和利用鼠標取值（3）利用zoom和ginput指令獲取零點的近似值。>>[t4,y4]=fzero(f,tt(4))t4=0.5993y4=0【說明】用戶可以試試將內(nèi)聯(lián)函數(shù)用字符串代替求解函數(shù)零點，此時的輸入宗量必須是x。當用戶執(zhí)行指令ginput()時，鼠標光標將變成十字形，如圖4.3所示。用戶單擊鼠標即可獲得圖中相應(yīng)點的坐標，并將其送給工作空間中的變量。（4）求靠近tt(4)的精確零點。MATLAB中函數(shù)fsolve()可以用來計算多元函數(shù)的零點。調(diào)用格式為：?[x,fval]=fsolve(fun,x0)?計算函數(shù)fun在x0附近的零點。輸出宗量x為零點橫坐標坐標，fval為零點縱坐標。2.多元函數(shù)的零點>>x0=[2,2,2];>>[x,fval]=fsolve('2*x(1)+3*x(2)^2+x(1)*x(3)',x0)Warning:……x=-0.0040-0.06501.1829fval=4.5057e-011或者執(zhí)行語句“fsolve('2*x(1)+3*x(2)^2+x(1)*x(3)',[222])”即可。【例4-30】計算方程在(2,2,2)附近的零點4.4.2函數(shù)的極值點

數(shù)學上，函數(shù)的極值點是通過確定函數(shù)導(dǎo)數(shù)為零的點，解析求出這些極值點。有的函數(shù)很難找到導(dǎo)數(shù)為零的點，因此利用解析法求極值點難度很大，有時甚至是不可行的。MATLAB提供了在數(shù)值上找尋函數(shù)極值點的指令。對于函數(shù)的極值點存在極大值和極小值的區(qū)別，MATLAB僅僅提供了獲取函數(shù)極小值的指令，并且這種“極小”是一個“局部極小”，即是在給定范圍內(nèi)的“極小”。計算一個函數(shù)的極大值等價于計算該函數(shù)相反數(shù)的極小值。函數(shù)fminbnd()可以用來求解一定區(qū)間內(nèi)函數(shù)的極小值，調(diào)用格式：?[x,y]=fminbnd(fun,x1,x2)?計算區(qū)間(x1,x2)內(nèi)函數(shù)fun的極小值。輸出宗量x、y分別為極小值點的橫坐標和縱坐標。如果只有一個輸出宗量即為橫坐標?！纠?-31】計算函數(shù)在(0,3)附近的極小值點和極大值點。編寫文件名為exm4_31的腳本文件：clear,f=inline('2*exp(-x)*cos(x)');[x,y_min]=fminbnd(f,0,3)%輸出區(qū)間（0,3）內(nèi)函數(shù)的極小值坐標1，一元函數(shù)的極值點%計算區(qū)間（0,3）的極大值[x,y_min]=fminbnd('-2*exp(-x)*cos(x)',0,3);y_max=-y_min在指令窗中運行exm4_31，有：x=2.3562y_min=-0.1340y_max=1.9999【例4-31】計算函數(shù)在(0,3)附近的極小值點和極大值點。求多元函數(shù)極小值點常見的兩種方法：單純行下山法(Downhillsimplexmethods)和擬牛頓法(Qussi-Newtonmethods)。對應(yīng)的MATLAB指令為：x=fminsearch(fun,x0)[x,fval]=fminsearch(fun,x0)

單純行下山法求多元函數(shù)極小值最簡格式x=fminunc(fun,x0)[x,fval]=fminunc(fun,x0)擬牛頓法求多元函數(shù)極小值最簡格式其中x為極小值點的坐標。fval輸出對應(yīng)于極小值點x的極小值。其中x0為找尋極小值的起點，x0可以是標量、向量或者矩陣。2.多元函數(shù)的極值點>>clear,clc,>>x0=[-1.21.2];%猜測零點位置>>ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');>>[x,fval]=fminsearch(ff,x0)%單純行法求極小值點x=1.00001.0000fval=5.4009e-010>>[x,fval]=fminunc(ff,x0)%擬牛頓法求極小值點Warning:……x=1.00001.0000fval=

2.2044e-011【例4-22】計算“Banana”測試函數(shù)的極小值點。該測試函數(shù)有一片淺谷，許多算法難以逾越此谷。它的理論極小值點是(1,1)。4.5數(shù)值微積分當已知函數(shù)的表達式時，理論上可以通過公式進行微積分計算。但在實際應(yīng)用中，往往需要處理的對象并不能用公式來進行計算。因此，有必要介紹這些函數(shù)的微分和積分的數(shù)值算法。4.5.1差分和偏導(dǎo)函數(shù)f(x)在x點的偏導(dǎo)數(shù)定義為：設(shè)h>0，△f(x)=f(x+h)-f(x)稱為f(x)在x點的向前差分。在MATLAB中，計算差分的指令為diff()，調(diào)用格式為：?DX=diff(A,n,dim)?計算矩陣A的n階差分。dim=1時，按列計算差分；dim=2，按行計算差分。dim缺省時按列計算；n缺省時為計算1階差分。指令gradient()也可以計算數(shù)值偏導(dǎo)，格式如下：?[dfdx,dfdy]=gradient(f,dx,dy)?計算矩陣f分別對x和y方向的數(shù)值導(dǎo)數(shù)，即計算矩陣在x、y方向的梯度。clear,clc,dt=pi/20;t=0:dt:pi/2;x=exp(t).*cos(t);%創(chuàng)建數(shù)據(jù)dxdt_diff=diff(x)/dtdxdt_grad=gradient(x,dt)dxdt_diff=Columns1through90.99110.93210.79750.56720.2191-0.2701-0.9243-1.7666-2.8178Column10-4.0943dxdt_grad=Columns1through90.99110.96160.86480.68240.3931-0.0255-0.5972-1.3455-2.2922Columns10through11

-3.4561-4.0943【例4-23】已知，應(yīng)用指令diff()和gradient()計算該函數(shù)在區(qū)間的導(dǎo)函數(shù)。注意：函數(shù)diff()輸出結(jié)果比原數(shù)組長度少一.4.5.2數(shù)值積分

在數(shù)值計算中，數(shù)值積分的基本思想都是將整個積分區(qū)間[a,b]分成n個子區(qū)間,i=1,2,…,n,其中，。積分在各個子區(qū)間上進行，最后求和得到總的積分結(jié)果。數(shù)值積分有定積分和不定積分的差別，MATLAB中進行數(shù)值計算時調(diào)用的執(zhí)行函數(shù)是一致的，只是在輸入宗量上有差別：對于定積分，必須給出積分區(qū)間。當沒有積分區(qū)間時，則計算不定積分。常用的一重積分函數(shù)有兩個，即quad()和quadl().quad函數(shù)是基于變步長辛普生法的求定積分函數(shù)，調(diào)用格式如下：?[I,n]=quad(fun,a,b,tol)?計算被積函數(shù)fun在區(qū)間[a,b]上的定積分。輸出宗量I為定積分結(jié)果，n為被積函數(shù)調(diào)用次數(shù)。如果只有一個輸出宗量，即為定積分結(jié)果。tol缺省時為。quadl是基于洛巴托法求定積分函數(shù)，調(diào)用格式如下;?[I,n]=quadl(fun,a,b,tol)?功能同函數(shù)quad()。該函數(shù)可以更精確地求出定積分的值，且一般情況下函數(shù)調(diào)用的步數(shù)明顯小于quad函數(shù)，從而保證能以更高的效率求出所需的定積分值。tol缺省時為。以表格形式定義的函數(shù)，MATLAB計算定積分的函數(shù)為trapz(),格式為：?St=trapz(x,y)?采用梯形法求函數(shù)y關(guān)于自變量x的積分。1.一重積分編寫文件名為exm4_34的腳本文件：clear,clc,formatlong[s_quad,n1]=quad('sqrt(1-1/2*sin(x))',0,pi/2,1e-11)%設(shè)置絕對精度為10^（-11）[s_quadl,n2]=quadl('sqrt(1-1/2*sin(x))',0,pi/2,1e-11)x=linspace(0,pi/2,500);y=sqrt(1-1/2*sin(x));%構(gòu)建表格數(shù)據(jù)s_trapz=trapz(x,y)%采用trape函數(shù)計算積分【例4-34】求,其精確值1.288898992310102….在指令窗中執(zhí)行exm4_34，結(jié)果為：s_quad=1.288898992309974n1=121s_quadl=

1.288898992310102n2=48s_trapz=1.288899198751829【說明】采用trapz()函數(shù)計算積分時，用戶不能控制計算結(jié)果的精度quad()和quadl()可以事先指定數(shù)值積分的精度，默認精度為。函數(shù)quadl()計算數(shù)值積分時調(diào)用被積函數(shù)的次數(shù)遠小于quad()【例4-34】求,其精確值1.288898992310102….?S=dblquad(fun,xmin,xmax,ymin,ymax,tol)?S=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)?采用遞推自適應(yīng)辛普森法計算二重dblquad()和三重triplequad()定積分。tol缺省時，積分的絕對精度為【例4-35】求的結(jié)果。2.二重積分(1)編寫文件名為fin的函數(shù)文件。

functionf=fin(x,y)%FINFunctionfortheintegratedcalculation%x,y被積函數(shù)自變量%f被積函數(shù)%K統(tǒng)計被積函數(shù)被調(diào)用次數(shù)globalK;%定義K為全局變量K=K+1;f=exp(-x.^2/3).*sin(x.^2+y.^2);end