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

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

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

1.矩陣轉(zhuǎn)置的MATLAB運(yùn)算符為單引號(hào)“’”。?A'?求A的轉(zhuǎn)置，其中A可以是行向量、列向量和矩陣?！纠?-1】矩陣轉(zhuǎn)置的MATLAB實(shí)例。>>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.‘%的點(diǎn)轉(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)變換對(duì)于復(fù)數(shù)矩陣，MATLAB完成的是復(fù)數(shù)矩陣的共軛轉(zhuǎn)置。復(fù)數(shù)矩陣的非共軛轉(zhuǎn)置的運(yùn)算符為“.‘”。注意：矩陣轉(zhuǎn)置的操作符必須在英文狀態(tài)下輸入。4.1.1矩陣的結(jié)構(gòu)變換指令flipud()和fliplr()完成矩陣的對(duì)稱變換，調(diào)用格式如下：?B=flipud(A)?上下方向翻轉(zhuǎn)的矩陣。如果是列向量，返回相反順序的向量；如果是行向量，返回原向量。?B=fliplr(A)?水平方向翻轉(zhuǎn)的矩陣。如果是行向量，返回相反順序的向量；如果是列向量，返回原向量。另外，MATLAB還提供了一個(gè)以指定維翻轉(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.對(duì)稱變換>>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)變換實(shí)例。指令rot90()可以完成矩陣逆時(shí)針旋轉(zhuǎn)90度。

?rot90(A,k)

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

?rank(A)

?計(jì)算矩陣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還提供了一個(gè)化簡(jiǎn)函數(shù)rref()，其格式如下：

?rref(A)

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

