1、6/15/20221MATLAB2009MATLAB2009从入门到精通从入门到精通6/15/20222课程主要内容课程主要内容 第1章 MATLAB简介 第2章 数值运算 第3章 单元数组和结构 第4章 字符串 第5章 符号运算 第6章 MATLAB绘图基础 第7章 程序设计 第8章 计算方法的MATLAB实现 第9章 优化设计 第10章 SIMULINK仿真初探6/15/20223第第8章章 计算方法的计算方法的MATLAB实现实现 随着计算机的迅速发展与广泛运用，在众多的领域，科学计算方法的应用越来越广泛了，而MATLAB在进行科学计算方面有着无与伦比的优势。本章介绍MATLAB在计算方
2、法中的运用。6/15/202248.1 方程求根 roots见多项式求根；roots(多项式系数矩阵) fzero可求解非线性方程，它的格式为： fzero(function,x0) 其中function为求解的方程，x0为估计的根，x0可为标量或长度为2的向量，为向量时函数的两端的值应该符号相反，此时求区间上的解。只能求解x0附近的一个解。即使在某个区间内有多个解，但是区间端点符号相同的话仍然出错。6/15/20225程序实例 fzero(x3-3*x-1,2) ans = 1.8794 fzero(x3-3*x-1,1,4) ans = 1.8794 fzero(x3-3*x-1,2,4)
3、 ? Error using = fzero The function values at the interval endpoints must differ in sign.6/15/20226程序实例 fzero(x2-3*x+2,0) ans = 1.0000 fzero(x2-3*x+2,3) ans = 2.0000 fzero(x2-3*x+2,0,4) ? Error using = fzero The function values at the interval endpoints must differ in sign.6/15/202278.2 线性方程组数值解法 线性
4、方程组的求解不仅在工程技术领域涉及到，而且在其他的许多领域也经常碰到，因此这是一个应用相当广泛的课题。 关于线性方程组的数值解法一般分为两类：一类是直接法，就是在没有舍入误差的情况下，通过有限步四则运算求得方程组准确解的方法。另一类是迭代法，就是先给定一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法。6/15/202288.2.1 直接解法 关于线性方程组的直接解法有许多种，比如Gauss消去法、列主元消去法、平方根法等。而在MATLAB中，线性方程组的直接解法只需用符号“/”或“”就解决问题。还可以使用逆阵函数来求解： x=inv(A)*B。6/15/20229程序实例
5、求解下列方程组615318153312321321321xxxxxxxxx6/15/202210程序实例 a=12 -3 3;-18 3 -1;1 1 1; b=15;-15;6; x1=ab x1 = 1.0000 2.0000 3.0000 x2=inv(a)*b x2 = 1 2 36/15/2022118.2.2 线性方程组求解中的变换 上三角变换 U=triu(x)返回矩阵x的上三角部分； U=triu(x,k)返回第k条对角线以上部分的元素。6/15/202212程序实例 a=12 -3 3;-18 3 -1;1 1 1; triu(a) ans = 12 -3 3 0 3 -1
6、0 0 16/15/202213 triu(a,1) ans = 0 -3 3 0 0 -1 0 0 0 triu(a,-1) ans = 12 -3 3 -18 3 -1 0 1 1程序实例6/15/202214 下三角变换 U=tril(x)返回矩阵x的下三角部分； U=tril(x,k)返回第k条对角线以上下部分的元素。6/15/202215程序实例 a=12 -3 3;-18 3 -1;1 1 1; tril(a) ans = 12 0 0 -18 3 0 1 1 16/15/202216程序实例 tril(a,1) ans = 12 -3 0 -18 3 -1 1 1 1 tril(
7、a,-1) ans = 0 0 0 -18 0 0 1 1 06/15/202217 对角变换 U=diag(x)返回矩阵x主对角线上的元素，返回结果是一列向量形式； U=diag(x,k)返回第k条对角线上的元素值。 当x为向量时生成矩阵。6/15/202218程序实例 a=12 -3 3;-18 3 -1;1 1 1; diag(a) ans = 12 3 16/15/202219程序实例 a=12 -3 3;-18 3 -1;1 1 1; diag(a,1) ans = -3 -1 diag(a,-1) ans = -18 16/15/2022208.2.3 迭代解法 迭代解法非常适用于
8、求解大型稀疏系数矩阵的方程组，在线性方程组常用的迭代解法主要有Jacobi迭代法、Gauss-Seidel迭代法。6/15/202221Jacobi迭代法6/15/202222Jacobi.m function s=jacobi(a,b,x0,eps) %jacobi迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin=3 eps=1.0e-6; elseif nargin=eps x0=s; s=B*x0+f; end return6/15/202224程序实例 用上面编写的jacobi函数求解下列方程组10521510232103213
9、21321xxxxxxxxx6/15/202225程序实例 a=10 -2 -1;-2 10 -1;-1 -2 5; b=3 15 10; x=jacobi(a,b,0 0 0,eps) x = 1 2 36/15/202226Gauss-Saidel迭代法6/15/202227gauss.m function s=gauss(a,b,x0,eps) %gauss-seidel迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin=3 eps=1.0e-6; elseif nargin=eps x0=s; s=B*x0+f; end retu
10、rn6/15/202229程序实例 用上面编写的gauss函数求解下列方程组105210226210321321321xxxxxxxxx6/15/202230程序实例 a=10 -2 -1;-2 2 -1;-1 -2 5; b=6;10;10; x=gauss(a,b,0 0 0,eps) x = 4 13 86/15/2022318.3 非线性方程组数值解法 与线性方程组的求解一样，非线性方程组的求解也是应用很广泛的课题。一般情况，非线性方程组的数据值解法主要采用迭代法来求解。比较常用的迭代法主要有不动点迭代法、Newton迭代法、拟Newton迭代法等几种方法。6/15/202232不动点
11、迭代法6/15/202233staticiterate.m function s=staticiterate(x,eps) %不动点迭代法解非线性方程组，x为迭代初值，eps为允许误差 if nargin=1 eps=1.0e-6; elseif nargin=eps x=xx; xx=fx(x); end s=xx; return6/15/202235程序实例 用不动点迭代法求解下面的方程组081008102122122121xxxxxxx6/15/202236fx.m 首先编写上述非线性方程组的M文件fx.m function y=fx(x) y(1)=0.1*(x(1)*x(1)+x(2
12、)*x(2)+8); y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8); y=y(1) y(2);6/15/202237程序实例 x=staticiterate(0 0) x = 1.0000 1.00006/15/202238Newton迭代法6/15/202239newtoniterate.m function s=newtoniterate(x,eps) %newton迭代法解非线性方程组，x为迭代初值，eps为允许误差 if nargin=1 eps=1.0e-6; elseif nargin=eps x=x0+x; x1=fx1(x); x2=-dfx1(x); x3
13、=inv(x2); x0=x3*x1; end s=x0+x; return6/15/202241程序实例 用上面编写的newton迭代函数求解下列方程组081008102122122121xxxxxxx6/15/202242fx1.m和dfx1.m 首先，编写上述非线性方程组的M文件fx1.m function y=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.m function y=dfx1(x) y(1
14、)=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);6/15/202243程序实例 x=newtoniterate(0 0) x = 1 16/15/2022448.4 插值与拟合 在生产实践及科学实验中，插值与拟合的应用非常广泛。下面，就对如何用MATLAB来处理插值与拟合作一介绍。6/15/2022458.4.1 一维插值 yi=interp1(x,y,xi)返回在插值向量xi处的函数向量yi，它是根据向量x和y插值而来。如y是矩阵，则对y每一列进行插值，如xi中元素不
15、在x内，返回NaN。 yi=interp1(y,xi)省略x，表示x=1:N，此时N为向量y的长度或为矩阵y的行数。 yi=interp1(x,y,xi,method)表示用method指定的插值方法进行插值。6/15/202246 Method可取如下的值： nearest最近插值 linear线性插值 spline三次样条插值 cubic三次插值 Method默认值为线性插值，上述插值要求向量x单调。6/15/202247程序实例 x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp
16、1(x,y,x1) y1 = 5.4000 7.0500 7.50006/15/202248程序实例 x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,linear) y1 = 5.4000 7.0500 7.50006/15/202249程序实例 x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,nearest) y1 = 5 7 86/1
17、5/202250程序实例 x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,spline) y1 = 5.5520 7.1114 7.67856/15/202251程序实例 x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,cubic) y1 = 5.5006 7.0814 7.54766/15/202252程序实例 x=linspace(0
18、,2*pi,100); y=sin(x); x0=linspace(0,2*pi,6); y0=sin(x0); x1=0.62 1.9 3.2 4.34 5.55; y1=interp1(x0,y0,x1); plot(x,y,x0,y0,ro,x1,y1,+) title(默认插值)6/15/2022536/15/202254程序结果 y1 y1 = 0.4692 0.7651 -0.0546 -0.7526 -0.5549 y10=sin(x1) y10 = 0.5810 0.9463 -0.0584 -0.9315 -0.6692 deltay1=y1-y10 deltay1 = -0
19、.1118 -0.1812 0.0037 0.1789 0.11436/15/202255程序实例 x=linspace(0,2*pi,100); y=sin(x); x0=linspace(0,2*pi,6); y0=sin(x0); x1=0.62 1.9 3.2 4.34 5.55; y1=interp1(x0,y0,x1,linear); plot(x,y,x0,y0,ro,x1,y1,+) title(linear插值)6/15/2022566/15/202257程序结果 y1 y1 = 0.4692 0.7651 -0.0546 -0.7526 -0.5549 y10=sin(x1
20、) y10 = 0.5810 0.9463 -0.0584 -0.9315 -0.6692 deltay1=y1-y10 deltay1 = -0.1118 -0.1812 0.0037 0.1789 0.11436/15/202258程序实例 x=linspace(0,2*pi,100); y=sin(x); x0=linspace(0,2*pi,6); y0=sin(x0); x1=0.62 1.9 3.2 4.34 5.55; y1=interp1(x0,y0,x1,nearest); plot(x,y,x0,y0,ro,x1,y1,+) title(nearest插值)6/15/202
21、2596/15/202260程序结果 y1 y1 = 0 0.5878 -0.5878 -0.5878 -0.9511 y10=sin(x1) y10 = 0.5810 0.9463 -0.0584 -0.9315 -0.6692 deltay1=y1-y10 deltay1 = -0.5810 -0.3585 -0.5294 0.3437 -0.28186/15/202261程序实例 x=linspace(0,2*pi,100); y=sin(x); x0=linspace(0,2*pi,6); y0=sin(x0); x1=0.62 1.9 3.2 4.34 5.55; y1=interp
22、1(x0,y0,x1,spline); plot(x,y,x0,y0,ro,x1,y1,+) title(spline插值)6/15/2022626/15/202263程序结果 y1 y1 = 0.6415 0.9212 -0.0592 -0.9073 -0.7219 y10=sin(x1) y10 = 0.5810 0.9463 -0.0584 -0.9315 -0.6692 deltay1=y1-y10 deltay1 = 0.0605 -0.0251 -0.0008 0.0242 -0.05276/15/202264程序实例 x=linspace(0,2*pi,100); y=sin(x
23、); x0=linspace(0,2*pi,6); y0=sin(x0); x1=0.62 1.9 3.2 4.34 5.55; y1=interp1(x0,y0,x1,cubic); plot(x,y,x0,y0,ro,x1,y1,+) title(cubic插值)6/15/2022656/15/202266程序结果 y1 y1 = 0.6697 0.8339 -0.0689 -0.8194 -0.7563 y10=sin(x1) y10 = 0.5810 0.9463 -0.0584 -0.9315 -0.6692 deltay1=y1-y10 deltay1 = 0.0887 -0.11
24、24 -0.0106 0.1121 -0.08706/15/2022678.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)表
25、示用method指定的插值方法进行插值。Method取值同上，要求x与y都单调且具有相同格式。6/15/202268程序实例 z1=3 5 7 9 11 10 4 15;1 3 2 6 11 5 7 13; z=interp2(z1,1 3,1 2) z = 3 26/15/202269程序实例 z1=3 5 7 9 11 10 4 15;1 3 2 6 11 5 7 13; z=interp2(z1)或z=interp2(z1,1) z = Columns 1 through 12 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.000
26、0 11.0000 10.5000 10.0000 7.0000 2.0000 3.0000 4.0000 4.2500 4.5000 6.0000 7.5000 9.2500 11.0000 9.2500 7.5000 6.5000 1.0000 2.0000 3.0000 2.5000 2.0000 4.0000 6.0000 8.5000 11.0000 8.0000 5.0000 6.0000 Columns 13 through 15 4.0000 9.5000 15.0000 5.5000 9.7500 14.0000 7.0000 10.0000 13.00006/15/2022
27、70程序实例 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插值)6/15/2022716/15/202272程序实例 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插值)6/15/2022736/15/202274程序实例 x,y,z=p
28、eaks(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插值)6/15/2022756/15/202276程序实例 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插值)6/15/2022776/15/202278程序实例 x,y=meshgrid(-3:0.3:3,-3:
29、0.3:3); z=peaks(x,y); meshc(x,y,z) title(原始图)6/15/2022796/15/2022808.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)同前。6/15/202281程序实例 x,y,z,v=flow(10); xi,yi,zi=meshgrid(0.1:0.25:10,-3:0.25:3,-3:0.25:3
30、); vi=interp3(x,y,z,v,xi,yi,zi); slice(xi,yi,zi,vi,6 9.5,2,-2,0.2)6/15/2022826/15/2022838.4.4 Lagrange插值6/15/202284lagrange.m function s=lagrange(x,y,x0) %lagrange插值，x与y为已知的插值点及其函数值，x0为需要求的插值点的x值 nx=length(x); ny=length(y); if nx=ny warning(向量x与y的长度应该相同) return; end m=length(x0);6/15/202285 %按照公式，对需
31、要求的插值点向量x0的每个元素进行计算 for i=1:m t=0.0; for j=1:nx u=1.0; for k=1:nx if k=j u=u*(x0(i)-x(k)/(x(j)-x(k); end end t=t+u*y(j); end s(m)=t; end return6/15/202286程序实例 x=1 2 3 4 5; y=2 4 6 8 10; y1=lagrange(x,y,1.6) y1 = 3.2000 y2=lagrange(x,y,1.4 2.5 3.7) y2 = 0 0 7.40006/15/2022878.4.5 Newton插值6/15/202288n
32、ewton.m function s=newton(x,y,x0,nn) %newton插值，x与y为已知的插值点及其函数值 %x0为需要求的插值点的x坐标值。nn为newton插值多项式的次数 nx=length(x); ny=length(y); if nx=ny warning(向量x与y的长度应该相同) return; end m=length(x0); %按照公式，对需要求的插值点向量x0的每个元素进行计算6/15/202289 for i=1:m t=0.0; j=1; yy=y; kk=j; %求各级均差 while(kk=nn) kk=kk+1; for k=kk:nx yy(
33、k)=(yy(k)-yy(kk-1)/(x(k)-x(kk-1); end end %求插值结果6/15/202290 t=yy(1); for k=2:nn u=1.0; jj=1; while(jj x=0.4 0.55 0.65 0.8 0.9 1.05; y=0.41075 0.57815 0.69675 0.88811 1.02652 1.25382; y1=newton(x,y,0.596,4) y1 = 0.63196/15/2022928.4.6 三次样条插值 众所周知，使用高阶多项式的插值常常会产生病态的结果，而三次样条插值就能消除这种病态。 MATLAB提供的三次样条插值函
34、数有spline与interp1两个，下面重点来介绍spline函数。spline函数的调用格式如下： yy=spline(x,y,xx)利用三次样条插值法寻找在插值点xx处的插值函数值yy，插值函数根据输入参数x与y的关系得来。x与y为向量形式，而xx可以为向量形式，也可以是标量形式。此函数的作用等同于interp1(x,y,xx,spline)。6/15/202293程序实例 x=0:10; y=sin(x); xx=0:0.25:10; yy=spline(x,y,xx); plot(x,y,o,xx,yy)6/15/2022946/15/2022958.4.7 最小二乘法曲线拟合 在实
35、际工程应用与科学实践中，经常要得到一条光滑的曲线，而实际却只能测得一些分散的数据点。此时，就需要利用这些离散的点，运用各种拟合方法来生成一条连续的曲线。在这些拟合方法中，最常用的是最小二乘曲线拟合法。 已知一组数据(xi,yi)，从中找出自变量x与因变量y之间的函数关系y=f (x)。最小二乘法并不要求y=f (x)在每个点上都完全相等，它只要求在给定点xi 上使误差(f(xi)-yi)2 最小。 在MATLAB中，可以运用函数polyfit来进行最小二乘多项式拟合，以实现最小二乘法曲线拟合。6/15/202296 p=polyfit(x,y,n)返回拟合的多项式的系数，系数存储在向量p中。
36、p,s=polyfit(x,y,n)返回是拟合生成的多项式系数向量p及用polyval函数获得的误差预测结果s。6/15/202297程序实例 x=1 3 4 5 6 7 8 9 10; y=10 5 4 2 1 1 2 3 4 ; p,s=polyfit(x,y,4); y1=polyval(p,x); plot(x,y,go,x,y1,b-)6/15/2022986/15/2022998.5 数值积分 在工程实践与科学应用中，经常要计算函数的积分与导数（即微分）。在已知函数形式求函数的积分时，理论上可以利用牛顿-莱布尼兹公式来计算，但在实际应用中，经常接触到的许多函数都找不到其积分函数，或
37、者是有些函数形式非常复杂，在用牛顿-莱布尼兹公式求解时也非常复杂，有时甚至根本计算不出来。此时，就可以应用数值积分对函数进行求积。在微积分学中，求函数的导数一般来说比较容易，但若所给的函数是由表格形式给出的离散点拟合求得时，求函数的导数就不那么容易了，此时就要运用数值微分来求函数的导数。下面，对MATLAB如何实现数值积分与数值微分作一介绍。6/15/20221008.5.1 Newton-Cotes求积公式 Newton-Cotes求积公式适合于等间距节点的情形，因此也叫等距节点求积公式。6/15/20221011、梯形法数值积分 梯形法数值积分用函数trapz来实现。trapz函数的调用格
38、式如下： 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就从这维开始计算。6/15/2022102程序实例 x=0:0.1:pi; y=sin(x); s=trapz(x,y) s = 1.99756/15/2022103程序实例
39、 x=0:0.01:pi; y=sin(x); s=trapz(x,y) s = 2.00006/15/20221042、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。6/15/2022105程序实例 s1=quad(sin(x),0,pi) s1 = 2.0000
40、 s2=quad(sin(x),0,2*pi) s2 = 06/15/20221063、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。6/15/2022107程序实例 s1=quad(cos(x),0,pi/2) s1 = 1.0
41、000 s2=quad(cos(x),0,pi) s2 = -1.1102e-0166/15/20221088.5.2 Gauss求积公式 在Newton-Cotes求积公式中，节点是等距的，从而限制了精度的提高。而Gauss求积公式将取消这个限制条件，使求积公式的精度尽可能的高。Gauss求积公式如下： nkkkxfAxxf111d6/15/2022109 公式中的xk 称为Gauss求积点，Ak 称为求积系数。对于一般区间a b上的求积，如果用Gauss求积公式，那么必须作变量替换tabbax5 . 05 . 0 txtabbafabxxfbad21212d11 以使a b-1 1，并有
42、对于上式的右边，可以应用Gauss求积公式来进行数值积分。6/15/2022110gausslegendre.m function s=gausslegendre(a,b) %用4点gauss-legendre求积公式球数值积分，其中a与b分别为积分区间 x=0.8611363116 -0.8611363116 0.3399810436 -0.3399810436; u=0.3478548451 0.3478548451 0.6521451549 0.6521451549; t=0.0; for i=1:4 y=x(i)*(b-a)*0.5+(a+b)*0.5; t=t+u(i)*gaussf
43、f(y); end s=t*(b-a)*0.5; return6/15/2022111程序实例 用上面年写的函数求下面的积分102d xexIx 首先编写求积函数的M文件gaussff.m6/15/2022112gaussff.m function y=gaussff(x) y=x*x*exp(x);6/15/2022113程序实例 s=gausslegendre(0,1) s = 0.71836/15/20221148.5.3 Romberg（龙贝格）求积公式 梯形求积公式进行数值积分的算法简单，但收敛慢，精度低。因此人们关心的是如何发扬梯形法的优点，克服它的缺点，形成一个新的算法。这就是下
44、面引进的Romberg求积公式。 Romberg求积公式是由对近似值进行修正而得到更近似的公式，它能自动改变积分步长，以使其相邻两次值的绝对误差或相对误差小于预先设定的允许误差。Romberg求积公式，据此可以编写romberg求积函数的M文件。14411jkjkjjkjfTfTfT6/15/2022115romberg.m function s=romberg(a,b,eps) %romberg求积法进行数值积分，其中a与b为积分区间，eps为允许的误差值 if nargin=2 eps=1.0e-6; elseif nargin=eps area=0.0; %n=n+1; h=(b-a)/
45、2(n-1); for i=1:(2(n-1) area=area+0.5*h*(rombergff(a+h*(i-1)+rombergff(h*i+a); end t(n,1)=area;6/15/2022117 for j=2:n for i=1:(n-j+1) t(i,j)=(4(j-1)*t(i+1,j-1)-t(i,j-1)/(4(j-1)-1); end end t1=t(1,n); t2=t(1,n-1); n=n+1; end s=t1; return6/15/2022118程序实例 用上面编写的romberg积分函数计算如下的积分10dx)2cos(xI 首先编写积分函数ro
46、mbergff.m6/15/2022119rombergff.m function y=rombergff(x) y=cos(pi*x/2);6/15/2022120程序实例 s=romberg(0,1) s = 0.63666/15/20221218.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)6/15/2022122程序实例 q=dblquad(3*x.2+3*y.2,0,1,0,1)
47、q = 26/15/20221238.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)6/15/2022124程序实例 q=triplequad(x.2+y.2+z.2,0,1,0,1,0,1) q = 16/15/20221258.6 常微分方程的数值解法 在高等数学课程里讨论的常微分方程求解方法都是一些求典型方程的解析解方法。然而，在工程
48、实际与科学研究中遇到的微分方程往往比较复杂，在很多情况下，都不能给出解析表达式；有些虽然能给出解析表达式，但因计算量太大而不实用。以上说明用求解析解的基本方法来计算微分方程的数值解往往是不适宜的，甚至很难办到。所以，研究微分方程的数值解法就显得十分的必要了。下面，就讨论常微分方程初值问题在MATLAB中的解法。6/15/20221268.6.1 Euler方法 Euler方法是数值求解一阶常微分方程初值问题的最常用方法之一，按照计算精度的不同，有Euler折线法、梯形法、改进的Euler法等。下面，我们重点介绍计算精度较高的改进的Euler法。改进的Euler法实际上是Euler折线法与梯形法
49、联合使用而得来的。如下所示，即为改进的Euler法的公式*111*1,2,iiiiiiiiiiyxfyxfhyyyxhfyy6/15/2022127 为了便于编写程序，可将上式改写为下列形式，根据上述的公式，编写改进的Euler法的M函数。pcipiiciiipyyyyxhfyyyxhfyy,21,116/15/2022128euler.mfunction s=euler(fun,x0,xn,y0,n)%用euler法计算场微分方程初值问题%x0为初值条件y(x0)=y0中的x0，y0为初始条件中的y0%xn为x取值区间的最后一个节点的横坐标值，n为把区间分成的等份数if nargin5 er
50、ror returnendh=(xn-x0)/n;for i=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;return6/15/2022129程序实例 是用上面编写的euler函数求下面的初值问题 y=-2xy 0=x y=euler(eulerff,0,1.2,1,100) y = 0.23706/15/20221328.6.2 Runge-kutta法 Runge-Kutta方法是求解常微分方程初值问题的最重要的方法之一。MATLAB中专门提供了几个采用Runge
