matlab第2 數(shù)值數(shù)組及其運算(續(xù))

第二講數(shù)值數(shù)組及其運算

(續(xù))一、多項式運算1、多項式的創(chuàng)建（1）多項式系數(shù)向量的直接輸入法（2）利用指令：P=poly(AR)產(chǎn)生多項式系數(shù)向量若AR是方陣，則多項式P就是該方陣的特征多項式若AR=［ar1ar2ar3……arn］，則AR的元素被認為是多項式P的根【例7】求3階方陣A的特征多項式。A=[111213;141516;171819];PA=poly(A)PPA=poly2str(PA,'s')PA=1.0000-45.0000-18.00000.0000PPA=s^3-45s^2-18s+1.8303e-0142、多項式運算函數(shù)指令含義p=conv(p1,p2)P是多項式p1和p2的乘積多項式[q,r]=deconv(p1,p2)多項式p1被p2除的商多項式為q，而余多項式為rp=ploy(AR)求方陣AR的特征多項式p；或求向量AR指定根所對應的多項式dp=polyder(p)求多項式p的導數(shù)多項式dpdp=polyder(p1,p2)求多項式p1,p2乘積的導數(shù)多項式dp[Num,Den]=ployder(p1,p2)對有理分式(p1/p2)求導數(shù)所得的有理分式為(Num/Den)P=polyfit(x,y,n)求x,y向量給定的數(shù)據(jù)的n階擬合多項式ppA=polyval(p,S)按數(shù)組運算規(guī)則計算多項式值；p為多項式，S為矩陣PM=polyvalm(p,S)按矩陣運算規(guī)則計算多項式值；p為多項式，S為矩陣[r,p,k]=residue(b,a)部分分式展開，b,a分別是分子、分母多項式系數(shù)向量；r、p、k分別是留數(shù)、極點和直項r=roots(p)求多項式p的根【例1】求的“商”及“余”多項式。p1=conv([1,0,2],conv([1,4],[1,1]));

%計算分子多項式 p2=[1011];

%注意缺項補零[q,r]=deconv(p1,p2);cq='商多項式為';cr='余多項式為';disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])商多項式為s+5余多項式為5s^2+4s+3【例2】

求多項式的微分polyder(p)

polyder(a,b):求多項式a,b乘積的微分[p,q]=polyder(a,b):求多項式a,b商的微分a=[12345];poly2str(a,'x')ans=x^4+2x^3+3x^2+4x+5b=polyder(a)b=4664poly2str(b,'x')ans=4x^3+6x^2+6x+43．代數(shù)多項式求值polyval函數(shù)用來求代數(shù)多項式的值，其調(diào)用格式為：Y=polyval(P,x)若x為一數(shù)值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值?！纠?】已知多項式x4+8x3-10，分別取x=1.2和一個2×3矩陣為自變量計算該多項式的值。A=[1,8,0,0,-10];%4次多項式系數(shù)x=1.2;%取自變量為一數(shù)值y1=polyval(A,x)y1=5.8976x=[-1,1.2,-1.4;2,-1.8,1.6];%給出一個矩陣xy2=polyval(A,x)y2=-17.00005.8976-28.110470.0000-46.158429.32164．矩陣多項式求值polyvalm函數(shù)用來求矩陣多項式的值，其調(diào)用格式與polyval相同，但含義不同。polyvalm函數(shù)要求x為方陣,它以方陣為自變量求多項式的值。設A為方陣，P代表多項式x3-5x2+8，那么polyvalm(P,A)的含義是：A*A*A-5*A*A+8*eye(size(A))而polyval(P,A)的含義是：A.*A.*A-5*A.*A+8*ones(size(A))【例4】仍以多項式x4+8x3-10為例，取一個2×2矩陣為自變量分別用polyval和polyvalm計算該多項式的值。A=[1,8,0,0,-10];%多項式系數(shù)x=[-1,1.2;2,-1.8];%給出一個矩陣xy1=polyval(A,x)%計算代數(shù)多項式的值y1=-17.00005.897670.0000-46.1584y2=polyvalm(A,x)%計算矩陣多項式的值y2=-60.584050.649684.4160-94.35045多項式求根n次多項式具有n個根，當然這些根可能是實根，也可能含有若干對共軛復根。MATLAB提供的roots函數(shù)用于求多項式的全部根，其調(diào)用格式為：x=roots(P)其中P為多項式的系數(shù)向量，求得的根賦給向量x，即x(1),x(2),…,x(n)分別代表多項式的n個根。【例5】求多項式x4+8x3-10的根。命令如下：A=[1,8,0,0,-10];x=roots(A)x=-8.01941.0344-0.5075+0.9736i-0.5075-0.9736i二、數(shù)據(jù)統(tǒng)計處理(一)最大值和最小值

