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

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

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

%計(jì)算分子多項(xiàng)式 p2=[1011];

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

求多項(xiàng)式的微分polyder(p)

polyder(a,b):求多項(xiàng)式a,b乘積的微分[p,q]=polyder(a,b):求多項(xiàng)式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ù)多項(xiàng)式求值polyval函數(shù)用來求代數(shù)多項(xiàng)式的值，其調(diào)用格式為：Y=polyval(P,x)若x為一數(shù)值，則求多項(xiàng)式在該點(diǎn)的值；若x為向量或矩陣，則對(duì)向量或矩陣中的每個(gè)元素求其多項(xiàng)式的值?！纠?】已知多項(xiàng)式x4+8x3-10，分別取x=1.2和一個(gè)2×3矩陣為自變量計(jì)算該多項(xiàng)式的值。A=[1,8,0,0,-10];%4次多項(xiàng)式系數(shù)x=1.2;%取自變量為一數(shù)值y1=polyval(A,x)y1=5.8976x=[-1,1.2,-1.4;2,-1.8,1.6];%給出一個(gè)矩陣xy2=polyval(A,x)y2=-17.00005.8976-28.110470.0000-46.158429.32164．矩陣多項(xiàng)式求值polyvalm函數(shù)用來求矩陣多項(xiàng)式的值，其調(diào)用格式與polyval相同，但含義不同。polyvalm函數(shù)要求x為方陣,它以方陣為自變量求多項(xiàng)式的值。設(shè)A為方陣，P代表多項(xiàng)式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】仍以多項(xiàng)式x4+8x3-10為例，取一個(gè)2×2矩陣為自變量分別用polyval和polyvalm計(jì)算該多項(xiàng)式的值。A=[1,8,0,0,-10];%多項(xiàng)式系數(shù)x=[-1,1.2;2,-1.8];%給出一個(gè)矩陣xy1=polyval(A,x)%計(jì)算代數(shù)多項(xiàng)式的值y1=-17.00005.897670.0000-46.1584y2=polyvalm(A,x)%計(jì)算矩陣多項(xiàng)式的值y2=-60.584050.649684.4160-94.35045多項(xiàng)式求根n次多項(xiàng)式具有n個(gè)根，當(dāng)然這些根可能是實(shí)根，也可能含有若干對(duì)共軛復(fù)根。MATLAB提供的roots函數(shù)用于求多項(xiàng)式的全部根，其調(diào)用格式為：x=roots(P)其中P為多項(xiàng)式的系數(shù)向量，求得的根賦給向量x，即x(1),x(2),…,x(n)分別代表多項(xiàng)式的n個(gè)根。【例5】求多項(xiàng)式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)計(jì)處理(一)最大值和最小值

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

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

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

在MATLAB中，提供了計(jì)算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std。

std(X)：返回一個(gè)標(biāo)準(zhǔn)方差。

std(A)：返回一個(gè)行向量，每個(gè)元素是矩陣A各列或各行的標(biāo)準(zhǔn)方差。

Y=std(A,flag,dim)

缺省flag=0，dim=1當(dāng)dim=1時(shí)，求各列元素的標(biāo)準(zhǔn)方差；當(dāng)dim=2時(shí)，則求各行元素的標(biāo)準(zhǔn)方差。當(dāng)flag=0時(shí)，按S1所列公式計(jì)算標(biāo)準(zhǔn)方差；當(dāng)flag=1時(shí)，按S2所列公式計(jì)算標(biāo)準(zhǔn)方差。

【例6】對(duì)二維矩陣x，從不同維方向求出其標(biāo)準(zhǔn)方差。命令如下：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．相關(guān)系數(shù)

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

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

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

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

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

一維插值函數(shù)根據(jù)X,Y的值，計(jì)算函數(shù)在X1處的值。X,Y是兩個(gè)等長的已知向量，分別描述采樣點(diǎn)和樣本值，X1是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)，Y1是與X1等長的插值結(jié)果。method是插值方法，例如：‘linear’線性插值‘nearest’最近點(diǎn)插值‘cubic’三次多項(xiàng)式插值‘spline’三次樣條插值或Y1=spline(X,Y,X1)計(jì)算三次樣條插值的專門函數(shù)三、數(shù)據(jù)插值注意：X1的取值范圍不能超出X的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。【例9】某觀測站測得某日6:00時(shí)至18:00時(shí)之間每隔2小時(shí)的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時(shí)之間每隔2小時(shí)各點(diǎn)的近似溫度(℃)。解：設(shè)時(shí)間變量h為一行向量，溫度變量t為一個(gè)兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲(chǔ)存室外溫度。命令如下：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次樣條插值計(jì)算(二)二維數(shù)據(jù)插值

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

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

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

