第8章 計(jì)算方法的MATLAB實(shí)現(xiàn)

MATLAB7.0從入門到精通哈爾濱工業(yè)大學(xué)(威海)汽車工程學(xué)院2/1/20231課程主要內(nèi)容第1章MATLAB簡介第2章數(shù)值運(yùn)算第3章單元數(shù)組和結(jié)構(gòu)第4章字符串第5章符號運(yùn)算第6章MATLAB繪圖基礎(chǔ)第7章程序設(shè)計(jì)第8章計(jì)算方法的MATLAB實(shí)現(xiàn)第9章優(yōu)化設(shè)計(jì)第10章SIMULINK仿真初探2/1/20232第8章計(jì)算方法的MATLAB實(shí)現(xiàn)

隨著計(jì)算機(jī)的迅速發(fā)展與廣泛運(yùn)用，在眾多的領(lǐng)域，科學(xué)計(jì)算方法的應(yīng)用越來越廣泛了，而MATLAB在進(jìn)行科學(xué)計(jì)算方面有著無與倫比的優(yōu)勢。本章介紹MATLAB在計(jì)算方法中的運(yùn)用。2/1/202338.1方程求根roots見多項(xiàng)式求根；roots(多項(xiàng)式系數(shù)矩陣)fzero可求解非線性方程，它的格式為：fzero(‘function’,x0)其中function為求解的方程，x0為估計(jì)的根，x0可為標(biāo)量或長度為2的向量，為向量時(shí)函數(shù)的兩端的值應(yīng)該符號相反，此時(shí)求區(qū)間上的解。只能求解x0附近的一個(gè)解。即使在某個(gè)區(qū)間內(nèi)有多個(gè)解，但是區(qū)間端點(diǎn)符號相同的話仍然出錯(cuò)。2/1/20234程序?qū)嵗?gt;>fzero('x^3-3*x-1',2)ans=1.8794>>fzero('x^3-3*x-1',[1,4])ans=1.8794>>fzero('x^3-3*x-1',[2,4])???Errorusing==>fzeroThefunctionvaluesattheintervalendpointsmustdifferinsign.2/1/20235程序?qū)嵗?gt;>fzero('x^2-3*x+2',0)ans=1.0000>>fzero('x^2-3*x+2',3)ans=2.0000>>fzero('x^2-3*x+2',[0,4])???Errorusing==>fzeroThefunctionvaluesattheintervalendpointsmustdifferinsign.2/1/202368.2線性方程組數(shù)值解法線性方程組的求解不僅在工程技術(shù)領(lǐng)域涉及到，而且在其他的許多領(lǐng)域也經(jīng)常碰到，因此這是一個(gè)應(yīng)用相當(dāng)廣泛的課題。關(guān)于線性方程組的數(shù)值解法一般分為兩類：一類是直接法，就是在沒有舍入誤差的情況下，通過有限步四則運(yùn)算求得方程組準(zhǔn)確解的方法。另一類是迭代法，就是先給定一個(gè)解的初始值，然后按一定的法則逐步求出解的各個(gè)更準(zhǔn)確的近似值的方法。2/1/202378.2.1直接解法關(guān)于線性方程組的直接解法有許多種，比如Gauss消去法、列主元消去法、平方根法等。而在MATLAB中，線性方程組的直接解法只需用符號“/”或“\”就解決問題。還可以使用逆陣函數(shù)來求解：x=inv(A)*B。2/1/20238程序?qū)嵗蠼庀铝蟹匠探M2/1/20239程序?qū)嵗?gt;>a=[12-33;-183-1;111];>>b=[15;-15;6];>>x1=a\bx1=1.00002.00003.0000>>x2=inv(a)*bx2=1232/1/2023108.2.2線性方程組求解中的變換上三角變換U=triu(x)返回矩陣x的上三角部分；U=triu(x,k)返回第k條對角線以上部分的元素。2/1/202311程序?qū)嵗?gt;>a=[12-33;-183-1;111];>>triu(a)ans=12-3303-10012/1/202312>>triu(a,1)ans=0-3300-1000>>triu(a,-1)ans=12-33-183-1011程序?qū)嵗?/1/202313下三角變換U=tril(x)返回矩陣x的下三角部分；U=tril(x,k)返回第k條對角線以上下部分的元素。2/1/202314程序?qū)嵗?gt;>a=[12-33;-183-1;111];>>tril(a)ans=1200-18301112/1/202315程序?qū)嵗?gt;>tril(a,1)ans=12-30-183-1111>>tril(a,-1)ans1/202316對角變換U=diag(x)返回矩陣x主對角線上的元素，返回結(jié)果是一列向量形式；U=diag(x,k)返回第k條對角線上的元素值。當(dāng)x為向量時(shí)生成矩陣。2/1/202317程序?qū)嵗?gt;>a=[12-33;-183-1;111];>>diag(a)ans=12312/1/202318程序?qū)嵗?gt;>a=[12-33;-183-1;111];>>diag(a,1)ans=-3-1>>diag(a,-1)ans=-1812/1/2023198.2.3迭代解法迭代解法非常適用于求解大型稀疏系數(shù)矩陣的方程組，在線性方程組常用的迭代解法主要有Jacobi迭代法、Gauss-Seidel迭代法。2/1/202320Jacobi迭代法2/1/202321Jacobi.mfunctions=jacobi(a,b,x0,eps)%jacobi迭代法皆線性方程組%a為系數(shù)矩陣，b為方程組ax=b中的右邊的矩陣b，x0為迭代初值ifnargin==3