MATLAB提供的求數(shù)據(jù)序列的最大值和最小值的函數(shù)分別為max和min，兩個函數(shù)的調(diào)用格式和操作過程類似。1．求向量的最大值和最小值求一個向量X的最大值的函數(shù)有兩種調(diào)用格式，分別是：(1)y=max(X)：返回向量X的最大值存入y，如果X中包含復數(shù)元素，則按模取最大值。(2)[y,I]=max(X)：返回向量X的最大值存入y，最大值的序號存入I，如果X中包含復數(shù)元素，則按模取最大值。2．求矩陣的最大值和最小值

求矩陣A的最大值的函數(shù)有3種調(diào)用格式，分別是：(1)max(A)：返回一個行向量，向量的第i個元素是矩陣A的第i列上的最大值。(2)[Y,U]=max(A)：返回行向量Y和U，Y向量記錄A的每列的最大值，U向量記錄每列最大值的行號。(3)max(A,[],dim)：dim取1或2。dim取1時，該函數(shù)和max(A)完全相同；dim取2時，該函數(shù)返回一個列向量，其第i個元素是A矩陣的第i行上的最大值。3．兩個向量或矩陣對應元素的比較函數(shù)max和min還能對兩個同型的向量或矩陣進行比較，調(diào)用格式為：(1)U=max(A,B)：A,B是兩個同型的向量或矩陣，結(jié)果U是與A,B同型的向量或矩陣，U的每個元素等于A,B對應元素的較大者。(2)U=max(A,n)：n是一個標量，結(jié)果U是與A同型的向量或矩陣，U的每個元素等于A對應元素和n中的較大者。min函數(shù)的用法和max完全相同。(二)求和與求積數(shù)據(jù)序列求和與求積的函數(shù)是sum和prod，其使用方法類似。設X是一個向量，A是一個矩陣，函數(shù)的調(diào)用格式為：sum(X)：返回向量X各元素的和。prod(X)：返回向量X各元素的乘積。sum(A)：返回一個行向量，第i個元素是A的第i列的元素和。prod(A)：返回一個行向量，第i個元素是A的第i列的元素乘積。sum(A,dim)：dim=1時，該函數(shù)等同于sum(A)；dim=2時，返回一個列向量，其第i個元素是A的第i行的各元素之和。prod(A,dim)：dim=1時，該函數(shù)等同于prod(A)；dim=2時，返回一個列向量，其第i個元素是A的第i行的各元素乘積。(三)平均值和中值求數(shù)據(jù)序列平均值的函數(shù)是mean，求數(shù)據(jù)序列中值的函數(shù)是median。兩個函數(shù)的調(diào)用格式為：mean(X)：返回向量X的算術(shù)平均值。median(X)：返回向量X的中值。mean(A)：返回一個行向量，其第i個元素是A的第i列的算術(shù)平均值。median(A)：返回一個行向量，其第i個元素是A的第i列的中值。mean(A,dim)：dim=1時，該函數(shù)等同于mean(A)；dim=2時，返回一個列向量，其第i個元素是A的第i行的算術(shù)平均值。median(A,dim)：dim=1時，該函數(shù)等同于median(A)；dim=2時，返回一個列向量，其第i個元素是A的第i行的中值。(四)累加和與累乘積

在MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素的累加和與累乘積向量，函數(shù)的調(diào)用格式為：cumsum(X)：返回向量X累加和向量。cumprod(X)：返回向量X累乘積向量。cumsum(A)：返回一個矩陣，其第i列是A的第i列的累加和向量。cumprod(A)：返回一個矩陣，其第i列是A的第i列的累乘積向量。cumsum(A,dim)：dim=1時，函數(shù)等同于cumsum(A)；dim=2時，返回一個矩陣，其第i行是A的第i行的累加和向量。cumprod(A,dim)：dim=1時，函數(shù)等同于cumprod(A)；dim=2時，返回一個向量，其第i行是A的第i行的累乘積向量。(五)標準方差與相關系數(shù)1．求標準方差

在MATLAB中，提供了計算數(shù)據(jù)序列的標準方差的函數(shù)std。

std(X)：返回一個標準方差。

std(A)：返回一個行向量，每個元素是矩陣A各列或各行的標準方差。

Y=std(A,flag,dim)

缺省flag=0，dim=1當dim=1時，求各列元素的標準方差；當dim=2時，則求各行元素的標準方差。當flag=0時，按S1所列公式計算標準方差；當flag=1時，按S2所列公式計算標準方差。

【例6】對二維矩陣x，從不同維方向求出其標準方差。命令如下：x=[4,5,6;1,4,8];y1=std(x,0,1)y2=std(x,1,1)y3=std(x,0,2)y4=std(x,1,2)2．相關系數(shù)

MATLAB提供了corrcoef函數(shù)，可以求出數(shù)據(jù)的相關系數(shù)矩陣。corrcoef函數(shù)的調(diào)用格式為：corrcoef(X)：返回從矩陣X形成的一個相關系數(shù)矩陣。此相關系數(shù)矩陣的大小與矩陣X一樣。它把矩陣X的每列作為一個變量，然后求它們的相關系數(shù)。corrcoef(X,Y)：在這里，X,Y是向量，它們與corrcoef([X,Y])的作用一樣?！纠?】生成滿足正態(tài)分布的10000×5隨機矩陣，然后求各列元素的均值和標準方差，再求這5列隨機數(shù)據(jù)的相關系數(shù)矩陣。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)(六)排序

MATLAB中對向量X是排序函數(shù)是sort(X)，函數(shù)返回一個對X中的元素按升序排列的新向量。

[Y,I]=sort(A,dim)：對矩陣A的各列或各行重新排序若dim=1，則按列排；若dim=2，則按行排。

Y是排序后的矩陣，而I記錄Y中的元素在A中位置。【例8】對下列矩陣做各種排序。A=[1,-8,5;4,12,6;13,7,-13];sort(A)%對A的每列按升序排序-sort(-A,2)%對A的每行按降序排序[X,I]=sort(A)

%對A按列排序，并將每個元素所在的行號送矩陣I(一)一維數(shù)據(jù)插值數(shù)據(jù)插值的任務是根據(jù)采集到的離散數(shù)據(jù)構(gòu)造一個函數(shù)g(x)，既與真實函數(shù)f(x)接近，又有很好的性質(zhì)。Y1=interp1(X,Y,X1,'method')

一維插值函數(shù)根據(jù)X,Y的值，計算函數(shù)在X1處的值。X,Y是兩個等長的已知向量，分別描述采樣點和樣本值，X1是一個向量或標量，描述欲插值的點，Y1是與X1等長的插值結(jié)果。method是插值方法，例如：‘linear’線性插值‘nearest’最近點插值‘cubic’三次多項式插值‘spline’三次樣條插值或Y1=spline(X,Y,X1)計算三次樣條插值的專門函數(shù)三、數(shù)據(jù)插值注意：X1的取值范圍不能超出X的給定范圍，否則，會給出“NaN”錯誤?！纠?】某觀測站測得某日6:00時至18:00時之間每隔2小時的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時之間每隔2小時各點的近似溫度(℃)。解：設時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲存室外溫度。命令如下：h=6:2:18;t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';XI=6.5:2:17.5YI=interp1(h,t,XI,'spline')%用3次樣條插值計算(二)二維數(shù)據(jù)插值

