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

MATLAB基礎(chǔ)及其應(yīng)用教程第3章:MATLAB數(shù)值運算第3章提綱3.1多項式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多項式3.1.1多項式的表示和創(chuàng)建

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

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

%畫出函數(shù)計算值和3種插值比較圖例11：取余弦曲線上11個點的自變量和函數(shù)值點作為已知數(shù)據(jù)，再選取41個自變量點，分別用3種插值方法計算確定插值函數(shù)的值。3種插值方法比較圖2、多項式擬合（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：對向量X=[-2.8-10.22.15.26.8]和Y=[3.14.62.31.22.3-1.1]分別進行3、4、5階的多項式擬合。3～5階多項式擬合曲線圖插值擬合函數(shù)3.3.1

微分和差分（diff）返回x對預設(shè)獨立變量的一次微分值；返回x對獨立變量t的一次微分值；返回x對預設(shè)獨立變量的n次微分值；返回x對獨立變量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：計算方程式的微分：，，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：利用矩形法計算積分x為向量，則返回一個向量，其第i個元素為向量x的前i個元素之和；x為矩陣，則返回一個同維矩陣，矩陣中包含x各列的累積和；矩形積分公式為，利用MATLAB為cumsum（x）*h，h為子區(qū)間步長，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)，計算函數(shù)fun在區(qū)間[a,b]上的積分，其精確度為1e-6。q=quad('fun',a,b,tol)，指定允許誤差，指定的誤差tol需大于1e-6。q=quad('fun',a,b,tol,trace)，跟蹤迭代過程，分別為計算函數(shù)值的次數(shù)、當前積分區(qū)間的左邊界、步長和該區(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()來實現(xiàn)。首先介紹函數(shù)dblquad()，該函數(shù)的基本格式如下：q=dblquad(fun,xmin,xmax,ymin,ymax)，函數(shù)的參數(shù)分別為函數(shù)句柄、兩個自變量的積分限，返回積分結(jié)果。q=dblquad(fun,xmin,xmax,ymin,ymax,tol)，指定積分結(jié)果的精度。triplequad()函數(shù)的調(diào)用格式和dblquad()基本相同，在調(diào)用triplequad()函數(shù)時，需要六個參數(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;>>ima