eps=1.0e-6;elseif

nargin<3errorreturnend2/1/202322D=diag(diag(a));%求出對角矩陣D=inv(D);%求出對角矩陣的逆L=tril(a,-1);%求出嚴(yán)格下三角矩陣U=triu(a,1);%求出嚴(yán)格上三角矩陣B=-D*(L+U);f=D*b;s=B*x0+f;whilenorm(s-x0)>=epsx0=s;s=B*x0+f;endreturn2/1/202323程序?qū)嵗蒙厦婢帉懙膉acobi函數(shù)求解下列方程組2/1/202324程序?qū)嵗?gt;>a=[10-2-1;-210-1;-1-25];>>b=[31510]';>>x=jacobi(a,b,[000]',eps)x=1232/1/202325Gauss-Saidel迭代法2/1/202326gauss.mfunctions=gauss(a,b,x0,eps)%gauss-seidel迭代法皆線性方程組%a為系數(shù)矩陣，b為方程組ax=b中的右邊的矩陣b，x0為迭代初值ifnargin==3

eps=1.0e-6;elseif

nargin<3errorreturnend2/1/202327L=tril(a,-1);%求出嚴(yán)格下三角矩陣D=diag(diag(a));%求出對角矩陣U=triu(a,1);%求出嚴(yán)格上三角矩陣C=inv(D+L);B=-C*U;f=C*b;s=B*x0+f;whilenorm(s-x0)>=epsx0=s;s=B*x0+f;endreturn2/1/202328程序?qū)嵗蒙厦婢帉懙膅auss函數(shù)求解下列方程組2/1/202329程序?qū)嵗?gt;>a=[10-2-1;-22-1;-1-25];>>b=[6;10;10];>>x=gauss(a,b,[000]',eps)x=41382/1/2023308.3非線性方程組數(shù)值解法與線性方程組的求解一樣，非線性方程組的求解也是應(yīng)用很廣泛的課題。一般情況，非線性方程組的數(shù)值解法主要采用迭代法來求解。比較常用的迭代法主要有不動(dòng)點(diǎn)迭代法、Newton迭代法、擬Newton迭代法等幾種方法。2/1/202331不動(dòng)點(diǎn)迭代法2/1/202332staticiterate.mfunctions=staticiterate(x,eps)%不動(dòng)點(diǎn)迭代法解非線性方程組，x為迭代初值，eps為允許誤差ifnargin==1

eps=1.0e-6;elseif

nargin<1errorreturnend2/1/202333xx=fx(x);%第一次迭代whilenorm(xx-x)>=epsx=xx;xx=fx(x);ends=xx;return2/1/202334程序?qū)嵗貌粍?dòng)點(diǎn)迭代法求解下面的方程組2/1/202335fx.m首先編寫上述非線性方程組的M文件fx.mfunctiony=fx(x)y(1)=0.1*(x(1)*x(1)+x(2)*x(2)+8);y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8);y=[y(1)y(2)];2/1/202336程序?qū)嵗?gt;>x=staticiterate([00])x=1.00001.00002/1/202337Newton迭代法2/1/202338newtoniterate.mfunctions=newtoniterate(x,eps)%newton迭代法解非線性方程組，x為迭代初值，eps為允許誤差ifnargin==1