在MATLAB中，提供了解決二維插值問題的函數(shù)interp2，其調(diào)用格式為：Z1=interp2(X,Y,Z,X1,Y1,'method')

其中X,Y是兩個向量，分別描述兩個參數(shù)的采樣點，Z是與參數(shù)采樣點對應的函數(shù)值，X1,Y1是兩個向量或標量，描述欲插值的點。Z1是根據(jù)相應的插值方法得到的插值結(jié)果。method的取值與一維插值函數(shù)相同。X,Y,Z也可以是矩陣形式。同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會給出“NaN”錯誤。【例10】設z=x2+y2，對z函數(shù)在[0,1]×[0,2]區(qū)域內(nèi)進行插值。x=0:0.1:1;y=0:0.2:2;[X,Y]=meshgrid(x,y);%產(chǎn)生自變量網(wǎng)格坐標Z=X.^2+Y.^2;

%求對應的函數(shù)值interp2(x,y,Z,0.5,0.5)

%在(0.5,0.5)點插值【例11】某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度(℃)。試用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi);mesh(xi,hi,TI)在MATLAB中，用polyfit函數(shù)來求得最小二乘擬合多項式的系數(shù)，再用polyval函數(shù)按所得的多項式計算所給出的點上的函數(shù)近似值。

[P,S]=polyfit(X,Y,m)函數(shù)根據(jù)采樣點X和采樣點函數(shù)值Y，產(chǎn)生一個m次多項式P及其在采樣點的誤差向量S。其中X,Y是兩個等長的向量，P是一個長度為m+1的向量，P的元素為多項式系數(shù)。

