1、第第6章章 MATLAB數值計算數值計算6.1 數據處理與多項式計算數據處理與多項式計算6.2 數值微積分數值微積分6.3 離散傅立葉變換離散傅立葉變換6.4 線性方程組求解線性方程組求解6.5 非線性方程與最優化問題求解非線性方程與最優化問題求解6.6 常微分方程的數值求解常微分方程的數值求解6.7 稀疏矩陣稀疏矩陣6.1 數據處理與多項式計算數據處理與多項式計算1. 求矩陣最大元素和最小元素求矩陣最大元素和最小元素MATLAB提供的求數據序列的最大值和最小提供的求數據序列的最大值和最小值的函數分別為值的函數分別為max和和min，兩個函數的調，兩個函數的調用格式和操作過程類似。用格式和操作

2、過程類似。（1）求向量的最大值和最小值）求向量的最大值和最小值 y=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，如，如果果X中包含復數元素，則按模取最大值。中包含復數元素，則按模取最大值。 y,I=max(X)：返回向量：返回向量X的最大值存入的最大值存入y，最大，最大值的序號存入值的序號存入I，如果，如果X中包含復數元素，則按模中包含復數元素，則按模取最大值。取最大值。求向量求向量X的最小值的函數是的最小值的函數是min(X)，用法和，用法和max(X)完全相同。完全相同。例例 求向量求向量x的最大值。的最大值。命令如下：命令如下：x=-43,72,9,16,23,47;y

3、=max(x) %求向量求向量x中的最大值中的最大值y,l=max(x) %求向量求向量x中的最大值及其該元素的位置中的最大值及其該元素的位置（2）求矩陣的最大值和最小值）求矩陣的最大值和最小值求矩陣求矩陣A的最大值的函數有的最大值的函數有3種調用格式，分種調用格式，分別是：別是： max(A)：返回一個行向量，向量的第：返回一個行向量，向量的第i個元個元素是矩陣素是矩陣A的第的第i列上的最大值。列上的最大值。 Y,U=max(A)：返回行向量：返回行向量Y和和U，Y向量向量記錄記錄A的每列的最大值，的每列的最大值，U向量記錄每列最向量記錄每列最大值的行號。大值的行號。 max(A,dim)：

4、dim取取1或或2。dim取取1時，時，該函數和該函數和max(A)完全相同；完全相同；dim取取2時，該時，該函數返回一個列向量，其第函數返回一個列向量，其第i個元素是個元素是A矩矩陣的第陣的第i行上的最大值。行上的最大值。 求最小值的函數是求最小值的函數是min，其用法和，其用法和max完全完全相同。相同。例例6.1 分別矩陣分別矩陣A中各列和各行元素中的最大中各列和各行元素中的最大值，并求整個矩陣的最大值和最小值。值，并求整個矩陣的最大值和最小值。（3）兩個向量或矩陣對應元素的比較）兩個向量或矩陣對應元素的比較 函數函數max和和min還能對兩個同型的向量或矩陣進還能對兩個同型的向量或矩

5、陣進行比較，調用格式為：行比較，調用格式為： U=max(A,B)：A,B是兩個同型的向量或矩陣，結是兩個同型的向量或矩陣，結果果U是與是與A,B同型的向量或矩陣，同型的向量或矩陣，U的每個元素等的每個元素等于于A,B對應元素的較大者。對應元素的較大者。 U=max(A,n)：n是一個標量，結果是一個標量，結果U是與是與A同型同型的向量或矩陣，的向量或矩陣，U的每個元素等于的每個元素等于A對應元素和對應元素和n中的較大者。中的較大者。min函數的用法和函數的用法和max完全相同。完全相同。 例例 求兩個求兩個23矩陣矩陣x, y所有同一位置上的較大元所有同一位置上的較大元素構成的新矩陣素構成的

6、新矩陣p。2. 求矩陣的平均值和中值求矩陣的平均值和中值求數據序列平均值的函數是求數據序列平均值的函數是mean，求數據序列中值的函數，求數據序列中值的函數是是median。兩個函數的調用格式為：。兩個函數的調用格式為：mean(X)：返回向量：返回向量X的算術平均值。的算術平均值。median(X)：返回向量：返回向量X的中值。的中值。mean(A)：返回一個行向量，其第：返回一個行向量，其第i個元素是個元素是A的第的第i列的算術列的算術平均值。平均值。median(A)：返回一個行向量，其第：返回一個行向量，其第i個元素是個元素是A的第的第i列的中列的中值。值。mean(A,dim)：當：