eps=1.0e-6;elseif

nargin<1errorreturnend2/1/202339x1=fx1(x);%非線性方程組x2=-dfx1(x);%非線性方程組導(dǎo)數(shù)x3=inv(x2);x0=x3*x1';whilenorm(x0)>=epsx=x0'+x;x1=fx1(x);x2=-dfx1(x);x3=inv(x2);x0=x3*x1';ends=x0'+x;return2/1/202340程序?qū)嵗蒙厦婢帉懙膎ewton迭代函數(shù)求解下列方程組2/1/202341fx1.m和dfx1.m首先，編寫上述非線性方程組的M文件fx1.mfunctiony=fx1(x)y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;y(2)=x(1)*x(2)*x(2)-10*x(2)+x(1)+8;y=[y(1)y(2)];然后，編寫上述非線性方程組導(dǎo)數(shù)的M文件dfx1.mfunctiony=dfx1(x)y(1)=2*x(1)-10;y(2)=2*x(2);y(3)=x(2)*x(2)+1;y(4)=2*x(1)*x(2)-10;y=[y(1)y(2);y(3)y(4)];2/1/202342程序?qū)嵗?gt;>x=newtoniterate([00])x=112/1/2023438.4插值與擬合在生產(chǎn)實(shí)踐及科學(xué)實(shí)驗(yàn)中，插值與擬合的應(yīng)用非常廣泛。下面，就對如何用MATLAB來處理插值與擬合作一介紹。2/1/2023448.4.1一維插值yi=interp1(x,y,xi)返回在插值向量xi處的函數(shù)向量yi，它是根據(jù)向量x和y插值而來。如y是矩陣，則對y每一列進(jìn)行插值，如xi中元素不在x內(nèi)，返回NaN。yi=interp1(y,xi)省略x，表示x=1:N，此時(shí)N為向量y的長度或?yàn)榫仃噛的行數(shù)。yi=interp1(x,y,xi,’method’)表示用method指定的插值方法進(jìn)行插值。2/1/202345Method可取如下的值：‘nearest’最近插值‘linear’線性插值‘spline’三次樣條插值‘cubic’三次插值Method默認(rèn)值為線性插值，上述插值要求向量x單調(diào)。2/1/202346程序?qū)嵗?gt;>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1)y1=5.40007.05007.50002/1/202347程序?qū)嵗?gt;>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'linear')y1=5.40007.05007.50002/1/202348程序?qū)嵗?gt;>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'nearest')y1=5782/1/202349程序?qū)嵗?gt;>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'spline')y1=5.55207.11147.67852/1/202350程序?qū)嵗?gt;>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'cubic')y1=5.50067.08147.54762/1/202351程序?qū)嵗?gt;>x=linspace(0,2*pi,100);>>y=sin(x);>>x0=linspace(0,2*pi,6);>>y0=sin(x0);>>x1=[0.621.93.24.345.55];>>y1=interp1(x0,y0,x1);>>plot(x,y,x0,y0,'ro',x1,y1,'+')>>title('默認(rèn)插值')2/1/2023522/1/202353程序結(jié)果>>y1y1=0.46920.7651-0.0546-0.7526-0.5549>>y10=sin(x1)y10=0.58100.9463-0.0584-0.9315-0.6692>>deltay1=y1-y10deltay1=-0.1118-0.18120.00370.17890.11432/1/202354程序?qū)嵗?gt;>x=linspace(0,2*pi,100);>>y=sin(x);>>x0=linspace(0,2*pi,6);>>y0=sin(x0);>>x1=[0.621.93.24.345.55];>>y1=interp1(x0,y0,x1,'linear');>>plot(x,y,x0,y0,'ro',x1,y1,'+')>>title('linear插值')2/1/2023552/1/202356程序結(jié)果>>y1y1=0.46920.7651-0.0546-0.7526-0.5549>>y10=sin(x1)y10=0.58100.9463-0.0584-0.9315-0.6692>>deltay1=y1-y10deltay1=-0.1118-0.18120.00370.17890.11432/1/202357程序?qū)嵗?gt;>x=linspace(0,2*pi,100);>>y=sin(x);>>x0=linspace(0,2*pi,6);>>y0=sin(x0);>>x1=[0.621.93.24.345.55];>>y1=interp1(x0,y0,x1,'nearest');>>plot(x,y,x0,y0,'ro',x1,y1,'+')>>title('nearest插值')2/1/2023582/1/202359程序結(jié)果>>y1y1=00.5878-0.5878-0.5878-0.9511>>y10=sin(x1)y10=0.58100.9463-0.0584-0.9315-0.6692>>deltay1=y1-y10deltay1=-0.5810-0.3585-0.52940.3437-0.28182/1/202360程序?qū)嵗?gt;>x=linspace(0,2*pi,100);>>y=sin(x);>>x0=linspace(0,2*pi,6);>>y0=sin(x0);>>x1=[0.621.93.24.345.55];>>y1=interp1(x0,y0,x1,'spline');>>plot(x,y,x0,y0,'ro',x1,y1,'+')>>title('spline插值')2/1/2023612/1/202362程序結(jié)果>>y1y1=0.64150.9212-0.0592-0.9073-0.7219>>y10=sin(x1)y10=0.58100.9463-0.0584-0.9315-0.6692>>deltay1=y1-y10deltay1=0.0605-0.0251-0.00080.0242-0.05272/1/202363程序?qū)嵗?gt;>x=linspace(0,2*pi,100);>>y=sin(x);>>x0=linspace(0,2*pi,6);>>y0=sin(x0);>>x1=[0.621.93.24.345.55];>>y1=interp1(x0,y0,x1,'cubic');>>plot(x,y,x0,y0,'ro',x1,y1,'+')>>title('cubic插值')2/1/2023642/1/202365程序結(jié)果>>y1y1=0.66970.8339-0.0689-0.8194-0.7563>>y10=sin(x1)y10=0.58100.9463-0.0584-0.9315-0.6692>>deltay1=y1-y10deltay1=0.0887-0.1124-0.01060.1121-0.08702/1/2023668.4.2二維插值zi=interp2(x,y,z,xi,yi)返回在插值向量xi、yi處的函數(shù)值向量，它是根據(jù)向量x、y與z插值而來，這里的x、y與z也可以是矩陣形式。如果xi、yi有元素不在x、y范圍內(nèi)，則返回NaN。zi=interp2(z,xi,yi)省略x、y，表示x=1:N，y=1:M。此處[M,N]=size(z)。zi=interp2(z,ntimes)表示在z的各點(diǎn)之間插入數(shù)據(jù)點(diǎn)來擴(kuò)展z，依次執(zhí)行ntimes次迭代，ntimes默認(rèn)為1。zi=interp2(x,y,z,xi,yi,’method’)表示用method指定的插值方法進(jìn)行插值。Method取值同上，要求x與y都單調(diào)且具有相同格式。2/1/202367程序?qū)嵗?gt;>z1=[35791110415;1326115713];>>z=interp2(z1,[13],[12])z=322/1/202368程序?qū)嵗?gt;>z1=[35791110415;1326115713];>>z=interp2(z1)或z=interp2(z1,1)z=Columns1through123.00004.00005.00006.00007.00008.00009.000010.000011.000010.500010.00007.00002.00003.00004.00004.25004.50006.00007.50009.250011.00009.25007.50006.50001.00002.00003.00002.50002.00004.00006.00008.500011.00008.00005.00006.0000Columns13through154.00009.500015.00005.50009.750014.00007.000010.000013.00002/1/202369程序?qū)嵗?gt;>[x,y,z]=peaks(10);>>[xi,yi]=meshgrid(-3:0.3:3,-3:0.3:3);>>zi=interp2(x,y,z,xi,yi,'neareat');>>meshc(xi,yi,zi)>>title('nearest插值')2/1/2023702/1/202371程序?qū)嵗?gt;>[x,y,z]=peaks(10);>>[xi,yi]=meshgrid(-3:0.3:3,-3:0.3:3);>>zi=interp2(x,y,z,xi,yi,'linear');>>meshc(xi,yi,zi)>>title('linear插值')2/1/2023722/1/202373程序?qū)嵗?gt;>[x,y,z]=peaks(10);>>[xi,yi]=meshgrid(-3:0.3:3,-3:0.3:3);>>zi=interp2(x,y,z,xi,yi,'spline');>>meshc(xi,yi,zi)>>title('spline插值')2/1/2023742/1/202375程序?qū)嵗?gt;>[x,y,z]=peaks(10);>>[xi,yi]=meshgrid(-3:0.3:3,-3:0.3:3);>>zi=interp2(x,y,z,xi,yi,'cubic');>>meshc(xi,yi,zi)>>title('cubic插值')2/1/2023762/1/202377程序?qū)嵗?gt;>[x,y]=meshgrid(-3:0.3:3,-3:0.3:3);>>z=peaks(x,y);>>meshc(x,y,z)>>title('原始圖')2/1/2023782/1/2023798.4.3三維插值vi=interp3(x,y,z,v,xi,yi,zi)返回三維函數(shù)v在插值向量xi，yi，zi處的函數(shù)向量vi，它們大小形式相同。vi=interp3(v,xi,yi,zi)省略x、y、z。同前。vi=interp3(x,y,z,v,xi,yi,zi,’method’)同前。2/1/202380程序?qū)嵗?gt;>[x,y,z,v]=flow(10);>>[xi,yi,zi]=meshgrid(0.1:0.25:10,-3:0.25:3,-3:0.25:3);>>vi=interp3(x,y,z,v,xi,yi,zi);>>slice(xi,yi,zi,vi,[69.5],2,[-2,0.2])2/1/2023812/1/2023828.4.4Lagrange插值2/1/202383lagrange.mfunctions=lagrange(x,y,x0)%lagrange插值，x與y為已知的插值點(diǎn)及其函數(shù)值，x0為需要求的插值點(diǎn)的x值nx=length(x);ny=length(y);ifnx~=nywarning('向量x與y的長度應(yīng)該相同')return;endm=length(x0);2/1/202384%按照公式，對需要求的插值點(diǎn)向量x0的每個(gè)元素進(jìn)行計(jì)算fori=1:mt=0.0;forj=1:nxu=1.0;fork=1:nxifk~=ju=u*(x0(i)-x(k))/(x(j)-x(k));endendt=t+u*y(j);end