根據(jù)MATLAB的輸出結(jié)果，用戶可以得到矩陣的秩，即行階梯陣中非零行的行數(shù)【例4-8】（續(xù)例4-7）行階梯陣化簡(jiǎn)函數(shù)rref（）應(yīng)用實(shí)例。>>rref(a)%將矩陣a化簡(jiǎn)為行階梯陣ans=100010001>>rref(b)%將矩陣b化簡(jiǎn)為行階梯陣ans=10-2012000由化簡(jiǎn)結(jié)果可知，矩陣a的秩為3，矩陣b的秩為2.。1秩把方陣A看做行列式|A|，對(duì)其按照行列式的規(guī)則求值，稱該值為行列式的值。MATLAB的實(shí)現(xiàn)指令為det()，其格式為：?det(A)?計(jì)算方陣A對(duì)應(yīng)的行列式的值行列式的值為0時(shí)，相應(yīng)的矩陣成為奇異矩陣，否則稱為非奇異矩陣（或滿秩矩陣）。對(duì)于非奇異矩陣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)%計(jì)算矩陣對(duì)應(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）矩陣對(duì)應(yīng)行列式的值、逆矩陣和廣義逆矩陣運(yùn)算實(shí)例。對(duì)于奇異矩陣，調(diào)用MATLAB指令inv()也可以獲得計(jì)算結(jié)果.此時(shí)MATLAB會(huì)給出警告：Warning:Matrixisclosetosingularorbadlyscaled.用戶不要將輸出矩陣認(rèn)為是該奇異矩陣的逆矩陣，必須通過逆矩陣定義加以驗(yàn)證?！纠?-10】逆矩陣的應(yīng)用：設(shè)A=,B=,C=,求矩陣X，使?jié)M足:?！菊f明】由線性代數(shù)可知：，或者。本題給出兩種計(jì)算程序，并比較計(jì)算結(jié)果和耗時(shí)。編寫文件名為exm4_10的腳本文件：clear,clc,A=[123;234;343];B=[2,5;43];C=[13;40;21];%顯示利用逆矩陣求解時(shí)的耗時(shí)tic%開啟計(jì)時(shí)器X1=inv(A)*C*inv(B)toc%結(jié)束計(jì)時(shí)，并輸出時(shí)間%顯示利用矩陣除法求解時(shí)的耗時(shí)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)用利用逆矩陣求解方程組時(shí)需要兩步計(jì)算，即求逆矩陣和矩陣乘法，因此較左除法求解方程組速度要慢。程序在不同計(jì)算機(jī)上執(zhí)行時(shí)耗時(shí)不唯一，但時(shí)長(zhǎng)順序不會(huì)改變。如果矩陣A不是方陣或是奇異矩陣時(shí)，該矩陣存在與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ù)解。矩陣的跡：矩陣的對(duì)角線元素之和，即矩陣的特征值之和。MATLAB中求跡函數(shù)是trace()，格式如下：?trace(A)?求矩陣A的跡。向量種度量形式。使用范數(shù)可以測(cè)量?jī)蓚€(gè)向量或矩陣之間的距離。范數(shù)有多種定義形式。MATLAB計(jì)算向量或矩陣的范數(shù)指令為norm()，格式如下：?norm(A,k)?計(jì)算向量或矩陣A的k—范數(shù):k可以取1、2、Inf、'fro',分別為計(jì)算矩陣的1-范數(shù)、2-范數(shù)、無窮大-范數(shù)、Frobenius-范數(shù),k缺省時(shí)為計(jì)算2-范數(shù)。3.矩陣的跡和范數(shù)>>clear,A=[86-11;43-12;53-34;34-7-2];>>trace(A)%計(jì)算矩陣的跡ans=6>>norm(A,1)%計(jì)算矩陣的1-范數(shù)ans=20>>norm(A)%計(jì)算矩陣的2-范數(shù)ans=14.6918【例4-12】計(jì)算矩陣的跡和范數(shù)。>>norm(A,inf)%計(jì)算矩陣的無窮大-范數(shù)ans=16>>norm(A,'fro')%計(jì)算矩陣的Frobenius-范數(shù)ans=16.4012>>norm(A(:,2),1)%計(jì)算向量的1-范數(shù)ans=16【例4-12】計(jì)算矩陣的跡和范數(shù)。矩陣的條件數(shù)是解線性方程組是判斷系數(shù)矩陣的變化對(duì)解的影響程度的一個(gè)參數(shù)。矩陣的條件數(shù)的定義依賴于范數(shù)。矩陣的條件數(shù)接近1時(shí)表示方程為良性的，否則為病態(tài)。MATLAB計(jì)算矩陣條件數(shù)的指令為cond(),格式如下：?cond(A,p)?計(jì)算矩陣A的p—范數(shù):p可以取1、2、Inf、'fro',分別為計(jì)算矩陣的1-范數(shù)條件數(shù)、2-范數(shù)條件數(shù)、無窮大-范數(shù)條件數(shù)、Frobenius-范數(shù)條件數(shù),p缺省時(shí)為計(jì)算2-范數(shù)條件數(shù)。4.條件數(shù)>>cond(A,1)%計(jì)算矩陣的1-范數(shù)條件數(shù)ans=108.1720>>cond(A)%計(jì)算矩陣的2-范數(shù)條件數(shù)ans=51.7487>>cond(A,inf)%計(jì)算矩陣的無窮大-范數(shù)條件數(shù)ans=67.0968>>cond(A,'fro')%計(jì)算矩陣的Frobenius-范數(shù)條件數(shù)ans=58.0285【例4-13】（續(xù)例4-12）計(jì)算矩陣A的條件數(shù)。4.1.3矩陣的特征值分析設(shè)n階方陣A，如果數(shù)和n維非零列向量，使關(guān)系式成立，則數(shù)稱為方陣A的特征值，非零向量稱為矩陣A的對(duì)應(yīng)于特征值的特征向量，MATLAB計(jì)算矩陣特征值和特征向量的函數(shù)為eig()，格式為：?[V,D]=eig(A)?計(jì)算矩陣A的特征值D和對(duì)應(yīng)的特征向量V，使得AV=VD。如果函數(shù)只有一個(gè)輸出宗量，則只給出特征值。>>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】計(jì)算3階魔方矩陣的特征值和特征向量。4.1.4矩陣的分解矩陣分解是線性代數(shù)中比較重要的內(nèi)容，針對(duì)不同的分解方法，MATLAB提供了多個(gè)關(guān)于矩陣分解的函數(shù)，具體見表4-1所示。分解類別調(diào)用格式功能LU分解[L,U]=lu(A)計(jì)算上三角矩陣U和交換下三角矩陣L，使LU=A[L,U,P]=lu(A)計(jì)算上三角矩陣U、有單位對(duì)角線的下三角矩陣L和交換矩陣P，滿足LU=PACholesky因式分解R=chol(A)計(jì)算正定矩陣A的Cholesky因子，是一個(gè)上三角矩陣，滿足R'*R=A。如果A不是一個(gè)正定矩陣，則給出一個(gè)錯(cuò)誤信息。[G,p]=chol(A)計(jì)算矩陣A的Cholesky因子G。如果A不是一個(gè)正定矩陣，則不給出錯(cuò)誤信息，而是將p設(shè)為正整數(shù)。qr因式分解[Q,R]=qr(A)計(jì)算m×n矩陣A分解得到的m×n酉矩陣Q和m×n的上三角矩陣R。Q的列形成了一個(gè)正交基。Q和R滿足A=Q*R[Q,R,P]=qr(A)計(jì)算矩陣Q、上三角矩陣R和交換矩陣P。Q的列形成一個(gè)正交基，R的對(duì)角線元素按大小降序排列，滿足A*P=Q*RSVD分解s=svd(A)計(jì)算一個(gè)包含矩陣A的奇異值的向量s[U,S,V]=svd(A)計(jì)算奇異值為對(duì)角線且與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分解法求欠定方程組的一個(gè)特解。4.1.5線性方程組的求解線性方程組分為齊次方程組和非齊次方程組兩類。對(duì)于系數(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】解方程組齊次線性方程組求解對(duì)于系數(shù)陣A為奇異陣的齊次線性方程組，可以利用函數(shù)rref(A)將系數(shù)陣化為行階梯陣，對(duì)照化簡(jiǎn)結(jié)果可以直接寫出方程組的解?！纠?-17】解方程組對(duì)于系數(shù)矩陣為非滿秩矩陣，即奇異矩陣，MATLAB提供了針對(duì)齊次線性方程組和非齊次線性方程組的不同求解函數(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ù)：當(dāng)有r時(shí)，輸出有理基；當(dāng)無r時(shí)，輸出正交規(guī)范基?！纠?-18】（續(xù)例4-17）利用函數(shù)null（）求解齊次線性方程組示例。>>x=null(A,'r')x=4/36-23/31則方程組的解為。結(jié)果與【例4-14】一致。【例4-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)的秩且等于未知變量的個(gè)數(shù)m（3）無窮多解：系數(shù)陣A的秩n等于增廣陣(A,b)的秩且小于未知變量的個(gè)數(shù)m，則方程組有無窮多解。對(duì)于情況（3），MATLAB必須要給出該方程的一個(gè)特解和對(duì)應(yīng)齊次方程的(n-m)個(gè)通解

