第3章 MATLAB數(shù)值運(yùn)算

MATLAB基礎(chǔ)及其應(yīng)用教程第3章:MATLAB數(shù)值運(yùn)算第3章提綱3.1多項(xiàng)式3.2插值和擬合3.3數(shù)值微積分3.4線性方程組解3.5稀疏矩陣3.6常微分方程的數(shù)值解①polyadd(x,y)②conv(a,b)③[q,r]=deconv(a,b)④poly(r)⑤

polyval(p,x)⑥

poly2sym(p)⑦roots(a)①

interp1(x,y,xi,method)②polyfit(x,y,n)①

diff②cumsum③

trapz④

quad⑤

quadl①

rref②jacobi③

gseidel④

sor

3.1多項(xiàng)式3.1.1多項(xiàng)式的表示和創(chuàng)建

在高等數(shù)學(xué)中，多項(xiàng)式一般可表示為以下形式：在MATLAB中多項(xiàng)式用一個(gè)行向量表示，向量中的元素為該多項(xiàng)式的系數(shù)，按照降序排列。如多項(xiàng)式可以表示為向量p=[9743]。若多項(xiàng)式某些系數(shù)為0，則在向量對(duì)應(yīng)位置補(bǔ)0。3.1.2多項(xiàng)式的四則運(yùn)算由于多項(xiàng)式是利用向量來(lái)表示，多項(xiàng)式的四則運(yùn)算可以轉(zhuǎn)化為向量的運(yùn)算。多項(xiàng)式的加減為對(duì)應(yīng)項(xiàng)系數(shù)的加減，因此可以通過(guò)向量的加減來(lái)實(shí)現(xiàn)。多項(xiàng)式的乘法實(shí)際上是多項(xiàng)式系數(shù)向量之間的卷積運(yùn)算，可以通過(guò)MATLAB中的卷積函數(shù)conv來(lái)完成。多項(xiàng)式的除法為乘法的逆運(yùn)算，可以通過(guò)反卷積函數(shù)deconv來(lái)實(shí)現(xiàn)。1、多項(xiàng)式相加、減（polyadd）>>a=[1234];>>b=[14916];>>c=polyadd(a,b)例1：完成兩個(gè)同階次多項(xiàng)式相加：例2：完成兩個(gè)同階次多項(xiàng)式相減：>>a=[1234];>>b=[14916];>>c=polyadd(a,-b)2、多項(xiàng)式乘法（conv）例3：完成兩個(gè)同階次多項(xiàng)式相乘：>>a=[1234];>>b=[14916];>>e=conv(a,b)3、多項(xiàng)式除法（deconv）例4：求>>a=[162050758464];>>b=[14916];>>[f,r]=deconv(a,b)例5：求

的“商”及“余”多項(xiàng)式。p1=conv([1,0,2],conv([1,4],[1,1])); p2=[1011]; [q,r]=deconv(p1,p2);3.1.3多項(xiàng)式求值、求根運(yùn)算和多項(xiàng)式構(gòu)造>>p=[12-12-17];>>z=polyval(p,3)例6：求函數(shù)在處的值例7：估計(jì)矩陣多項(xiàng)式在已知矩陣處的值，其中X=[121;-102;412]>>X=[121;-102;412];>>P=[1,0,-2,-1];>>Y=polyvalm(P,X)1、多項(xiàng)式求值（polyval、polyvalm）2、多項(xiàng)式求根（roots）例8：求代數(shù)方程的根>>p=[1-36546-453622449-67284118124-10958440320];>>roots(p)3、多項(xiàng)式構(gòu)造（poly2sym、poly）例9：構(gòu)造多項(xiàng)式>>T=[13-15-29];>>poly2sym(T)例10：用多項(xiàng)式的根構(gòu)造多項(xiàng)式>>r=roots(p)；>>poly(r)多項(xiàng)式函數(shù)3.2插值和擬合插值：根據(jù)給定有限個(gè)已知數(shù)據(jù)（樣本），構(gòu)造一簡(jiǎn)單函數(shù)，使得函數(shù)通過(guò)全部的數(shù)據(jù)點(diǎn)。擬合：根據(jù)給定的數(shù)據(jù)，選取適當(dāng)階數(shù)的多項(xiàng)式逼近給定的數(shù)據(jù)，不要求多項(xiàng)式通過(guò)全部的給定數(shù)據(jù)。區(qū)別：擬合要找出一個(gè)曲線函數(shù)，插值僅求出插值數(shù)值即可。根據(jù)實(shí)驗(yàn)數(shù)據(jù)來(lái)確定函數(shù)的方法1、多項(xiàng)式插值（interp1）在MATLAB中，一維插值由函數(shù)interp1實(shí)現(xiàn)。該函數(shù)的調(diào)用格式為 yi=interp1(x,y,xi,method)x、y：采用數(shù)據(jù)的x坐標(biāo)和y坐標(biāo)xi：待插值的位置method：采用的插值方法該語(yǔ)句返回函數(shù)在點(diǎn)xi處的插值結(jié)果。該語(yǔ)句中的參數(shù)method可以選擇的內(nèi)容如表所示。>>x=0:10;y=cos(x);>>xi=0:.25:10;>>y0=cos(xi);