s(m)=t;endreturn2/1/202385程序?qū)嵗?gt;>x=[12345];>>y=[246810];>>y1=lagrange(x,y,1.6)y1=3.2000>>y2=lagrange(x,y,[1.42.53.7])y2=007.40002/1/2023868.4.5Newton插值2/1/202387newton.mfunctions=newton(x,y,x0,nn)%newton插值，x與y為已知的插值點(diǎn)及其函數(shù)值%x0為需要求的插值點(diǎn)的x坐標(biāo)值。nn為newton插值多項(xiàng)式的次數(shù)nx=length(x);ny=length(y);ifnx~=nywarning('向量x與y的長度應(yīng)該相同')return;endm=length(x0);%按照公式，對需要求的插值點(diǎn)向量x0的每個(gè)元素進(jìn)行計(jì)算2/1/202388fori=1:mt=0.0;j=1;

yy=y;

kk=j;%求各級均差

while(kk<=nn)

kk=kk+1;fork=kk:nx

yy(k)=(yy(k)-yy(kk-1))/(x(k)-x(kk-1));endend%求插值結(jié)果2/1/202389

t=yy(1);fork=2:nnu=1.0;

jj=1;

while(jj<k)u=u*(x0(i)-x(jj));

jj=jj+1;endt=t+yy(k)*u;end

