第8讲计算方法MATLAB实现课件

MATLAB2009从入门到精通7/23/20231课程主要内容第1章MATLAB简介第2章数值运算第3章单元数组和结构第4章字符串第5章符号运算第6章MATLAB绘图基础第7章程序设计第8章计算方法的MATLAB实现第9章优化设计第10章SIMULINK仿真初探7/23/20232第8章计算方法的MATLAB实现

随着计算机的迅速发展与广泛运用，在众多的领域，科学计算方法的应用越来越广泛了，而MATLAB在进行科学计算方面有着无与伦比的优势。本章介绍MATLAB在计算方法中的运用。7/23/202338.1方程求根roots见多项式求根；roots(多项式系数矩阵)fzero可求解非线性方程，它的格式为：fzero(‘function’,x0)其中function为求解的方程，x0为估计的根，x0可为标量或长度为2的向量，为向量时函数的两端的值应该符号相反，此时求区间上的解。只能求解x0附近的一个解。即使在某个区间内有多个解，但是区间端点符号相同的话仍然出错。7/23/20234程序实例>>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.7/23/20235程序实例>>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.7/23/202368.2线性方程组数值解法线性方程组的求解不仅在工程技术领域涉及到，而且在其他的许多领域也经常碰到，因此这是一个应用相当广泛的课题。关于线性方程组的数值解法一般分为两类：一类是直接法，就是在没有舍入误差的情况下，通过有限步四则运算求得方程组准确解的方法。另一类是迭代法，就是先给定一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法。7/23/202378.2.1直接解法关于线性方程组的直接解法有许多种，比如Gauss消去法、列主元消去法、平方根法等。而在MATLAB中，线性方程组的直接解法只需用符号“/”或“\”就解决问题。还可以使用逆阵函数来求解：x=inv(A)*B。7/23/20238程序实例求解下列方程组7/23/20239程序实例>>a=[12-33;-183-1;111];>>b=[15;-15;6];>>x1=a\bx1=1.00002.00003.0000>>x2=inv(a)*bx2=1237/23/2023108.2.2线性方程组求解中的变换上三角变换U=triu(x)返回矩阵x的上三角部分；U=triu(x,k)返回第k条对角线以上部分的元素。7/23/202311程序实例>>a=[12-33;-183-1;111];>>triu(a)ans=12-3303-10017/23/202312>>triu(a,1)ans=0-3300-1000>>triu(a,-1)ans=12-33-183-1011程序实例7/23/202313下三角变换U=tril(x)返回矩阵x的下三角部分；U=tril(x,k)返回第k条对角线以上下部分的元素。7/23/202314程序实例>>a=[12-33;-183-1;111];>>tril(a)ans=1200-18301117/23/202315程序实例>>tril(a,1)ans=12-30-183-1111>>tril(a,-1)ans23/202316对角变换U=diag(x)返回矩阵x主对角线上的元素，返回结果是一列向量形式；U=diag(x,k)返回第k条对角线上的元素值。当x为向量时生成矩阵。7/23/202317程序实例>>a=[12-33;-183-1;111];>>diag(a)ans=12317/23/202318程序实例>>a=[12-33;-183-1;111];>>diag(a,1)ans=-3-1>>diag(a,-1)ans=-1817/23/2023198.2.3迭代解法迭代解法非常适用于求解大型稀疏系数矩阵的方程组，在线性方程组常用的迭代解法主要有Jacobi迭代法、Gauss-Seidel迭代法。7/23/202320Jacobi迭代法7/23/202321Jacobi.mfunctions=jacobi(a,b,x0,eps)%jacobi迭代法皆线性方程组%a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值ifnargin==3

eps=1.0e-6;elseif

nargin<3errorreturnend7/23/202322D=diag(diag(a));%求出对角矩阵D=inv(D);%求出对角矩阵的逆L=tril(a,-1);%求出严格下三角矩阵U=triu(a,1);%求出严格上三角矩阵B=-D*(L+U);f=D*b;s=B*x0+f;whilenorm(s-x0)>=epsx0=s;s=B*x0+f;endreturn7/23/202323程序实例用上面编写的jacobi函数求解下列方程组7/23/202324程序实例>>a=[10-2-1;-210-1;-1-25];>>b=[31510]';>>x=jacobi(a,b,[000]',eps)x=1237/23/202325Gauss-Saidel迭代法7/23/202326gauss.mfunctions=gauss(a,b,x0,eps)%gauss-seidel迭代法皆线性方程组%a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值ifnargin==3