polyval函數(shù)的功能是按多項式的系數(shù)計算x點多項式的值。四、曲線擬合【例12】用一個三次多項式在區(qū)間[0,2π]內(nèi)逼近函數(shù)sin(x)。X=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,3)%得到3次多項式的系數(shù)和誤差五、數(shù)值積分與微分(一)數(shù)值積分1、變步長辛普生法基本思想都是將整個積分區(qū)間[a,b]分成n個子區(qū)間[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。這樣求定積分問題就分解為求和問題。MATLAB給出了quad和quadl函數(shù)來求定積分。[I,n]=quad('filename',a,b,tol,trace)[I,n]=quadl('filename',a,b,tol,trace)其中filename是被積函數(shù)名。a和b分別是定積分的下限和上限。tol用來控制積分精度，缺省時取tol=10-6。trace控制是否展現(xiàn)積分過程，若取非0則展現(xiàn)積分過程，取0則不展現(xiàn)，缺省時取trace=0。返回參數(shù)I即定積分值，n為被積函數(shù)的調(diào)用次數(shù)?！纠?3】求定積分解：(1)建立被積函數(shù)文件fesin.m：functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)調(diào)用數(shù)值積分函數(shù)quad求定積分：[S,n]=quad(‘fesin’,0,3*pi)或[S,n]=quadl(‘fesin’,0,3*pi)或[S,n]=quad(＠fesin,0,3*pi)S=0.9008n=77quad函數(shù)采用自適應變步長方法求取定積分quadl函數(shù)則采用Lobbato方法求取定積分，計算精度和速度要高于quad函數(shù)方法。利用quadgk函數(shù)求解廣義積分問題。2．被積函數(shù)由一個表格定義在MATLAB中，對由表格形式定義的函數(shù)關系的求定積分問題用trapz(X,Y)函數(shù)。其中向量X,Y定義函數(shù)關系Y=f(X)?！纠?4】用trapz函數(shù)計算定積分命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函數(shù)關系數(shù)據(jù)向量trapz(X,Y)ans=0.28583、二重積分的數(shù)值求解使用MATLAB提供的dblquad函數(shù)就可以直接求出二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為：I=dblquad('f(x,y)',a,b,c,d,tol,trace)該函數(shù)求f(x,y)在[a,b]×[c,d]區(qū)域上的二重定積分。參數(shù)tol，trace的用法與函數(shù)quad完全相同?！纠?5】計算二重定積分解：(1)建立一個函數(shù)文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于統(tǒng)計被積函數(shù)的調(diào)用次數(shù)f=exp(-x.^2/2).*sin(x.^2+y);(2)調(diào)用dblquad函數(shù)globalki;ki=0;I=dblquad('fxy',-2,2,-1,1)I=1.57454、一般二重積分的數(shù)值求解MATLAB中沒有提供求解一般的雙重積分問題的函數(shù)在Mathworks公司的網(wǎng)站上有一個數(shù)值積分工具箱NIT（NumericalIntegrationToolbox）可以解決此問題。利用工具箱中的函數(shù)quad2dggen()可以直接求解此類問題。J=quad2dggen(Fun,Flower,Fupper,xm,xM,tol)默認tol=10-3，F(xiàn)un描述二元被積函數(shù)，需要注意是：該函數(shù)指定的積分順序是先x后y，若要求先y后x的積分，需要交換被積函數(shù)中x和y的順序?！纠?6】計算二重定積分解：fh=@(x)sqrt(1-x.^2/2);fl=@(x)-sqrt(1-x.^2/2);f=@(y,x)exp(-x.^2/2).*sin(x.^2+y);y=quad2dggen(f,fl,fh,-1/2,1,eps)y=0.41195、三重及多重積分的數(shù)值求解長方體區(qū)域的三重積分問題的函數(shù)采用triplequad()函數(shù)可以求解此類問題。J=triplequad(Fun,xm,xM,ym,yM,zm,zM,tol,@quadl)默認tol=10-6。NIT工具箱也可以解決多重超長方體邊界的定積分問題。使用quadndg()函數(shù)。該問題的求解語句為：J=quadndg(Fun[x1m,x2m,…,xpm],[x1M,x2M,…,xpM],tol)【例17】計算三重定積分解：f=@(x,y,z)4*x.*z.*exp(-x.*x.*y-z.*z);triplequad(f,0,2,0,pi,0,pi,1e-7,@equadl)ans=3.1081該工具箱中還有單重積分函數(shù)quadg()的調(diào)用格式同quad()，其效率也高于quadl()。(二)數(shù)值微分1、數(shù)值差分與差商高等數(shù)學關心的是導函數(shù)的形式和性質(zhì)，而數(shù)值分析關心的問題是怎樣計算導函數(shù)：

f'(x)=g(x)

在一串離散點X=(x1,x2,…,xn)的近似值G=(g1,g2,…,gn)以及所計算的近似值有多大誤差。引進記號：稱分別為函數(shù)在x點處以h(h>0)為步長的向前差分、向后差分和中心差分。稱分別為函數(shù)在x點處以h(h>0)為步長的向前差商、向后差商和中心差商。當步長h(h>0)充分小時，函數(shù)f在點x的微分接近于函數(shù)在該點的任意種差分。f在點x的導數(shù)接近于函數(shù)在該點的任意種差商。2、數(shù)值微分的實現(xiàn)在MATLAB中，沒有直接提供求數(shù)值導數(shù)的函數(shù)，只有計算向前差分的函數(shù)diff，其調(diào)用格式為：

DX=diff(X)：計算X的向前差分計算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。