7、當dim為為1時，該函數等同于時，該函數等同于mean(A)；當；當dim為為2時，返回一個列向量，其第時，返回一個列向量，其第i個元素是個元素是A的第的第i行的算術行的算術平均值。平均值。median(A,dim)：當：當dim為為1時，該函數等同于時，該函數等同于median(A)；當；當dim為為2時，返回一個列向量，其第時，返回一個列向量，其第i個元素是個元素是A的第的第i行的行的中值。中值。3. 矩陣元素求和與求積矩陣元素求和與求積數據序列求和與求積的函數是數據序列求和與求積的函數是sum和和prod，其使用方法類似。設其使用方法類似。設X是一個向量，是一個向量，A是一是一個矩陣，函

8、數的調用格式為：個矩陣，函數的調用格式為：sum(X)：返回向量：返回向量X各元素的和。各元素的和。prod(X)：返回向量：返回向量X各元素的乘積。各元素的乘積。sum(A)：返回一個行向量，其第：返回一個行向量，其第i個元素是個元素是A的第的第i列的元素和。列的元素和。prod(A)：返回一個行向量，其第：返回一個行向量，其第i個元素是個元素是A的第的第i列的元素乘積。列的元素乘積。sum(A,dim)：當：當dim為為1時，該函數等同于時，該函數等同于sum(A)；當；當dim為為2時，返回一個列向量，時，返回一個列向量，其第其第i個元素是個元素是A的第的第i行的各元素之和。行的各元素之

9、和。prod(A,dim)：當：當dim為為1時，該函數等同于時，該函數等同于prod(A)；當；當dim為為2時，返回一個列向量，時，返回一個列向量，其第其第i個元素是個元素是A的第的第i行的各元素乘積。行的各元素乘積。例例6.2 求矩陣求矩陣A的每行元素的乘積和全部元素的每行元素的乘積和全部元素的乘積。的乘積。4. 矩陣元素累加和與累乘積矩陣元素累加和與累乘積在在MATLAB中，使用中，使用cumsum和和cumprod函數能方便地求得函數能方便地求得向量和矩陣元素的累加和與累乘積向量，函數的調用格式向量和矩陣元素的累加和與累乘積向量，函數的調用格式為：為：cumsum(X)：返回向量：返

10、回向量X累加和向量。累加和向量。cumprod(X)：返回向量：返回向量X累乘積向量。累乘積向量。cumsum(A)：返回一個矩陣，其第：返回一個矩陣，其第i列是列是A的第的第i列的累加和向列的累加和向量。量。cumprod(A)：返回一個矩陣，其第：返回一個矩陣，其第i列是列是A的第的第i列的累乘積列的累乘積向量。向量。cumsum(A,dim)：當：當dim為為1時，該函數等同于時，該函數等同于cumsum(A)；當當dim為為2時，返回一個矩陣，其第時，返回一個矩陣，其第i行是行是A的第的第i行的累加行的累加和向量。和向量。cumprod(A,dim)：當：當dim為為1時，該函數等同于

11、時，該函數等同于cumprod(A)；當當dim為為2時，返回一個向量，其第時，返回一個向量，其第i行是行是A的第的第i行的累乘行的累乘積向量。積向量。5求標準方差求標準方差在在MATLAB中，提供了計算數據序列的標準方差的中，提供了計算數據序列的標準方差的函數函數std。對于向量。對于向量X，std(X)返回一個標準方差。返回一個標準方差。對于矩陣對于矩陣A，std(A)返回一個行向量，它的各個元返回一個行向量，它的各個元素便是矩陣素便是矩陣A各列或各行的標準方差。各列或各行的標準方差。std函數的函數的一般調用格式為：一般調用格式為：Y=std(A,flag,dim)其中其中dim取取1或

12、或2。當。當dim=1時，求各列元素的標準方時，求各列元素的標準方差；當差；當dim=2時，則求各行元素的標準方差。時，則求各行元素的標準方差。flag取取0或或1，當，當flag=0時，按時，按1所列公式計算標準方所列公式計算標準方差，當差，當flag=1時，按時，按2所列公式計算標準方差。所列公式計算標準方差。缺省缺省flag=0，dim=1。例例6.4 對二維矩陣對二維矩陣x，從不同維方向求出其標準方差。，從不同維方向求出其標準方差。6相關系數相關系數MATLAB提供了提供了corrcoef函數，可以求出數函數，可以求出數據的相關系數矩陣。據的相關系數矩陣。corrcoef函數的調用格函

13、數的調用格式為：式為：corrcoef(X)：返回從矩陣：返回從矩陣X形成的一個相關系形成的一個相關系數矩陣。此相關系數矩陣的大小與矩陣數矩陣。此相關系數矩陣的大小與矩陣X一一樣。它把矩陣樣。它把矩陣X的每列作為一個變量，然后的每列作為一個變量，然后求它們的相關系數。求它們的相關系數。corrcoef(X,Y)：在這里，：在這里，X,Y是向量，它們與是向量，它們與corrcoef(X,Y)的作用一樣。的作用一樣。例例6.5 生成滿足正態分布的生成滿足正態分布的100005隨機矩隨機矩陣，然后求各列元素的均值和標準方差，陣，然后求各列元素的均值和標準方差，再求這再求這5列隨機數據的相關系數矩陣。

