《数值计算方法》

1、数值计算方法上机实验指导书邹昌文编2009年10月“数值计算方法”上机实验指导书实验一 误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式显然该多项式的全部根为1,2,20共计20个，且每个根都是单重的。现

2、考虑该多项式的一个扰动其中是一个非常小的数。这相当于是对（1.1）中的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个MATLAB函数：“roots”和“poly”。其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。设a的元素依次为，则输出u的各分量是多项式方程的全部根；而函数 的输出b是一个n+1维向量，它是以n维向量v的各分量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。上述简单的MATLAB程序便得到（1.2）的全部根，程序中的“ess”即是（1.2

3、）中的。实验要求：（1） 选择充分小的ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？（2） 将方程（1.2）中的扰动项改成或其它形式，实验中又有怎样的现象出现？（3） （选作部分）请从理论上分析产生这一问题的根源。注意我们可以将方程（1.2）写成展开的形式， 同时将方程的解x看成是系数的函数，考察方程的某个解关于的扰动是否敏感，与研究它关于的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于的变化更敏感？思考题一：（上述实验的改进）在上述实验中

4、我们会发现用roots函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考MATLAB的帮助。思考题二：（二进制产生的误差）用MATLAB计算。结果居然有误差！因为从十进制数角度分析，这一计算应该是准确的。实验反映了计算机内部的二进制本质。思考题三：（一个简单公式中产生巨大舍入误差的例子）可以用下列式子计算自然对数的底数 这个极限表明随着n的增加，计算e值的精度是不确定的。现编程计算与exp(1)值的差。n大到什么程度的时候误差最大？你能解释其中的原因吗？相关MATLAB函数提示：p

5、oly(a) 求给定的根向量a生成其对应的多项式系数（降序）向量roots(p) 求解以向量p为系数的多项式（降序）的所有根poly2sym(p) 将多项式向量p表示成为符号多项式（降序）sym(arg) 将数字、字符串或表达式arg转换为符号对象syms arg1 arg2 argk 将字符arg1,arg2,argk定义为基本符号对象solve(eq1) 求符号多项式方程eq1的符号解实验二 插值法实验2.1（多项式插值的振荡现象）问题提出：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，是否也更加靠近被

6、逼近的函数。龙格（Runge）给出一个例子是极著名并富有启发性的。设区间-1,1上函数 实验内容：考虑区间-1,1的一个等距划分，分点为 则拉格朗日插值多项式为 其中的是n次拉格朗日插值基函数。实验要求：（1） 选择不断增大的分点数目n=2,3.，画出原函数f(x)及插值多项式函数在-1，1上的图像，比较并分析实验结果。（2）选择其他的函数，例如定义在区间-5，5上的函数重复上述的实验看其结果如何。（3）区间a,b上切比雪夫点的定义为 以为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。实验2.2（样条插值的收敛性）问题提出：多项式插值是不收敛的，即插值的节点多，效果不一定就好。对样条

7、函数插值又如何呢？理论上证明样条插值的收敛性是比较困难的，但通过本实验可以验证这一理论结果。实验内容：请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验2.1中的函数或选择其他你有兴趣的函数，可以用MATLAB的函数“spline”作此函数的三次样条插值。实验要求：（1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较。（2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下：xk012345678910yk0.00.791.532.19