s(i)=t;endreturn2/1/202390程序?qū)嵗?gt;>x=[0.40.550.650.80.91.05];>>y=[0.410750.578150.696750.888111.026521.25382];>>y1=newton(x,y,0.596,4)y1=0.63192/1/2023918.4.6三次樣條插值眾所周知，使用高階多項(xiàng)式的插值常常會(huì)產(chǎn)生病態(tài)的結(jié)果，而三次樣條插值就能消除這種病態(tài)。MATLAB提供的三次樣條插值函數(shù)有spline與interp1兩個(gè)，下面重點(diǎn)來介紹spline函數(shù)。spline函數(shù)的調(diào)用格式如下：yy=spline(x,y,xx)利用三次樣條插值法尋找在插值點(diǎn)xx處的插值函數(shù)值yy，插值函數(shù)根據(jù)輸入?yún)?shù)x與y的關(guān)系得來。x與y為向量形式，而xx可以為向量形式，也可以是標(biāo)量形式。此函數(shù)的作用等同于interp1(x,y,xx,’spline’)。2/1/202392程序?qū)嵗?gt;>x=0:10;>>y=sin(x);>>xx=0:0.25:10;>>yy=spline(x,y,xx);>>plot(x,y,'o',xx,yy)2/1/2023932/1/2023948.4.7最小二乘法曲線擬合在實(shí)際工程應(yīng)用與科學(xué)實(shí)踐中，經(jīng)常要得到一條光滑的曲線，而實(shí)際卻只能測得一些分散的數(shù)據(jù)點(diǎn)。此時(shí)，就需要利用這些離散的點(diǎn)，運(yùn)用各種擬合方法來生成一條連續(xù)的曲線。在這些擬合方法中，最常用的是最小二乘曲線擬合法。已知一組數(shù)據(jù)(xi,yi)，從中找出自變量x與因變量y之間的函數(shù)關(guān)系y=f(x)。最小二乘法并不要求y=f(x)在每個(gè)點(diǎn)上都完全相等，它只要求在給定點(diǎn)xi上使誤差Σ(f(xi)-yi)2最小。在MATLAB中，可以運(yùn)用函數(shù)polyfit來進(jìn)行最小二乘多項(xiàng)式擬合，以實(shí)現(xiàn)最小二乘法曲線擬合。2/1/202395p=polyfit(x,y,n)返回?cái)M合的多項(xiàng)式的系數(shù)，系數(shù)存儲(chǔ)在向量p中。[p,s]=polyfit(x,y,n)返回是擬合生成的多項(xiàng)式系數(shù)向量p及用polyval函數(shù)獲得的誤差預(yù)測結(jié)果s。2/1/202396程序?qū)嵗?gt;>x=[1345678910];>>y=[1054211234];>>[p,s]=polyfit(x,y,4);>>y1=polyval(p,x);>>plot(x,y,'go',x,y1,'b--')2/1/2023972/1/2023988.5數(shù)值積分在工程實(shí)踐與科學(xué)應(yīng)用中，經(jīng)常要計(jì)算函數(shù)的積分與導(dǎo)數(shù)（即微分）。在已知函數(shù)形式求函數(shù)的積分時(shí)，理論上可以利用牛頓-萊布尼茲公式來計(jì)算，但在實(shí)際應(yīng)用中，經(jīng)常接觸到的許多函數(shù)都找不到其積分函數(shù)，或者是有些函數(shù)形式非常復(fù)雜，在用牛頓-萊布尼茲公式求解時(shí)也非常復(fù)雜，有時(shí)甚至根本計(jì)算不出來。此時(shí)，就可以應(yīng)用數(shù)值積分對函數(shù)進(jìn)行求積。在微積分學(xué)中，求函數(shù)的導(dǎo)數(shù)一般來說比較容易，但若所給的函數(shù)是由表格形式給出的離散點(diǎn)擬合求得時(shí)，求函數(shù)的導(dǎo)數(shù)就不那么容易了，此時(shí)就要運(yùn)用數(shù)值微分來求函數(shù)的導(dǎo)數(shù)。下面，對MATLAB如何實(shí)現(xiàn)數(shù)值積分與數(shù)值微分作一介紹。2/1/2023998.5.1Newton-Cotes求積公式Newton-Cotes求積公式適合于等間距節(jié)點(diǎn)的情形，因此也叫等距節(jié)點(diǎn)求積公式。2/1/20231001、梯形法數(shù)值積分梯形法數(shù)值積分用函數(shù)trapz來實(shí)現(xiàn)。trapz函數(shù)的調(diào)用格式如下：z=trapz(y)表示通過梯形求積法計(jì)算y的數(shù)值積分。對于向量，trapz(y)返回y的積分；對于矩陣，trapz(y)返回一行向量，向量中的元素分別對應(yīng)矩陣中每列對y進(jìn)行積分后的結(jié)果；對于N-D數(shù)組，trapz(y)從第一個(gè)非獨(dú)立維進(jìn)行計(jì)算；z=trapz(x,y)表示通過梯形求積法計(jì)算y對x的積分值。x與y必須是長度相等的向量，或者x必須是一個(gè)列向量，而y是第一個(gè)非獨(dú)立維長度與x等長的數(shù)組，trapz就從這維開始計(jì)算。2/1/2023101程序?qū)嵗?gt;>x=0:0.1:pi;>>y=sin(x);>>s=trapz(x,y)s=1.99752/1/2023102程序?qū)嵗?gt;>x=0:0.01:pi;>>y=sin(x);>>s=trapz(x,y)s=2.00002/1/20231032、simpson法數(shù)值積分此方法的數(shù)值積分用函數(shù)quad來實(shí)現(xiàn)。q=quad(‘f’,a,b)表示使用自適應(yīng)遞歸的simpson方法從區(qū)間a到b對函數(shù)f(x)積分。求積相對誤差在1e-3范圍內(nèi)。q=quad(‘f’,a,b,tol)表示使用自適應(yīng)遞歸的simpson方法從區(qū)間a到b對函數(shù)f(x)積分。求積相對誤差在tol范圍內(nèi)。此函數(shù)最大遞歸調(diào)用次數(shù)為10次，如超出則返回Inf。2/1/2023104程序?qū)嵗?gt;>s1=quad('sin(x)',0,pi)s1=2.0000>>s2=quad('sin(x)',0,2*pi)s2=02/1/20231053、cotes(科特斯)法數(shù)值積分此方法的數(shù)值積分用函數(shù)quad8來實(shí)現(xiàn)。q=quad8(‘f’,a,b)表示使用自適應(yīng)遞歸的newton-cotes8方法從區(qū)間a到b對函數(shù)f(x)積分。求積相對誤差在1e-3范圍內(nèi)。q=quad8(‘f’,a,b,tol)表示使用自適應(yīng)遞歸的newton-cotes8方法從區(qū)間a到b對函數(shù)f(x)積分。求積相對誤差在tol范圍內(nèi)。此函數(shù)最大遞歸調(diào)用次數(shù)為10次，如超出則返回Inf。2/1/2023106程序?qū)嵗?gt;>s1=quad('cos(x)',0,pi/2)s1=1.0000>>s2=quad('cos(x)',0,pi)s2=-1.1102e-0162/1/20231078.5.2Gauss求積公式在Newton-Cotes求積公式中，節(jié)點(diǎn)是等距的，從而限制了精度的提高。而Gauss求積公式將取消這個(gè)限制條件，使求積公式的精度盡可能的高。Gauss求積公式如下：2/1/2023108公式中的xk