14、列隨機數據的相關系數矩陣。命令如下：命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)7. 排序排序MATLAB中對向量中對向量X是排序函數是是排序函數是sort(X)，函數返，函數返回一個對回一個對X中的元素按升序排列的新向量。中的元素按升序排列的新向量。sort函數也可以對矩陣函數也可以對矩陣A的各列或各行重新排序，其的各列或各行重新排序，其調用格式為：調用格式為：Y,I=sort(A,dim)其中其中dim指明對指明對A的列還是行進行排序。若的列還是行進行排序。若dim=1，則按列排；若則按列排；若dim=2，則按行排。，則按行排。Y

15、是排序后的矩是排序后的矩陣，而陣，而I記錄記錄Y中的元素在中的元素在A中位置。中位置。1. 一維數據插值一維數據插值在在MATLAB中，實現這些插值的函數是中，實現這些插值的函數是interp1，其調用格，其調用格式為：式為：Y1=interp1(X,Y,X1,method)函數根據函數根據X,Y的值，計算函數在的值，計算函數在X1處的值。處的值。X,Y是兩個等長是兩個等長的已知向量，分別描述采樣點和樣本值，的已知向量，分別描述采樣點和樣本值，X1是一個向量是一個向量或標量，描述欲插值的點，或標量，描述欲插值的點，Y1是一個與是一個與X1等長的插值結等長的插值結果。果。method是插值方法，

16、允許的取值有是插值方法，允許的取值有linear、nearest、cubic、spline。注意：注意：X1的取值范圍不能超出的取值范圍不能超出X的給定范圍，否則，的給定范圍，否則，會給出會給出“NaN”錯誤。錯誤。 例例6.7 給出概率積分的數據表如表給出概率積分的數據表如表6.1所示，用不同的插值所示，用不同的插值方法計算方法計算f(0.472)。 例例6.8 某檢測參數某檢測參數f隨時間隨時間t的采樣結果如表的采樣結果如表5.1，用數據插，用數據插值法計算值法計算t=2，7，12，17，22，17，32，37，42，47，52，57時的時的f值。值。2. 二維數據插值二維數據插值在在MA

17、TLAB中，提供了解決二維插值問題的函數中，提供了解決二維插值問題的函數interp2，其調用格式為：，其調用格式為：Z1=interp2(X,Y,Z,X1,Y1,method)其中其中X,Y是兩個向量，分別描述兩個參數的采樣點，是兩個向量，分別描述兩個參數的采樣點，Z是與參數采樣點對應的函數值，是與參數采樣點對應的函數值，X1,Y1是兩個向是兩個向量或標量，描述欲插值的點。量或標量，描述欲插值的點。Z1是根據相應的插是根據相應的插值方法得到的插值結果。值方法得到的插值結果。 method的取值與一維插的取值與一維插值函數相同。值函數相同。X,Y,Z也可以是矩陣形式。也可以是矩陣形式。同樣，同

18、樣，X1,Y1的取值范圍不能超出的取值范圍不能超出X,Y的給定范圍，的給定范圍，否則，會給出否則，會給出“NaN”錯誤。錯誤。例例6.9 設設z=x2+y2，對，對z函數在函數在0,10,2區域內進行區域內進行插值。插值。例例6.10 某實驗對一根長某實驗對一根長10米的鋼軌進行熱源的溫度米的鋼軌進行熱源的溫度傳播測試。用傳播測試。用x表示測量點表示測量點0:2.5:10(米米)，用，用h表示表示測量時間測量時間0:30:60(秒秒)，用，用T表示測試所得各點的溫表示測試所得各點的溫度度()。試用線性插值求出在一分鐘內每隔。試用線性插值求出在一分鐘內每隔10秒、秒、鋼軌每隔鋼軌每隔0.5米處的

19、溫度。米處的溫度。在在MATLAB中，用中，用polyfit函數來求得最小二乘擬合多項式的函數來求得最小二乘擬合多項式的系數，再用系數，再用polyval函數按所得的多項式計算所給出的點上函數按所得的多項式計算所給出的點上的函數近似值。的函數近似值。polyfit函數的調用格式為：函數的調用格式為：P,S=polyfit(X,Y,m)函數根據采樣點函數根據采樣點X和采樣點函數值和采樣點函數值Y，產生一個，產生一個m次多項式次多項式P及其在采樣點的誤差向量及其在采樣點的誤差向量S。其中。其中X,Y是兩個等長的向量，是兩個等長的向量，P是一個長度為是一個長度為m+1的向量，的向量，P的元素為多項式

