第8章-计算方法的MATLAB实现(MATLAB课件).ppt

1、MATLAB 7.0MATLAB 7.0从入从入门到精通门到精通哈尔滨工业大学哈尔滨工业大学(威海威海)汽车工程学院汽车工程学院2023-2-112课程主要内容课程主要内容 第第1章章 MATLAB简介简介 第第2章章 数值运算数值运算 第第3章章 单元数组和结构单元数组和结构 第第4章章 字符串字符串 第第5章章 符号运算符号运算 第第6章章 MATLAB绘图基础绘图基础 第第7章章 程序设计程序设计 第第8章章 计算方法的计算方法的MATLAB实现实现 第第9章章 优化设计优化设计 第第10章章 SIMULINK仿真初探仿真初探2023-2-113第第8章章 计算方法的计算方法的MATLA

2、B实现实现 随着计算机的迅速发展与广泛运用，在众随着计算机的迅速发展与广泛运用，在众多的领域，科学计算方法的应用越来越广多的领域，科学计算方法的应用越来越广泛了，而泛了，而MATLAB在进行科学计算方面有在进行科学计算方面有着无与伦比的优势。本章介绍着无与伦比的优势。本章介绍MATLAB在在计算方法中的运用。计算方法中的运用。2023-2-1148.1 方程求根方程求根 roots见多项式求根；见多项式求根；roots(多项式系数矩阵多项式系数矩阵)fzero可求解非线性方程，它的格式为：可求解非线性方程，它的格式为：fzero(function,x0)其中其中function为求解的方程，为

3、求解的方程，x0为估计的根，为估计的根，x0可为标量或长度为可为标量或长度为2的向量，为向量时函的向量，为向量时函数的两端的值应该符号相反，此时求区间数的两端的值应该符号相反，此时求区间上的解。只能求解上的解。只能求解x0附近的一个解。即使附近的一个解。即使在某个区间内有多个解，但是区间端点符在某个区间内有多个解，但是区间端点符号相同的话仍然出错。号相同的话仍然出错。2023-2-115程序实例程序实例 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)?Error using=fzero

4、The function values at the interval endpoints must differ in sign.2023-2-116程序实例程序实例 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.2023-2-1178.2 线性方程组数值解法线性方程组数值解法 线性方程组的求解不仅在工程技术领域涉线性方程

5、组的求解不仅在工程技术领域涉及到，而且在其他的许多领域也经常碰到，及到，而且在其他的许多领域也经常碰到，因此这是一个应用相当广泛的课题。因此这是一个应用相当广泛的课题。关于线性方程组的数值解法一般分为两类：关于线性方程组的数值解法一般分为两类：一类是直接法，就是在没有舍入误差的情一类是直接法，就是在没有舍入误差的情况下，通过有限步四则运算求得方程组准况下，通过有限步四则运算求得方程组准确解的方法。另一类是迭代法，就是先给确解的方法。另一类是迭代法，就是先给定一个解的初始值，然后按一定的法则逐定一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法。步求出解的各个更准确的近似值的方

6、法。2023-2-1188.2.1 直接解法直接解法 关于线性方程组的直接解法有许多种，比关于线性方程组的直接解法有许多种，比如如Gauss消去法、列主元消去法、平方根法消去法、列主元消去法、平方根法等。而在等。而在MATLAB中，线性方程组的直接中，线性方程组的直接解法只需用符号解法只需用符号“/”或或“”就解决问题。还就解决问题。还可以使用逆阵函数来求解：可以使用逆阵函数来求解：x=inv(A)*B。2023-2-119程序实例程序实例 求解下列方程组求解下列方程组615318153312321321321xxxxxxxxx2023-2-1110程序实例程序实例 a=12-3 3;-18

7、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 32023-2-11118.2.2 线性方程组求解中的变换线性方程组求解中的变换 上三角变换上三角变换 U=triu(x)返回矩阵返回矩阵x的上三角部分；的上三角部分；U=triu(x,k)返回第返回第k条对角线以上部分的元条对角线以上部分的元素。素。2023-2-1112程序实例程序实例 a=12-3 3;-18 3-1;1 1 1;triu(a)ans=12 -3 3 0 3 -1 0 0 12023-2-1113 triu(a,1)ans=0 -3