eps=1.0e-6;elseif

nargin<3errorreturnend7/23/202327L=tril(a,-1);%求出严格下三角矩阵D=diag(diag(a));%求出对角矩阵U=triu(a,1);%求出严格上三角矩阵C=inv(D+L);B=-C*U;f=C*b;s=B*x0+f;whilenorm(s-x0)>=epsx0=s;s=B*x0+f;endreturn7/23/202328程序实例用上面编写的gauss函数求解下列方程组7/23/202329程序实例>>a=[10-2-1;-22-1;-1-25];>>b=[6;10;10];>>x=gauss(a,b,[000]',eps)x=41387/23/2023308.3非线性方程组数值解法与线性方程组的求解一样，非线性方程组的求解也是应用很广泛的课题。一般情况，非线性方程组的数据值解法主要采用迭代法来求解。比较常用的迭代法主要有不动点迭代法、Newton迭代法、拟Newton迭代法等几种方法。7/23/202331不动点迭代法7/23/202332staticiterate.mfunctions=staticiterate(x,eps)%不动点迭代法解非线性方程组，x为迭代初值，eps为允许误差ifnargin==1

eps=1.0e-6;elseif

nargin<1errorreturnend7/23/202333xx=fx(x);%第一次迭代whilenorm(xx-x)>=epsx=xx;xx=fx(x);ends=xx;return7/23/202334程序实例用不动点迭代法求解下面的方程组7/23/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)];7/23/202336程序实例>>x=staticiterate([00])x=1.00001.00007/23/202337Newton迭代法7/23/202338newtoniterate.mfunctions=newtoniterate(x,eps)%newton迭代法解非线性方程组，x为迭代初值，eps为允许误差ifnargin==1

eps=1.0e-6;elseif