20、系數。的元素為多項式系數。polyval函數的功能是按多項式的系數計算函數的功能是按多項式的系數計算x點多項式的值。點多項式的值。例例6.11 用一個用一個3次多項式在區間次多項式在區間0，2內逼近函數。內逼近函數。命令如下：命令如下：X=linspace(0,2*pi,50);Y=sin(X);P=polyfit(X,Y,3) %得到得到3次多項式的系數和誤差次多項式的系數和誤差1. 多項式的四則運算多項式的四則運算（1）多項式的加減運算）多項式的加減運算（2）多項式乘法運算）多項式乘法運算函數函數conv(P1,P2)用于求多項式用于求多項式P1和和P2的乘積。的乘積。這里，這里，P1、P

21、2是兩個多項式系數向量。是兩個多項式系數向量。（3）多項式除法）多項式除法函數函數Q,r=deconv(P1,P2)用于對多項式用于對多項式P1和和P2作除作除法運算。其中法運算。其中Q返回多項式返回多項式P1除以除以P2的商式，的商式，r返返回回P1除以除以P2的余式。這里，的余式。這里，Q和和r仍是多項式系數仍是多項式系數向量。向量。deconv是是conv的逆函數，即有的逆函數，即有P1=conv(P2,Q)+r。2. 多項式的導函數多項式的導函數對多項式求導數的函數是：對多項式求導數的函數是：p=polyder(P)：求多項式：求多項式P的導函數的導函數p=polyder(P,Q)：求

22、：求PQ的導函數的導函數p,q=polyder(P,Q)：求：求P/Q的導函數，導函數的分的導函數，導函數的分子存入子存入p，分母存入，分母存入q。上述函數中，參數上述函數中，參數P,Q是多項式的向量表示，結果是多項式的向量表示，結果p,q也是多項式的向量表示。也是多項式的向量表示。3. 多項式求值多項式求值MATLAB提供了兩種求多項式值的函數：提供了兩種求多項式值的函數：polyval與與polyvalm，它們的輸入參數均為多項式系數向量，它們的輸入參數均為多項式系數向量P和自變量和自變量x。兩者的區別在于前者是代數多項式求。兩者的區別在于前者是代數多項式求值，而后者是矩陣多項式求值。值，

23、而后者是矩陣多項式求值。(1)代數多項式求值代數多項式求值polyval函數用來求代數多項式的值，其調用函數用來求代數多項式的值，其調用格式為：格式為：Y=polyval(P,x)若若x為一數值，則求多項式在該點的值；若為一數值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值。元素求其多項式的值。例例6.14 已知多項式已知多項式x4+8x3-10，分別取，分別取x=1.2和和一個一個23矩陣為自變量計算該多項式的值。矩陣為自變量計算該多項式的值。(2)矩陣多項式求值矩陣多項式求值polyvalm函數用來求矩陣多項式的值，其