8、3 0 0 -1 0 0 0 triu(a,-1)ans=12 -3 3 -18 3 -1 0 1 1程序实例程序实例2023-2-1114 下三角变换下三角变换 U=tril(x)返回矩阵返回矩阵x的下三角部分；的下三角部分；U=tril(x,k)返回第返回第k条对角线以上下部分的条对角线以上下部分的元素。元素。2023-2-1115程序实例程序实例 a=12-3 3;-18 3-1;1 1 1;tril(a)ans=12 0 0 -18 3 0 1 1 12023-2-1116程序实例程序实例 tril(a,1)ans=12 -3 0 -18 3 -1 1 1 1 tril(a,-1)an

9、s=0 0 0 -18 0 0 1 1 02023-2-1117 对角变换对角变换 U=diag(x)返回矩阵返回矩阵x主对角线上的元素，返主对角线上的元素，返回结果是一列向量形式；回结果是一列向量形式；U=diag(x,k)返回第返回第k条对角线上的元素值。条对角线上的元素值。当当x为向量时生成矩阵。为向量时生成矩阵。2023-2-1118程序实例程序实例 a=12-3 3;-18 3-1;1 1 1;diag(a)ans=12 3 12023-2-1119程序实例程序实例 a=12-3 3;-18 3-1;1 1 1;diag(a,1)ans=-3 -1 diag(a,-1)ans=-18

10、 12023-2-11208.2.3 迭代解法迭代解法 迭代解法非常适用于求解大型稀疏系数矩阵迭代解法非常适用于求解大型稀疏系数矩阵的方程组，在线性方程组常用的迭代解法主的方程组，在线性方程组常用的迭代解法主要有要有Jacobi迭代法、迭代法、Gauss-Seidel迭代法。迭代法。2023-2-1121Jacobi迭代法迭代法2023-2-1122Jacobi.m function s=jacobi(a,b,x0,eps)%jacobi迭代法皆线性方程组迭代法皆线性方程组%a为系数矩阵，为系数矩阵，b为方程组为方程组ax=b中的右边中的右边的矩阵的矩阵b，x0为迭代初值为迭代初值 if na

11、rgin=3 eps=1.0e-6;elseif nargin=eps x0=s;s=B*x0+f;end return2023-2-1124程序实例程序实例 用上面编写的用上面编写的jacobi函数求解下列方程组函数求解下列方程组1052151023210321321321xxxxxxxxx2023-2-1125程序实例程序实例 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 32023-2-1126Gauss-Saidel迭代法迭代法2023-2-1127gauss.m function s=gauss(a,b,

12、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 return2023-2-1129程序实例程序实例 用上面编写的用上面编写的gauss函数求解下列方程组函数求解下列方程组105210226210321321321xxxxxxxxx2023-2-1130程序实例程序实例 a=10-2-1;-2 2-1;-1-2 5;b=6;10;10;x=g

13、auss(a,b,0 0 0,eps)x=4 13 82023-2-11318.3 非线性方程组数值解法非线性方程组数值解法 与线性方程组的求解一样，非线性方程组与线性方程组的求解一样，非线性方程组的求解也是应用很广泛的课题。一般情况，的求解也是应用很广泛的课题。一般情况，非线性方程组的数值解法主要采用迭代法非线性方程组的数值解法主要采用迭代法来求解。比较常用的迭代法主要有不动点来求解。比较常用的迭代法主要有不动点迭代法、迭代法、Newton迭代法、拟迭代法、拟Newton迭代法迭代法等几种方法。等几种方法。2023-2-1132不动点迭代法不动点迭代法2023-2-1133staticite

14、rate.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;return2023-2-1135程序实例程序实例 用不动点迭代法求解下面的方程组用不动点迭代法求解下面的方程组081008102122122121xxxxxxx2023-2-1136fx.mfx.m 首先编写上述非线性方程组的首先编写上述非线性方程组的M文件文件fx.m funct

15、ion y=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);2023-2-1137程序实例程序实例 x=staticiterate(0 0)x=1.0000 1.00002023-2-1138Newton迭代法迭代法2023-2-1139newtoniterate.m function s=newtoniterate(x,eps)%newton迭代法解非线性方程组，迭代法解非线性方程组，x为迭代为迭代初值，初值，eps为允许误差为允许误差 if nargin=1 eps=1.0e

16、-6;elseif nargin=eps x=x0+x;x1=fx1(x);x2=-dfx1(x);x3=inv(x2);x0=x3*x1;end s=x0+x;return2023-2-1141程序实例程序实例 用上面编写的用上面编写的newton迭代函数求解下列迭代函数求解下列方程组方程组081008102122122121xxxxxxx2023-2-1142fx1.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(

17、2)*x(2)-10*x(2)+x(1)+8;y=y(1)y(2);然后，编写上述非线性方程组导数的然后，编写上述非线性方程组导数的M文件文件dfx1.m function y=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);2023-2-1143程序实例程序实例 x=newtoniterate(0 0)x=1 12023-2-11448.4 插值与拟合插值与拟合 在生产实践及科学实验中，插值与拟合的在生产实践及科学实验中，插值与拟合的应用非常广泛。下面，就对如