DX=diff(X,n)：計算向量X的向前n階差分計算X的n階向前差分。例如，diff(X,2)=diff(diff(X))。DX=diff(A,n,dim)計算矩陣A的n階差分，dim=1時(缺省狀態(tài))，按列計算差分；dim=2，按行計算差分?！纠?8】設用不同的方法求函數(shù)f(x)的數(shù)值導數(shù)，并在同一個坐標系中做出f'(x)的圖像。可以考慮用三種方法：1、用一個5次多項式p(x)擬合函數(shù)f(x)，并對p(x)求一般意義下的導數(shù)dp(x)，求出dp(x)在假設點的值；2、直接求f(x)在假設點的數(shù)值導數(shù)；3、求出g(x)=f‘(x)的表達式，然后求f’(x)在假設點的數(shù)值導數(shù)。程序如下：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);%1。用5次多項式p擬合f(x)dp=polyder(p);%對擬合多項式p求導數(shù)dpdpx=polyval(dp,x);%求dp在假設點的函數(shù)值dx=diff(f([x,3.01]))/0.01;%2。直接對f(x)求數(shù)值導數(shù)gx=g(x);%3。求函數(shù)f的導函數(shù)g在假設點的導數(shù)plot(x,dpx,x,dx,'.',x,gx,'-');%作圖六、方程組求解(一)線性方程組求解1、直接求解法利用左除運算符的直接解法對于線性方程組Ax=b，可以利用左除運算符“\”求解：x=A\b2、利用矩陣的分解求解線性方程組矩陣分解是指根據(jù)一定的原理用某種算法將一個矩陣分解成若干個矩陣的乘積。常見的矩陣分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇異分解等?！纠?9】用直接求解法求解下列線性方程組命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b(1)LU分解矩陣的LU分解就是將一個矩陣表示為一個下三角矩陣L和一個上三角矩陣U的乘積形式。線性代數(shù)中已經(jīng)證明，只要方陣A是非奇異的，LU分解總是可以進行的。

MATLAB提供的lu函數(shù)用于對矩陣進行LU分解，其調(diào)用格式為：

[L,U]=lu(A)：產(chǎn)生一個上三角陣U和一個變換形式的下三角陣L(行交換)，使之滿足A=LU。注意，這里的矩陣A必須是方陣。

[L,U,P]=lu(A)：產(chǎn)生一個上三角陣U和一個下三角陣L以及一個置換矩陣P，使之滿足PA=LU。當然矩陣A同樣必須是方陣。實現(xiàn)LU分解后，線性方程組Ax=b的解x=U\(L\b)或x=U\(L\P＊b)，這樣可以大大提高運算速度。

(2)QR分解對矩陣A進行QR分解，就是把A分解為一個正交矩陣Q和一個上三角矩陣R的乘積形式。QR分解只能對方陣進行。MATLAB的qr函數(shù)可用于對矩陣進行QR分解，其調(diào)用格式為：

[Q,R]=qr(A)：產(chǎn)生一個一個正交矩陣Q和一個上三角矩陣R，使之滿足A=QR。

[Q,R,E]=qr(A)：產(chǎn)生一個一個正交矩陣Q、一個上三角矩陣R以及一個置換矩陣E，使之滿足AE=QR。實現(xiàn)QR分解后，線性方程組Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。

(3)Cholesky分解如果矩陣A是對稱正定的，則Cholesky分解將矩陣A分解成一個下三角矩陣和上三角矩陣的乘積。設上三角矩陣為R，則下三角矩陣為其轉(zhuǎn)置，即A=R'R。MATLABchol函數(shù)用于對矩陣A進行Cholesky分解，其調(diào)用格式為：

R=chol(A)：產(chǎn)生一個上三角陣R，使R'R=A。若X為非對稱正定，則輸出一個出錯信息。