24、調用格式與函數用來求矩陣多項式的值，其調用格式與polyval相同，但含義不同。相同，但含義不同。polyvalm函數要求函數要求x為方陣，它以方為方陣，它以方陣為自變量求多項式的值。設陣為自變量求多項式的值。設A為方陣，為方陣，P代表多項式代表多項式x3-5x2+8，那么，那么polyvalm(P,A)的含義是：的含義是：A*A*A-5*A*A+8*eye(size(A)而而polyval(P,A)的含義是：的含義是：例例6.15 仍以多項式仍以多項式x4+8x3-10為例，取一個為例，取一個22矩陣為自變量矩陣為自變量分別用分別用polyval和和polyvalm計算該多項式的值。計算該多

25、項式的值。4. 多項式求根多項式求根n次多項式具有次多項式具有n個根，當然這些根可能是實個根，當然這些根可能是實根，也可能含有若干對共軛復根。根，也可能含有若干對共軛復根。MATLAB提供的提供的roots函數用于求多項式的函數用于求多項式的全部根，其調用格式為：全部根，其調用格式為：x=roots(P)其中其中P為多項式的系數向量，求得的根賦給向為多項式的系數向量，求得的根賦給向量量x，即，即x(1),x(2),x(n)分別代表多項式的分別代表多項式的n個根。個根。例例6.16 求多項式求多項式x4+8x3-10的根。的根。命令如下：命令如下：A=1,8,0,0,-10;x=roots(A)

26、若已知多項式的全部根，則可以用若已知多項式的全部根，則可以用poly函數建立起函數建立起該多項式，其調用格式為：該多項式，其調用格式為：P=poly(x)若若x為具有為具有n個元素的向量，則個元素的向量，則poly(x)建立以建立以x為其為其根的多項式，且將該多項式的系數賦給向量根的多項式，且將該多項式的系數賦給向量P。例例6.17 已知已知 f(x)(1) 計算計算f(x)=0 的全部根。的全部根。(2) 由方程由方程f(x)=0的根構造一個多項式的根構造一個多項式g(x)，并，并與與f(x)進行對比。進行對比。命令如下：命令如下：P=3,0,4,-5,-7.2,5;X=roots(P) %

27、求方程求方程f(x)=0的根的根G=poly(X) %求多項式求多項式g(x)6.2 數值微積分數值微積分1. 數值差分與差商數值差分與差商2. 數值微分的實現數值微分的實現在在MATLAB中，沒有直接提供求數值導數的函數，只有計中，沒有直接提供求數值導數的函數，只有計算向前差分的函數算向前差分的函數diff，其調用格式為：，其調用格式為：DX=diff(X)：計算向量：計算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：計算：計算X的的n階向前差分。例如，階向前差分。例如，diff(X,2)=diff(diff(X)。DX=di

28、ff(A,n,dim)：計算矩陣：計算矩陣A的的n階差分，階差分，dim=1時時(缺省狀缺省狀態態)，按列計算差分；，按列計算差分；dim=2，按行計算差分。，按行計算差分。例例6.18 設設x由由0,2間均勻分布的間均勻分布的10個點組成，求個點組成，求sinx的的13階差分。階差分。命令如下：命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y); %計算計算Y的一階差分的一階差分D2Y=diff(Y,2); %計算計算Y的二階差分的二階差分，也可用命令，也可用命令diff(DY)計算計算D3Y=diff(Y,3); %計算計算Y的三階差分的三階差分，

29、也可用，也可用diff(D2Y)或或diff(DY,2)例例6.19 用不同的方法求函數用不同的方法求函數f(x)的數值導數，并在同一個坐的數值導數，并在同一個坐標系中做出標系中做出f(x)的圖像。的圖像。程序如下：程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多項式次多項式p擬合擬合f(x)dp=polyder(p); %對擬

30、合多項式對擬合多項式p求導數求導數dpdpx=polyval(dp,x); %求求dp在假設點的函數值在假設點的函數值dx=diff(f(x,3.01)/0.01; %直接對直接對f(x)求數值導數求數值導數gx=g(x); %求函數求函數f的導函數的導函數g在假設點的導數在假設點的導數plot(x,dpx,x,dx,.,x,gx,-); %作圖作圖1. 數值積分基本原理數值積分基本原理 求解定積分的數值方法多種多樣，如簡單求解定積分的數值方法多種多樣，如簡單的梯形法、辛普生的梯形法、辛普生(Simpson) 法、牛頓法、牛頓柯特斯柯特斯(Newton-Cotes)法等都是經常采用的法等都是經

31、常采用的方法。它們的基本思想都是將整個積分區方法。它們的基本思想都是將整個積分區間間a,b分成分成n個子區間個子區間xi,xi+1，i=1,2,n，其中其中x1=a，xn+1=b。這樣求定積分問題就分。這樣求定積分問題就分解為求和問題。解為求和問題。2. 數值積分的實現數值積分的實現(1)被積函數是一個解析式被積函數是一個解析式MATLAB提供了提供了quad函數和函數和quadl函數來求函數來求定積分。它們的調用格式為：定積分。它們的調用格式為： quad(filename,a,b,tol,trace) quadl(filename,a,b,tol,trace)例例6.20 用兩種不同的方法

32、求定積分。用兩種不同的方法求定積分。先建立一個函數文件先建立一個函數文件ex.m：function ex=ex(x)ex=exp(-x.2);然后在然后在MATLAB命令窗口命令窗口，輸入命令：，輸入命令：format longI=quad(ex,0,1) %注意函數名應加字符引號注意函數名應加字符引號I = 0.74682418072642I=quadl(ex,0,1)I = 0.74682413398845也可不建立關于被積函數的函數文件，而使用語句函數也可不建立關于被積函數的函數文件，而使用語句函數(內聯函數內聯函數)求解，命令如求解，命令如下：下：g=inline(exp(-x.2);

33、 %定義一個語句函數定義一個語句函數g(x)=exp(-x2)I=quadl(g,0,1) %注意函數名不加注意函數名不加號號I = 0.74682413398845format short(2)被積函數由一個表格定義被積函數由一個表格定義在科學實驗和工程應用中，函數關系往往是不知道在科學實驗和工程應用中，函數關系往往是不知道的，只有實驗測定的一組樣本點和樣本值，這時，的，只有實驗測定的一組樣本點和樣本值，這時，就無法使用就無法使用quad函數計算其定積分。在函數計算其定積分。在MATLAB中，對由表格形式定義的函數關系的求定積分問中，對由表格形式定義的函數關系的求定積分問題用題用trapz(

34、X,Y)函數。其中向量函數。其中向量X、Y定義函數關定義函數關系系Y=f(X)。X、Y是兩個等長的向量：是兩個等長的向量：X=(x1,x2,xn)，Y=(y1,y2,yn)，并且，并且x1x2 cholMatrix must be positive definite命令執行時，出現錯誤信息，說明命令執行時，出現錯誤信息，說明A為非正定矩陣。為非正定矩陣。迭代解法非常適合求解大型系數矩陣的方程組。在數值分析迭代解法非常適合求解大型系數矩陣的方程組。在數值分析中，迭代解法主要包括中，迭代解法主要包括 Jacobi迭代法、迭代法、Gauss-Serdel迭代迭代法、超松弛迭代法和兩步迭代法。法、超松

35、弛迭代法和兩步迭代法。1Jacobi迭代法迭代法對于線性方程組對于線性方程組Ax=b，如果，如果A為非奇異方陣，即為非奇異方陣，即aii0(i=1,2,n)，則可將，則可將A分解為分解為A=D-L-U，其中，其中D為對為對角陣，其元素為角陣，其元素為A的對角元素，的對角元素，L與與U為為A的下三角陣和上的下三角陣和上三角陣，于是三角陣，于是Ax=b化為：化為：x=D-1(L+U)x+D-1b與之對應的迭代公式為：與之對應的迭代公式為：x(k+1)=D-1(L+U)x(k)+D-1b這就是這就是Jacobi迭代公式。如果序列迭代公式。如果序列x(k+1)收斂于收斂于x，則，則x必必是方程是方程A

36、x=b的解。的解。Jacobi迭代法的迭代法的MATLAB函數文件函數文件Jacobi.m如下：如下：function y,n=jacobi(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=B*x0+f; n=n+1;end例例6.28 用用Jacobi迭代法求解線性方程組。設迭代初值為迭代法求解線性方程組。設迭代初值為0，迭代精度為迭代精度為10-6。在命令中調用函數文件在命令中調用函數文件Jacobi.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=jacobi(

37、A,b,0,0,0,1.0e-6)2Gauss-Serdel迭代法迭代法在在Jacobi迭代過程中，計算時，已經得到，不必再用，即原迭代過程中，計算時，已經得到，不必再用，即原來的迭代公式來的迭代公式Dx(k+1)=(L+U)x(k)+b可以改進為可以改進為Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：，于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b該式即為該式即為Gauss-Serdel迭代公式。和迭代公式。和Jacobi迭代相比，迭代相比，Gauss-Serdel迭代用新分量代替舊分量，精度會高些。迭代用新分量代替舊分量，精度會高些。Gauss-Serdel

38、迭代法的迭代法的MATLAB函數文件函數文件gauseidel.m如下：如下：function y,n=gauseidel(A,b,x0,eps)if nargin=3 eps=1.0e-6;elseif nargin=eps x0=y; y=G*x0+f; n=n+1;end例例6.29 用用Gauss-Serdel迭代法求解下列線性方程組。設迭代迭代法求解下列線性方程組。設迭代初值為初值為0，迭代精度為，迭代精度為10-6。在命令中調用函數文件在命令中調用函數文件gauseidel.m，命令如下：，命令如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=ga

39、useidel(A,b,0,0,0,1.0e-6)例例6.30 分別用分別用Jacobi迭代和迭代和Gauss-Serdel迭代法求解下列線迭代法求解下列線性方程組，看是否收斂。性方程組，看是否收斂。命令如下：命令如下：a=1,2,-2;1,1,1;2,2,1;b=9;7;6;x,n=jacobi(a,b,0;0;0)x,n=gauseidel(a,b,0;0;0)線性方程組的求解分為兩類：一類是求方程組的惟一解即特線性方程組的求解分為兩類：一類是求方程組的惟一解即特解，另一類是求方程組的無窮解即通解。這里對線性方程解，另一類是求方程組的無窮解即通解。這里對線性方程組組 Ax=b的求解理論作一

40、個歸納。的求解理論作一個歸納。(1)當系數矩陣當系數矩陣A是一個滿秩方陣時，方程是一個滿秩方陣時，方程Ax=b稱為恰定方程，稱為恰定方程，方程有惟一解方程有惟一解x=A-1b，這是最基本的一種情況。一般用，這是最基本的一種情況。一般用x=Ab求解速度更快。求解速度更快。(2)當方程組右端向量當方程組右端向量b=0時，方程稱為齊次方程組。齊次方時，方程稱為齊次方程組。齊次方程組總有零解，因此稱解程組總有零解，因此稱解x=0為平凡解。當系數矩陣為平凡解。當系數矩陣A的的秩小于秩小于n(n為方程組中未知變量的個數為方程組中未知變量的個數)時，齊次方程組有時，齊次方程組有無窮多個非平凡解，其通解中包含

41、無窮多個非平凡解，其通解中包含n-rank(A)個線性無關個線性無關的解向量，用的解向量，用MATLAB的函數的函數null(A,r)可求得基礎解系。可求得基礎解系。(3)當方程組右端向量當方程組右端向量b0時，系數矩陣的秩時，系數矩陣的秩rank(A)與其增廣矩陣的秩與其增廣矩陣的秩rank(A,b)是判斷其是否有解是判斷其是否有解的基本條件：的基本條件：當當rank(A)=rank(A,b)=n時，方程組有惟一解：時，方程組有惟一解：x=Ab 或或 x=pinv(A)*b。當當rank(A)=rank(A,b)n時，方程組有無窮多個時，方程組有無窮多個解，其通解解，其通解=方程組的一個特解

42、方程組的一個特解+對應的齊次方程對應的齊次方程組組Ax=0的通解。可以用的通解。可以用Ab求得方程組的一個特求得方程組的一個特解，用解，用null(A,r)求得該方程組所對應的齊次方程求得該方程組所對應的齊次方程組的基礎解系，基礎解系中包含組的基礎解系，基礎解系中包含n-rank(A)個線性個線性無關的解向量。無關的解向量。當當rank(A)rank(A,b)時，方程組無解。時，方程組無解。 有了上面這些討論，可以設計一個求解線有了上面這些討論，可以設計一個求解線性方程組的函數文件性方程組的函數文件line_solution.m。在例。在例中可以調用中可以調用line_solution.m文件

43、來解線性方文件來解線性方程組。程組。 6.5 非線性方程與最優化問題求解非線性方程與最優化問題求解1. 單變量非線性方程求解單變量非線性方程求解 在在MATLAB中提供了一個中提供了一個fzero函數，可以用來求單變量函數，可以用來求單變量非線性方程的根。該函數的調用格式為：非線性方程的根。該函數的調用格式為： z=fzero(fname,x0,tol,trace)其中其中fname是待求根的函數文件名，是待求根的函數文件名，x0為搜索的起點。一個為搜索的起點。一個函數可能有多個根，但函數可能有多個根，但fzero函數只給出離函數只給出離x0最近的那個根。最近的那個根。tol控制結果的相對精度

44、，缺省時取控制結果的相對精度，缺省時取tol=eps，trace 指定指定迭代信息是否在運算中顯示，為迭代信息是否在運算中顯示，為1時顯示，為時顯示，為0時不顯示，時不顯示，缺省時取缺省時取trace=0。例例6.33 求求f(x)在在x0=-5和和x0=1作為迭代初值時的零點。作為迭代初值時的零點。先建立函數文件先建立函數文件fz.m：function f=fz(x)f=x-1/x+5;然后調用然后調用fzero函數求根。：函數求根。：fzero(fz,-5) %以以-5作為迭代初值作為迭代初值ans = -5.1926fzero(fz,1) %以以1作為迭代初值作為迭代初值ans = 0.

45、1926 2. 非線性方程組的求解非線性方程組的求解 對于非線性方程組對于非線性方程組F(X)=0，用，用fsolve函數求其數值解。函數求其數值解。fsolve函數的調用格式為：函數的調用格式為： X=fsolve(fun,X0,option)其中其中X為返回的解，為返回的解，fun是用于定義需求解的非線性方程組的是用于定義需求解的非線性方程組的函數文件名，函數文件名，X0是求根過程的初值，是求根過程的初值，option為最優化工具為最優化工具箱的選項設定。最優化工具箱提供了箱的選項設定。最優化工具箱提供了20多個選項，用戶可多個選項，用戶可以使用以使用optimset命令將它們顯示出來。如

46、果想改變其中某命令將它們顯示出來。如果想改變其中某個選項，則可以調用個選項，則可以調用optimset()函數來完成。例如，函數來完成。例如，Display選項決定函數調用時中間結果的顯示方式，其中選項決定函數調用時中間結果的顯示方式，其中off為不顯示，為不顯示，iter表示每步都顯示，表示每步都顯示，final只顯示只顯示最終結果。最終結果。optimset(Display,off)將設定將設定Display選項為選項為off。例例6.34 求下列方程組在求下列方程組在(1，1，1)附近的解并對結果進行驗附近的解并對結果進行驗證。證。首先建立函數文件首先建立函數文件myfun.m。func

47、tion F=myfun (X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z2*exp(x);F(2)=x+y+z;F(3)=x*y*z;在給定的初值在給定的初值x0=1，y0=1，z0=1下，調用下，調用fsolve函數求方程的根。函數求方程的根。X=fsolve(myfun,1,1,1,optimset(Display, off)X = 0.0224 -0.0224 -0.0000 在實際應用中，許多科學研究和工程計算問題都可以歸結為在實際應用中，許多科學研究和工程計算問題都可以歸結為一個最小化問題，如能量最小、時間最短等。一個最小化問題，如能量最小、時間最短等

48、。MATLAB提供了提供了3個求最小值的函數，它們的調用格式為：個求最小值的函數，它們的調用格式為：(1)x,fval=fminbnd(filename,x1,x2,option)：求一元函數在：求一元函數在(xl,x2)區間中的極小值點區間中的極小值點x和最小值和最小值fval。(2)x,fval=fminsearch(filename,x0,option)：基于單純形算法：基于單純形算法求多元函數的極小值點求多元函數的極小值點x和最小值和最小值fval。(3)x,fval=fminunc(filename,x0,option)：基于擬牛頓法求多：基于擬牛頓法求多元函數的極小值點元函數的極小

49、值點x和最小值和最小值fval。MATLAB沒有專門提供求函數最大值的函數，但只要注意沒有專門提供求函數最大值的函數，但只要注意到到-f(x)在區間在區間(a,b)上的最小值就是上的最小值就是f(x)在在(a,b)的最大值，的最大值，所以所以fminbnd(-f,x1,x2)返回函數返回函數f(x)在區間在區間(x1,x2)上的最大上的最大值。值。例例6.36 求函數在區間求函數在區間(-10，1)和和(1，10)上的最小值點。上的最小值點。首先建立函數文件首先建立函數文件fx.m：function f=f(x)f=x-1/x+5;上述函數文件也可用一個語句函數代替：上述函數文件也可用一個語句

50、函數代替：f=inline(x-1/x+5)再在再在MATLAB命令窗口，輸入命令：命令窗口，輸入命令：fminbnd(fx,-10,-1) %求函數在求函數在(-10，-1)內的最小值點和最小值內的最小值點和最小值fminbnd(f,1,10) %求函數在求函數在(1，10)內的最小值點。注意函數名內的最小值點。注意函數名f不用不用加加例例6.37 求函數求函數f在在(0.5,0.5,0.5)附近的最小值。附近的最小值。建立函數文件建立函數文件fxyz.m：function f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;在在MA

51、LAB命令窗口，輸入命令：命令窗口，輸入命令：U,fmin=fminsearch(fxyz,0.5,0.5,0.5) %求函數的最小值點和最小值求函數的最小值點和最小值MATLAB最優化工具箱提供了一個最優化工具箱提供了一個fmincon函數，專門用于求解各種約函數，專門用于求解各種約束下的最優化問題。該函數的調用格式為：束下的最優化問題。該函數的調用格式為：x,fval=fmincon(filename,x0,A,b, Aeq,beq,Lbnd,Ubnd, NonF,option)其中其中x、fval、filename、x0和和option的含義與求最小值函數相同。其余參的含義與求最小值函數

52、相同。其余參數為約束條件，參數數為約束條件，參數NonF為非線性約束函數的為非線性約束函數的M文件名。如果某個文件名。如果某個約束不存在，則用空矩陣來表示。約束不存在，則用空矩陣來表示。例例6.38 求解有約束最優化問題。求解有約束最優化問題。首先編寫目標函數首先編寫目標函數M文件文件fop.m。function f=fop(x)f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3;再設定約束條件，并調用再設定約束條件，并調用fmincon函數求解此約束最優化問題。函數求解此約束最優化問題。x0=0.5;0.5;A=-1,-0.5;-0.5,-1;b=-0.4

53、;-0.5;lb=0;0;option=optimset; option.LargeScale=off; option.Display =off;x,f=fmincon(fop ,x0,A,b,lb,option)6.6 常微分方程的數值求解常微分方程的數值求解 基于龍格庫塔法，基于龍格庫塔法，MATLAB提供了求常微分方程數值提供了求常微分方程數值解的函數，一般調用格式為：解的函數，一般調用格式為： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中其中fname是定義是定義f(t,y)的函數文件名，該函數文件必須返回的函數文件名，該函數文件必須返回一個列向量。一個列向量。tspan形式為形式為t0,tf,表示求解區間。表示求解區間。y0是初始是初始狀態列向量。狀態列向量。t和和y分別給出時間向量和相應的狀態向量。分別給出時間向量和相應的狀態向量。例例6.39 設有初值問題，試求其數值解，并與精確解相比較設有初值問題，試求其數值解，并與精確解相比較(精確解為精確解為y(t)=)。 (1) 建立函數文件建立函數文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。求解微分方程。t0=0;tf=10;y0=2