18、何用应用非常广泛。下面，就对如何用MATLAB来处理插值与拟合作一介绍。来处理插值与拟合作一介绍。2023-2-11458.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,met

19、hod)表示用表示用method指指定的插值方法进行插值。定的插值方法进行插值。2023-2-1146 Method可取如下的值：可取如下的值：nearest最近插值最近插值 linear线性插值线性插值 spline三次样条插值三次样条插值 cubic三次插值三次插值 Method默认值为线性插值，上述插值要求默认值为线性插值，上述插值要求向量向量x单调。单调。2023-2-1147程序实例程序实例 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)y1=5.4000 7

20、.0500 7.50002023-2-1148程序实例程序实例 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.50002023-2-1149程序实例程序实例 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 82023-2-1150程序实例程序实例 x=1 2 4 6

21、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.67852023-2-1151程序实例程序实例 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.54762023-2-1152程序实例程序实例 x=linspace(0,2*pi,100);y=sin(x);x0=lins

22、pace(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(默认插值默认插值)2023-2-11532023-2-1154程序结果程序结果 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.1118 -0.1812 0.0037 0.1789 0.114

23、32023-2-1155程序实例程序实例 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插值插值)2023-2-11562023-2-1157程序结果程序结果 y1 y1=0.4692 0.7651 -0.0546 -0.7526 -0.5549 y10=sin(x1)y10=0.5810 0.9463 -0.0584 -0.93

24、15 -0.6692 deltay1=y1-y10 deltay1=-0.1118 -0.1812 0.0037 0.1789 0.11432023-2-1158程序实例程序实例 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插值插值)2023-2-11592023-2-1160程序结果程序结果 y1 y1=0 0.5878

25、-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.28182023-2-1161程序实例程序实例 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,spline);plot(x,y,x0,y0,ro,x1,y1

26、,+)title(spline插值插值)2023-2-11622023-2-1163程序结果程序结果 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.05272023-2-1164程序实例程序实例 x=linspace(0,2*pi,100);y=sin(x);x0=linspace(0,2*pi,6);y0=sin(x0);x1=0.62

27、 1.9 3.2 4.34 5.55;y1=interp1(x0,y0,x1,cubic);plot(x,y,x0,y0,ro,x1,y1,+)title(cubic插值插值)2023-2-11652023-2-1166程序结果程序结果 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.1124 -0.0106 0.1121 -0.08702023-2-11678.4.2 二维插值二维插

28、值 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次迭代，次迭代，nti

29、mes默认为默认为1。zi=interp2(x,y,z,xi,yi,method)表示用表示用method指定的插值方法进行插值。指定的插值方法进行插值。Method取值同上，要求取值同上，要求x与与y都单调且具有相同格都单调且具有相同格式。式。2023-2-1168程序实例程序实例 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 22023-2-1169程序实例程序实例 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

30、=Columns 1 through 12 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 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 1

31、5 4.0000 9.5000 15.0000 5.5000 9.7500 14.0000 7.0000 10.0000 13.00002023-2-1170程序实例程序实例 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插值插值)2023-2-11712023-2-1172程序实例程序实例 x,y,z=peaks(10);xi,yi=meshgrid(-3:0.3:3,-3:0.3:3);zi=interp2(x,y,z

32、,xi,yi,linear);meshc(xi,yi,zi)title(linear插值插值)2023-2-11732023-2-1174程序实例程序实例 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插值插值)2023-2-11752023-2-1176程序实例程序实例 x,y,z=peaks(10);xi,yi=meshgrid(-3:0.3:3,-3:0.3:3);zi=interp2(x,y,z,xi,yi,cubic

33、);meshc(xi,yi,zi)title(cubic插值插值)2023-2-11772023-2-1178程序实例程序实例 x,y=meshgrid(-3:0.3:3,-3:0.3:3);z=peaks(x,y);meshc(x,y,z)title(原始图原始图)2023-2-11792023-2-11808.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。同前。同

34、前。vi=interp3(x,y,z,v,xi,yi,zi,method)同前。同前。2023-2-1181程序实例程序实例 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,6 9.5,2,-2,0.2)2023-2-11822023-2-11838.4.4 Lagrange插值插值2023-2-1184lagrange.m function s=lagrange(x,y,x0)%lagrange插值，插值，x与与