nargin<1errorreturnend7/23/202339x1=fx1(x);%非线性方程组x2=-dfx1(x);%非线性方程组导数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;return7/23/202340程序实例用上面编写的newton迭代函数求解下列方程组7/23/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)];然后，编写上述非线性方程组导数的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)];7/23/202342程序实例>>x=newtoniterate([00])x=117/23/2023438.4插值与拟合在生产实践及科学实验中，插值与拟合的应用非常广泛。下面，就对如何用MATLAB来处理插值与拟合作一介绍。7/23/2023448.4.1一维插值yi=interp1(x,y,xi)返回在插值向量xi处的函数向量yi，它是根据向量x和y插值而来。如y是矩阵，则对y每一列进行插值，如xi中元素不在x内，返回NaN。yi=interp1(y,xi)省略x，表示x=1:N，此时N为向量y的长度或为矩阵y的行数。yi=interp1(x,y,xi,’method’)表示用method指定的插值方法进行插值。7/23/202345Method可取如下的值：‘nearest’最近插值‘linear’线性插值‘spline’三次样条插值‘cubic’三次插值Method默认值为线性插值，上述插值要求向量x单调。7/23/202346程序实例>>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1)y1=5.40007.05007.50007/23/202347程序实例>>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'linear')y1=5.40007.05007.50007/23/202348程序实例>>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'nearest')y1=5787/23/202349程序实例>>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'spline')y1=5.55207.11147.67857/23/202350程序实例>>x=[12468910131516];>>y=[57810131415171920];>>x1=[1.22.13];>>y1=interp1(x,y,x1,'cubic')y1=5.50067.08147.54767/23/202351程序实例>>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('默认插值')7/23/2023527/23/202353程序结果>>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.11437/23/202354程序实例>>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插值')7/23/2023557/23/202356程序结果>>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.11437/23/202357程序实例>>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插值')7/23/2023587/23/202359程序结果>>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.28187/23/202360程序实例>>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插值')7/23/2023617/23/202362程序结果>>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.05277/23/202363程序实例>>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插值')7/23/2023647/23/202365程序结果>>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.08707/23/2023668.4.2二维插值zi=interp2(x,y,z,xi,yi)返回在插值向量xi、yi处的函数值向量，它是根据向量x、y与z插值而来，这里的x、y与z也可以是矩阵形式。如果xi、yi有元素不在x、y范围内，则返回NaN。zi=interp2(z,xi,yi)省略x、y，表示x=1:N，y=1:M。此处[M,N]=size(z)。zi=interp2(z,ntimes)表示在z的各点之间插入数据点来扩展z，依次执行ntimes次迭代，ntimes默认为1。zi=interp2(x,y,z,xi,yi,’method’)表示用method指定的插值方法进行插值。Method取值同上，要求x与y都单调且具有相同格式。7/23/202367程序实例>>z1=[35791110415;1326115713];>>z=interp2(z1,[13],[12])z=327/23/202368程序实例>>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.00007/23/202369程序实例>>[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插值')7/23/2023707/23/202371程序实例>>[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插值')7/23/2023727/23/202373程序实例>>[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插值')7/23/2023747/23/202375程序实例>>[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插值')7/23/2023767/23/202377程序实例>>[x,y]=meshgrid(-3:0.3:3,-3:0.3:3);>>z=peaks(x,y);>>meshc(x,y,z)>>title('原始图')7/23/2023787/23/2023798.4.3三维插值vi=interp3(x,y,z,v,xi,yi,zi)返回三维函数v在插值向量xi，yi，zi处的函数向量vi，它们大小形式相同。vi=interp3(v,xi,yi,zi)省略x、y、z。同前。vi=interp3(x,y,z,v,xi,yi,zi,’method’)同前。7/23/202380程序实例>>[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])7/23/2023817/23/2023828.4.4Lagrange插值7/23/202383lagrange.mfunctions=lagrange(x,y,x0)%lagrange插值，x与y为已知的插值点及其函数值，x0为需要求的插值点的x值nx=length(x);ny=length(y);ifnx~=nywarning('向量x与y的长度应该相同')return;endm=length(x0);7/23/202384%按照公式，对需要求的插值点向量x0的每个元素进行计算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;endreturn7/23/202385程序实例>>x=[12345];>>y=[246810];>>y1=lagrange(x,y,1.6)y1=3.2000>>y2=lagrange(x,y,[1.42.53.7])y2=007.40007/23/2023868.4.5Newton插值7/23/202387newton.mfunctions=newton(x,y,x0,nn)%newton插值，x与y为已知的插值点及其函数值%x0为需要求的插值点的x坐标值。nn为newton插值多项式的次数nx=length(x);ny=length(y);ifnx~=nywarning('向量x与y的长度应该相同')return;endm=length(x0);%按照公式，对需要求的插值点向量x0的每个元素进行计算7/23/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%求插值结果7/23/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;endreturn7/23/202390程序实例>>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.63197/23/2023918.4.6三次样条插值众所周知，使用高阶多项式的插值常常会产生病态的结果，而三次样条插值就能消除这种病态。MATLAB提供的三次样条插值函数有spline与interp1两个，下面重点来介绍spline函数。spline函数的调用格式如下：yy=spline(x,y,xx)利用三次样条插值法寻找在插值点xx处的插值函数值yy，插值函数根据输入参数x与y的关系得来。x与y为向量形式，而xx可以为向量形式，也可以是标量形式。此函数的作用等同于interp1(x,y,xx,’spline’)。7/23/202392程序实例>>x=0:10;>>y=sin(x);>>xx=0:0.25:10;>>yy=spline(x,y,xx);>>plot(x,y,'o',xx,yy)7/23/2023937/23/2023948.4.7最小二乘法曲线拟合在实际工程应用与科学实践中，经常要得到一条光滑的曲线，而实际却只能测得一些分散的数据点。此时，就需要利用这些离散的点，运用各种拟合方法来生成一条连续的曲线。在这些拟合方法中，最常用的是最小二乘曲线拟合法。已知一组数据(xi,yi)，从中找出自变量x与因变量y之间的函数关系y=f(x)。最小二乘法并不要求y=f(x)在每个点上都完全相等，它只要求在给定点xi上使误差Σ(f(xi)-yi)2最小。在MATLAB中，可以运用函数polyfit来进行最小二乘多项式拟合，以实现最小二乘法曲线拟合。7/23/202395p=polyfit(x,y,n)返回拟合的多项式的系数，系数存储在向量p中。[p,s]=polyfit(x,y,n)返回是拟合生成的多项式系数向量p及用polyval函数获得的误差预测结果s。7/23/202396程序实例>>x=[1345678910];>>y=[1054211234];>>[p,s]=polyfit(x,y,4);>>y1=polyval(p,x);>>plot(x,y,'go',x,y1,'b--')7/23/2023977/23/2023988.5数值积分在工程实践与科学应用中，经常要计算函数的积分与导数（即微分）。在已知函数形式求函数的积分时，理论上可以利用牛顿-莱布尼兹公式来计算，但在实际应用中，经常接触到的许多函数都找不到其积分函数，或者是有些函数形式非常复杂，在用牛顿-莱布尼兹公式求解时也非常复杂，有时甚至根本计算不出来。此时，就可以应用数值积分对函数进行求积。在微积分学中，求函数的导数一般来说比较容易，但若所给的函数是由表格形式给出的离散点拟合求得时，求函数的导数就不那么容易了，此时就要运用数值微分来求函数的导数。下面，对MATLAB如何实现数值积分与数值微分作一介绍。7/23/2023998.5.1Newton-Cotes求积公式Newton-Cotes求积公式适合于等间距节点的情形，因此也叫等距节点求积公式。7/23/20231001、梯形法数值积分梯形法数值积分用函数trapz来实现。trapz函数的调用格式如下：z=trapz(y)表示通过梯形求积法计算y的数值积分。对于向量，trapz(y)返回y的积分；对于矩阵，trapz(y)返回一行向量，向量中的元素分别对应矩阵中每列对y进行积分后的结果；对于N-D数组，trapz(y)从第一个非独立维进行计算；z=trapz(x,y)表示通过梯形求积法计算y对x的积分值。X与y必须是长度相等的向量，或者x必须是一个列向量，而y是第一个非独立维长度与x等长的数组，trapz就从这维开始计算。7/23/2023101程序实例>>x=0:0.1:pi;>>y=sin(x);>>s=trapz(x,y)s=1.99757/23/2023102程序实例>>x=0:0.01:pi;>>y=sin(x);>>s=trapz(x,y)s=2.00007/23/20231032、simpson法数值积分此方法的数值积分用函数quad来实现。q=quad(‘f’,a,b)表示使用自适应递归的simpson方法从区间a到b对函数f(x)积分。求积相对误差在1e-3范围内。q=quad(‘f’,a,b,tol)表示使用自适应递归的simpson方法从区间a到b对函数f(x)积分。求积相对误差在tol范围内。此函数最大递归调用次数为10次，如超出则返回Inf。7/23/2023104程序实例>>s1=quad('sin(x)',0,pi)s1=2.0000>>s2=quad('sin(x)',0,2*pi)s2=07/23/20231053、cotes(科特斯)法数值积分此方法的数值积分用函数quad8来实现。q=quad8(‘f’,a,b)表示使用自适应递归的newton-cotes8方法从区间a到b对函数f(x)积分。求积相对误差在1e-3范围内。q=quad8(‘f’,a,b,tol)表示使用自适应递归的newton-cotes8方法从区间a到b对函数f(x)积分。求积相对误差在tol范围内。此函数最大递归调用次数为10次，如超出则返回Inf。7/23/2023106程序实例>>s1=quad('cos(x)',0,pi/2)s1=1.0000>>s2=quad('cos(x)',0,pi)s2=-1.1102e-0167/23/20231078.5.2Gauss求积公式在Newton-Cotes求积公式中，节点是等距的，从而限制了精度的提高。而Gauss求积公式将取消这个限制条件，使求积公式的精度尽可能的高。Gauss求积公式如下：7/23/2023108公式中的xk