%精確值>>y1=interp1(x,y,xi);

%線性插值結(jié)果>>y2=interp1(x,y,xi,'cubic');

%三次方程式插值結(jié)果>>y3=interp1(x,y,xi,'spline');

%三次樣條式插值結(jié)果>>plot(xi,y0,'o',xi,y1,xi,y2,'-.',xi,y3)

%畫(huà)出函數(shù)計(jì)算值和3種插值比較圖例11：取余弦曲線上11個(gè)點(diǎn)的自變量和函數(shù)值點(diǎn)作為已知數(shù)據(jù)，再選取41個(gè)自變量點(diǎn)，分別用3種插值方法計(jì)算確定插值函數(shù)的值。3種插值方法比較圖2、多項(xiàng)式擬合（polyfit）>>x=[-2.8-10.22.15.26.8];>>y=[3.14.62.31.22.3-1.1];>>p3=polyfit(x,y,3);>>p4=polyfit(x,y,4);>>p5=polyfit(x,y,5);>>xcurve=-3.5:0.1:7.2;>>p3curve=polyval(p3,xcurve);>>p4curve=polyval(p4,xcurve);>>p5curve=polyval(p5,xcurve);>>plot(xcurve,p3curve,'--',xcurve,p4curve,'-.',xcurve,p5curve,'-',x,y,'*')>>gridon例12：對(duì)向量X=[-2.8-10.22.15.26.8]和Y=[3.14.62.31.22.3-1.1]分別進(jìn)行3、4、5階的多項(xiàng)式擬合。3～5階多項(xiàng)式擬合曲線圖插值擬合函數(shù)3.3.1

微分和差分（diff）返回x對(duì)預(yù)設(shè)獨(dú)立變量的一次微分值；返回x對(duì)獨(dú)立變量t的一次微分值；返回x對(duì)預(yù)設(shè)獨(dú)立變量的n次微分值；返回x對(duì)獨(dú)立變量t的n次微分值；3.3數(shù)值微積分>>S1='6*x^3-4*x^2+b*x-5';>>S2='sin(a)';>>S3='(1-t^3)/(1+t^4)';>>diff(S1);>>diff(S1,2);>>diff(S1,’b’);>>diff(S2);>>diff(S3);例13：計(jì)算方程式的微分：，，3.3.2數(shù)值積分（cumsum、trapz、quad、quadl）矩形法數(shù)值積分>>x=linspace(0,pi,100);>>y=sin(x);>>T=cumsum(y)*pi/(100-1);>>I=T(100)例14：利用矩形法計(jì)算積分x為向量，則返回一個(gè)向量，其第i個(gè)元素為向量x的前i個(gè)元素之和；x為矩陣，則返回一個(gè)同維矩陣，矩陣中包含x各列的累積和；矩形積分公式為，利用MATLAB為cumsum（x）*h，h為子區(qū)間步長(zhǎng)，cumsum（x）為梯形法數(shù)值積分>>x=linspace(0,pi,100);>>y=sin(x);>>t=trapz(x,y)辛普森數(shù)值積分>>q=quad('sin',0,pi)>>quad('1./(x.^3-2.*x-5)',0,2)>>F='1./(x.^3-2.*x-5)';>>quad(F,0,2):采用低階的自適應(yīng)遞歸Simpson方法求積分q=quad('fun',a,b)，計(jì)算函數(shù)fun在區(qū)間[a,b]上的積分，其精確度為1e-6。q=quad('fun',a,b,tol)，指定允許誤差，指定的誤差tol需大于1e-6。q=quad('fun',a,b,tol,trace)，跟蹤迭代過(guò)程，分別為計(jì)算函數(shù)值的次數(shù)、當(dāng)前積分區(qū)間的左邊界、步長(zhǎng)和該區(qū)間內(nèi)的積分值。科茨數(shù)值積分:采用高階自適應(yīng)Lobatto方法求積分>>z=quadl('exp(-x.^2)',-1,1)>>z=quadl('1./(x.^3-2*x-5)',0,2)q=quadl('fun',a,b)q=quadl('fun',a,b,tol)q=quadl('fun',a,b,trace)MATLAB中二重積分和三重積分分別由函數(shù)dblquad()和函數(shù)triplequad()來(lái)實(shí)現(xiàn)。首先介紹函數(shù)dblquad()，該函數(shù)的基本格式如下：q=dblquad(fun,xmin,xmax,ymin,ymax)，函數(shù)的參數(shù)分別為函數(shù)句柄、兩個(gè)自變量的積分限，返回積分結(jié)果。q=dblquad(fun,xmin,xmax,ymin,ymax,tol)，指定積分結(jié)果的精度。triplequad()函數(shù)的調(diào)用格式和dblquad()基本相同，在調(diào)用triplequad()函數(shù)時(shí)，需要六個(gè)參數(shù)指定積分限。