[R,p]=chol(A)：這個命令格式將不輸出出錯信息。當A為對稱正定的，則p=0，R與上述格式得到的結(jié)果相同；否則p為一個正整數(shù)。如果A為滿秩矩陣，則R為一個階數(shù)為q=p-1的上三角陣，且滿足R'R=A(1:q,1:q)。實現(xiàn)Cholesky分解后，線性方程組Ax=b變成R'Rx=b，所以x=R\(R'\b)。(二)非線性方程數(shù)值求解1、單變量非線性方程求解

在MATLAB中提供了一個fzero函數(shù)，可以用來求單變量非線性方程的根。該函數(shù)的調(diào)用格式為：

z=fzero('fname',x0,tol,trace)

其中fname是待求根的函數(shù)文件名，x0為搜索的起點。一個函數(shù)可能有多個根，但fzero函數(shù)只給出離x0最近的那個根。tol控制結(jié)果的相對精度，缺省時取tol=eps，trace指定迭代信息是否在運算中顯示，為1時顯示，為0時不顯示，缺省時取trace=0?！纠?0】求f(x)=x-10x+2=0在x0=0.5附近的根。步驟如下：

(1)建立函數(shù)文件funx.m。

functionfx=funx(x)

fx=x-10.^x+2;(2)調(diào)用fzero函數(shù)求根。

z=fzero(‘funx’,0.5)或z=fzero(＠funx,0.5)z=0.37582、非線性方程組的求解

對于非線性方程組F(X)=0，用fsolve函數(shù)求其數(shù)值解。fsolve函數(shù)的調(diào)用格式為：

X=fsolve('fun',X0,option)

其中X為返回的解，fun是用于定義需求解的非線性方程組的函數(shù)文件名，X0是求根過程的初值，option為最優(yōu)化工具箱的選項設定。最優(yōu)化工具箱提供了20多個選項，用戶可以使用optimset命令將它們顯示出來。如果想改變其中某個選項，則可以調(diào)用optimset()函數(shù)來完成。例如，‘Display’選項決定函數(shù)調(diào)用時中間結(jié)果的顯示方式，其中‘off’為不顯示，‘iter’表示每步都顯示，‘final’只顯示最終結(jié)果。optimset('Display','off')將設定Display選項為‘off’?！纠?1】求下列非線性方程組在(0.5,0.5)附近的數(shù)值解。

(1)建立函數(shù)文件myfun.m。functionq=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y);(2)在給定的初值x0=0.5,y0=0.5下，調(diào)用fsolve函數(shù)求方程的根。x=fsolve('myfun',[0.5,0.5]',optimset('Display','off'))x=0.63540.3734

將求得的解代回原方程，可以檢驗結(jié)果是否正確，命令如下：q=myfun(x)q=1.0e-009*0.23750.2957

可見得到了較高精度的結(jié)果。(三)

微分方程

MATLAB能夠求解的微分方程類型包括:常微分方程初值問題常微分方程邊值問題時滯微分方程初值問題偏微分方程任何高階常微分方程都可以變換成一階微分方程，即表示成右函數(shù)形式，這是利用龍格-庫塔法求解微分方程的前提。

1．常微分方程初值問題的求解格式：[T,Y]=ode45(‘odefun’,tspan,y0)功能：返回由文件‘odefun’所定義的具有初始條件為y0、時間t變化范圍為[t0,tfinal]的微分方程y‘=f(t,y)的解，其中tspan=[t0,tfinal]。向量T中的每一列對應著矩陣Y的每一列。ode45：此方法被推薦為首選方法ode23：這是一個比ode45低階的方法ode113：用于更高階或大的標量計算ode23t：用于解決難度適中的問題ode23s：用于解決難度較大的微分方程組ode15s：與ode23相同，但要求的精度更高ode23tb：用于解決難度較大的問題

【例22】設有初值問題

求其數(shù)值解，并與精確解比較(精確解為)。

(1)建立函數(shù)文件funt.m。functionyp=funt(t,y)yp=(y^2-t-2)/4/(t+1);(2)求解微分方程。t0=0;tf=10;y0=2;[t,y]=ode23('funt',[t0,tf],y0);%求數(shù)值解y1=sqrt(t+1)+1;%求精確解plot(t,y,'b.',t,y1,'r-');%通過圖形來比較2．常微分方程邊值條件的求解對于如下的微分方程：

用于邊界條件的常微分方程求解問題:函數(shù)bvp4c格式：sol=bvp4c(‘odefun’,bcfun,solinit)其中‘odefun’為常微分方程函數(shù)，bcfun為邊界條件函數(shù)，solinit為求解的初始值。輸出sol是一個結(jié)構(gòu)，它有4個屬性?！il.x

bvp4v選擇的網(wǎng)格節(jié)點·sil.y

網(wǎng)格點sil.x處y(x)的近似值·sil.yp

網(wǎng)格點sil.x處y’(x)的近似值·sil.parameters

未知參數(shù)的值。如果存在未知參數(shù)，則求解函數(shù)會自動求得?！纠?3】求某微分方程初值問題在[1，3]區(qū)間內(nèi)的數(shù)值解，并將結(jié)果與解析解進行比較。先建立一個該函數(shù)的m文件fxy1.m：functionf=f(x,y)f=-2.*y./x+4*x%注意使用點運算符return再輸入命令：[X,Y]=ode45('fxy1',[1,3],2);X'%顯示自變量的一組采樣點Y'%顯示求解函數(shù)與采樣點對應的一組數(shù)值解(X.^2+1./X.^2)'%顯示求解函數(shù)與采樣點對應的一組解析解3．常微分方程的解析解(即符號計算解微分方程)該函數(shù)的具體調(diào)用格式為r=dsolve(‘eq1,eq2,...’,‘cond1,cond2,...’,‘v’)r=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')其中eq1、eq2等表示待求解的方程，默認的自變量為x，也可以自己規(guī)定如v，y等。微分方程中用D表示各階導數(shù)項，如Dy表示或；如果在D后面帶有數(shù)字，則表示多階導數(shù)，如D2y表示或。cond1、cond2等表示初始值，通常表示為y(a)=b或者Dy(a)=b。如果不指定初始值，或者初始值方程的個數(shù)少于因變量的個數(shù)，則最后得到的結(jié)果中會包含常數(shù)項，表示為C1、C2等。dsolve