稱為Gauss求積點(diǎn)，Ak

稱為求積系數(shù)。對于一般區(qū)間[ab]上的求積，如果用Gauss求積公式，那么必須作變量替換以使[ab]→[-11]，并有對于上式的右邊，可以應(yīng)用Gauss求積公式來進(jìn)行數(shù)值積分。2/1/2023109gausslegendre.mfunctions=gausslegendre(a,b)%用4點(diǎn)gauss-legendre求積公式球數(shù)值積分，其中a與b分別為積分區(qū)間x=[0.8611363116-0.86113631160.3399810436-0.3399810436];u=[0.34785484510.34785484510.65214515490.6521451549];t=0.0;fori=1:4y=x(i)*(b-a)*0.5+(a+b)*0.5;t=t+u(i)*gaussff(y);ends=t*(b-a)*0.5;return2/1/2023110程序?qū)嵗蒙厦婺陮懙暮瘮?shù)求下面的積分首先編寫求積函數(shù)的M文件gaussff.m2/1/2023111gaussff.mfunctiony=gaussff(x)y=x*x*exp(x);2/1/2023112程序?qū)嵗?gt;>s=gausslegendre(0,1)s=0.71832/1/20231138.5.3Romberg（龍貝格）求積公式梯形求積公式進(jìn)行數(shù)值積分的算法簡單，但收斂慢，精度低。因此人們關(guān)心的是如何發(fā)揚(yáng)梯形法的優(yōu)點(diǎn)，克服它的缺點(diǎn)，形成一個(gè)新的算法。這就是下面引進(jìn)的Romberg求積公式。

