下载文档原格式
MATLAB的基本使用教程MATLAB是一种强大的数学计算软件，广泛应用于科学、工程和技术领域。
它提供了丰富的功能和工具，能够快速、有效地处理和分析各种数学问题。
本文将介绍MATLAB的基本使用方法，帮助初学者快速入门。
一、MATLAB的安装与启动1、下载和安装MATLAB软件：在MathWorks官方网站上下载适合自己操作系统的MATLAB软件，并根据安装提示进行安装。
安装完成后，会生成一个MATLAB的启动图标。
2、启动MATLAB：双击MATLAB的启动图标，或者在命令行中输入"matlab"命令，即可启动MATLAB。
二、MATLAB的基本操作1、工作环境：MATLAB提供了一个强大的集成开发环境（IDE），可以在其中编写和运行代码。
在MATLAB的界面中，包括主窗口、命令窗口、变量窗口、编辑器等。
2、命令窗口：在命令窗口中可以输入和执行MATLAB命令。
可以直接在命令窗口中输入简单的计算，例如输入"2+3"并按下回车键，即可输出计算结果。
3、脚本文件：MATLAB可以编写和运行脚本文件，将一系列命令组织起来，并按顺序执行。
在编辑器中编写MATLAB代码，并将文件保存为.m扩展名的脚本文件。
然后在命令窗口中输入脚本文件的文件名（不带扩展名），按下回车键即可执行脚本文件中的代码。
4、变量和赋值：在MATLAB中，可以创建和操作各种类型的变量。
例如，可以使用"="符号将一个值赋给一个变量，例如"A=5"。
在后续的计算和分析中，可以使用这个变量，例如输入"B=A+3"，结果B 将被赋值为8。
5、矩阵和向量：MATLAB中的基本数据结构是矩阵和向量。
可以使用方括号[]来创建矩阵和向量，并使用逗号或空格来分隔不同的元素。
例如，"[1,2,3]"表示一个包含3个元素的行向量。
6、矩阵运算：MATLAB提供了丰富的矩阵运算符和函数，可以对矩阵进行各种运算。
matlab_简明实例教程MATLAB是一种强大的科学计算工具，广泛应用于科学研究、数据分析和工程计算等领域。
它具有简单易用的语法和丰富的函数库，可以快速实现复杂的计算任务。
本教程将为你提供一些简单实例，帮助你入门MATLAB。
1.计算圆的面积和周长```matlabradius = input('请输入圆的半径：');area = pi * radius^2;circumference = 2 * pi * radius;disp(['圆的面积为：', num2str(area)]);disp(['圆的周长为：', num2str(circumference)]);```2.计算两个向量的点积```matlabv1 = input('请输入向量1（用逗号分隔元素）：');v2 = input('请输入向量2（用逗号分隔元素）：');dot_product = dot(v1, v2);disp(['两个向量的点积为：', num2str(dot_product)]);```3.绘制正弦曲线```matlabx = 0:0.1:2*pi;y = sin(x);plot(x, y);xlabel('x');ylabel('sin(x)');title('正弦曲线');```4.求解方程```matlabsyms x;eqn = x^2 - 4 == 0;sol = solve(eqn, x);disp(['方程的解为：', char(sol)]); ```5.读取和写入文件```matlabfilename = 'data.txt';data = importdata(filename);disp('文件中的数据：');disp(data);output = [1 2 3; 4 5 6; 7 8 9];dlmwrite('result.txt', output, 'delimiter', '\t', 'precision', 4);disp('结果已保存到result.txt文件中。
实验一特殊函数与图形一、问题背景与实验目的二、相关函数（命令）及简介三、实验内容四、自己动手一、问题背景与实验目的著名的Riemann函数大家都很熟悉了，但是关于它的图像你是否清楚呢除了最上面那几点，其他都很难画吧你想不想看看下面那些“挤在一起”的点是怎样分布的呢还有几何中的马鞍面、单叶双曲面等是怎样由直线生成的，是不是也想目睹一下呢这些，都离不开绘图．实际上绘图一直是数学中的一种重要手段，借助图形，往往可以化繁为简，使抽象的对象得到明白直观的体现．比如函数的基本性质，一个图形常可以使之一目了然，非常有效．它虽不能代替严格的分析与证明，但在问题的研究过程中，可以帮助研究人员节约相当一部分精力．此外，它还可以使计算、证明、建模等的结果得到更明白易懂的表现，有时，这比科学论证更有说服力．同时，数学的教学与学习过程也离不开绘图．借助直观的图形，常可以使初学者更容易接受新知识．如数学分析中有不少函数，其解析式着实让人望而生畏，即使对其性质作了详尽的分析，还是感到难明就里；但如果能看到它的图形，再配合理论分析，则问题可以迎刃而解．又如在几何的学习中，会遇到大量的曲线与曲面，也离不开图形的配合．传统的手工作图，往往费力耗时，效果也不尽理想．计算机恰恰弥补了这个不足，使你可以方便地指定各种视角、比例、明暗，从各个角度进行观察．本实验通过对函数的图形表示和几个曲面（线）图形的介绍，一方面展示它们的特点，另一方面，也将就Matlab软件的作图功能作一个简单介绍．大家将会看到，Matlab 的作图功能非常强大．二、相关函数（命令）及简介1．平面作图函数：plot，其基本调用形式：plot(x,y,s)以x作为横坐标，y作为纵坐标．s是图形显示属性的设置选项．例如：x=-pi:pi/10:pi;y=sin(x);plot(x,y,'--rh','linewidth',2,'markeredgecolor','b','markerfacecolor','g')图1在使用函数plot时，应当注意到当两个输入量同为向量时，向量x与y必须维数相同，而且必须同是行向量或者同是列向量．绘图时，可以制定标记的颜色和大小，也可以用图形属性制定其他线条特征，这些属性包括：linewidth 指定线条的粗细．markeredgecolor 指定标记的边缘色markerfacecolor 指定标记表面的颜色．markersize 指定标记的大小．若在一个坐标系中画几个函数，则plot的调用格式如下：plot(x1,y1,s1,x2,y2,s2,……)2．空间曲线作图函数：plot3，它与plot相比，只是多了一个维数而已．其调用格式如下：plot3(x,y,z,s)．例如：x=0:pi/30:20*pi;y=sin(x);z=cos(x);plot3(x,y,z)得到三维螺旋线：图23．空间曲面作图函数：（1）mesh函数．绘制彩色网格面图形．调用格式：mesh(z)，mesh(x,y,z)和mesh(x,y,z,c)．其中，mesh(x,y,z,c)画出颜色由c指定的三维网格图．若x、y均为向量，则length(x)=n，length(y)=m，[m,n]=size(z)．（2）surf在矩形区域内显示三维带阴影曲面图．调用格式与mesh类似．（3）ezmesh用符号函数作三维曲面网格图．调用格式：ezmesh(x,y,z)其中x = x(s,t), y = y(s,t),z = z(s,t)．画图区域默认为： -2*pi < s < 2*pi 且-2*pi < t < 2*pi.或者用格式：ezmesh(x,y,z,[smin,smax,tmin,tmax])（4）ezsurf用符号函数作三维曲面图．调用格式与ezmesh类似．（5）sphere画球体命令．4．meshgrid，调用格式：[x,y]=meshgrid(m,n)，这里的m，n为给定的向量，可以定义网格划分区域和划分方法．矩阵x和矩阵y是网格划分后的数据矩阵．5．图像的修饰与其他函数：（1）axis equal 控制各个坐标轴的分度，使其相等；（2）colormap设置绘图颜色．调用格式：colormap([r g b])其中r,g,b都是0-1之间的数．或者用格式：colormap(s)s颜色映像相应的颜色系颜色映像相应的颜色系autumn红黄色系hsv色调饱和色系gray线性灰色系hot黑红黄白色系cool青和洋红色系pink柔和色系（3（4）find找出符合条件的元素在数组中的位置．调用格式：y=find(条件)例如：输入：a=[4 5 78 121 4 665 225 4 1];b=find(a>7)输出： b =3 4 6 7三、实验内容数学分析中，特别是积分部分，我们接触了不少有趣的函数，由于其中有的不是一一对应的，用上面的方法无法画出它们的图像，这时就只能用参数了．此外还有些图形只能用参数来画，比如空间曲线，在计算机上不接受“两个曲面的交线”这种表示，所以也只能用参数来实现．用参数方式作图的关键在于找出合适的参数表示，尤其是不能有奇点，最好也不要用到开方．所以要找的参数最好是有几何意义的．当然这也不可一概而论，需要多积累经验．1．利用函数plot在一个坐标系中画以下几个函数图像，要求采用不同颜色、不同线形、不同的符号标记．函数为：．程序如下：t=0:pi/20:2*pi;x=sin(t);y=cos(t);z=sin(2*t);plot(t, x, '--k*', t, y, '-rs', t, z, ':bo')图像如下：图32．绘制类似田螺线的一条三维螺线(方程自己设计)．程序如下：t=0:.1:30;x=2*(cos(t)+t.*sin(t));y=2*(sin(t)-t.*cos(t));z=*t;plot3(x,y,-z) %取–z 主要是为了画图看起来更清楚axis equal图像如下：图43．利用函数，绘制一个墨西哥帽子的图形．程序如下：[a,b]=meshgrid(-8:.5:8); %先生成一个网格c=sqrt(a.^2+b.^2)+eps;z=sin(c)./c;mesh(a,b,z)axis square图像如下：图5思考：这里的 eps 是什么其作用是什么4．利用surf绘制马鞍面图形（函数为：）．程序如下：[x,y]=meshgrid(-25:1:25,-25:1:25);z=x.^2/9-y.^2/4;surf(x,y,z)title('马鞍面')grid off图像如下：图65．分别用ezmesh和ezsurf各绘制一个圆环面，尝试将两个圆环面放在一个图形界面内，观察它们有什么不同之处．提示：圆环面的方程为：，而圆环面的参数方程为：程序参见附录1．图像如下：图76．绘制黎曼函数图形，加深对黎曼函数的理解．说明：黎曼函数的定义为程序参见附录2．图像如下：图8四、自己动手1．作出下图所示的三维图形：图9提示：图形为圆环面和球面的组合.2．作出下图所示的墨西哥帽子及其剪裁图形：图103．画出球面、椭球面、双叶双曲面、单叶双曲面．4．若要求田螺线的一条轴截面的曲边是一条抛物线：时．试重新设计田螺线的参数方程，并画出该田螺线．5．作出下图所示的马鞍面（颜色为灰色，并有一个标题：“马鞍面”）：图116．绘制图8所示的黎曼函数图形，要求分母的最大值的数值由键盘输入（提示：使用input语句）．回目录下一页实验二定积分的近似计算一、问题背景与实验目的二、相关函数（命令）及简介三、实验内容1．矩形法2．梯形法3．抛物线法4. 直接应用Matlab命令计算结果四、自己动手一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，Matlab有专门函数可用．二、相关函数（命令）及简介1．sum(a)：求数组a的和．2．format long：长格式，即屏幕显示15位有效数字．（注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）．3．double()：若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值．4．quad()：抛物线法求数值积分．格式： quad(fun,a,b) ，注意此处的fun是函数，并且为数值形式的，所以使用*、/、^等运算时要在其前加上小数点，即 .*、./、.^等．例：Q = quad('1./(x.^3-2*x-5)',0,2);5．trapz()：梯形法求数值积分．格式：trapz(x,y)其中x为带有步长的积分区间；y为数值形式的运算（相当于上面介绍的函数fun）例：计算x=0:pi/100:pi;y=sin(x);trapz(x,y)6．dblquad()：抛物线法求二重数值积分．格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun可以用inline定义，也可以通过某个函数文件的句柄传递．例1：Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2，通过计算，比较Q1 与Q2结果（或加上手工验算），找出积分变量x、y的上下限的函数代入方法．Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi)例2：Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件：function z = integrnd(x, y)z = y*sin(x);7．fprintf（文件地址，格式，写入的变量）：把数据写入指定文件．例：x = 0:.1:1;y = [x; exp(x)];fid = fopen('','w'); %打开文件fprintf(fid,'% %\n',y); %写入fclose(fid) %关闭文件8．syms 变量1 变量2 …：定义变量为符号．9．sym('表达式')：将表达式定义为符号．解释：Matlab中的符号运算事实上是借用了Maple的软件包，所以当在Matlab中要对符号进行运算时，必须先把要用到的变量定义为符号．10．int(f,v,a,b)：求f关于v积分，积分区间由a到b．11．subs(f，'x'，a)：将 a 的值赋给符号表达式 f 中的 x，并计算出值．若简单地使用subs(f)，则将f的所有符号变量用可能的数值代入，并计算出值．三、实验内容1．矩形法根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度．针对不同的取法，计算结果会有不同，我们以为例（取），（1）左点法：对等分区间，在区间上取左端点，即取，，理论值，此时计算的相对误差（2）右点法：同（1）中划分区间，在区间上取右端点，即取，，理论值，此时计算的相对误差（3）中点法：同（1）中划分区间，在区间上取中点，即取，，理论值，此时计算的相对误差如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物线法就是这一指导思想的产物．2．梯形法等分区间，相应函数值为（）．曲线上相应的点为（）将曲线的每一段弧用过点，的弦（线性函数）来代替，这使得每个上的曲边梯形成为真正的梯形，其面积为，．于是各个小梯形面积之和就是曲边梯形面积的近似值，，即，称此式为梯形公式．仍用的近似计算为例，取，，理论值，此时计算的相对误差很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多．3．抛物线法由梯形法求近似值，当为凹曲线时，它就偏小；当为凸曲线时，它就偏大．若每段改用与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法．将积分区间作等分，分点依次为，，对应函数值为（），曲线上相应点为（）．现把区间上的曲线段用通过三点，，的抛物线来近似代替，然后求函数从到的定积分：由于，代入上式整理后得同样也有……将这个积分相加即得原来所要计算的定积分的近似值：，即这就是抛物线法公式，也称为辛卜生（Simpson）公式．仍用的近似计算为例，取，=，理论值，此时计算的相对误差4. 直接应用Matlab命令计算结果（1）数值计算方法1：int('1/(1+x^2)','x',0,1) (符号求积分)方法2：quad('1./(1+x.^2)',0,1) （抛物线法求数值积分）方法3：x=0::1;y=1./(1+x.^2);trapz(x,y) （梯形法求数值积分）（2）数值计算方法1：int(int('x+y^2','y',-1,1),'x',0,2) （符号求积分）方法2：dblquad(inline('x+y^2'),0,2,-1,1) （抛物线法二重数值积分）四、自己动手1．实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算，取，并比较三种方法的精确程度．2．分别用梯形法与抛物线法，计算，取．并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异．3．试计算定积分．（注意：可以运用trapz()、quad()或附录程序求解吗为什么）4．将的近似计算结果与Matlab中各命令的计算结果相比较，试猜测Matlab中的数值积分命令最可能采用了哪一种近似计算方法并找出其他例子支持你的观点．5．通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型的函数（具有某种单调特性或凹凸特性），用某种近似计算方法所得结果更接近于实际值6．学习的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环．上一页回目录下一页实验三求代数方程的近似根（解）一、问题背景和实验目的二、相关函数（命令）及简介三、实验内容四、自己动手一、问题背景和实验目的求代数方程的根是最常见的数学问题之一（这里称为代数方程，主要是想和后面的微分方程区别开．为简明起见，在本实验的以下叙述中，把代数方程简称为方程），当是一次多项式时，称为线性方程，否则称之为非线性方程．当是非线性方程时，由于的多样性，尚无一般的解析解法可使用，但如果对任意的精度要求，能求出方程的近似根，则可以认为求根的计算问题已经解决，至少能满足实际要求．本实验介绍一些求方程实根的近似值的有效方法，要求在使用这些方法前先确定求根区间，或给出某根的近似值．在实际问题抽象出的数学模型中，可以根据物理背景确定；也可根据的草图等方法确定，还可用对分法、迭代法以及牛顿切线法大致确定根的分布情况．通过本实验希望你能：1. 了解对分法、迭代法、牛顿切线法求方程近似根的基本过程；2. 求代数方程（组）的解．二、相关函数（命令）及简介1．abs( )：求绝对值函数．2．diff(f)：对独立变量求微分，f 为符号表达式．diff(f, 'a')：对变量a求微分，f 为符号表达式．diff(f, 'a', n)：对变量 a 求 n 次微分，f 为符号表达式．例如：syms x tdiff(sin(x^2)*t^6, 't', 6)ans=720*sin(x^2)3．roots([c(1), c(2), …, c(n+1)])：求解多项式的所有根．例如：求解：．p = [1 -6 -72 -27];r = roots(p)r =4．solve('表达式')：求表达式的解．solve('2*sin(x)=1')ans =1/6*pi5．linsolve(A, b)：求线性方程组 A*x=b 的解．例如：A= [9 0; -1 8]; b=[1; 2];linsolve(A, b)ans=[ 1/9][19/72]6．fzero(fun, x0)：在x0附近求fun 的解．其中fun为一个定义的函数，用“@函数名”方式进行调用．例如：fzero(@sin, 3)ans=7．subs(f, 'x ', a)：将 a 的值赋给符号表达式 f 中的 x，并计算出值．例如：subs('x^2 ', 'x ', 2)ans = 4三、实验内容首先，我们介绍几种与求根有关的方法：1．对分法对分法思想：将区域不断对分，判断根在某个分段内，再对该段对分，依此类推，直到满足精度为止．对分法适用于求有根区间内的单实根或奇重实根．设在上连续，，即，或，．则根据连续函数的介值定理，在内至少存在一点，使．下面的方法可以求出该根：(1)令，计算；(2)若，则是的根，停止计算，输出结果．若，则令，，若，则令，；．……，有、以及相应的．(3) 若 (为预先给定的精度要求)，退出计算，输出结果；反之，返回(1)，重复(1)，(2)，(3)．以上方法可得到每次缩小一半的区间序列，在中含有方程的根．当区间长很小时，取其中点为根的近似值，显然有以上公式可用于估计对分次数．分析以上过程不难知道，对分法的收敛速度与公比为的等比级数相同．由于，可知大约对分10次，近似根的精度可提高三位小数．对分法的收敛速度较慢，它常用来试探实根的分布区间，或求根的近似值．2. 迭代法1)迭代法的基本思想：由方程构造一个等价方程从某个近似根出发，令，可得序列，这种方法称为迭代法．若收敛，即，只要连续，有即可知，的极限是的根，也就是的根．当然，若发散，迭代法就失败．以下给出迭代过程收敛的一些判别方法：定义：如果根的某个邻域中，使对任意的，迭代过程，收敛，则称迭代过程在附近局部收敛．定理1：设，在的某个邻域内连续，并且，，则对任何，由迭代决定的序列收敛于．定理2：条件同定理 1，则定理3：已知方程，且(1) 对任意的，有．(2) 对任意的，有，则对任意的，迭代生成的序列收敛于的根，且．以上给出的收敛定理中的条件要严格验证都较困难，实用时常用以下不严格的标准：当根区间较小，且对某一，明显小于1时，则迭代收敛(参见附录3)．2) 迭代法的加速：a) 松弛法：若与同是的近似值，则是两个近似值的加权平均，其中称为权重，现通过确定看能否得到加速．迭代方程是：其中，令，试确定：当时，有，即当，时，可望获得较好的加速效果，于是有松弛法：，松弛法的加速效果是明显的 (见附录4)，甚至不收敛的迭代函数经加速后也能获得收敛．b) Altken方法：松弛法要先计算，在使用中有时不方便，为此发展出以下的 Altken 公式：，是它的根，是其近似根．设，，因为，用差商近似代替，有，解出，得由此得出公式；；，这就是Altken 公式，它的加速效果也是十分明显的，它同样可使不收敛的迭代格式获得收敛(见附录5)．3. 牛顿(Newton)法（牛顿切线法）1) 牛顿法的基本思想：是非线性方程，一般较难解决，多采用线性化方法．记：是一次多项式，用作为的近似方程．的解为记为，一般地，记即为牛顿法公式.2) 牛顿法的收敛速度：对牛顿法，迭代形式为：注意分子上的，所以当时，，牛顿法至少是二阶收敛的，而在重根附近，牛顿法是线性收敛的．牛顿法的缺点是：(1)对重根收敛很慢；(2)对初值要求较严，要求相当接近真值．因此，常用其他方法确定初值，再用牛顿法提高精度．4. 求方程根(解)的其它方法(1) solve('x^3-3*x+1=0')(2) roots([1 0 -3 1])(3) fzero('x^3-3*x+1', -2)(4) fzero('x^3-3*x+1',(5) fzero('x^3-3*x+1',(6) linsolve([1, 2, 3; 4, 5, 6; 7, 8, 0], [1, 2, 3]')体会一下，(2)(5) 用了上述 1 3 中的哪一种方法以下是本实验中的几个具体的实验，详细的程序清单参见附录．具体实验1：对分法先作图观察方程：的实根的分布区间，再利用对分法在这些区间上分别求出根的近似值．输入以下命令，可得的图象：f='x^3-3*x+1';g='0';ezplot(f, [-4, 4]);hold on;ezplot(g, [-4, 4]); %目的是画出直线 y=0，即 x 轴grid on;axis([-4 4 -5 5]);hold off请填写下表：实根的分布区间该区间上根的近似值在某区间上求根的近似值的对分法程序参见附录1．具体实验2：普通迭代法采用迭代过程：求方程在附近的根，精确到第 4 位小数．构造等价方程：用迭代公式：，用 Matlab 编写的程序参见附录2．请利用上述程序填写下表：分析：将附录2第4行中的分别改为以及，问运行的结果是什么你能分析得到其中的原因吗看看下面的“具体实验3”是想向你表达一个什么意思．用 Matlab 编写的程序参见附录3．具体实验3：收敛/发散判断设方程的三个根近似地取，和，这些近似值可以用上面的对分法求得．迭代形式一：收敛 (很可能收敛，下同)不收敛 (很可能不收敛，下同)不收敛迭代形式二：收敛不收敛不收敛迭代形式三：不收敛收敛收敛具体实验4：迭代法的加速1——松弛迭代法，，迭代公式为程序参见附录4．具体实验5：迭代法的加速2——Altken迭代法迭代公式为：，，程序参见附录5．具体实验6：牛顿法用牛顿法计算方程在-2到2之间的三个根．提示：，迭代公式：程序参见附录6 （牛顿法程序）．具体实验7：其他方法求下列代数方程（组）的解：（1）命令：solve('x^5-x+1=0')（2）命令：[x, y]=solve('2*x+3*y=0', '4*x^2+3*y=1')(3) 求线性方程组的解，已知，命令：for i=1:5for j=1:5m(i, j)=i+j-1;endendm(5, 5)=0;b=[1:5]'linsolve(m, b)思考：若，或是类似的但阶数更大的稀疏方阵，则应如何得到四、自己动手1．对分法可以用来求偶重根附近的近似解吗为什么2．对照具体实验2、4、5，你可以得出什么结论3．选择适当的迭代过程，分别使用：（1）普通迭代法；（2）与之相应的松弛迭代法和 Altken 迭代法．求解方程在附近的根，精确到4位小数，请注意迭代次数的变化．4．分别用对分法、普通迭代法、松弛迭代法、Altken 迭代法、牛顿切法线等5种方法，求方程的正的近似根，．(建议取．时间许可的话，可进一步考虑的情况．)上一页回目录下一页。
MATLAB是一种用于数学计算、数据可视化和编程的高级技术计算语言和交互式环境。
它是许多工程和科学领域中的首选工具之一，能够帮助用户快速解决各种小问题。
本文将通过例子和代码，介绍MATLAB是如何解决小问题的。
1. 读取和绘制数据假设我们有一组实验数据，保存在一个名为"data.csv"的文件中。
我们可以使用MATLAB的csvread函数读取数据，然后使用plot函数绘制图形，如下所示：```matlabdata = csvread('data.csv'); % 读取数据plot(data(:,1), data(:,2)); % 绘制数据xlabel('x轴'); % 添加x轴标签ylabel('y轴'); % 添加y轴标签title('数据可视化'); % 添加标题```2. 拟合曲线现在我们想对这组数据进行曲线拟合，以便更好地理解数据的特征。
我们可以使用MATLAB的polyfit函数来进行多项式拟合，然后使用polyval函数绘制拟合曲线，如下所示：```matlabp = polyfit(data(:,1), data(:,2), 2); % 二次多项式拟合y_fit = polyval(p, data(:,1)); % 计算拟合曲线的值plot(data(:,1), data(:,2)); % 绘制原始数据hold on;plot(data(:,1), y_fit, 'r--'); % 绘制拟合曲线xlabel('x轴'); % 添加x轴标签ylabel('y轴'); % 添加y轴标签title('数据拟合'); % 添加标题legend('原始数据', '拟合曲线'); % 添加图例```3. 解方程假设我们需要解一个简单的方程，例如x^2-5x+6=0。
MATLAB 数学工具软件实例简明教程王正盛编写南京航空航天大学第一章 MATLAB简介MALAB 译于矩阵实验室MATrix LABoratory是用来提供通往LINPACK 和EISPACK 矩阵软件包接口的后来它渐渐发展成了通用科技计算图视交互系统和程序语言MATLAB 的基本数据单位是矩阵它的指令表达与数学工程中常用的习惯形式十分相似比如矩阵方程Ax=b 在MATLAB 中被写成A*x=b 而若要通过A,b求x 那么只要写x=A\b 即可完全不需要对矩阵的乘法和求逆进行编程因此用MATLAB 解算问题要比用C Fortran 等语言简捷得多MATLAB 发展到现在已经成为一个系列产品 MATLAB 主包和各种可选的toolbox 工具包主包中有数百个核心内部函数迄今所有的三十几个工具包又可分为两类功能性工具包和学科性工具包功能性工具包主要用来扩充MATLAB 的符号计算功能图视建模仿真功能文字处理功能以及硬件实时交互功能这种功能性工具包用于多种学科而学科性工具包是专业性比较强的如控制工具包Control Toolbox信号处理工具包( SignalProcessing Toolbox) 通信工具包(Communication Toolbox)等都属此类开放性也许是MATLAB 最重要最受人欢迎的特点除内部函数外所有MATLAB 主包文件和各工具包文件都是可读可改的源文件用户可通过对源文件的修改或加入自己编写文件去构成新的专用工具包MATLAB 已经受了用户的多年考验在欧美发达国家MATLAB 已经成为应用线性代数自动控制理论数理统计数字信号处理时间序列分析动态系统仿真等高级课程的基本教学工具成为攻读学位的大学生硕士生博士生必须掌握的基本技能在设计研究单位和工业部门MATLAB 被广泛地用于研究和解决各种具体工程问题第二章 MATLAB 入门2.1工作窗和指令行的操作除了通过菜单选项对工作窗进行控制外 MATLAB 还提供了许多通过键盘输入的控制指令如下表MATLAB工作窗中的部分通用指令在运行中显示变量和文字内容显示所有指定文件的全部内容启动MATLAB 后就可以利用它工作了由于MATLAB 是一种交互式语言随时输入指令即时给出运算结果是它的主要工作方式当然更可以编制程序详见第七章2sin( 0.3π)比如要计算的值只要在光标位置处键入1+ 52*sin(0.3*pi)/(1+sqrt(5))然后按[Enter]键,该指令便被执行并给出结果ans = 0.5000下面介绍控制光标对指令进行编辑的一些常用操作键常用操作键在MATLAB 中矩阵输入的方法有多种此处只简单介绍矩阵的直接输入法详细介绍见第三章在MATLAB 中不必对矩阵维数做任何说明存储将自动配置在直接输入矩阵时矩阵元素用空格或逗号分隔矩阵行用隔离整个矩阵放在方括号[ ] 里[例1]A=[1,2,3;4,5,6;7,8,9;10,11,12]A=1 2 34 5 67 8 910 11 12说明指令执行后矩阵A 被保存在MATLAB 的工作间Workspace中以备后用如果用户不用clear 指令清除它或对它重新定义该矩阵会一直保存在工作间中直到本MATLAB 指令窗被关闭为止[例2]矩阵分行输入A=[1 2 3 45 6 7 80 1 2 3]A = 12 3 45 6 7 80 1 2 3[例3]矩阵元素输入B(1,2)=3;B(4,4)=6;B(4,2)=11B= 0 3 0 00 0 0 00 0 0 00 11 0 62.3语句与变量MATLAB 采用表达式语句用户输入语句由MATLAB 系统结实运行MATLAB语句有两种常见的形式1 表达式 2 变量=表达式说明1 表达式由算符函数变量名和数字构成2在第一种形式中表达式被执行后产生的矩阵将被自动赋给名为ans 的变量并显示在屏幕上ans 是一个缺省变量名它会被以后类似的操作刷新3在第二种形式中等号右边的表达式是被演绎后产生的矩阵将被赋给等号左边的变量存入内存并显示在屏幕上4书写表达式时运算符号= + 以及* 等两侧允许有空格以增加可读性但在复数或符号表达式中要尽量避免装饰性空格以防出错5变量名函数名以一个字母打头后面最多可接19 个字母或数字注意MATLAB 是区分字母的大小写的[例1] 表达式的计算结果2001/81ans =24.7037[例2]运算结果的赋值s=1-1/2+1/3-1/4+1/5-1/6+1/7-1/8;说明结尾的分号作用是指令执行结果将不会显示在屏幕上但变量s 仍将驻留在内存中如想看s 的值只要键入s s = 0.63452.4Who Whos 和永久变量Who 和Whos 这两个指令的作用都是列出在MATLAB 工作间中已经驻留的变量名清单不过Whos 在给出变量名的同时还给出它们的维数及性质[例1]用 who 检查内存变量whoYour variables are:s[例2]用whos 检查驻留变量的详细情况whosName Size Bytes Class s1x1 8 double array Grandtotal is 1 elements using 8 bytes在MATLAB 工作内存中还驻留几个由系统本身在启动时定义的变量如下值表称为永久变量Permanent variables或称为预定义变量Predefinedvariables系统预定义的变量2它们不会被清除内存变量指令clear 所清除 3 他们可以重新定义为其他值但用clear 可清除重定义值恢复预定义[例1]无穷大s=1/0Warning: Divide by zero.s = Inf无穷大a=Inf/inf a= NaN2.5数与表达式MATLAB 的数值采用习惯的十进制表示可以带小数点或负号如下是合法的3-99 0.0013 9.2445154 1.2434e-6 4.673e33在采用IEEE 浮点算法的计算机上数值的相对精度是eps 即大约保持16 位有效数字数值范围大致为1×10−308 ~ 1×10308表达式由下列算符构成并按习惯的优先次序进行运算+ 加法减法 * 乘法 / 右除 \ 左除 ^ 乘方注意设置两种除法是为了方便矩阵的运算对标量而言两者作用相同[例1]x=2*pi/3+2^3/5-0.3e-3x =3.69412.6复数和复矩阵MATLAB 认识复数并定义i 和j 作为虚数单位矩阵元素允许是复数复变量和由它们组成的表达式[例1]z1=3+4*i,z2=2*exp(i*pi/6)z=z1*z2z1 =3.0000 +4.0000iz2 =1.7321 + 1.0000iz =1.1962 + 9.9282i[例2]A=[1,3;2,4]-i*[5,8;6,9]B=[1+5*i,2+6*i;3+8*i,4+9*i]C=A*BA=1.0000 - 5.0000i 3.0000 - 8.0000i2.0000 - 6.0000i 4.0000 - 9.0000iB=1.0000 + 5.0000i2.0000 + 6.0000i3.0000 + 8.0000i4.0000 + 9.0000iC=1.0e+002 *0.9900 1.1600 - 0.0900i1.1600 + 0.0900i 1.37002.7函数MATLAB 的强大功能可函数中略见一斑本质上讲分为三类1内部函数2系统附带各种工具包中的M 文件所提供的大量函数3用户自己增加的函数这一特点是其他许多软件平台无法比拟的MATLAB 提供的通用数理类函数包括基本数学函数 特殊函数 基本矩阵函数 特殊矩阵函数 矩阵分解和分析函数 数据分析函数 微分方程求解 多项式函数 非线性方程及其优化函数 数值积分函数 信号处理函数[例] z=1233.344x=sqrt(log(z))z=1.233344000000000e+003x =2.667861401680282.8显示格式在缺省的状态下MATLAB 以短格式short 格式显示计算结果可以用MATLAB 命令窗口中format 指令来改变数字的显示格式由于MATLAB 以双精度执行所有运算显示格式的设置仅影响矩阵的显示不影响矩阵的计算与存储如果矩阵的所有元素都是整数则矩阵以不带小数点的格式显示如果有一个元素不是整数则有几种输出格式默认格式为short 格式只显示5 位有效数字其他的显示格式可显示更多的有效数字还可用科学表示法[例] x=[4/31.2345e-6] 默认short格式x =1.3333 0.0000format short e 短格式科学表示xx =1.3333e+000 1.2345e-006format long 长格式xx =1.33333333333333 0.00000123450000format long e 长格式科学表示xx =1.333333333333333e+000 1.234500000000000e-006format bank 银行格式xx = 1.330.00format hex 十六进制格式xx = 3ff5555555555555 3eb4b6231abfd271format + +格式用于显示大矩阵的紧凑格式+ 空格分别表示正数负数和零x x =++另外还有一种命令为format compact(紧凑格式) 它消去了矩阵之间的间隔行这样可在一屏中显示更多的信息2.9变量的存储与调用quit 和exit 指令都可退出MATLAB 结束MATLAB 任务会删除工作间中的变量在退出前可以保存工作空间以备再次调出使用这些变量保存的指令格式1save 工作间中的所有变量保存在磁盘上名为matlab.mat 的文件中2save [文件名] [变量名] 将指定的变量保存在指定文件中如save temp x y z 把x,y,z 这三个变量保存在文件temp.mat 中在下次加载MATLAB 时可以利用load 指令将保存在文件中的变量恢复到工作间中其格式有load 将保存在matlab.mat 中的变量装入到MATLAB 工作间中load [文件名] [变量名] 从指定的文件中将指定的变量装入MATLAB 工作间如load temp x 从文件temp.mat 中只将变量x 装入到MATLAB 工作间中2.10图形图形是MATLAB 的主要特色之一MATLAB 图形指令具有自然简洁灵活及易扩充的特点MATLAB 的指令很多这里仅介绍几个简单的绘图指令详见第六章[例1]作多条曲线t=0:pi/50:4*pi;y0=exp(-t/3); y=exp(-t/3).*sin(3*t);plot(t,y,t,y0,t,-y0)grid[例2]三维曲面x=-8:0.5:8;y=x';X=ones(size(y))*x;Y=y*ones(size(x));R=sqrt(X.^2+Y.^2)+eps;Z=sin(R)./R; mesh(Z);colormap([1,0,0])2.11lp 指令lookfor 指令及其他帮助指令MATLAB 的在线帮助系统相当完备就查询系统的调用方式而言可分为两种1从MATLAB 指令窗的 help 菜单选项中寻求帮助此与一般windows 的求助方法一样2在MATLAB 指令窗中直接键入求助指令(i)help 不带任何参数显示出MATLAB 的目录项产生清单信息helpHELP topics:12matlab\general - General purpose commands. matlab\ops - Operators and special characters. matlab\lang - Programming language constructs. matlab\elmat - Elementary matrices and matrix manipulation. matlab\elfun - Elementary math functions. matlab\specfun - Specialized math functions. matlab\matfun - Matrix functions - numerical linear algebra. matlab\datafun - Data analysis and Fourier transforms.matlab\polyfun - Interpolation and polynomials. matlab\funfun - Function functions and ODEsolvers.matlab\sparfun - Sparse matrices.matlab\graph2d - Two dimensional graphs.matlab\graph3d - Three dimensional graphs.matlab\specgraph - Specialized graphs.matlab\graphics - Handle Graphics.matlab\uitools - Graphical user interface tools. matlab\strfun - Character strings. matlab\iofun - File input/output. matlab\timefun - Time and dates.matlab\datatypes - Data types and structures.matlab\winfun - Windows Operating System Interface Files (DDE/ActiveX) matlab\demos - Examples and demonstrations.toolbox\runtime - MATLAB Runtime Server Development Kit rtw\windows - Real Time Windows Target.daq\daq - Data Acquisition Toolbox daq\daqdemos - Data Acquisition Toolbox - Data Acquisition Demos. toolbox\dials - Dials & Gauges Blockset toolbox\rptgenext - Simulink Report Generator toolbox\rptgen - MATLAB Report Generator database\database - Database Toolbox. database\dbdemos - Database Toolbox Demonstration Functions.powersys\powerdemo - Power System Blockset Demos. powersys\powersys - Power System Blocksettoolbox\compiler - MATLAB Compiler (and Compiler1.2.1) comm\comm - Communications Toolbox. comm\commmasks - Communications Toolbox mask helper functions. comm\commsfun - Communications Toolbox S-functions. comm\commsim - Communications Toolbox Simulink files.toolbox\symbolic - Symbolic Math Toolbox.nag\nag - NAG Foundation Toolbox - Numerical & Statistical Librarynag\examples - NAG Foundation Toolbox - Numerical & Statistical Librarymap\map - Mapping Toolboxmap\mapdisp - Mapping Toolbox Map Definition and Display.map\mapproj - Mapping Toolbox Projections.wavelet\wavelet - Wavelet Toolbox.wavelet\wavedemo - Wavelet Toolbox Demos.toolbox\pde - Partial Differential EquationToolbox.finance\finance - Financial Toolbox.finance\calendar - Financial Toolbox calendarfunctions.finance\findemos - Financial Toolbox demonstrationfunctions. lmi\lmictrl - LMI Control Toolbox:Control Applicationslmi\lmilab - LMI Control Toolboxqft\qft - QFT Control DesignToolbox.qft\qftdemos - QFT Control Design ToolboxDemos toolbox\fixpoint - Fixed-Point Blocksetfixpoint\fxpdemos - Fixed-Point Blockset Demosfixpoint\obsolete - Obsolete Fixed-Point Blocksetdspblks\dspblks - DSP Blockset.dspblks\dspmex - (No table of contents file)dspblks\dspdemos - DSP Blockset demonstrations andexamples. dspblks\dspmasks - DSP Blockset maskhelper functions.fuzzy\fuzzy - Fuzzy Logic Toolbox.fuzzy\fuzdemos - Fuzzy Logic Toolbox Demos.mpc\mpccmds - Model Predictive Control Toolbox.mpc\mpcdemos - Model Predictive Control Toolboxfdident\fdident - Frequency Domain IdentificationToolbox. fdident\fddemos - Demonstrations for theFDIDENT Toolboxhosa\hosa - Higher-Order Spectral AnalysisToolbox. hosa\hosademo - Higher-Order SpectralAnalysis Toolbox - Demo suite toolbox\stats -Statistics Toolbox. toolbox\ncd - NonlinearControl Design Blockset images\images - ImageProcessing Toolbox. images\imdemos - ImageProcessing Toolbox --- demos and sample images nnet\nnet - Neural Network Toolbox. nnet\nndemos - NeuralNetwork Demonstrations. nnet\nnutils - (No tableof contents file) nnet\nnobsolete - (No table ofcontents file) mutools\commands - Mu-Analysis andSynthesis Toolbox. mutools\subs - Mu-Analysisand Synthesis Toolbox.signal\signal - Signal Processing Toolbox.signal\siggui - Signal Processing Toolbox GUIsignal\sigdemos - Signal Processing Toolbox Demonstrations toolbox\splines - Spline Toolbox.toolbox\optim - Optimization Toolbox.toolbox\robust - Robust Control Toolbox.toolbox\ident - System Identification Toolbox.toolbox\control - Control System Toolbox.control\ctrlguis - Control System Toolbox -- GUIsupport functions. control\obsolete - Control SystemToolbox -- obsolete commands. toolbox\rtw -Real-Time Workshop rtw\rtwdemos - (No table ofcontents file) stateflow\sfdemos - Stateflow demonstrations and samples.toolbox\sb2sl - SystemBuild to Simulink Translator stateflow\stateflow - Stateflow simulink\simulink - Simulink simulink\blocks - Simulink block library.simulink\simdemos - Simulink 3 demonstrations and samples.simulink\dee - Differential Equation EditorMATLAB53\work - (No table of contents file)toolbox\local - Preferences.For more help on directory/topic, type "help topic".(ii) help 目录名显示指定目录中的所有命令及其函数help lang 将列出与MATLAB 编程语言的所有命令及其函数help matfun 将列出与数值线性代数有关的所有矩阵函数help elfun 列出所有基本函数iiihelp 命令名/函数名/符号显示指定的命令名/函数名/符号的详细信息[例]help eig 显示计算矩阵特征值和特征向量的函数eig 的说明help eigEIG Eigenvalues and eigenvectors.E = EIG(X) is a vector containing the eigenvalues ofa square matrix X.[V,D] = EIG(X) produces a diagonal matrix D ofeigenvalues and afull matrix V whose columns are thecorresponding eigenvectors so that X*V = V*D.[V,D] = EIG(X,'nobalance') performs the computationwith balancingdisabled, which sometimes gives more accurateresults for certain problems with unusual scaling.E = EIG(A,B) is a vector containing thegeneralized eigenvalues of square matrices A and B.[V,D] = EIG(A,B) produces a diagonal matrix D of generalizedeigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that A*V = B*V*D.See also CONDEIG, EIGS.Overloaded methodshelp sym/eig.m helplti/eig.m注意help 的工作机理是是把指定名字的那个M 文件的第一段注释内容显示出来以构成自己文件的再线求助lookfor 指令可以根据用户提供的完整或不完整的关键词去搜索出一组与之有关的指令[例1]查找有关积分的指令lookfor integral[例2]ELLIPKE Complete elliptic integral.[例3]EXPINT Exponential integral function.[例4]DBLQUAD Numerically evaluate double integral.[例5]INNERLP Used with DBLQUAD to evaluate inner loop of integral.[例6]QUAD Numerically evaluate integral, low order method.[例7]QUAD8 Numerically evaluate integral, higher order method.[例8]COSINT Cosine integral function.[例9]SININT Sine integral function.[例10]ASSEMA Assembles area integral contributions in a PDE problem.[例11]COSINT Cosine integral function.[例12]FOURIER Fourier integral transform.[例13]IFOURIER Inverse Fourier integral transform.[例14]SININT Sine integral function.BLKPIDCON The output of the block is the sum ofproportional, integral and[例2]查找有关傅里叶变换的指令lookfor fourierFFT Discrete Fourier transform.FFT2 Two-dimensional discrete Fourier Transform.FFTN N-dimensional discrete Fourier Transform. IFFTInverse discrete Fourier transform. IFFT2 Two-dimensional inverse discrete Fourier transform.IFFTN N-dimensional inverse discrete Fourier transform. XFOURIER Graphics demo of Fourier series expansion. INSTDFFT Inverse non-standard 1-D fast Fourier transform. NSTDFFT Non-standard 1-D fast Fourier transform.EXPFOU Write data to a Fourier vector or a (maybe existing) file (for ELIS).IMPFOU Read complex amplitudes from a Fourier vector orfile (used by ELIS).MODIFYFV Modify Fourier (maybe also variance) data bygiven transfer function.SIMFOU Generate simulated Fourier amplitudes.PLOTFOU Plot contents of Fourier filesDFTMTX Discrete Fourier transform matrix.FOURIER Fourier integral transform.IFOURIER Inverse Fourier integral transform.注意 lookfor 指令的机制对MATLAB 中的每个M 文件注释区的第一行进行扫描一旦发现包含要查询的字符串就显示出来提示用户也可利用此机理建立自己文件的在线帮助其他帮助指令按扩展名分类列出在搜索路径中指定目录上的文件名列出指定名字文件所在的目录2.12为了保护MATLAB 目录结构的严整为了用户自己用MATLAB 所创建修改的M 文件和其他文件的方便用户应建立自己的工作目录MATLAB 启动后的默认目录是C:\MATLAB\BIN 若不建子目录则MATLAB环境产生的数据文件就登陆在这个缺省目录上建立工作目录两种方法1 在DOS 环境中建立 ;(2)在windows 环境下建立MATLAB 只能在启动时由mathabrc.m设定的路径上搜索不能与原定路径以外的其他目录交换信息可用以下三种方法扩充(1 在MATLAB 指令窗口中键入 CD C:\MYDIR2利用path 指令扩展搜索路径 path(path, 'c:\mydir')3在MATLAB 环境下键入pathtool 或者在MATLAB 指令窗口菜单上File 中的Set path 项设置第三章 MATLAB 的数值计算功能3.1数值矩阵的创建保存和数据格式3.1.1创建矩阵的直接输入法前面已述此出不赘述[例]x=14 ;y=4.32;A=[x,2*x-y,0;sin(pi/4),3*y+x,sqrt(y)]A =14.0000 23.6800 00.7071 26.9600 2.07853.1.2利用MATLAB 函数和语句创建数值矩阵[ 例1]利用指令reshape 创建数值矩阵av=1:12 %产生12个元素的行向量av以%开头的是注释行 bm=reshape(av,3,4) %利用向量av 创建3×4矩阵bmav =1 2 3 4 5 6 7 8 9 10 1112bm =1 4 7 102 5 8 113 6 9 12[例2]利用指令diag 产生对角阵ar=rand(4,4) %产生4×4的0-1均匀分布随即矩阵ar d=diag(ar) %用矩阵的主对角线元素形成向量d D=diag(d) %用向量d构成对角矩阵Dar =0.9501 0.8913 0.8214 0.92180.2311 0.7621 0.4447 0.73820.6068 0.4565 0.6154 0.17630.4860 0.0185 0.79190.4057 d = 0.95010.76210.61540.4057D =0.9501 0 0 00 0.7621 0 00 0 0.6154 00 0 0 0.40571.0.0利用M 文件创建和保存矩阵本节方法既适用于数值矩阵又适用于符号矩阵[例1]创建和保存矩阵AM 的matrix.m 文件生成过程步骤1 使用DOS 的编辑器edit ,Windows 的书写器(write)记事本notepad 或其他字处理软件如Word 等编辑如下AM=[1 2 3;3 4 5]步骤2 把此内容以纯文本方式ASCII 保存在用户自己的目录下名为matrix.m 的文件中步骤3 在MATLAB 指令窗中只要键入matrix 矩阵AM 就会自动生成于MATLAB 工作内存中即产生一个名为AM 的变量供显示和调用3.1.4通过MAT 文件保存和获取矩阵MAT 文件是MATLAB 保存数据的一种标准格式二进制文件MAT 文件的生成和调用由指令save 和load 进行[例1]把矩阵AR 保存到文件大他data.mat步骤1 在矩阵AR 存在于MATLAB 内存空间的前提下键入save data AR步骤2在下次进入MATLAB 后需要矩阵AR 时键入如下边可将data.mat 中的内容读入MATLAB 内存空间load data说明MATLAB 默认扩展名为.mat默认路径为matlab\bin 子目录用户如把data.mat 登陆在指定目录可用如下命令保存或调入save c:\mydir\data AR load c:\mydir\data AR3.1.5利用外部数据文件装入到指定矩阵假如磁盘中已有名为c:\mydir\data.dat 的ASCII 数据文件利用指令loadc:\mydir\data.dat 可在MATLAB 工作空间产生一个名为 data 的矩阵即变量当然也可以用指令fopen""fread"及其他MATLAB 底层数据输入输出I/O 指令实现可查看帮助如help fopen3.2矩阵的标识矩阵的元素子矩阵可以通过标识向量冒号的标识来援引和赋值[例1]b=[1 2 3 4 5; 6 7 8 9 10 ;11 12 13 14 15]b23=b(2,3) b1=b(1:2,[1 3 5]) b2=b([31],:)b([1 3],[2 4])=zeros(2)b =1 2 3 4 56 7 8 910 11 12 1314 15 b23 = 8b1 =1 35 6 810 b2 =11 12 13 1415 1 2 34 5 b =1 0 3 0 56 7 8 9 1011 0 13 0 15[例2]x=[1 2 3 4 5] %产生1×5向量x =1 2 3 4 5l=x<=3 %标出小于等于3 的元素的位置l=1 1 1 0 0x=x(l) %获得元素不超过3 的子向量x =1 2 33.3矩阵运算和数组运算矩阵运算和数组运算是Matlab 的数值运算中的两大类运算矩阵运算是按矩阵运算法则进行的运算数组运算无论是何种运算操作都是对元素逐个进行矩阵运算和数组运算指令对照汇总例a=[1 2 3; 4 5 6; 7 8 9];b=[1 2 3; 3 2 1;1 4 5]; c=[1 1 1;2 3 1;1 0 2];d=a*c^2+bd =32 31 3682 79 82128 129 1343.4矩阵函数和数组函数3.4.1基本数组函数数组函数是对各个元素的函数设计的f(.)基本函数表基本矩阵函数指令矩阵的条件数最大奇异范数全部奇平方根矩阵A矩阵的秩非零奇异值的a=magic(5);s=svd(a)'d=det(a),t=trace(a),rk=rank(a),c=cond(a)n1=norm(a,1),n2=norm(a),ninf=norm(a,inf),nf=norm(a,'fro') s =65.0000 22.5471 21.6874 13.403611.9008 d = 5070000 t = 65 rk = 5c =5.4618 n1= 65n2 =65.0000ninf =65 nf =74.33033.5 线性方程组的直接解法线性方程组Ax=b A 是n× m 的系数矩阵1)当n=m 且非奇异时此方程称为恰定方程Properly DeterminedEquation2)当n>m 时此方程称为超定方程Overdetermined Equation3)当n<m 时此方程称为欠定方程Underdetermined Equation3.5.1矩阵逆和除法解恰定方程组1采用求逆运算x=inv(A)*b2采用左除运算x=A\b说明1 由于MATLAB 遵循IEEE 算法所以即使A 阵奇异该运算也照样进行但在运算结束时一方面给出警告Warning:Matrix is singular to working precision 另一方面所得逆阵的元素都是Inf无穷大1当A 为病态时也给出警告信息2在MATLAB 中inv 指令不很有用MATLAB 推荐尽量使用除运算少用逆运算例1求逆法和左除法解恰定方程组的性能对比为对比两种方法的性能先用以下指令构造一个条件数很大的高阶恰定方程组rand('seed',12);%选定随机种子目的是可重复产生随机矩阵AA=rand(500)+1.e8;%rand(500)生成500×500均匀分布的随机矩阵 %每个随机矩阵元素加一个数的目的是使A的条件数变大x=ones(500,1); %令解向量x为全1的500元列向量b=A*x; %为使Ax=b方程一致用A和x生成b向量cond(A) %计算矩阵A的条件数ans =1.7608e+013过程 1 求逆法解恰定方程组的误差残差和所用计算时间tic %启动记时器Stopwatch Timerxi=inv(A)*b; %xi是用求逆法解恰定方程组所得的解toc %关闭计时器并显示解方程所用的时间eri=norm(x-xi) %解向量xi与真解向量的2-范误差rei=norm(A*xi-b)/norm(b) %方程的2-范相对残差elapsed_time =7.2500 eri = 0.0066 rei = 1.5488e-006过程2 左除法解恰定方程组的误差残差和所用计算时间ticxd=A\b;% xd是用左除法解恰定方程组所得的解toc erd=norm(x-xd)red=norm(A*xd-b)/norm(b)elapsed_time =3.3500 erd =0.0021 red =1.2695e-015说明1计算结果表明除法求解比求逆求解速度明显快精度相当但除法的相对残差几乎是机器零而逆阵法的相对残差高得多2MATLAB在设计求逆函数inv 时采用的是Gauss 消去法3MATLAB在设计左除运算解恰定方程时并不求逆而是直接采用Gauss 消去法求解有效地减少了残差所以即便在高条件数下也能得到较好的结果3.5.2矩阵除法解超定方程组1求正则方程Normal equations( A T A)x = A T b 的解2用Householder 变换Householder transformation 直接求原超定方程的最小二乘解由于第二种方程法采用的是正交变换据最小二乘理论可知第二种方法所得的解的准确性可靠性都比第一种方法好得多MATLAB 解超定方程组用的就是第二种方法例除运算解超定方程的简单算例a=[1 2 3;4 5 -6;7 8 9;10 11 12];b=[1:4]'; x=a\b x = -0.33330.66670.00003.5.3矩阵除法解欠定方程组欠定方程的解是不唯一的用除法运算所得的解有两个重要特征在解中至多有Rank(A)个非零元素 2 它是这类型解中范数最小的一个例除运算解欠定方程的简单算例a=[1 2 3;4 5 -6;7 8 9;10 11 12]; b=a'; c=[1 3 3]'; x=b\c x = 2.00000.1667 0-0.16673.6 矩阵分解函数 3.6.1 LU 分解[L,U]=lu(X)产生一个上三角矩阵U 和一个心理上下三角矩阵L 即由下三角矩阵和置换矩阵构成使得X=L*U[L,U,P]=lu(X) 产生一个上三角矩阵U 和一个下三角矩阵L 以及一个置换矩阵 P 使得P*X=L*U 注意X 必须是方阵例1 用lu 分解的两种指令格式对矩阵A 进行 LU 分解解 A=[1 –1 1; 5 –4 3; 2 1 1] [L,U]=lu(A) [L,U,P]=lu(A) A = 1 -1 1 5 -4 3 2 1 1 [ L,U]=lu(a ) L =0.2000 -0.0769 1.0000 1.0000 0 0 0.4000 1.0000 0 U =5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 [ L,U,P]=lu(a ) L =1.0000 0 0 0.4000 1.0000 0 0.2000 -0.0769 1.0000 U =5.0000 -4.0000 3.0000 0 2.6000 -0.2000 0 0 0.3846 P =0 1 00 0 1110 0例2 A为严格对角占优矩阵A=[4 1 2 ;2 5 –1; 5 3 11] [L,U]=lu(a) A =4 1 22 5 -15 3 11[ L,U]=lu(A )L=0.8000 -0.3684 1.00000.4000 1.0000 01.0000 0 0U =5.0000 3.0000 11.00000 3.8000 -5.40000 0 -8.7895例3 利用 LU 分解求解线性方程组AX=b A=[1,-1,1;5,-4,3;2,1,1]A = 1 -1 15 -43 2 11 b=[2;-3;1] b = 2-31[ L,U,P]=lu(A )L =1.0000 0 00.4000 1.0000 00.2000 -0.0769 1.0000U =5.0000 -4.0000 3.00000 2.6000 -0.20000 0 0.3846P =0 1 00 0 1 1 00 y=L\(P*b) y = -3.00002.20002.7692x=U\y x =-3.80001.40007.20001.0.0QR 分解[Q,R]=qr(A) 产生一个与矩阵A 相同维数的上三角矩阵R 和一个酉阵 Q 使得A=Q*R[Q,R,E]=qr(A)产生一个置换矩阵E 上三角矩阵R 和一个酉阵 Q使得A*E=Q*R例1a=[1 –1 1;5 –4 3;2 1 1][ Q,R]=qr(a )[Q,R,E]=qr(a)a =1-1 15 -4 32 1 1[ Q,R]=qr(a )Q= -0.1826 -0.1501 -0.9717-0.9129 -0.3412 0.2242-0.3651 0.9279 -0.0747R=-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737[ Q,R,E]=qr(a )Q=-0.1826 -0.1501 -0.9717-0.9129 -0.3412 0.2242-0.3651 0.9279 -0.0747R=-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737E =1 0 00 1 00 0 1例2 a为严格对角占优矩阵a=[10 –1 1;5 –14 3;2 1 21][ Q,R]=qr(a )Q=-0.1826 -0.1501 -0.9717-0.9129 -0.3412 0.2242-0.3651 0.9279 -0.0747R=-5.4772 3.4689 -3.28630 2.4427 -0.24560 0 -0.3737[ Q,R,E]=qr(a )Q=-0.0471 -0.0678 -0.9966-0.1413 -0.9872 0.0738-0.9889 0.1443 0.0369R=-21.2368 1.0359 -3.15490 14.0331 -5.32540 0 -9.5230E =0 0 10 1 010 03.6.3 Choleshy 分解Cholesky 分解要求被分解矩阵A 为对称正定矩阵R=chol(A) 产生一个上三角矩阵R 使得A = R T * R 例a=[2,1,1;1,2,-1;1,-1,3]R=chol(a)a =2 1 11 2 -11 -1 3R=chol(a)R =1.4142 0.7071 0.70710 1.2247 -1.22470 0 1.00003.7 多项式3.7.1 多项式的表达和创建多项式表达方式的约定多项式p( x ) = a0x n + a1x n−1 + + a n−1x + a n用以下系数行向量表示p = [ a0 ,a1, ,a n−1,a n ]多项式行向量的创建方法1多项式系数向量的直接输入法2利用指令p=poly(AR) 产生多项式系数向量若AR 是方阵则多项式为特征多项式若AR 是向量即AR = [ ar1,ar2 , ,ar n ] 则所得的多项式满足( x −ar1 )( x −ar2 ) ( x −ar n ) = a0x n + a1x n−1 + + a n−1x + a n例求3 阶方阵A 的特征多项式a=[11 12 13;14 15 16;17 18 19]; pa=poly(a)ppa=poly2str(pa,'s')pa =1.0000 -45.0000 -18.0000 0.0000 ppa =s^3 - 45 s^2 - 18 s + 5.3518e-015 注意1 n 阶方阵的特征多项式系数向量一定是(n+1)维的 2 特征多项式向量的第一的元素必是1 例由给定根向量求多项式系数向量r=[-0.5,-0.3+0.4*i,-0.3-0.4*i]; p=poly(r) pr=real(p) ppr=poly2str(pr,'x')p = 1.0000 1.1000 0.55000.1250 pr = 1.0000 1.1000 0.5500 0.1250 ppr = x^3 + 1.1 x^2 + 0.55 x +0.125 1 要形成实数多项式根向量中的复数根必须共轭成对2 含复数根的根向量所生成的多项式系数向量如p 的系数有可能带在截断误差数量级的虚部此时可采用取实部的指令'real'把虚部去掉3 poly2str 是一个函数文件它存在于MATLAB 的控制工具包( Control Toolbox)中3.6.2 常用多项式运算指令R=roots(p) 求多项式向量p 的根PA=polyval(p,S) 按数组运算规则计算多项式值p 为多项式S 为矩阵 PM=polyvalm(p;S) 按矩阵运算规则计算多项式值p 为多项式S 为矩 阵[r,p,k]=residue(b,a) 部分分式展开b,a 分别是分子分母多项式系数 向量r,p,k 分别是留数极点直项向 量P=polyfit(x,y,n) 用n 阶多项式拟合x,y 向量给定的数据 例1 求多项式x 3 −6x 2 −72x −27 的根解 r=roots([1,-6,-72,-27]) r = 12.1229 -5.7345说明-0.3884注意MATLAB 约定多项式系数用行向量表示一组根用列向量表示例2 两种多项式求值指令的差别s=pascal(4)p=poly(s);pp=poly2str(p,'x')pa=polyval(p,s)pm=polyvalm(p,s)s=1 1 1 11 2 3 41 3 610 1 410 20 pp =x^4 - 29 x^3 + 72 x^2 - 29 x + 1pa =1.0e+004 *0.0016 0.0016 0.0016 0.00160.0016 0.0015 -0.0140 -0.05630.0016 -0.0140 -0.2549 -1.20890.0016 -0.0563 -1.2089 -4.3779 pm = 1.0e-011 *-0.0077 0.0053 -0.0096 0.0430-0.0068 0.0481 -0.0110 0.12220.0075 0.1400 -0.0095 0.26080.0430 0.2920 -0.0007 0.4737说明: pm 中的元素都很小它是运算误差造成的理论上它应该是零这就是著名的Caylay-Hamilton 定理任何一个矩阵满足它自己的特征多项式2 x −t2 dt 进行最小二例3 用6 阶多项式对区间[0 2.5]上的误差函数y( x ) =∫0 eπ乘拟合解x=0:0.1:2.5;y=erf(x); %计算误差函数在[0,2.5]内的数据点p=polyfit(x,y,6) px=poly2str(p,'x')p =0.0084 -0.0983 0.4217 -0.7435 0.14711.1064 0.0004 px =0.0084194 x^6 - 0.0983 x^5 + 0.42174 x^4 - 0.74346 x^3+ 0.1471 x^ 2+ 1.1064 x + 0.00044117例4 有效拟合的区间性图示用[0,2.5]区间数据拟合曲线拟合[0,5]区间数据x=0:0.1:5; x1=0:0.1:2.5 y=erf(x); y1=erf(x1);p=polyfit(x1,y1,6) f=polyval(p,x); plot(x,y,'bo',x,f,'r-') axis([0,5,0,2]) legend('拟合曲线','原数据线')x1 =Columns 1 through 700.1000 0.2000 0.3000 0.4000 0.50000.6000Columns 8 through 140.7000 0.8000 0.9000 1.0000 1.10001.20001.3000Columns 15 through 211.4000 1.5000 1.6000 1.7000 1.80001.90002.0000Columns 22 through 262.1000 2.2000 2.3000 2.40002.5000 p =0.0084 -0.0983 0.4217 -0.7435 0.14711.1064 0.0004说明在[0,2.5]区间之外两条曲线的表现完全不同3.8 数值积分求函数数值积分指令S=quad('fname',a,b,tol,trace) 自适用递推Simpson 数值积分法S=quad8('fname',a,b,tol,trace) 自适用递推Newton -Cotes 数值积分法说明1输入参数'fname'是被积分函数表达式字符串或函数文件名2输入参数a,b 分别是积分上下限3输入参数tol 用来控制积分精度缺省时取tol=0.0014输入参数trace 若取1 则用图形展先积分过程取0 无图形缺省时不画图5quad8 比quad 有更高的积分精度但无论是quad8 还是quad,都不能处理可积的软奇异点例1设f ( x ) = e −0.5x sin( x + 6 ) 求S=∫03π f ( x )dx 解(1)建立被积函数文件fesin.mfunction f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6);(2)把函数文件fesin.m 存入自己的工作目录如d:\wzs 并在matlab 命令窗口中用下指令使该目录成为当前目录 cd d:\wzs(3)调用数值积分函数quad 求定积分 S=quad('fesin',0,3*pi) S=0.9008建立被积函数文件时注意所写程序应允许向量作为输入参数所以在该函数文件中采用了.*运算符号2.5 −xdx 的近似值并在同一积分精度例2分别用 quad 函数和 quad8 函数求∫1 e 下比较被积函数被调用的次数解为了统计被积函数被调用的次数可在被积函数文件中定义一个全局变量 num 每调用一次 num 加1 被积分函数为ftest.m function fx=ftest(x) global num; num=num+1; fx=exp(-x);调用函数quad 求定积分 global num; num=0;format long;I=quad('ftest',1,2.5,1e-6,1),num I= 0.28579444302661 global num;num=0; format long;πI=quad8('ftest',1,2.5,1e-6,1),numI = 0.28579444254754num = 4可见 quad8 32 次而且quad8的计算精度和效率也优于quad3.9 优化和解非线性方程非线性方程的求根方法很多常用的有牛顿迭代法弦截法等MATLAB 提供了有关指令函数用于非线性方程求根3.9.1多项式非线性函数求根 r=roots(p) 例p=[1 2 3 4]; r=roots(p)r =-1.6506-0.1747 + 1.5469i-0.1747 - 1.5469i3.9.2单变量非线性方程求解单变量函数求零点z=fzero('fname',x0,tol,trace)其中fname 是待求根的函数名x0 为搜索的起点一个函数可能有多个根但fzero 函数只给出离x0 最近的那个根tol 控制结果的相对精度缺省时取tol=eps;trace 指定迭代信息是否在运算中显示为1 时显示为0 时不显示缺省时trace=0例求函数y(t ) = 0.2t −e−0.5t sin(t +π)在t=2 附近的零点6步骤一建立函数文件f2.mfunction y=f2(t)y=0.2*t-exp(-0.5*t).*sin(t+pi/6); 步骤二求方程的根函数的零点z=fzero('f2',2)z= 1.69933.9.3一般非线性方程组求解此部分内容已超出MATLAB 的主包范围它需要用到MATLAB 的优化Optimization 工具包但由于非线性方程组是常用遇到的一类数学问题因此在此一并简单介绍对于一般非线性方程组F(X)=0 其数值解X 的求解指令X=fsolve('fname',X0)其中fname 是待求根的函数名x0 是对解的初始猜测值sin( x ) + y^2 + ln( z ) = 7求 3x + 2^ y −z^3+1 = 0 的数值解例x + y + z = 5 步骤一建立函数文件myxyz.mfunction q=myxyz(p) x=p(1);y=p(2);z=p(3);q=zeros(3,1) % 此指令可省略q(1)=sin(x)+y^2+log(z)-7; q(2)=3*x+2^y-z^3+1; q(3)=x+y+z-5;步骤二在给定的初值下调用fsolve 函数求方程的根x=fsolve('myxyz',[1 1 1]) 或xyz0=[1 1 1];x=fsolve('myxyz',xyz0)x= 0.59902.39592.0050说明从原理上讲非线性方程组也可用符号工具包中的solve 和vpa 指令求数值解详细可用help 指查询但计算时间难以接受3.10 微分方程的数值解[ t,x]=ode23('xprime',to,tf,x0,tol,trace )[ t,x]=ode45('xprime',to,tf,x0,tol,trace )或[ t,x]=ode23('xprime',[t0,tf],x0,tol,trace )[ t,x]=ode45('xprime',[t0,tf],x0,tol,trace )说明1两个指令的调用格式完全相同均为用Runge-Kutta 法2该指令是针对一阶常微分设计的因此假如待解的是高阶微分方程那么它必须先演化为形如x = f ( x,t ) 的一阶微分方程组即状态方程3'xprime'是定义f(x,t)的函数文件名该函数文件必须以x 为一个列向量输出以t,x 为输入参量注意变量的次序不可颠倒一定先时间变量后状态变量4输入参量t0 和tf 分别是积分的起始值和终止值5输入参量x0 为初始状态列向量6输出参量t 和x 分别给出时间向量和相应的状态向量7tol 控制解的精度可缺省缺省时ode23 默认tol=1.e-3 ode45 默认tol=1.e-68输入参量trace 控制求解的中间结果是否显示可缺省缺省时默认tol=0 不显示中间结果9一般地两者分别采用自适应变步长即当解的变化较慢时采用较大的步长从而使得计算速度快当解的变化速度较快时步长会自动地变小从而使得计算精度高的二三阶Runge-Kutta 算法和四五阶Runge_Kutta 算法ode45 比ode23 的积分分段少而运算速度快y2 −t −2例1 求初值问题 y' = 4(t +1) ,0 ≤ t ≤10 的数值解并与解析解y( 0 ) = 2y(t ) = t +1+1相比较解(1)建立函数文件funt.mfunction yp=funt(t,y) yp=(y^2-t-2)/4/(t+1);(2)求解微分方程 t0=0;tf=10;y0=2;[t,y]=ode23('funt',[t0,tf],y0);y1=sqrt(t+1)+1; t',y',y1'plot(t,y,'-r',t,y1,':b') legend('数值解 ','解析解')。