函數(shù)最多接受12個輸入?yún)?shù)。例：某一階微分方程dsolve('Dx=y','Dy=x','x(0)=0','y(0)=1')ans=x(t)=sin(t),y(t)=cos(t)例：某二階微分方程dsolve('D2y=-a^2*y','y(0)=1','Dy(pi/a)=0')ans=cos(a*x)例y=dsolve('D2y+2*Dy+2*y=0','y(0)=1','Dy(0)=0')ans=exp(-x)*cos(x)+exp(-x)*sin(x)ezplot(y)——方程解y(t)的時間曲線圖求該方程的解MATLAB可以求解一般的偏微分方程，也可以利用偏微分方程工具箱（PDE

Toolbox）中給出的相當函數(shù)求解一些偏微分方程。1、偏微分方程組的求解考慮如下的偏微分方程：(四)

偏微分方程該偏微分方程可以編寫下面的函數(shù)表示，其入口為：其中pdefun為函數(shù)名。由給定的輸入變量即可計算出c,f,s這三個函數(shù)。邊界條件可以用以下函數(shù)描述：邊值函數(shù)可以編寫為一個MATLAB函數(shù)：初始條件函數(shù)可以表示為：即：MATLAB提供了pdepe()函數(shù)來求解該問題的數(shù)值解。其基本調(diào)用格式為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)其中有t0<t<tf

，a<x<b,，m=0、1或2。如果m>0，則必須a>0?！纠?4】試求解下面的偏微分方程。其中：初始條件：邊界條件：原方程改寫為：由此可見m=0，且：描述偏微分方程的函數(shù)：function[c,f,s]=c7mpde(x,t,u,du)c=[1;1];y=u(1)-u(2);F=exp(5.73*y)-exp(-11.46*y);s=F*[-1;1];f=[0.024*du(1);0.17*du(2)];寫出邊值方程：左邊界：右邊界：描述邊界條件的MATLAB函數(shù)：function[pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];描述初值的MATLAB函數(shù)：u0=@(x)[1;0];求解此微分方程，并得出解.x=0:.05:1;t=0:.05:2;m=0;sol=pdepe(m,@c7mpde,u0,@c7mpbc,x,t);surf(x,t