Romberg求積公式是由對近似值進(jìn)行修正而得到更近似的公式，它能自動(dòng)改變積分步長，以使其相鄰兩次值的絕對誤差或相對誤差小于預(yù)先設(shè)定的允許誤差。Romberg求積公式，據(jù)此可以編寫romberg求積函數(shù)的M文件。2/1/2023114romberg.mfunctions=romberg(a,b,eps)%romberg求積法進(jìn)行數(shù)值積分，其中a與b為積分區(qū)間，eps為允許的誤差值ifnargin==2

eps=1.0e-6;elseif

nargin<2errorreturnendt1=10000;t2=-10000;n=2;t(1,1)=0.5*(b-a)*(rombergff(a)+rombergff(b));2/1/2023115whileabs(t2-t1)>=epsarea=0.0;%n=n+1;h=(b-a)/2^(n-1);fori=1:(2^(n-1))area=area+0.5*h*(rombergff(a+h*(i-1))+rombergff(h*i+a));endt(n,1)=area;2/1/2023116

forj=2:nfori=1:(n-j+1)

t(i,j)=(4^(j-1)*t(i+1,j-1)-t(i,j-1))/(4^(j-1)-1);endendt1=t(1,n);t2=t(1,n-1);n=n+1;ends=t1;return2/1/2023117程序?qū)嵗蒙厦婢帉懙膔omberg積分函數(shù)計(jì)算如下的積分首先編寫積分函數(shù)rombergff.m2/1/2023118rombergff.mfunctiony=rombergff(x)y=cos(pi*x/2);2/1/2023119程序?qū)嵗?gt;>s=romberg(0,1)s=0.63662/1/20231208.5.4二重積分二重積分函數(shù)的使用格式如下q=dblquad(fun,x1,x2,y1,y2)q=dblquad(fun,x1,x2,y1,y2,tol)tol指定精度q=dblquad(fun,x1,x2,y1,y2,tol,method)2/1/2023121程序?qū)嵗?gt;>q=dblquad('3*x.^2+3*y.^2',0,1,0,1)q=22/1/20231228.5.5三重積分二重積分函數(shù)的使用格式如下q=triplequad(fun,x1,x2,y1,y2,z1,z2)q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol)tol指定精度q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol,method)2/1/2023123程序?qū)嵗?gt;>q=triplequad('x.^2+y.^2+z.^2',0,1,0,1,0,1)q=12/1/20231248.6常微分方程的數(shù)值解法在高等數(shù)學(xué)課程里討論的常微分方程求解方法都是一些求典型方程的解析解方法。然而，在工程實(shí)際與科學(xué)研究中遇到的微分方程往往比較復(fù)雜，在很多情況下，都不能給出解析表達(dá)式；有些雖然能給出解析表達(dá)式，但因計(jì)算量太大而不實(shí)用。以上說明用求解析解的基本方法來計(jì)算微分方程的數(shù)值解往往是不適宜的，甚至很難辦到。所以，研究微分方程的數(shù)值解法就顯得十分的必要了。下面，就討論常微分方程初值問題在MATLAB中的解法。2/1/20231258.6.1Euler方法Euler方法是數(shù)值求解一階常微分方程初值問題的最常用方法之一，按照計(jì)算精度的不同，有Euler折線法、梯形法、改進(jìn)的Euler法等。下面，我們重點(diǎn)介紹計(jì)算精度較高的改進(jìn)的Euler法。改進(jìn)的Euler法實(shí)際上是Euler折線法與梯形法聯(lián)合使用而得來的。如下所示，即為改進(jìn)的Euler法的公式2/1/2023126為了便于編寫程序，可將上式改寫為下列形式，根據(jù)上述的公式，編寫改進(jìn)的Euler法的M函數(shù)。2/1/2023127euler.mfunctions=euler(fun,x0,xn,y0,n)%用euler法計(jì)算場微分方程初值問題%x0為初值條件y(x0)=y0中的x0，y0為初始條件中的y0%xn為x取值區(qū)間的最后一個(gè)節(jié)點(diǎn)的橫坐標(biāo)值，n為把區(qū)間分成的等份數(shù)ifnargin<5errorreturnendh=(xn-x0)/n;fori=1:n

yp=y0+h*feval(fun,x0,y0);x0=x0+h;

yc=y0+h*feval(fun,x0,yp);y0=(yp+yc)/2;ends=y0;return2/1/2023128程序?qū)嵗怯蒙厦婢帉懙膃uler函數(shù)求下面的初值問題y’=-2xy0<=x<=1.2y(0)=1首先編寫所求函數(shù)文件eulerff.m2/1/2023129eulerff.mfunctiony=eulerff(x,y)y=-2*x*y;2/1/2023130程序?qū)嵗?gt;>y=euler('eulerff',0,1.2,1,100)y=0.23702/1/20231318.6.2Runge-kutta法Runge-Kutta方法是求解常微分方程初值問題的最重要的方法之一。MATLAB中專門提供了幾個(gè)采用Runge-Kutta方法來求解常微分方程的函數(shù)。重點(diǎn)介紹最常用的ode23