称为Gauss求积点，Ak

称为求积系数。对于一般区间[ab]上的求积，如果用Gauss求积公式，那么必须作变量替换以使[ab]→[-11]，并有对于上式的右边，可以应用Gauss求积公式来进行数值积分。7/23/2023109gausslegendre.mfunctions=gausslegendre(a,b)%用4点gauss-legendre求积公式球数值积分，其中a与b分别为积分区间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;return7/23/2023110程序实例用上面年写的函数求下面的积分首先编写求积函数的M文件gaussff.m7/23/2023111gaussff.mfunctiony=gaussff(x)y=x*x*exp(x);7/23/2023112程序实例>>s=gausslegendre(0,1)s=0.71837/23/20231138.5.3Romberg（龙贝格）求积公式梯形求积公式进行数值积分的算法简单，但收敛慢，精度低。因此人们关心的是如何发扬梯形法的优点，克服它的缺点，形成一个新的算法。这就是下面引进的Romberg求积公式。

Romberg求积公式是由对近似值进行修正而得到更近似的公式，它能自动改变积分步长，以使其相邻两次值的绝对误差或相对误差小于预先设定的允许误差。Romberg求积公式，据此可以编写romberg求积函数的M文件。7/23/2023114romberg.mfunctions=romberg(a,b,eps)%romberg求积法进行数值积分，其中a与b为积分区间，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));7/23/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;7/23/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;return7/23/2023117程序实例用上面编写的romberg积分函数计算如下的积分首先编写积分函数rombergff.m7/23/2023118rombergff.mfunctiony=rombergff(x)y=cos(pi*x/2);7/23/2023119程序实例>>s=romberg(0,1)s=0.63667/23/20231208.5.4二重积分二重积分函数的使用格式如下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)7/23/2023121程序实例>>q=dblquad('3*x.^2+3*y.^2',0,1,0,1)q=27/23/20231228.5.5三重积分二重积分函数的使用格式如下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)7/23/2023123程序实例>>q=triplequad('x.^2+y.^2+z.^2',0,1,0,1,0,1)q=17/23/20231248.6常微分方程的数值解法在高等数学课程里讨论的常微分方程求解方法都是一些求典型方程的解析解方法。然而，在工程实际与科学研究中遇到的微分方程往往比较复杂，在很多情况下，都不能给出解析表达式；有些虽然能给出解析表达式，但因计算量太大而不实用。以上说明用求解析解的基本方法来计算微分方程的数值解往往是不适宜的，甚至很难办到。所以，研究微分方程的数值解法就显得十分的必要了。下面，就讨论常微分方程初值问题在MATLAB中的解法。7/23/20231258.6.1Euler方法Euler方法是数值求解一阶常微分方程初值问题的最常用方法之一，按照计算精度的不同，有Euler折线法、梯形法、改进的Euler法等。下面，我们重点介绍计算精度较高的改进的Euler法。改进的Euler法实际上是Euler折线法与梯形法联合使用而得来的。如下所示，即为改进的Euler法的公式7/23/2023126为了便于编写程序，可将上式改写为下列形式，根据上述的公式，编写改进的Euler法的M函数。7/23/2023127euler.mfunctions=euler(fun,x0,xn,y0,n)%用euler法计算场微分方程初值问题%x0为初值条件y(x0)=y0中的x0，y0为初始条件中的y0%xn为x取值区间的最后一个节点的横坐标值，n为把区间分成的等份数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;return7/23/2023128程序实例是用上面编写的euler函数求下面的初值问题y’=-2xy0<=x<=1.2y(0)=1首先编写所求函数文件eulerff.m7/23/2023129eulerff.mfunctiony=eulerff(x,y)y=-2*x*y;7/23/2023130程序实例>>y=euler('eulerff',0,1.2,1,100)y=0.23707/23/20231318.6.2Runge-kutta法Runge-Kutta方法是求解常微分方程初值问题的最重要的方法之一。MATLAB中专门提供了几个采用Runge-Kutta方法来求解常微分方程的函数。重点介绍最常用的ode23,ode45这两个函数。二三阶runge