求解通解用函數(shù)null()，特解可以利用矩陣LU分解法求得。2.非齊次線性方程組求解functionX=equs(A,b,n)%EQUSSolvingthehomogeneousorinhomogeneouslinearequations%A方程組的系數(shù)矩陣%b常數(shù)項(xiàng)對(duì)應(yīng)的列向量%n方程組中變量個(gè)數(shù)%X線性方程組的解ifnargin>3error('輸入宗量太多，請(qǐng)重新輸入')endB=[A,b];%增廣矩陣R_A=rank(A);%求系數(shù)矩陣的秩R_B=rank(B);%求增廣矩陣的秩ifR_A>n|R_B>nerror('輸入變量個(gè)數(shù)有誤，請(qǐng)重新輸入')end編寫一個(gè)文件名為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ù)文件，可以用來求解任意類型的齊次或者非齊次線性方程組。編寫一個(gè)文件名為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多項(xiàng)式多項(xiàng)式運(yùn)算是數(shù)學(xué)運(yùn)算中最基本運(yùn)算方式之一，也是線性代數(shù)分析和設(shè)計(jì)中的重要內(nèi)容。MATLAB提供了多個(gè)函數(shù)可以幫助用戶完成多項(xiàng)式的運(yùn)算。4.2.1多項(xiàng)式的表達(dá)和創(chuàng)建MATLAB中將n階多項(xiàng)式p(x)存儲(chǔ)在長(zhǎng)度為n+1的行向量p中，行向量p的元素為多項(xiàng)式的系數(shù)，并按x的降冪排列。行向量代表的多項(xiàng)式為：注意：多項(xiàng)式中缺少某個(gè)冪次項(xiàng)，則在創(chuàng)建行向量時(shí)將該冪次項(xiàng)的系數(shù)寫為0.【例4-21】創(chuàng)建多項(xiàng)式>>p=[4,-3,0,0,1,0]%三次項(xiàng)、二次項(xiàng)、常數(shù)項(xiàng)系數(shù)均為0p=4-300104.2.2多項(xiàng)式的運(yùn)算

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

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

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

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

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

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

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

4.4.1函數(shù)的零點(diǎn)

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

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

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

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