%在(0.5,0.5)點(diǎn)插值【例11】某實(shí)驗(yàn)對(duì)一根長10米的鋼軌進(jìn)行熱源的溫度傳播測試。用x表示測量點(diǎn)0:2.5:10(米)，用h表示測量時(shí)間0:30:60(秒)，用T表示測試所得各點(diǎn)的溫度(℃)。試用線性插值求出在一分鐘內(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ù)來求得最小二乘擬合多項(xiàng)式的系數(shù)，再用polyval函數(shù)按所得的多項(xiàng)式計(jì)算所給出的點(diǎn)上的函數(shù)近似值。

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

polyval函數(shù)的功能是按多項(xiàng)式的系數(shù)計(jì)算x點(diǎn)多項(xiàng)式的值。四、曲線擬合【例12】用一個(gè)三次多項(xiàng)式在區(qū)間[0,2π]內(nèi)逼近函數(shù)sin(x)。X=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,3)%得到3次多項(xiàng)式的系數(shù)和誤差五、數(shù)值積分與微分(一)數(shù)值積分1、變步長辛普生法基本思想都是將整個(gè)積分區(qū)間[a,b]分成n個(gè)子區(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用來控制積分精度，缺省時(shí)取tol=10-6。trace控制是否展現(xiàn)積分過程，若取非0則展現(xiàn)積分過程，取0則不展現(xiàn)，缺省時(shí)取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ù)采用自適應(yīng)變步長方法求取定積分quadl函數(shù)則采用Lobbato方法求取定積分，計(jì)算精度和速度要高于quad函數(shù)方法。利用quadgk函數(shù)求解廣義積分問題。2．被積函數(shù)由一個(gè)表格定義在MATLAB中，對(duì)由表格形式定義的函數(shù)關(guān)系的求定積分問題用trapz(X,Y)函數(shù)。其中向量X,Y定義函數(shù)關(guān)系Y=f(X)。【例14】用trapz函數(shù)計(jì)算定積分命令如下：X=1:0.01:2.5;Y=exp(-X);%生成函數(shù)關(guān)系數(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】計(jì)算二重定積分解：(1)建立一個(gè)函數(shù)文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于統(tǒng)計(jì)被積函數(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)站上有一個(gè)數(shù)值積分工具箱NIT（NumericalIntegrationToolbox）可以解決此問題。利用工具箱中的函數(shù)quad2dggen()可以直接求解此類問題。J=quad2dggen(Fun,Flower,Fupper,xm,xM,tol)默認(rèn)tol=10-3，F(xiàn)un描述二元被積函數(shù)，需要注意是：該函數(shù)指定的積分順序是先x后y，若要求先y后x的積分，需要交換被積函數(shù)中x和y的順序?！纠?6】計(jì)算二重定積分解：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)默認(rèn)tol=10-6。NIT工具箱也可以解決多重超長方體邊界的定積分問題。使用quadndg()函數(shù)。該問題的求解語句為：J=quadndg(Fun[x1m,x2m,…,xpm],[x1M,x2M,…,xpM],tol)【例17】計(jì)算三重定積分解：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ù)學(xué)關(guān)心的是導(dǎo)函數(shù)的形式和性質(zhì)，而數(shù)值分析關(guān)心的問題是怎樣計(jì)算導(dǎo)函數(shù)：

f'(x)=g(x)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中fname是待求根的函數(shù)文件名，x0為搜索的起點(diǎn)。一個(gè)函數(shù)可能有多個(gè)根，但fzero函數(shù)只給出離x0最近的那個(gè)根。tol控制結(jié)果的相對(duì)精度，缺省時(shí)取tol=eps，trace指定迭代信息是否在運(yùn)算中顯示，為1時(shí)顯示，為0時(shí)不顯示，缺省時(shí)取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、非線性方程組的求解

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

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

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

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

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

微分方程

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

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

【例22】設(shè)有初值問題

求其數(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．常微分方程邊值條件的求解對(duì)于如下的微分方程：

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

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

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

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

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

函數(shù)最多接受12個(gè)輸入?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)的時(shí)間曲線圖求該方程的解MATLAB可以求解一般的偏微分方程，也可以利用偏微分方程工具箱（PDE

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

偏微分方程該偏微分方程可以編寫下面的函數(shù)表示，其入口為：其中pdefun為函數(shù)名。由給定的輸入變量即可計(jì)算出c,f,s這三個(gè)函數(shù)。邊界條件可以用以下函數(shù)描述：邊值函數(shù)可以編寫為一個(gè)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