二重積分和三重積分3.4線性方程組的數(shù)值解1、直接法矩陣相除法例14：求解線性方程組>>a=[1/21/31;15/33;24/35];>>b=[1;3;2];>>c=a\b例15：求解線性方程組>>a=[1-11-1;1-1-11;1-1-22];>>b=[1;0;-0.5];>>c=a\b例16：求解線性方程組>>a=[1/21/31;15/33;24/35;12/31];>>b=[1;3;2;2];>>c=a\b消去法例17：求解線性方程組>>a=[1-11-11;1-1-110;1-1-22-0.5];>>rref(a)2、迭代法Jacobi迭代法functiontx=jacobi(A,b,imax,x0,tol)del=10^-10;tx=[x0];n=length(x0);fori=1:ndg=A(i,i);ifabs(dg)<deldisp('diagonalelementistoosmall');returnendendfork=1:imaxfori=1:nsm=b(i);forj=1:nifj~=ism=sm-A(i,j)*x0(j);endendx(i)=sm/A(i,i);endtx=[tx;x];ifnorm(x-x0)<tolreturnelsex0=x;endend例18：求解線性方程組>>A=[10-120;-111-13;2-110-1;03-18];>>b=[625-1115]';>>tol=1.0*10^-6;>>imax=10;>>x0=zeros(1,4);>>tx=jacobi(A,b,imax,x0,tol);>>forj=1:size(tx,1)fprintf('%4d%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4))endGauss-Seidel迭代法functiontx=gseidel(A,b,imax,x0,tol)del=10^-10;tx=[x0];n=length(x0);fori=1:ndg=A(i,i);ifabs(dg)<deldisp('diagonalelementistoosmall');returnendendfork=1:imax

x=x0;fori=1:nsm=b(i);forj=1:nifj~=i

sm=sm-A(i,j)*x(j);endendx(i)=sm/A(i,i);endtx=[tx;x];ifnorm(x-x0)<tolreturnelsex0=x;endend例18：求解線性方程組>>A=[10-120;-111-13;2-110-1;03-18];>>b=[625-1115]';>>tol=1.0*10^-6;>>imax=10;>>x0=zeros(1,4);>>tx=gseidel(A,b,imax,x0,tol);>>forj=1:size(tx,1)fprintf('%4d%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4))endSOR迭代法functiontx=gseidel(A,b,imax,x0,tol)del=10^-10;tx=[x0];n=length(x0);fori=1:ndg=A(i,i);ifabs(dg)<deldisp('diagonalelementistoosmall');returnendendfork=1:imax

x=x0;fori=1:nsm=b(i);forj=1:nifj~=i

sm=sm-A(i,j)*x(j);endendx(i)=sm/A(i,i);

x(i)=w*x(i)+(1-w)*x0(i);endtx=[tx;x];ifnorm(x-x0)<tolreturnelsex0=x;endend例19：求解線性方程組>>A=[10-120;-111-13;2-110-1;03-18];>>b=[625-1115]';>>tol=1.0*10^-6;>>imax=10;>>x0=zeros(1,4);>>w=1.02;>>tx=sor(A,b,imax,x0,tol);>>forj=1:size(tx,1)fprintf('%4d%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4))end3.5稀疏矩陣1、稀疏矩陣的創(chuàng)建例20：創(chuàng)建4×5階矩陣的稀疏矩陣>>i=[124];>>j=[135];>>s=[678];>>A=sparse(i,j,s)例21：創(chuàng)建矩陣A的對(duì)角稀疏矩陣>>A=[050100000601103007012140080025009];>>[B,d]=spdiags(A)>>s=spdiags(B,d,A)2、稀疏矩陣的存儲(chǔ)例22：比較一個(gè)10×10單位矩陣A及對(duì)應(yīng)的稀疏矩陣B=sparse(A)所需的存儲(chǔ)量。>>

A=eye(10);>>B=sparse(A);>>whosAB稀疏矩陣的存儲(chǔ)格式：（1）第一個(gè)數(shù)組以浮點(diǎn)格式存放數(shù)組中所有的非零元素；（2）第二個(gè)數(shù)組存放非零元素對(duì)應(yīng)的行號(hào)，行號(hào)為整數(shù)；（3）第三個(gè)數(shù)組存放n+1個(gè)整形指針。3.6常微分方程的數(shù)值解常微分方程求解思路：對(duì)求解區(qū)間進(jìn)行剖分，然后把微分方程離散成在節(jié)點(diǎn)上的近似公式或近似方程，最后結(jié)合定解條件求出近似解。常微分方程求解方法：泰勒展開(kāi)法、數(shù)值積分法（Runge-Kutta）、數(shù)值微分法（Adams）1、泰勒（Taylor）公式2、歐拉（Euler）法前向歐拉公式后向歐拉公