MATLAB基础与实例进阶ppt课件

第10章 MATLAB科学计算,【学习目标】 掌握常用的插值函数和拟合函数的用法 会用MATLAB求各类微积分 掌握各种线性方程组和非线性方程组的求 解方法，并用MATLAB实现 会用PDETOOL求解偏微分方程 了解最优化问题及常用的函数,Page 1,第10章 MATLAB科学计算,Page 2,10.1 数据插值,数据插值是函数逼近的一种重要方法，根据数据的分布规律，找到一个函数表达式可以连接已知的各点，并用这一函数表达式预测两点之间任意位置上的函数值，经常用于科学与工程研究领域。最简单的插值法是多项式插值法，采用MATLAB，既可以利用其提供的插值函数，也可以利用编程来实现数据的插值。,Page 3,1.一维插值,1）一维插值函数interp1 MATLAB中通常用于一维数据插值（表格查找）的函数是interp1，常用的语法形式及描述如下：,yi = interp1(x,Y,xi) 或yi = interp1(Y,xi) yi = interp1(x,Y,xi,method) yi = interp1(x,Y,xi,method,extrap) yi = interp1(x,Y,xi,method,extrapval) pp = interp1(x,Y,method,pp),Page 4,2019/10/12,机械工业出版社,Page 5,method：用特定的方法插值，可取的值有nearest（最邻近插值）、linear（分段线性插值，默认取值）、spline（样条插值）、pchip（分段三次厄密多项式插值，可用pchip函数取代）、cubic（同pchip）和v5cubic（MATLAB 5中的三次多项式插值，不能用于外插）；如果x非等间距划分，则用spline方法；,Page 5,2019/10/12,机械工业出版社,Page 6, extrap：对超出范围的插值点进行外插, extrapval：对超出范围的插值点返回标量 extrapval，常取NaN和0, pp：用指定的方法产生Y的分段多项式形式（ ppform）。其中，method可取3）中除 v5cubic外的任一种；pp可由函数 ppval估计,函数interp1的功能及其各参数之间的关系如下图所示。,Page 6,2019/10/12,机械工业出版社,Page 7,Page 7,2019/10/12,机械工业出版社,Page 8,2）快速一维线性插值函数interp1q interp1q在对非均匀划分的数据进行插值时要比interp1快，因为它对输入不做检查。interp1q用线性插值方法对数据进行插值，其语法形式为：,yi = interp1q(x,Y,xi),x必须是单调增加的列向量,Y必须是一个列向量或一个行数为x的长度的矩阵,xi必须是一个列向量,Page 8,2019/10/12,机械工业出版社,Page 9,3）用FFT方法的一维插值函数interpft interpft用FFT（快速傅里叶变换）方法进行一维数据插值，即先将x用FFT方法转换到频域，然后再以更多的点转换回时域。语法形式为：,y = interpft(x,n)； y = interpft(x,n,dim),x为长度为m的向量，且采样间隔为dx时，则y的 采样间隔为dy=dx*m/n，此时，n必须大于等于m,x为一个矩阵时，interpft对x按列进行数据插值 ，返回矩阵y的行数为n，列数与x的列数相等,dim：interpft对x按指定维数进行操作,Page 9,2019/10/12,机械工业出版社,Page 10,2.二维插值和三维插值,1）二维插值函数interp2 如果已知点集是三维空间中的点，则相应的插值问题即是二维数据插值（表格查找）。MATLAB中用于二维数据插值的函数是interp2，其常用的语法形式及描述如下：,ZI = interp2(X,Y,Z,XI,YI) 或ZI = interp2(Z,XI,YI) ZI = interp2(Z,ntimes) ZI = interp2(X,Y,Z,XI,YI,method) ZI = interp2(.,method, extrapval),Page 10,2019/10/12,机械工业出版社,Page 11,X、Y：必须是严格单调的和具有栅格格式，一般由函数meshgrid产生，假定Z的大小为m, n，X、Y省略时，默认为X=1:n，Y=1:m,XI、YI：若为矩阵，返回【XI(i, j),YI(i, j)】处的Z值，对超出范围的插值点返回NaN；若为行向量xi和列向量yi，默认执行命令“meshgrid(xi, yi)”,ntimes：在任意两元素之间插入插值点以扩展Z，反复ntimes次，省略时默认为1,method：指定插值方法，取值有nearest、linear和spline（cubic和spline相同），要求X和Y必须单调，并且通常是由meshgrid函数产生的栅格格式,Page 11,2019/10/12,机械工业出版社,Page 12,extrapval：对超出已知点集(X, Y)的插值点执行外插操作，对任意不是由X或Y给出的XI或YI值，返回的ZI=extrapval。此时，必须给method参数赋值，默认为linear,函数interp2的功能及其各参数之间的关系如右图所示,Page 12,2019/10/12,机械工业出版社,Page 13,2）三维插值函数interp3 interp3的语法形式、实现的功能及参数的说明和二维插值函数interp2类似，读者可通过MATLAB的帮助系统了解。interp3的常用语法形式如下：,VI = interp3(X,Y,Z,V,XI,YI,ZI) VI = interp3(V,ntimes) VI = interp3(.,method,extrapval),Page 13,2019/10/12,机械工业出版社,Page 14,3.数据插值实例,用函数interp1（方法分别采用linear、nearest、spline和cubic）对正弦函数进行一维插值。,clc;clear all;close all; t=0:10;y=sin(t);tt=0:.25:10; y1=interp1(t,y,tt);y2=interp1(t,y,tt,nearest); y3=interp1(t,y,tt,spline);y4=interp1(t,y,tt,cubic); plot(t,y,o,tt,y1,-,tt,y2,-,tt,y3,:,tt,y4,-.); xlabel(t); ylabel(sin(t);title(对正弦信号进行一维插值); axis(-0.5 10.5 -1.2 1.2); legend(original data,linear,nearest,spline,cubic);,Page 14,2019/10/12,机械工业出版社,Page 15,Page 15,2019/10/12,机械工业出版社,Page 16,用函数interp1（方法分别采用linear、nearest、spline和cubic）对正弦函数进行一维外插运算。,clc;clear all;close all; t=0:10;y=sin(t); tt=0:.25:15; y1=interp1(t,y,tt,linear,extrap); y2=interp1(t,y,tt,nearest,extrap); y3=interp1(t,y,tt,spline,extrap); y4=interp1(t,y,tt,cubic,extrap); plot(t,y,o,tt,y1,-,tt,y2,-,tt,y3,:,tt,y4,-.); xlabel(t); ylabel(sin(t); title(对正弦信号进行一维外插); legend(original data,linear,nearest,spline,cubic);,Page 16,2019/10/12,机械工业出版社,Page 17,Page 17,2019/10/12,机械工业出版社,Page 18,用函数interpft对一个类三角信号进行5倍一维插值。,clc;clear all;close all; y=0 .5 1 1.5 2 1.5 1 .5 0 -.5 -1 -1.5 -2 -1.5 -1 -.5 0; N=length(y); L=5; M=N*L; x=0:L:L*N-1; xi=0:M-1; yi=interpft(y,M); plot(x,y,o,xi,yi,*); legend(Original data,Interpolated data);,Page 18,2019/10/12,机械工业出版社,Page 19,Page 19,2019/10/12,机械工业出版社,Page 20,二维插值函数interp2的应用：对peaks函数进行二维插值。,clc;clear all;close all; X,Y=meshgrid(-3:.25:3); Z=peaks(X,Y); XI,YI=meshgrid(-3:.125:3); ZI=interp2(X,Y,Z,XI,YI); mesh(X,Y,Z);hold on; mesh(XI,YI,ZI+15);hold off;axis(-3 3 -3 3 -5 20);,Page 20,2019/10/12,机械工业出版社,Page 21,Page 21,2019/10/12,机械工业出版社,Page 22,三维插值函数interp3的应用,clc;clear all;close all; x,y,z,v=flow(10); xi,yi,zi=meshgrid(.1:.25:10, -3:.25:3, -3:.25:3); vi=interp3(x,y,z,v,xi,yi,zi); slice(xi,yi,zi,vi,6 9.5,2,-2 .2); shading flat;,Page 22,2019/10/12,机械工业出版社,Page 23,Page 23,2019/10/12,机械工业出版社,Page 24,10.2 曲线拟合,与数据插值类似，曲线拟合也是数据分析中的常用方法，它依据一个区间或一个区域上的有限采样点值，构造一个近似函数（曲线）S来描述这组数据的内在规律。与数据插值要求逼近函数在采样点与被逼近函数相等不同，曲线拟合不要求拟合函数严格通过所有数据点，而使拟合函数与数据的总体误差最小。曲线拟合常采用的最优标准是最小二乘原理，所构造的函数是一个次数小于插值节点个数的多项式。曲线拟合在工程应用和科学研究中都有非常广泛的应用。,Page 24,2019/10/12,机械工业出版社,Page 25,1.多项式基础,在MATLAB中，多项式用由多项式系数（从最高阶次到最低阶次排列）构成的行向量来表示。例如：MATLAB中的向量p=1 2 3 4 5代表的多项式为,特殊地，如果多项式中某些项的系数为0，则在向量中用0表示，不能省略。例如：多项式,在MATLAB中可以用向量p=5 0 -2 0 0 7表示,Page 25,2019/10/12,机械工业出版社,Page 26,1）求多项式的根 已知多项式的系数向量求解多项式的根由函数roots完成，其语法形式为：,r=roots(c),c：一个行向量； r：一个列向量，其元素是c代表的多项式的根,2）由多项式的根创建多项式 多项式的创建由函数poly完成，其语法形式为,p = poly(A) 或p = poly(r),A：一个n行n列的矩阵，此时p是一个长为n+1的行向量 ，其元素为A的特征多项式的系数； r：一个列向量，此时p是一个行向量，其元素是以r的元 素为根的多项式的系数,Page 26,2019/10/12,机械工业出版社,Page 27, A =1 2 3;4 5 6;7 8 0 A = 1 2 3 4 5 6 7 8 0 p=poly(A) p = 1.0000 -6.0000 -72.0000 -27.0000 r=roots(p) r = 12.1229 -5.7345 -0.3884 pp=poly(r) pp = 1.0000 -6.0000 -72.0000 -27.0000,Page 27,2019/10/12,机械工业出版社,Page 28,X必须为一个方阵，函数将X作为一个整体代入多项式求值,Page 28,2019/10/12,机械工业出版社,Page 29, X = pascal(4) X = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 p = poly(X) p = 1.0000 -29.0000 72.0000 -29.0000 1.0000 r1=polyval(p,X) r1 = 1.0e+004 * 0.0016 0.0016 0.0016 0.0016 0.0016 0.0015 -0.0140 -0.0563 0.0016 -0.0140 -0.2549 -1.2089 0.0016 -0.0563 -1.2089 -4.3779 r2=polyvalm(p,X) r2 = 1.0e-010 * -0.0013 -0.0063 -0.0104 -0.0242 -0.0048 -0.0218 -0.0360 -0.0798 -0.0115 -0.0512 -0.0822 -0.1812 -0.0229 -0.0973 -0.1560 -0.3410,Page 29,2019/10/12,机械工业出版社,Page 30,4）多项式的运算 多项式的运算有乘法、除法、微分、部分分式展开等，分别由函数conv、deconv、polyder和residue完成，常用的语法形式为：,w = conv(u,v) q,r = deconv(v,u) k = polyder(p) 或k = polyder(a,b) r,p,k = residue(b,a) 或 b,a = residue(r,p,k),Page 30,2019/10/12,机械工业出版社,Page 31,功能：conv计算向量u和v的卷积，代数学上等同于u和v代表的两个多项式的乘积；deconv用长除法从v中解卷积出u，代数学上即为多项式相除，返回的向量q和r分别表示商和剩余多项式，关系式为v=conv(u,q)+r；polyder计算向量p或向量a和b代表的多项式的导数；r,p,k = residue(b,a)求多项式的部分分式展开形式，b,a = residue(r,p,k)将由r、p和k描述的部分分式形式转换为整体分式形式。其中：,Page 31,2019/10/12,机械工业出版社,Page 32,b、a：代表下式中分子多项式的系数和分母多项式的系数,r、p、k：代表下述部分分式的留数向量、极点向量和直接项系数向量,Page 32,2019/10/12,机械工业出版社,Page 33,在MATLAB命令窗口按下面步骤输入指令并运行 1）输入向量b和a代表两个多项式；, b = 5 3 -2 7; a = -4 0 8 3;,2）计算两个多项式的乘积；,q,r = deconv(c,b) q = -4 0 8 3 r = 1.0e-014 * 0 0 0 0 0.2220 0.4441 0.4441, c=conv(b,a) c = -20 -12 48 11 -7 50 21,3）用多项式c除以多项式b；,Page 33,2019/10/12,机械工业出版社,Page 34,4）求两个多项式乘积的导数；, d=polyder(b,a) d = -120 -60 192 33 -14 50,5）将二多项式构成的分式转化为部分分式形式；, r, p, k = residue(b,a) r = -1.4167 -0.6653 1.3320 p = 1.5737 -1.1644 -0.4093 k = -1.2500,Page 34,2019/10/12,机械工业出版社,Page 35,6）将5）的部分分式形式转换为整体分式形式。, bb,aa = residue(r,p,k) bb = -1.2500 -0.7500 0.5000 -1.7500 aa = 1.0000 -0.0000 -2.0000 -0.7500,向量bb和aa是以向量a的第一个元素为基准做归一化的结果。,Page 35,2019/10/12,机械工业出版社,Page 36,2.多项式曲线拟合应用实例,当曲线拟合所构造的拟合函数是多项式形式的函数时，即为多项式曲线拟合。MATLAB中用于多项式曲线拟合的函数是polyfit，常用的语法形式为：,p = polyfit(x,y,n),p是长为n+1的行向量，代表如下所示的n阶多项式p(x)，p(x(i)能在最小均方意义上拟合y(i)。,Page 36,2019/10/12,机械工业出版社,Page 37,多项式曲线拟合函数polyfit的应用。,本例是对误差函数的拟合，误差函数的定义如下,误差函数是一个界限函数，而多项式则是极大函数。因此，对误差函数的拟合并不一定能达到理想的效果，如下所示。,Page 37,2019/10/12,机械工业出版社,Page 38,clc;clear all;close all; x=(0:0.1:2.5); % 产生离散数据点并构造拟合函数 % y=erf(x); p=polyfit(x,y,6); f=polyval(p,x); % 计算拟合函数在原始数据点处的值并作比较 % table=x y f y-f x=(0:0.1:5); % 重新产生原始数据并求拟合函数在这些 数据点处的值 % y=erf(x); f=polyval(p,x); plot(x,y,o,x,f,-); % 画图观察拟合效果 % legend(原始数据,拟合曲线); axis(0 5 0 2);,Page 38,2019/10/12,机械工业出版社,Page 39,MATLAB命令窗口里的运行结果如下：,table = 0 0 0.0004 -0.0004 0.1000 0.1125 0.1119 0.0006 0.2000 0.2227 0.2223 0.0004 0.3000 0.3286 0.3287 -0.0001 0.4000 0.4284 0.4288 -0.0004 0.5000 0.5205 0.5209 -0.0004 0.6000 0.6039 0.6041 -0.0002 0.7000 0.6778 0.6778 0.0000 0.8000 0.7421 0.7418 0.0003 0.9000 0.7969 0.7965 0.0004 1.0000 0.8427 0.8424 0.0003 1.1000 0.8802 0.8800 0.0002,Page 39,2019/10/12,机械工业出版社,Page 40,1.2000 0.9103 0.9104 -0.0000 1.3000 0.9340 0.9342 -0.0002 1.4000 0.9523 0.9526 -0.0003 1.5000 0.9661 0.9664 -0.0003 1.6000 0.9763 0.9765 -0.0002 1.7000 0.9838 0.9838 0.0000 1.8000 0.9891 0.9889 0.0002 1.9000 0.9928 0.9925 0.0003 2.0000 0.9953 0.9951 0.0002 2.1000 0.9970 0.9969 0.0001 2.2000 0.9981 0.9982 -0.0001 2.3000 0.9989 0.9991 -0.0003 2.4000 0.9993 0.9995 -0.0002 2.5000 0.9996 0.9994 0.0002,Page 40,2019/10/12,机械工业出版社,Page 41,Page 41,2019/10/12,机械工业出版社,Page 42,10.3 微积分,微积分是高等数学中研究函数的微分、积分以及有关概念和应用的数学分支。它是数学的一个基础学科，内容包括积分、函数极限、级数求和、Taylor幂级数展开、Fourier级数展开、常微分方程等问题的直接求解。数值积分可用于计算解析定义的函数的积分，也可以计算以列表形式给出的函数的积分，是一种十分基础而且重要的运算；而很多和方程有关的数值计算，如插值与拟合、方程求根和微分方程求解都涉及微分的计算。MATLAB的符号数学工具箱提供了求各类积分和微分的函数，能够帮助工程师们解决工程应用中的各种求导和积分问题。,Page 42,2019/10/12,机械工业出版社,Page 43,1.积分,功能：计算expr（符号变量，由symvar定义）的不定积分或关于符号标量v的不定积分。, syms x; int(-2*x/(1 + x2)2) ans = 1/(x2 + 1) syms x z; int(x/(1+z2),z) ans = x*atan(z),Page 43,2019/10/12,机械工业出版社,Page 44,2）定积分 函数int用于计算定积分的语法形式如下：,int(expr,a,b) 或int(expr,v,a,b),功能：计算expr在a,b上的定积分或在a,b上关于v的定积分，a和b为符号标量或双精度标量。, syms x; int(x*log(1+x),0,1) ans = 1/4 int(1/(x3),1,inf) ans = 1/2 syms x t; int(2*x,sin(t),1) ans = cos(t)2,Page 44,2019/10/12,机械工业出版社,Page 45,3）二重积分 MATLAB提供了计算二重积分的函数dblquad，常用的语法形式为：,q=dblquad(fun,xmin,xmax,ymin,ymax) q=dblquad(fun,xmin,xmax,ymin,ymax,tol) q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method),功能：调用函数quad计算函数fun(x,y)在由xmin,xmax,ymin,ymax确定的矩形上的二重积分。,tol：用指定的精度tol代替默认精度,method：代替默认方法quad，取值为quadl或自定义的与quad和quadl有相同调用次序的函数。,Page 45,2019/10/12,机械工业出版社,Page 46,4）三重积分 MATLAB提供了计算三重积分的函数triplequad，常用的语法形式为：,triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax) triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol) triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol, method),功能：调用函数triplequad计算函数fun(x,y)在由xmin,xmax,ymin,ymax,zmin,zmax确定的矩形上的三重积分。,tol：用指定的精度tol代替默认精度,method：代替默认方法triplequad，取值为quadl或自定义的与quad和quadl有相同调用次序的函数。,Page 46,2019/10/12,机械工业出版社,Page 47,2.导数,1）导数 MATLAB的符号数学工具箱提供了求各类导数的函数diff，常用的语法形式为：,diff(expr) 或diff(expr, v) diff(expr, n) 或diff(expr, v, n),功能：求expr的导数、关于v的导数、n（正整数）阶导数和关于v的n（正整数）阶导数。, diff(sin(x2) ans = 2*x*cos(x2) diff(t6,6) ans = 720,Page 47,2019/10/12,机械工业出版社,Page 48,2）多元函数的偏导数 多元函数的偏导数也可以由函数diff实现。已知函数f(x,y)，若要求 ，则可以用下面的代码实现：,diff(diff(f,x,m),y,n) 或diff(diff(f,y,n),x,m), diff(sin(x*t2),t) ans = 2*t*x*cos(t2*x) diff(diff(3*x3*y2+sin(x*y),x),y) ans = cos(x*y) + 18*x2*y - x*y*sin(x*y),Page 48,2019/10/12,机械工业出版社,Page 49,-diff(f,xj)/diff(f,xi),已知隐函数 ，下面的代码用以计算导数, syms x y; f=exp(y)+y*sin(x)-exp(x); -diff(f,x)/diff(f,y) ans = (exp(x) - y*cos(x)/(exp(y) + sin(x),Page 49,2019/10/12,机械工业出版社,Page 50,diff(f,t)/diff(g,t), syms t; x=t*sin(t); y=t*(1-cos(t); diff(y,t)/diff(x,t) ans = (t*sin(t) - cos(t) + 1)/ (sin(t) + t*cos(t),Page 50,2019/10/12,机械工业出版社,Page 51,3.极限,MATLAB符号数学工具箱提供了计算函数极限的函数limit，常用的语法形式为：,limit(expr,x,a) 或limit(expr,a) 或limit(expr) limit(expr,x,a,left) 或limit(expr,x,a,right),功能：计算函数expr当xa时的双向极限、当默认变量趋向于a时的双向极限、当默认变量趋向于0时的双向极限、当x从左边趋向于a时的极限和当x从右边趋向于a时的极限。, syms x h a; limit(sin(x)/x) ans = 1 limit(1/x,x,0,right) ans = Inf limit(1/x,x,0,left) ans = -Inf v=(1+a/x)x,exp(-x); limit(v,x,inf) ans = exp(a), 0,Page 51,2019/10/12,机械工业出版社,Page 52,4.级数求和,MATLAB符号数学工具箱提供了求级数和的函数symsum，常用的语法形式为：,r=symsum(expr) 或r=symsum(expr,v) r=symsum(expr,a,b) 或r=symsum(expr,v,a,b),功能：计算符号表达式expr的级数，表达式的默认变量为由symvar定义的defaultVar，其变化范围为0defaultVar-1。,v：计算expr的级数，变量v的变化范围为0v-1； a、b：计算expr的级数，默认变量defaultVar的变化范围 为ab； v、a、b：计算expr的级数，变量v的变化范围为ab。,Page 52,2019/10/12,机械工业出版社,Page 53, syms k x; symsum(k) ans = k2/2 - k/2 symsum(k2) ans = k3/3 - k2/2 + k/6 symsum(k2,0,10) ans = 385 symsum(xk/sym(k!),k,0,inf) ans = exp(x),下面的代码用以求级数,Page 53,2019/10/12,机械工业出版社,Page 54,5.泰勒级数展开,MATLAB符号数学工具箱提供了用于泰勒级数展开的函数taylor，语法形式为：,taylor(f) 或taylor(f,n) 或taylor(f,a) 或taylor(f,n,v) 或taylor(f,n,v,a),功能：将f展开为5阶迈克劳林（Maclaurin）多项式、n-1阶的迈克劳林（Maclaurin）多项式（n为正整数）、关于a的5阶泰勒级数（a是一个实数）、关于v（表达式中的独立变量）的n-1阶迈克劳林多项式（v可以是一个字符串或符号变量）和关于v-a的n-1阶的泰勒级数（a可以是数值、字符或字符串，当a是字符或字符串时，v不能省略）。,Page 54,2019/10/12,机械工业出版社,Page 55,解析函数f(x)关于x-a的泰勒级数表达式如下：,则，函数taylor的不同语法形式所得出的结果如下表所示,Page 55,2019/10/12,机械工业出版社,Page 56,若函数f是一个二元或多元函数（f=f(x, y, )），则函数taylor的语法形式及得到的结果如下表所示。,Page 56,2019/10/12,机械工业出版社,Page 57,下面的代码用以求 在采用函数taylor的不同语法形式时的级数表达式。, syms x a; f=exp(x2); taylor(f) ans = x4/2 + x2 + 1 taylor(f,3) ans = x2 + 1 taylor(f,a) ans = exp(a2) + exp(a2)*(a - x)4*(2*a*(a/3 + 2*a*(a2/6 + 1/6) + (2*a2)/3 + 1/2) - xp(a2)*(a - x)5*(a/3 + 2*a*(a2/6 + 1/6) + 2*a*(2*a*(a/12 + 2*a*(a2/30 + 1/24) + a2/6 + 1/6) -exp(a2)*(a + 2*a*(2*a2)/3 + 1/2)*(a - x)3 - 2*a*exp(a2)*(a - x) + exp(a2)*(2*a2+ 1)*(a - x)2 taylor(f,3,1) ans = exp(1) + 2*exp(1)*(x - 1) + 3*exp(1)*(x - 1)2 taylor(f,3,x,a) ans = exp(a2) - 2*a*exp(a2)*(a - x) + exp(a2)*(2*a2 + 1)*(a - x)2,Page 57,2019/10/12,机械工业出版社,Page 58,10.4 线性方程组求解,在自然科学和工程技术中，很多问题可归结为求解线性方程组。采用MATLAB，不仅可以利用其提供的相关函数直接解决一些简单的线性方程组，也可以通过简单的编程来解决一些复杂的线性方程组。,Page 58,2019/10/12,机械工业出版社,Page 59,1.直接法,下面的代码用直接法求解如下方程组。, A=3 1 -4 1;1 -5 0 7; 0 2 1 1;1 6 -1 -4; b=13 -9 6 0; x=Ab x = 43.6667 -17.0000 20.3333 -19.6667,Page 59,2019/10/12,机械工业出版社,Page 60,2.矩阵求逆法,MATLAB中用于求逆矩阵的函数是inv。,下面的代码用矩阵求逆法求解上面的方程组（a）。, b=13 -9 6 0; x=inv(A)*b x = 43.6667 -17.0000 20.3333 -19.6667,Page 60,2019/10/12,机械工业出版社,Page 61,3.矩阵分解法,矩阵分解法是指根据一定的原理，用某种算法将方程组的系数矩阵分解为若干个矩阵的乘积，分解后的矩阵一般是某种特殊矩阵。常见的矩阵分解方法有：LU分解、QR分解、Cholesky分解、Schur分解、Hessenberg分解和奇异分解等。这里只对LU分解、QR分解以及Cholesky稍加介绍，请读者通过MATLAB的帮助命令查看其他几种矩阵分解法的用法。,Page 61,2019/10/12,机械工业出版社,Page 62,U、L（行交换）P、：上三角矩阵、变换下三角矩阵和置 换矩阵。,Page 62,2019/10/12,机械工业出版社,Page 63,下面的代码用LU分解法求解方程组（a）。, A=3 1 -4 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4; b=13 -9 6 0; L,U=lu(A); x=U(Lb) x = 43.6667 -17.0000 20.3333 -19.6667 L,U,P=lu(A); x=U(LP*b) x = 43.6667 -17.0000 20.3333 -19.6667,Page 63,2019/10/12,机械工业出版社,Page 64,2）QR分解 矩阵的QR分解是将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。MATLAB提供了用于对矩阵进行QR分解的函数qr，其常用的语法形式为：,Q,R = qr(A) 或 Q,R,E = qr(A),若矩阵A为m行n列，则函数返回m行n列的上三角矩阵R和m行m列的正交矩阵Q，A=Q*R；第二种语法形式返回正交矩阵Q、上三角矩阵R和置换矩阵E，A*E=Q*R。,Page 64,2019/10/12,机械工业出版社,Page 65,下面的代码用QR分解法求解方程组（a）。, A=3 1 -4 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4; b=13 -9 6 0; Q,R=qr(A); x=R(Qb) x = 43.6667 -17.0000 20.3333 -19.6667 Q,R,E=qr(A); x=E*(R(Qb) x = 43.6667 -17.0000 20.3333 -19.6667,Page 65,2019/10/12,机械工业出版社,Page 66,3）Cholesky分解 如果A是正定对称矩阵，Cholesky分解将其分解为一个上三角矩阵R和一个下三角矩阵R的乘积，R是R的转置，A= RR。MATLAB提供了用于对矩阵进行Cholesky分解的函数chol，其常用的语法形式为：,R = chol(A) 或 R,p = chol(A),产生上三角矩阵R和下三角矩阵R，R*R=A，矩阵A必须是正定的，否则，MATLAB将返回错误信息。,Page 66,2019/10/12,机械工业出版社,Page 67,下面的代码用Cholesky分解法求解下面方程组（b）。, A=4 -1 1;-1 4.25 2.75;1 2.75 3.5; b=13 -9 6; R=chol(A); x=R(Rb) x = -0.1992 -6.7344 7.0625,Page 67,2019/10/12,机械工业出版社,Page 68,若对方程组（a）用Cholesky分解法求解，由于系数矩阵非正定，MATLAB会提示错误信息，代码如下：, A=3 1 -4 1;1 -5 0 7;0 2 1 -1;1 6 -1 -4; b=13 -9 6 0; R=chol(A); ? Error using = chol Matrix must be positive definite.,Page 68,2019/10/12,机械工业出版社,Page 69,4.迭代法,迭代法是将求一组解转换为求一个近似解序列的过程，并用最终的近似解来逼近真实解。迭代法速度快，非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Seidel迭代法、超松弛迭代法和两步迭代法等。,Page 69,2019/10/12,机械工业出版社,Page 70,1）Jacobi迭代法 若系数矩阵A的主对角线元素全不为0，则可将A分解为,其中D为对角阵，其元素为A的对角线元素，L与U为A的下三角阵和上三角阵。,则线性方程组 化为,与之相应的迭代公式为,此即Jacobi迭代公式，如果序列x(k+1)收敛于x，则x必是方程组的解。,Page 70,2019/10/12,机械工业出版社,Page 71,按照上述原理，编写用迭代法解线性方程组的函数Jacobi。,function x,n=Jacobi(A,b,x0,eps) if nargin=3 eps=1.0e-6; elseif nargin=eps x0=x; x=B*x0+f; n=n+1; end,Page 71,2019/10/12,机械工业出版社,Page 72,下面的代码用Jacobi迭代法求解方程组（b）。, A=4 -1 1;-1 4.25 2.75;1 2.75 3.5; b=13 -9 6; x,n=Jacobi(A,b,zeros(3,1),1.0e-6) x = -0.1992 -6.7344 7.0625 n = 97,Page 72,2019/10/12,机械工业出版社,Page 73,2）Gauss-Seidel迭代法 将Jacobi的迭代公式改为,即有,此即Gauss-Serdel迭代公式。与Jacobi迭代相比，Gauss-Seidel迭代用新分量代替旧分量，精度高。,Page 73,2019/10/12,机械工业出版社,Page 74,按照上述原理，编写用Gauss-Seidel迭代法求解线性方程组的函数GaussSeidel。,function x,n=GaussSeidel(A,b,x0,eps) if nargin=3 eps=1.0e-6; elseif nargin=eps x0=x; x=G*x0+f; n=n+1; end,Page 74,2019/10/12,机械工业出版社,Page 75,下面的代码用Gauss-Seidel迭代法求解方程组（b）。, A=4 -1 1;-1 4.25 2.75;1 2.75 3.5; b=13 -9 6; x,n=GaussSeidel(A,b,ones(3,1),1.0e-6) x = -0.1992 -6.7344 7.0625 n = 52,Page 75,2019/10/12,机械工业出版社,Page 76,10.5 非线性方程组求解,实际中遇到的多是线性方程（组），非线性方程（组）的理论远不如线性方程（组）成熟，特别是非线性方程组的解存在的唯一性还未完全解决，判断其存在性和解的个数几乎没有可行的方法。利用MATLAB提供的函数，可以求解一些简单的非线性方程，通过编程也可以解决一些较为复杂的非线性方程。,Page 76,2019/10/12,机械工业出版社,Page 77,1.非线性方程的求根,1）fzero函数 fzero用于求单变量非线性方程，常用的语法形式为,x = fzero(fun,x0) x = fzero(fun,x0,options) x,fval = fzero(.),功能：x = fzero(fun,x0)求函数fun在x0附近的根；x = fzero(fun,x0,options)使用优化选项求函数fun在x0附近的根，options的取值可通过MATLAB的帮助命令查看；x,fval = fzero(.)返回目标函数fun在解x处的值，fval=f(x)。,Page 77,2019/10/12,机械工业出版社,Page 78,下面的代码用函数fzero求非线性方程 在x=0.5附近的根。, fzero(x2+x-1,0.5) ans = 0.6180,fzero还可以求函数在某个区间上的根。下面的代码用于求方程 在区间0, 2上的根。, fzero(sin(x)-0.5,0,2) ans = 0.5236,当某个多元函数中的其他元是确定的时，fzero也可以用于求函数在未知变量的某个点附近的根。下面的代码用于求函数 当y=2时在x=2附近的根。, u=inline(xy-3,x,y); fzero(u,2,2) ans = 1.7321,Page 78,2019/10/12,机械工业出版社,Page 79,2） fsolve函数 fsolve用于求解多变量非线性方程组，其常用的语法形式为：,x = fsolve(fun,x0) x = fsolve(fun,x0,options) x,fval = fsolve(fun,x0),功能：x=fsolve(fun,x0)求解非线性方程fun 初始值为x0；x=fsolve(fun,x0,options)使用优化选项解方程；x,fval=fsolve(fun,x0)返回目标函数fun在x处的值，fval=F(x)。,Page 79,2019/10/12,机械工业出版社,Page 80,首先，创建函数myfun10_7.m；,function F=myfun10_7(p) x=p(1);y=p(2); F(1)=x-0.4*cos(x)-0.3*sin(y); F(2)=y-0.6*cos(x)+0.5*sin(y);,然后，在MATLAB的命令窗口输入以下指令, fsolve(myfun10_7,0.5,0.4,optimset(Display,off) ans = 0.4636 0.3604,Page 80,2019/10/12,机械工业出版社,Page 81,2.非线性方程组的数值解法,1）不动点迭代法 由不动点迭代法求非线性方程组数值解的原理编写函数stablepoint，代码如下：,% stablepoint.m % function r,n=stablepoint(x0,eps) if nargin=1 eps=1.0e-4; end r=x0+myfunSP(x0); n=1; tol=1; while toleps x0=r; r=x0+myfunSP(x0); tol=norm(r-x0); n=n+1; if(n100000) disp(迭代步数太多，可能不收敛！); return; end end,Page 81,2019/10/12,机械工业出版社,Page 82,用不动点迭代法求解非线性方程组,初始迭代值为（0, 0）。,首先，建立函数myfunSP,function f=myfunSP(x) f(1)=0.5*sin(x(1)+0.1*cos(x(2)*x(1)-x(1); f(2)=0.5*cos(x(1)-0.1*sin(x(2)-x(2); f=f(1) f(2);,然后，在MATLAB的命令窗口输入下面的指令。, r,n=stablepoint(0,0) r = 0.1979 0.4470 n = 11,Page 82,2019/10/12,机械工业出版社,Page 83,2）牛顿法 由牛顿法求非线性方程组数值解的原理编写函数newton，代码如下：,% newton.m % function r,n=newton(x0,eps) if nargin=1 eps=1.0e-4; end r=x0-myfunNT(x0)/dmyfunNT(x0); n=1; tol=1; while toleps x0=r; r=x0-myfunNT(x0)/dmyfunNT(x0); tol=norm(r-x0); n=n+1; if(n100000) disp(迭代步数太多，可能不收敛！); return; end end,Page 83,2019/10/12,机械工业出版社,Page 84,用牛顿法求解非线性方程组（c），初始迭代值为（0, 0）。,首先，建立函数myfunNT,function f=myfunNT(x) f(1)=0.5*sin(x(1)+0.1*cos(x(2)*x(1)-x(1); f(2)=0.5*cos(x(1)-0.1*sin(x(2)-x(2); f=f(1) f(2);,再建立函数dmyfunNT,function df=dmyfunNT(x) df=0.5*cos(x(1)-0.1*x(2)*sin(x(2)*x(1)-1,- 0.1*x(1)*sin(x(2)*x(1);-0.5*sin(x(1),. -0.1*cos(x(2)-1;,在MATLAB的命令窗口输入, r,n=newton(0,0) r = 0.1979 0.4470 n = 5,Page 84,2019/10/12,机械工业出版社,Page 85,3）牛顿下降法 由牛顿下降法求非线性方程组数值解的原理编写函数newtonD，代码如下：,% newtonD.m % function r,n=newtonD(x0,eps) if nargin=1 eps=1.0e-4; end r=x0-myfunNTD(x0)/dmyfunNTD(x0); n=1; tol=1; while toleps x0=r; ttol=1; w=1; F1=norm(myfunNTD(x0); while ttol=0 r=x0-w*myfunNTD(x0)/dmyfunNTD(x0); ttol=norm(myfunNTD(r)-F1; w=w/2; end tol=norm(r-x0); n=n+1; if(n100000) disp(迭代步数太多，可能不收敛！); return; end end,Page 85,2019/10/12,机械工业出版社,Page 86,用牛顿下降法求解非线性方程组(c)，初始迭代值为(0, 0)。,首先，建立函数myfunNTD,function f=myfunNT(x) f(1)=0.5*sin(x(1)+0.1*cos(x(2)*x(1)-x(1); f(2)=0.5*cos(x(1)-0.1*sin(x(2)-x(2); f=f(1) f(2);,再建立函数dmyfunNTD,function df=dmyfunNT(x) df=0.5*cos(x(1)-0.1*x(2)*sin(x(2)*x(1)-1,- 0.1*x(1)*sin(x(2)*x(1);-0.5*sin(x(1),. -0.1*cos(x(2)-1;,在MATLAB的命令窗口输入, r,n=newtonD(0,0) r = 0.1