35、y为已知的插值点及其为已知的插值点及其函数值，函数值，x0为需要求的插值点的为需要求的插值点的x值值 nx=length(x);ny=length(y);if nx=ny warning(向量向量x与与y的长度应该相同的长度应该相同)return;end m=length(x0);2023-2-1185%按照公式，对需要求的插值点向量按照公式，对需要求的插值点向量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)

36、;end s(m)=t;end return2023-2-1186程序实例程序实例 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.40002023-2-11878.4.5 Newton插值插值2023-2-1188newton.m function s=newton(x,y,x0,nn)%newton插值，插值，x与与y为已知的插值点及其函数值为已知的插值点及其函数值%x0为需要求的插值点的为需要求的插值点的x坐标值。坐标值。nn为为newton插值多项式的

37、次数插值多项式的次数 nx=length(x);ny=length(y);if nx=ny warning(向量向量x与与y的长度应该相同的长度应该相同)return;end m=length(x0);%按照公式，对需要求的插值点向量按照公式，对需要求的插值点向量x0的每个元的每个元素进行计算素进行计算2023-2-1189 for i=1:m t=0.0;j=1;yy=y;kk=j;%求各级均差求各级均差 while(kk=nn)kk=kk+1;for k=kk:nx yy(k)=(yy(k)-yy(kk-1)/(x(k)-x(kk-1);end end%求插值结果求插值结果2023-2-1

38、190 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.63192023-2-11928.4.6 三次样条插值三次样条插值 众所周知，使用高阶多项式的插值常常会产生病众所周知，使用高阶多项式的插值常常会产生病态的结果，而三次样条插值就能消除这种病态。态的结果，而三次样条插值就能消除这种病态。MATLAB提供的三次样条插值函数有提供的三次样条插值函数有spl

39、ine与与interp1两个，下面重点来介绍两个，下面重点来介绍spline函数。函数。spline函数的调用格式如下：函数的调用格式如下：yy=spline(x,y,xx)利用三次样条插值法寻找在插值利用三次样条插值法寻找在插值点点xx处的插值函数值处的插值函数值yy，插值函数根据输入参数，插值函数根据输入参数x与与y的关系得来。的关系得来。x与与y为向量形式，而为向量形式，而xx可以为向可以为向量形式，也可以是标量形式。此函数的作用等同量形式，也可以是标量形式。此函数的作用等同于于interp1(x,y,xx,spline)。2023-2-1193程序实例程序实例 x=0:10;y=sin

40、(x);xx=0:0.25:10;yy=spline(x,y,xx);plot(x,y,o,xx,yy)2023-2-11942023-2-11958.4.7 最小二乘法曲线拟合最小二乘法曲线拟合 在实际工程应用与科学实践中，经常要得到一条在实际工程应用与科学实践中，经常要得到一条光滑的曲线，而实际却只能测得一些分散的数据光滑的曲线，而实际却只能测得一些分散的数据点。此时，就需要利用这些离散的点，运用各种点。此时，就需要利用这些离散的点，运用各种拟合方法来生成一条连续的曲线。在这些拟合方拟合方法来生成一条连续的曲线。在这些拟合方法中，最常用的是最小二乘曲线拟合法。法中，最常用的是最小二乘曲线拟

41、合法。已知一组数据已知一组数据(xi,yi)，从中找出自变量，从中找出自变量x与因变量与因变量y之间的函数关系之间的函数关系y=f(x)。最小二乘法并不要求。最小二乘法并不要求y=f(x)在每个点上都完全相等，它只要求在给定点在每个点上都完全相等，它只要求在给定点xi 上使误差上使误差(f(xi)-yi)2 最小。最小。在在MATLAB中，可以运用函数中，可以运用函数polyfit来进行最小来进行最小二乘多项式拟合，以实现最小二乘法曲线拟合。二乘多项式拟合，以实现最小二乘法曲线拟合。2023-2-1196 p=polyfit(x,y,n)返回拟合的多项式的系数，返回拟合的多项式的系数，系数存储

42、在向量系数存储在向量p中。中。p,s=polyfit(x,y,n)返回是拟合生成的多项式返回是拟合生成的多项式系数向量系数向量p及用及用polyval函数获得的误差预测函数获得的误差预测结果结果s。2023-2-1197程序实例程序实例 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-)2023-2-11982023-2-11998.5 数值积分数值积分 在工程实践与科学应用中，经常要计算函数的积在工程实践与科学应用中，经常要计算函数的积分与导数（即微

43、分）。在已知函数形式求函数的分与导数（即微分）。在已知函数形式求函数的积分时，理论上可以利用牛顿积分时，理论上可以利用牛顿-莱布尼兹公式来计莱布尼兹公式来计算，但在实际应用中，经常接触到的许多函数都算，但在实际应用中，经常接触到的许多函数都找不到其积分函数，或者是有些函数形式非常复找不到其积分函数，或者是有些函数形式非常复杂，在用牛顿杂，在用牛顿-莱布尼兹公式求解时也非常复杂，莱布尼兹公式求解时也非常复杂，有时甚至根本计算不出来。此时，就可以应用数有时甚至根本计算不出来。此时，就可以应用数值积分对函数进行求积。在微积分学中，求函数值积分对函数进行求积。在微积分学中，求函数的导数一般来说比较容易

44、，但若所给的函数是由的导数一般来说比较容易，但若所给的函数是由表格形式给出的离散点拟合求得时，求函数的导表格形式给出的离散点拟合求得时，求函数的导数就不那么容易了，此时就要运用数值微分来求数就不那么容易了，此时就要运用数值微分来求函数的导数。下面，对函数的导数。下面，对MATLAB如何实现数值积如何实现数值积分与数值微分作一介绍。分与数值微分作一介绍。2023-2-111008.5.1 Newton-Cotes求积公式求积公式 Newton-Cotes求积公式适合于等间距节点求积公式适合于等间距节点的情形，因此也叫等距节点求积公式。的情形，因此也叫等距节点求积公式。2023-2-111011、

45、梯形法数值积分、梯形法数值积分 梯形法数值积分用函数梯形法数值积分用函数trapz来实现。来实现。trapz函数的函数的调用格式如下：调用格式如下：z=trapz(y)表示通过梯形求积法计算表示通过梯形求积法计算y的数值积分。的数值积分。对于向量，对于向量，trapz(y)返回返回y的积分；对于矩阵，的积分；对于矩阵，trapz(y)返回一行向量，向量中的元素分别对应矩返回一行向量，向量中的元素分别对应矩阵中每列对阵中每列对y进行积分后的结果；对于进行积分后的结果；对于N-D数组，数组，trapz(y)从第一个非独立维进行计算；从第一个非独立维进行计算；z=trapz(x,y)表示通过梯形求积

46、法计算表示通过梯形求积法计算y对对x的积分的积分值。值。x与与y必须是长度相等的向量，或者必须是长度相等的向量，或者x必须是一必须是一个列向量，而个列向量，而y是第一个非独立维长度与是第一个非独立维长度与x等长的等长的数组，数组，trapz就从这维开始计算。就从这维开始计算。2023-2-11102程序实例程序实例 x=0:0.1:pi;y=sin(x);s=trapz(x,y)s=1.99752023-2-11103程序实例程序实例 x=0:0.01:pi;y=sin(x);s=trapz(x,y)s=2.00002023-2-111042、simpson法数值积分法数值积分 此方法的数值积

47、分用函数此方法的数值积分用函数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。2023-2-11105程序实例程序实例 s1=quad(sin

48、(x),0,pi)s1=2.0000 s2=quad(sin(x),0,2*pi)s2=02023-2-111063、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)积积分。