8、2.713.033.272.893.063.193.29yk0.80.2 思考题一：（二维插值）在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值，并由此找出最高点和该点的高程。yx100200300400100636697624478200698712630478300680674598412400662626552334相关MATLAB函数提示：plot(x,y) 作出以数据(x(i),y(i)为节点的折线图，其中x,y为同长度的向量subplot(m,n,k) 将图形窗口分为m*n个子图，并指向第k幅图yi=in

9、terp1(x,y,xi) 根据数据(x,y)给出在xi的分段线性插值结果yipp=spline(x,y) 返回样条插值的分段多项式(pp)形式结构pp=csape(x,y,边界类型,边界值) 生成各种边界条件的三次样条插值yi=ppval(pp,xi) pp样条在xi的函数值ZI=interp2(x,y,z,xi,yi) x,xi为行向量，y,yi为列向量，z为矩阵的双线性二维插值ZI=interp2(,spline) 使用二元三次样条插值ZI=griddata(x,y,z,xi,yi) x,y,z均为向量（不必单调）表示数据，xi,yi为网格向量的三角形线性插值（不规则数据的二维插值）实验

10、三 函数逼近与曲线拟合实验3.1（曲线逼近方法的比较）问题提出：曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。实验内容：考虑实验2.1中的著名问题。下面的MATLAB程序给出了该函数的二次和三次拟合多项式。x=-1:0.2:1;y=1/(1+25*x.*x);xx=-1:0.02:1;p2=polyfit(x,y,2);yy=polyval(p2,xx);plot(x,y,o,xx,yy);xlabel(x);ylabel(y);hold on;p3=polyfit(x,y,3);yy=polyval(p3,xx);plot(x

11、,y,o,xx,yy);hold off;适当修改上述MATLAB程序，也可以拟合其他你有兴趣的函数。实验要求：（1）将拟合的结果与拉格朗日插值及样条插值的结果比较。（2）归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。思考题一：（病态）考虑将0,130等分节点，用多项式生成数据，再用polyfit求其3次、5次、10次、15次拟合多项式，并分析误差产生的原因。相关MATLAB函数提示：p=polyfit(x,y,k) 用k次多项式拟合向量数据(x,y)，返回多项式的降幂系数，当k=n-1时，实现多项式插值，这里n是向量维数实验四 数值积分与数值

12、微分实验4.1 （高斯数值积分方法用于积分方程求解）问题提出：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。实验内容：求解第二类Fredholm积分方程首先将积分区间a,b等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。实验要求：分别使用如下方法，离散积分方程中的积分1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。求解如下的积分方程。1），方程的准确解为；2），方程的准确解为；比较各算法的计算量和误

13、差以分析它们的优劣。实验4.2（高维积分数值计算的蒙特卡罗方法）问题提出：高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。实验内容：对于一般的区域，计算其测度（只要理解为平面上的面积或空间中的体积）的一般方法是：先找一个规则的区域A包含，且A的测度是已知的。生成区域A中m个均匀分布的随机点，如果其中有n个落在区域中，则区域的测度m()为n/m。函数f(x)在区域上的积分可以近似为：区域的测度与函数f(x)在中n个随机点上平均值的乘积，即实验要求

14、：假设冰琪淋的下部为一锥体而上面为一半球，考虑冰琪淋体积问题：计算锥面上方和球面内部区域的体积。如果使用球面坐标，该区域可以表示为如下的积分：用蒙特卡罗方法可以计算该积分。另一方面，显然这样的冰琪淋可以装在如下立方体的盒子里而该立方体的体积为8。只要生成这个盒子里均匀分布的随机点，落入冰琪淋锥点的个数与总点数之比再乘以8就是冰琪淋锥的体积。比较两种方法所得到的结果。类似的办法可以计算复杂区域的测度（面积或体积）。试求由下列关系所界定区域的测度：相关MATLAB函数提示：diff(x) 如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分diff(x,k) k阶差分q=polyder(p

15、) 求得由向量p表示的多项式导函数的向量表示qFx=gradient(F,x) 返回向量F表示的一元函数沿x方向的导函数F(x)，其中x是与F同维数的向量z=trapz(x,y) x表示积分区间的离散化向量；y是与x同维数的向量，表示被积函数；z返回积分的近似值z=guad(fun,a,b,tol) 自适应步长Simpson积分法求得Fun在区间a,b上的定积分，Fun为M文件函数句柄，tol为积分精度z=dblquad(fun,a,b,c,d,tol,method) 求得二元函数Fun(x,y)的重积分z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元

16、函数Fun(x,y,z)的重积分实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵，则方程有解。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择

17、程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）将上述矩阵A中的主元改为0.00006再重新作一次数值实验看看。（5）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组的摄动满足 矩阵的条件

18、数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算通常要比求解方程还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动，使得充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程，以1-范数，给出的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结

19、果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和的估计，马上就可以给出的估计。（4）估计著名的Hilbert矩阵的条件数。思考题一：（Vadermonde矩阵）设 ，其中，（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛

20、顿插值法的原因吗？相关MATLAB函数提示：zeros(m,n) 生成m行，n列的零矩阵ones(m,n) 生成m行，n列的元素全为1的矩阵eye(n) 生成n阶单位矩阵rand(m,n) 生成m行,n列(0,1)上均匀分布的随机矩阵diag(x) 返回由向量x的元素构成的对角矩阵tril(A) 提取矩阵A的下三角部分生成下三角矩阵triu(A) 提取矩阵A的上三角部分生成上三角矩阵rank(A) 返回矩阵A的秩det(A) 返回方阵A的行列式inv(A) 返回可逆方阵A的逆矩阵V，D=eig(A) 返回方阵A的特征值和特征向量norm(A,p) 矩阵或向量的p范数cond(A,p) 矩阵的条

21、件数L，U，P=lu(A) 选列主元LU分解R=chol(X) 平方根分解Hi=hilb(n) 生成n阶Hilbert矩阵实验六 解线性方程组的迭代法实验6.1（病态的线性方程组的求解）问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵， 这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比

22、较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？思考题一：讨论病态问题求解的算法Jacobi迭代法与Gauss-seidel迭代法的比较：用Jacobi迭代法与Gauss-seidel迭代法解下列方程组：（1）（2）取，你能得出什么结论？思考题二：SOR（超松驰）迭代法松驰因子对收敛性及速度的影响试用SOR（超松驰）迭代法求解下列方程组：取，选择松驰因子0.8,0.9,1,1.1,1.2等，试看对算法收敛性的影响，并找出所选用的松驰因子的最佳者。要求编制矩阵迭代求解的函数文件。实验6.2方程组性态讨论（1） 求的解向量，其中： （2） 求

23、系数矩阵A的条件数；（3） 将改为，改为，求解向量; （4） 令，求解，并求系数矩阵PA的条件数；（5） 对PA中的和Pb中的给以的扰动，求解向量结合上述结果，讨论两个方程组的性态实验6.2大型稀疏方程组的数值解法设n阶方阵：B为A的各元素之和，显然Ax=b的解为，用下面三种方法对于阶数n=100,200,500，误差限为的各种组合求解，分析收敛速度（1） 选取列主元Gauss消元法；（2） Jacobi迭代法（3） Gauss-seidel迭代法算法的改进：由于是稀疏矩阵，请考虑改进矩阵的存储方式，减少存储空间的使用来提高计算效率。也就是用稀疏存储方式。实验七 非线性方程求根实验7.1（迭代

24、法、初始值与收敛性）实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程 针对上述方程，可以构造多种迭代法，如 在实轴上取初始值x0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。实验要求：（1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用MATLAB的图形功

25、能），分析三种迭代法的收敛性与初值选取的关系。（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。思考题一：用Newton法求方程在区间-3,3上误差不大于的根，分别取初值进行计算，比较它们的迭代次数。相关MATLAB函数提示：x=fzero(fun,x0) 返回一元函数fun的一个零点，其中fun为函数句柄，x0为标量时，返回在x0附近的零点；x0为向量a,b时，返回函数在a,b中的零点x,f,h=fsolve(fun,x0) 返回一元或多元函数x

26、0附近fun的一个零点，其中fun为函数句柄，x0为迭代初值；f返回fun在x的函数值，应该接近0； h返回值如果大于0，说明计算结果可靠，否则不可靠实验八 常微分方程初值问题数值解法实验8.1（Lorenz问题与混沌）问题提出：考虑著名的Lorenz方程 其中s,r,b为变化区域有一定限制的实参数。该方程形式简单，表面上看并无惊人之处，但由该方程提示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。实验内容：先取定初值y0=(0,0,0),参数s=10,r=28,b=8/3,用MATLAB的数值求常微分方程函数ods45编程对（8.1）进行求解实验要求：（1）

27、对目前取定的参数值s,r和b，选取不同的初值y0进行运算，观察计算的结果有什么特点？解的曲线是否有界？解的曲线是不是周期的或趋于某个固定点？（2）在问题允许的范围内适当改变其中的参数值s,r,b，再选取不同的初始值y0进行运算，观察并记录计算的结果有什么特点？是否发现什么不同的现象？思考题一：考虑如下常微分方程 先固定其中参数：r=0.9,b=3,a=1,k0=0,m=1,初始值为原点，t的范围为0-500分别取不同的k1,如1，1.5，2，2.5，3，3.5等，绘图观察相位图（x,x）的变化情况，能找到其中的规律吗？试试改变其它参数呢？相关MATLAB函数提示：t,y=ode45(odefu

28、n,tspan,y0) odefun表示f(t,y)的函数句柄，t是标量，y是标量或向量；tspan是二维向量t0,tf，表示自变量初值t0和终值tf；y0表示初值向量，若无输出参数，则作出图形附录MATLAB简介MATLAB作为线性系统的一种分析和仿真工具，是理工科大学生应该掌握的技术工具，它作为一种编程语言和可视化工具，可解决工程、科学计算和数学学科中许多问题。MATLAB建立在向量、数组和矩阵的基础上，使用方便，人机界面直观，输出结果可视化。其中，矩阵是MATLAB的核心！下面将从基本规则和操作，编程和作图以及文件的操作等等多个方面，来讲解数学建模中MATLAB的一些常用方法。一、 变量

29、与函数变量、函数与编程所形成的m文件是MATLAB操作的基本，在介绍它们的具体使用方法之前，先给出一些必须了解的基本规则。1、 变量MATLAB和其他编程工具一样，变量是必须的基本元素，它也是以字母开头，后接字母、数字或下划线的字符序列，用法也基本一样。其具体的命名规则是：（1）变量名必须是不含空格的单个词；（2）变量名区分大小写；（3）变量名最多不超过19个字符；（4）变量名必须以字母打头，之后可以是任意字母、数字或下划线，变量名中不允许使用标点符号。特殊变量表注意，自定义的变量不能和表中变量名相同。2、数学运算符号及标点符号（1）MATLAB的每条命令后，若为逗号或无标点符号，则显示命令的

30、结果；若命令后为分号，则禁止显示结果；（2）“%”后面所有文字为注释；（3）“.”表示续行。3、数学函数基本上，常用的数学函数在MATLAB中都有相应的命令，部分如下：4、M文件MATLAB的内部函数是有限的，有时为了研究某一个函数的各种性态，需要为MATLAB定义新函数，为此必须编写函数文件. 函数文件是文件名后缀为M的文件，这类文件的第一行必须是一特殊字符function开始，格式为： function 因变量名=函数名（自变量名）函数值的获得必须通过具体的运算实现，并赋给因变量。M文件建立方法：1. 在Matlab中，点:File-New-M-file 2. 在编辑窗口中输入程序内容 3

31、. 点：File-Save，存盘，M文件名必须与函数名一致。Matlab的应用程序也以M文件保存。例：定义函数 f(x1,x2)=100(x2-x1)+(1-x1)1.建立M文件：fun.mfunction f=fun(x)f=100*(x(2)-x(1)2)2+(1-x(1)22. 可以直接使用函数fun.m例如：计算 f(1,2), 只需在Matlab命令窗口键入命令：x=1 2fun(x)二、数组与矩阵矩阵是MATLAB最基本的数据对象，MATLAB的大部分运算或命令都是在矩阵运算的意义下执行的。在MATLAB中，不需对矩阵的维数和类型进行说明，MATLAB会根据用户所输入的内容自动进行

32、配置。1、数 组 数组可以看作是只有一行或一列的简单矩阵，但作为常用的计算单元，matlab也专门为其设计了一系列命令。（1）创建简单的数组x=a b c d e f 创建包含指定元素的行向量x=first：last 创建从first开始，加1计数，到last结束的行向量x=first：increment：last 创建从first开始，加increment计数，last结束的行向量x=linspace(first，last，n） 创建从first开始，到last结束，有n个元素的行向量x=logspace(first，last，n） 创建从开始，到结束，有n个元素的对数分隔行向量。（2）数组

33、元素的访问（i）访问一个元素：x(i)表示访问数组x的第i个元素；（ii）访问一块元素：x(a ：b ：c)表示访问数组x的从第a个元素开始，以步长为b到第c个元素（但不超过c），b可以为负数，b缺损时为1；（iii）直接使用元素编址序号：x(a b c d) 表示提取数组x的第a、b、c、d个元素构成一个新的数组 x(a) x(b) x(c) x(d)。（3）数组的方向 前面例子中的数组都是一行数列，是行方向分布的，称之为行向量。数组也可以是列向量，它的数组操作和运算与行向量是一样的，唯一的区别是结果以列形式显示。产生列向量有两种方法： 直接产生 例 c=1；2；3；4 转置产生 例 b=1

34、 2 3 4; c=b 说明：以空格或逗号分隔的元素指定的是不同列的元素，而以分号分隔的元素指定了不同行的元素。（4）数组的运算（i）标量-数组运算 数组对标量的加、减、乘、除、乘方是数组的每个元素对该标量施加相应的加、减、乘、除、乘方运算。设：a=a1,a2,an, c=标量则：a+c=a1+c,a2+c,an+c a.*c=a1*c,a2*c,an*c a./c= a1/c,a2/c,an/c(右除） a.c= c/a1,c/a2,c/an (左除） a.c= a1c,a2c,anc c.a= ca1,ca2,can （ii）数组-数组运算 当两个数组有相同维数时，加、减、乘、除、幂运算可

35、按元素对元素方式进行的，不同大小或维数的数组是不能进行运算的。设：a=a1,a2,an, b=b1,b2,bn则：a+b= a1+b1,a2+b2,an+bn a.*b= a1*b1,a2*b2,an*bn a./b= a1/b1,a2/b2,an/bn a.b=b1/a1,b2/a2,bn/an a.b=a1b1,a2b2,anbn2、矩阵 通过下面的介绍，大家可以发现矩阵的相关命令和数组基本类似。（1）矩阵的建立逗号或空格用于分隔某一行的元素，分号用于区分不同的行。除了分号，在输入矩阵时，按Enter键也表示开始一新行。输入矩阵时，严格要求所有行有相同的列。例 m=1 2 3 4 ；5 6

36、 7 8；9 10 11 12 p=1 1 1 1 2 2 2 2 3 3 3 3特殊矩阵的建立：a= 产生一个空矩阵，当对一项操作无结果时，返回空矩阵，空矩阵的大小为零；b=zeros(m，n) 产生一个m行、n列的零矩阵c=ones(m，n) 产生一个m行、n列的元素全为1的矩阵d=eye(m，n) 产生一个m行、n列的单位矩阵（2）矩阵中元素的操作（a）矩阵A的第r行：A（r，：）（b）矩阵A的第r列：A（：，r）（c）依次提取矩阵A的每一列，将A拉伸为一个列向量：A（：）（d）取矩阵A的第i1i2行、第j1j2列构成新矩阵:A(i1:i2, j1:j2)（e）以逆序提取矩阵A的第i1i

37、2行，构成新矩阵:A(i2:-1：i1，：）（f）以逆序提取矩阵A的第j1j2列，构成新矩阵:A(:, j2:-1：j1 ）（g）删除A的第i1i2行，构成新矩阵:A(i1:i2，：)= （h）删除A的第j1j2列，构成新矩阵:A(：，j1:j2)= （i）将矩阵A和B拼接成新矩阵：A B；A；B（3）矩阵的运算（i）标量-矩阵运算：同标量-数组运算；（ii）矩阵-矩阵运算： 1 元素对元素的运算，同数组-数组运算；2 矩阵运算：矩阵加法：A+B矩阵乘法：A*B方阵的行列式：det（A）方阵的逆：inv（A）方阵的特征值与特征向量：V，D=eigA。三、 MATLAB编程Matlab虽然提供了

38、大量现成的函数供我们计算时使用，但当要真正解决一些实际问题时，这还远远不够，编程仍然是不可或缺的一步。Matlab也可以象C，FORTRAN等高级计算机语言一样，进行程序设计。下面简单介绍Matlab中一些重要的编程手段。1、关系与逻辑运算除了传统的数学运算，Matlab支持关系和逻辑运算。一个重要的应用是控制基于真/假命题的一系列Matlab命令的流程，或执行次序。（1）关系操作符 Matlab关系操作符包括所有常用的比较。能用来比较两个同样大小的矩阵，或用来比较一个矩阵和一个标量（矩阵中每一个元素都与标量相比，结果与矩阵大小一样）。（2）逻辑运算符 逻辑操作符提供了一中组合或否定关系表达式

39、，具体如下：（3）控制流MATLAB提供三种决策或控制流结构：for循环、while循环、if-else-end结构。 这些结构经常包含大量的MATLAB命令，故经常出现在MATLAB程序中，而不是直接加在MATLAB提示符下。i、for循环：允许一组命令以固定的和预定的次数重复 for x=array commands end 在for和end语句之间的命令串commands按数组（array）中的每一列执行一次。在每一次迭代中，x被指定为数组的下一列，即在第n次循环中，x=array(：，n)。ii、While循环 与for循环以固定次数求一组命令相反，while循环以不定的次数求一组语句

40、的值。while expression commands end 只要在表达式(expression)里的所有元素为真，就执行while和end语句之间的命令串commands。iii、If-Else-End结构（a）有一个选择的一般形式是： if expression commands end 如果在表达式(expression)里的所有元素为真，就执行if和end语句之间的命令串commands。先建立M文件fun1.m定义函数f（x），再在Matlab命令窗口输入fun1(2),fun1(-1)即可。（b） 有三个或更多的选择的一般形式是： if （expression1） comman

41、ds1 else if （expression2） commands2 else if （expression3） commands3 else if elsecommands endendend end先建立M文件fun2.m定义函数f（x），再在Matlab命令窗口输入fun2(2),fun2(0.5), fun2(-1)即可。四、MATLAB作图强大的图形功能是Matlab的优点之一，它能方便快速的出图，给我们的工作带来了巨大的便利。1. 曲线图Matlab作图是通过描点、连线来实现的，故在画一个曲线图形之前，必须先取得该图形上的一系列的点的坐标（即横坐标和纵坐标），然后将该点集的坐标传

42、给Matlab函数画图。命令为：PLOT(X,Y,S)其中，X,Y是向量，分别表示点集的横坐标和纵坐标，S指定曲线的颜色、线形等：y：黄色 c：蓝绿色 r：红色 m：洋红.：点 -：连线 o：圈 :：短虚线x：x-符号 -.：长短线 +：加号 -：长虚线PLOT(X,Y)画实线PLOT(X,Y1,S1,X,Y2,S2,X,Yn,Sn)将多条线画在一起例 在0,2*pi用红线画sin(x),用绿圈画cos(x)。解：x=linspace(0,2*pi,30);y=sin(x);z=cos(x);plot(x,y,r,x,z,g0)2.符号函数(显函数、隐函数和参数方程)画图(1) ezplote

43、zplot(f(x),a,b) 表示在axb绘制显函数f=f(x)的函数图ezplot(f(x,y),xmin,xmax,ymin,ymax) 表示在区间xminxxmax和 yminyymax绘制隐函数f(x,y)=0的函数图ezplot(x(t),y(t),tmin,tmax) 表示在区间tminttmax绘制参数方程x=x(t),y=y(t)的函数图例 在0,pi上画y=cos(x)的图形；解：输入命令 ezplot(sin(x),0,pi)例 在0,2*pi上画，星形图；解：输入命令 ezplot(cos(t)3,sin(t)3,0.2*pi)例 在-2，0.5，0，2上画隐函数的图。

44、解：输入命令ezplot(exp(x)+sin(x*y),-2,0.5,0,2)(2) fplotfplot(fun,lims) 表示绘制字符串fun指定的函数在lims=xmin,xmax的图形注意：1 fun必须是M文件的函数名或是独立变量为x的字符串；2 fplot函数不能画参数方程和隐函数图形，但在一个图上可以画多个图形。例 在-1，2上画的 图形解：先建M文件myfun1.m function Y=myfun1(x) Y=exp(2*x)+sin(3*x.2)再输入命令：fplot(myfun1,-1,2) 例 在-2,2范围内绘制函数tanh的图形解：fplot(tanh,-2,2

45、)例 x、y的取值范围都在-，画函数tanh(x),sin(x),cos(x)的图形解：输入命令 fplot(tanh(x),sin(x),cos(x),2*pi*-1 1 1 1)3、空间曲线(1)一条曲线：PLOT3(x,y,z,s)其中，X，Y，Z是n维向量，分别表示曲线上点集的横坐标、纵坐标、函数值S则指定颜色、线形等例 在区间0，10*pi画出参数曲线x=sin(t)，y=cos(t)，z=t解：t=0:pi/50:10*pi; plot3(sin(t),cos(t),t) rotate3d %旋转(2)多条曲线：PLOT3(x,y,z)其中x，y，z是都是m*n矩阵，其对应的每一列

46、表示一条曲线.例 画多条曲线观察函数Z=(X+Y).2. 解：x=-3:0.1:3;y=1:0.1:5; X,Y=meshgrid(x,y); Z=(X+Y).2; plot3(X,Y,Z)（这里meshgrid(x,y)的作用是产生一个以向量x为行、向量y为列的矩阵）4、空间曲面(1)surf(x,y,z) 画出数据点（x，y，z）表示的曲面数据矩阵。分别表示数据点的横坐标、纵坐标、函数值例 画函数Z=(X+Y).2的图形. 解：x=-3:0.1:3; y=1:0.1:5; X,Y=meshgrid(x,y); Z=(X+Y).2; surf(X,Y,Z) shading flat %将当前

47、图形变得平滑(2)Mesh(x,y,z) 画网格曲面例 画出曲面Z=(X+Y).2在不同视角的网格图. 解：x=-3:0.1:3; y=1:0.1:5; X,Y=meshgrid(x,y); Z=(X+Y).2; mesh(X,Y,Z) (3)meshz(X,Y,Z) 在网格周围画一个curtain图(如,参考平面)例 绘peaks的网格图解：输入命令: X,Y=meshgrid(-3:.125:3); Z=praks(X,Y); Meshz(X,Y,Z)5、特殊二、三维图形二维：1、极坐标图：polar (theta,rho,s)用角度theta（弧度表示）和极半径rho作极坐标图，用s指定

48、线型。2、散点图：scatter（X,Y,S,C） 在向量X和Y的指定位置显示彩色圈X和Y必须大小相同3、平面等值线图：contour (x,y,z,n) 绘制n个等值线的二维等值线图三维：1、空间等值线图：contour 3(x,y,z,n) 其中n表示等值线数。2、三维散点图：scatter3（X,Y,Z,S,C） 在向量X,Y和Z指定的位置上显示彩色圆圈； 向量X,Y和Z的大小必须相同。五、文件操作Matlab提供了对数据文件建立、打开、读、写以及关闭等一系列函数，数据文件一般存放在磁盘等介质上，用文件名标识，系统对文件名没有特殊要求。文件数据格式有二种形式，一是二进制格式文件，二是AS

49、CII文本文件，系统对这两类文件提供了不同的读写功能函数。1、文件的打开与关闭（1）fopen打开文件在读写文件之前，必须先用fopen命令打开一个文件，并指定允许对该文件进行的操作。文件操作结束后，应及时关闭文件，以免数据的丢失或误修改。fopen函数格式为：Fid= fopen(filename，permission)其中filename为文件名，permission为文件格式，可以是下列格式之一： r 打开文件，读数据，文件必须存在。 w 打开文件，写数据，若文件不存在，系统会自动建立。 a 打开文件，在文件末尾添加数据。 r+ 打开文件，可以读和写数据，文件必须存在。 w+ 打开文件，

50、供读与写数据用。 a+ 打开文件，供读与添加数据用。 W 打开文件供写数据用，无自动刷新功能。 A 打开文件供添加数据用，无自动刷新功能。例如：打开一个名为std.dat的数据文件并进行读操作，其命令格式为：Fid=fopen( std.dat, r )上述打开格式均为二进制格式，如果想用ASCII文本格式，则必须在格式字符串中加上字符t，例如用r t表示以ASCII格式打开供读操作的数据文件。（2）fclose关闭文件文件在进行完读、写等操作后，应及时关闭，以保证文件的安全可靠。关闭文件命令格式为：Sta=fclose(Fid) 关闭Fid所表示的文件Sta表示关闭文件操作的返回代码，若关闭成功，返回0，否则返回1。2、文件的读写操作（1）二进制数据文件fread 读二进制数据文件。格式为：A,COUNT=fread(Fid,size,precision)其中A为数据矩阵，COUNT返回所读取的数据元素个数。size为可选项，若不选用则读取整个文件内容，若选用它的值可以是下列值：N 读取 N个元素到一个列向量。inf 读取整个文件。M,N 读