2.2044e-011【例4-22】計(jì)算“Banana”測(cè)試函數(shù)的極小值點(diǎn)。該測(cè)試函數(shù)有一片淺谷，許多算法難以逾越此谷。它的理論極小值點(diǎn)是(1,1)。4.5數(shù)值微積分當(dāng)已知函數(shù)的表達(dá)式時(shí)，理論上可以通過公式進(jìn)行微積分計(jì)算。但在實(shí)際應(yīng)用中，往往需要處理的對(duì)象并不能用公式來進(jìn)行計(jì)算。因此，有必要介紹這些函數(shù)的微分和積分的數(shù)值算法。4.5.1差分和偏導(dǎo)函數(shù)f(x)在x點(diǎn)的偏導(dǎo)數(shù)定義為：設(shè)h>0，△f(x)=f(x+h)-f(x)稱為f(x)在x點(diǎn)的向前差分。在MATLAB中，計(jì)算差分的指令為diff()，調(diào)用格式為：?DX=diff(A,n,dim)?計(jì)算矩陣A的n階差分。dim=1時(shí)，按列計(jì)算差分；dim=2，按行計(jì)算差分。dim缺省時(shí)按列計(jì)算；n缺省時(shí)為計(jì)算1階差分。指令gradient()也可以計(jì)算數(shù)值偏導(dǎo)，格式如下：?[dfdx,dfdy]=gradient(f,dx,dy)?計(jì)算矩陣f分別對(duì)x和y方向的數(shù)值導(dǎo)數(shù)，即計(jì)算矩陣在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()計(jì)算該函數(shù)在區(qū)間的導(dǎo)函數(shù)。注意：函數(shù)diff()輸出結(jié)果比原數(shù)組長(zhǎng)度少一.4.5.2數(shù)值積分

在數(shù)值計(jì)算中，數(shù)值積分的基本思想都是將整個(gè)積分區(qū)間[a,b]分成n個(gè)子區(qū)間,i=1,2,…,n,其中，。積分在各個(gè)子區(qū)間上進(jìn)行，最后求和得到總的積分結(jié)果。數(shù)值積分有定積分和不定積分的差別，MATLAB中進(jìn)行數(shù)值計(jì)算時(shí)調(diào)用的執(zhí)行函數(shù)是一致的，只是在輸入宗量上有差別：對(duì)于定積分，必須給出積分區(qū)間。當(dāng)沒有積分區(qū)間時(shí)，則計(jì)算不定積分。常用的一重積分函數(shù)有兩個(gè)，即quad()和quadl().quad函數(shù)是基于變步長(zhǎng)辛普生法的求定積分函數(shù)，調(diào)用格式如下：?[I,n]=quad(fun,a,b,tol)?計(jì)算被積函數(shù)fun在區(qū)間[a,b]上的定積分。輸出宗量I為定積分結(jié)果，n為被積函數(shù)調(diào)用次數(shù)。如果只有一個(gè)輸出宗量，即為定積分結(jié)果。tol缺省時(shí)為。quadl是基于洛巴托法求定積分函數(shù)，調(diào)用格式如下;?[I,n]=quadl(fun,a,b,tol)?功能同函數(shù)quad()。該函數(shù)可以更精確地求出定積分的值，且一般情況下函數(shù)調(diào)用的步數(shù)明顯小于quad函數(shù)，從而保證能以更高的效率求出所需的定積分值。tol缺省時(shí)為。以表格形式定義的函數(shù)，MATLAB計(jì)算定積分的函數(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è)置絕對(duì)精度為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ù)計(jì)算積分【例4-34】求,其精確值1.288898992310102….在指令窗中執(zhí)行exm4_34，結(jié)果為：s_quad=1.288898992309974n1=121s_quadl=

1.288898992310102n2=48s_trapz=1.288899198751829【說明】采用trapz()函數(shù)計(jì)算積分時(shí)，用戶不能控制計(jì)算結(jié)果的精度quad()和quadl()可以事先指定數(shù)值積分的精度，默認(rèn)精度為。函數(shù)quadl()計(jì)算數(shù)值積分時(shí)調(diào)用被積函數(shù)的次數(shù)遠(yuǎn)小于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)辛普森法計(jì)算二重dblquad()和三重triplequad()定積分。tol缺省時(shí)，積分的絕對(duì)精度為【例4-35】求的結(jié)果。2.二重積分(1)編寫文件名為fin的函數(shù)文件。

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