49、求积相对误差在分。求积相对误差在tol范围内。范围内。此函数最大递归调用次数为此函数最大递归调用次数为10次，如超出次，如超出则返回则返回Inf。2023-2-11107程序实例程序实例 s1=quad(cos(x),0,pi/2)s1=1.0000 s2=quad(cos(x),0,pi)s2=-1.1102e-0162023-2-111088.5.2 Gauss求积公式求积公式 在在Newton-Cotes求积公式中，节点是等距求积公式中，节点是等距的，从而限制了精度的提高。而的，从而限制了精度的提高。而Gauss求积求积公式将取消这个限制条件，使求积公式的公式将取消这个限制条件，使求积公

50、式的精度尽可能的高。精度尽可能的高。Gauss求积公式如下：求积公式如下：nkkkxfAxxf111d2023-2-11109 公式中的公式中的xk 称为称为Gauss求积点，求积点，Ak 称为求积系数。称为求积系数。对于一般区间对于一般区间a b上的求积，如果用上的求积，如果用Gauss求积公求积公式，那么必须作变量替换式，那么必须作变量替换tabbax5.05.0 txtabbafabxxfbad21212d11 以使以使a b-1 1，并有，并有 对于上式的右边，可以应用对于上式的右边，可以应用Gauss求积公式来进行求积公式来进行数值积分。数值积分。2023-2-11110gaussl