怎么用MATLABphp 从数组中随机取值生成一个取值在(5，80)的四阶整数矩阵
点击联系发帖人
时间：2018-07-06 14:50
产生一个 3x5 随机整数矩阵,每个元素的范围在-5 到5 之间;给出每个元素的符号(+1 -1)表示，用MATLAB做_百度知道
产生一个 3x5 随机整数矩阵,每个元素的范围在-5 到5 之间;给出每个元素的符号(+1 -1)表示，用MATLAB做
A=round(rand(3,5)*10-5） ;//产生随机矩阵B=sign(A)//产生符号矩阵
采纳率：71%
来自团队：
为您推荐：
其他类似问题
换一换
回答问题，赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。怎样用MATLAB生成一个4行5列的间的随机矩阵_百度知道
怎样用MATLAB生成一个4行5列的间的随机矩阵
我有更好的答案
rand(4,5)&%&四行五列(0,1)之间的随机矩阵%&如果范围不在0-1之间，可以相应调整ceil(rand(4,5)*10)&%&[1,10]之间随机整数
采纳率：91%
为您推荐：
其他类似问题
随机矩阵的相关知识
换一换
回答问题，赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。当前位置： >>
Matlab 2014软件教程(完美版)
Matlab 软件目录1、Matlab 帮助的使用 ........................................................................ 3 2、Matlab 数据输入与类型 ................................................................ 4 3、Matlab 中的 M 文件及程序调试 ................................................ 12 4、Matlab 绘图命令 .......................................................................... 17 5、Matlab 在高等数学中的应用 ...................................................... 34 6、Matlab 在线性代数中的应用 ...................................................... 60 7、Matlab 数据处理 .......................................................................... 67 9、评价方法 ...................................................................................... 82 10、预测方法 .................................................................................... 97 11、蒙特卡洛方法 .......................................................................... 110 12、智能算法 .................................................................................. 122 13、分形 .......................................................................................... 129 14、Simulink 初步 .......................................................................... 134 15、Matlab 在概率统计中的应用 .................................................. 147 参考文献 .......................................................................................... 1511 作为和 Mathematica 、 Maple 并列的三大数学软件。 其强项就是其强大的矩阵计算以 及仿真能力。要知道 Matlab 的由来就是 Matrix + Laboratory = Matlab，所以这个软件在 国内也被称作《矩阵实验室》 。每次 MathWorks 发布 Matlab 的同时也会发布仿真工具 Simulink。在欧美很多大公司在将产品投入实际使用之前都会进行仿真试验，他们所主 要使用的仿真软件就是 Simulink。 Matlab 提供了自己的编译器： 全面兼容 C++以及 Fortran 两大语言。所以 Matlab 是工程师，科研工作者手上最好的语言，最好的工具和环境。 Matlab 已经成为广大科研人员的最值得信赖的助手和朋友！ 目前 MATLAB 产品族可以用来进行： － 数值分析 － 数值和符号计算 － 工程与科学绘图 － 控制系统的设计与方针 － 数字图像处理 － 数字信号处理 － 通讯系统设计与仿真 － 财务与金融工程... Simulink 是基于 MATLAB 的框图设计环境，可以用来对各种动态系统进行建模、 分析和仿真，它的建模范围广泛，可以针对任何能够用数学来描述的系统进行建模，例 如航空航天动力学系统、卫星控制制导系统、通讯系统、船舶及汽车等等，其中了包括 连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等等。 Simulink 提供 了利用鼠标拖放的方法建立系统框图模型的图形界面，而且 Simulink 还提供了丰富的 功能块以及不同的专业模块集合，利用 Simulink 几乎可以做到不书写一行代码完成整 个动态系统的建模工作。2 1、Matlab 帮助的使用1.1 helphelp ?????????????%帮助总览 help elfun ???%关于基本函数的帮助信息 help exp ????? %指数函数 exp 的详细信息1.2 lookfor 指令当要查找具有某种功能但又不知道准确名字的指令时， help 的能力就不够了， lookfor 可以根据用户提供的完整或不完整的关键词，去搜索出一组与之相关的指令。 lookfor integral ?? %查找有关积分的指令 lookfor fourier ??%查找能进行傅里叶变换的指令1.3 超文本格式的帮助文件在 Matlab 中，关于一个函数的帮助信息可以用 doc 命令以超文本的方式给出，如?doc ??? doc doc ??? doc eig ????%eig 求矩阵的特征值和特征向量1.4 pdf 帮助文件可从 MathWorks 网站上下载有关的 pdf 帮助文件。 网站地址：http://www.mathworks.com/3 2、Matlab 数据输入与类型2.1 Matlab 中的变量MATLAB 程序中的基本数据单元称为阵列(Array)，是一个分为行与列的数据集合。 变量被看作是只有一行一列的阵列。MATLAB 语言不需要对变量进行事先声明，也不 需要指定变量类型，它会自动根据所赋予变量的值或对变量所进行的操作来确定变量的 类型。其命名规则为： （1）变量名的大小写是敏感的。 （2）变量的第一个字符必须为英文字母，而且不能超过 31 个字符。 （3）变量名可以包含下划线、数字，但不能为空格符、标点。 （4）命名变量时可以取一个容易记忆并且能表达出其含义的名称，如汇率，可以 定义为 exchange_rate。 对于变量作用域，默认情况是局部变量，使用 global 定义全局变量，而且全局变 量常用大写的英文字母表示。 MATLAB 预定义的变量如下表所示： ans eps pi inf NaN i 或 j nargin nargout realmax realmin flops 预设的计算结果的变量名 MATLAB 定义的正的极小值=2.2204e-16 内建的π 值 ∞值，无限大 无法定义一个数目 虚数单位 i=j=√-1 函数输入参数个数 函数输出参数个数 最大的正实数 2 最小的正实数 2 浮点运算次数1023-10224 注 1：在定义变量时要尽量与避免与这些名字相同，以免改变它们的值，如果已经 改变，可以通过 clear 复这些值。 注 2：数字的输入输出格式。所有数据均按 IEEE 浮点标准的长型格式存储。输入 格式沿用了 C 语言的风格和规则；输出格式使用 format 数据格式 命令控制，只影响 在屏幕上的显示结果，不影响内部的存储和运算。 变量名 来恢复它的初始值，也可以通过重新启动 MATLAB 恢2.2 向量及其运算1、向量的生成 ①命令窗口直接输入，使用[ ]，元素之间用空格、逗号或者分号隔开。 ②使用冒号表达式，基本形式为 x=x0:step:xn，其中 xn 为尾元素数值限，而不一定 是尾元素的值。当 step=1 时可省略步长。 ③生成线性等分向量，使用 linspace 函数。Y=linspace(x1,x2,n) ④生成对数等分向量，使用 logspace 函数。Y=logspace(x1,x2,n) 2、向量的基本运算5 数加（减） 、数乘、点积(dot 函数)、叉积 (cross 函数)、混合积 dot(a,cross(b,c))2.3 矩阵及其运算1、简单矩阵的输入 （1）要直接输入矩阵时，矩阵一行中的元素用空格或逗号分隔；矩阵行与行之间 用分号“;”隔离，整个矩阵放在方括号“[ ]”里。 A=[1,2,3;4,5,6;7,8,9] ??? 说明：指令执行后，矩阵 A 被保存在 Matlab 的工作空间中，以备后用。如果用户 不用 clear 指令清除它，或对它进行重新赋值，那么该矩阵会一直保存在工作空间中， 直到本次指令窗关闭为止。 （2）矩阵的分行输入，此时回车键作为分行标志， A=[1,2,3 4,5,6 7,8,9] （3）使用 M 文件创建大矩阵，当矩阵维数非常大时，可以创建 m 文件，在 m 文 件中输入数据或者导出数据文件。 2、矩阵的基本运算 ①矩阵的四则运算。其中乘法运算要注意相乘的双方有相邻公共维，除法分为左除 “\” (A\B=inv(A)*B)和右除”/”（ A/B=A*inv (B)）(需要计算逆矩阵 ) ②矩阵的逆运算。inv 函数。 ③矩阵的幂运算。^ 。 ④矩阵的指数运算。exp(返回每个元素的指数值)，expm([V,D] = EIG(X) and EXPM(X) = V*diag(exp(diag(D)))/V)，expm1 (exp(x)-1) ⑤矩阵的对数运算。logm ⑥矩阵的特征值函数。eig 和 eigs(适合于大型稀疏方阵) ⑦矩阵的奇异值函数。svd（[U,S,V] = SVD(X)， X = U*S*V'）和 svds ⑧矩阵的条件数函数。cond （矩阵 A 的条件数等于 A 的范数与 A 的逆的范数的 乘积，c = cond(A,p)等价于 norm(A,p) * norm(inv(A),p)） ，condest(1 范数的条件数的估计 值)， rcond ⑨特征值的条件数函数。codeig（[V,D,s] = condeig(A) 等价于 [V,D] = eig(A); s6 =condeig(A);） ⑩范数函数。norm （1- 范数：即列范数，矩阵的各列绝对值之和的最大值；2- 范数： 所有元素的平方和开根号(默认)；无穷范数：即行范数，矩阵各行的绝对值之和的最大 值） ， normest (矩阵的 2 范数的估计值) 其他还有秩函数 rank， 迹函数 trace ， 零空间函数 null （又称为核空间， X=null(A)， 则 A*X=0 ，X? *X=I ） ，正 交 空 间 函 数 orth(B = orth(A),则 B'*B = eye(rank(A)))，伪逆函数 pinv 等。 3、特殊向量和特殊矩阵 （1）特殊向量 t=[0:0.1:10] %产生从 0 到 10 的行向量，元素之间间隔为 0.1 t=linspace(n1,n2,n) %产生 n1 和 n2 之间线性均匀分布的 n 个数 (缺省 n 时,产生 100 个数) t=logspace(n1,n2,n) ( 缺省 n 时,产生 50 个数) %在和之间按照对数距离等间距产生 n 个数。 （2）特殊矩阵 i)单位矩阵 eye(m), eye(m,n) 可得到一个可允许的最大单位矩阵而其余处补 0, eye(size(a)) 可以得到与矩阵 a 同样大小的单位矩阵。 ii) 所有元素为 1 的矩阵 ones(n)，ones(size(a))，ones(m， n)。 iii）所有元素为 0 的矩阵 zeros(n)， zeros(m,n)。 iv）空矩阵是 一个特殊矩阵，这在线性代数中是不存在的。例如 q=[ ] 矩阵 q 在工作空间之中，但它的大小为零。通过空矩阵的办法可以删除矩阵的行与 列。例如 a(:,3)=[] 表示删除矩阵 a 的第 3 列。 v）随机数矩阵7 rand(m,n) 产生 m× n 矩阵，其中的元素是服从[0,1]上均匀分布的随机数。 randint(m,n,[min,max]) 产生 m× n 矩阵，其中的元素是[min,max]上的随机整数。 normrnd(mu,sigma,m,n)产生 m× n 矩阵，其中的元素是服从均值为 mu，标准差 为 sigma 的正态分布的随机数。 exprnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的指数分布的随 机 数。 poissrnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的泊松 （Poisson） 分布的随机数。 unifrnd(a,b,m,n) 产生 m× n 矩阵， 其中的元素是服从区间[a,b]上均匀分布的随机数。 r = mvnrnd(MU,SIGMA,cases) 产生 cases 对均值向量为 MU，协方差阵为 SIGMA 的多维正态分布的随机数。 vi）随机置换 randperm(n)产生 1 到 n 的一个随机全排列。 perms([1:n])产生 1 到 n 的所有全排列。 vii) 稀疏矩阵 稀疏矩阵是指矩阵中零元素很多，非零元素很少的矩阵。对于稀疏矩阵，只要存放 非零元素的行标、列标、非零元素的值即可，可以按如下方式存储 （非零元素的行地址，非零元素的列地址），非零元素的值。 在 Matlab 中无向图 和有向图邻接矩阵的使用上有很大差异。 对于有向图，只要写出邻接矩阵，直接 使用 Matlab 的命令 s parse 命令，就可以把 邻接矩阵转化为稀疏矩阵的表示方式。 对于无向图，由于邻接矩阵是对称阵，Matlab 中只需使用邻接矩阵的下三角元素， 即 Matlab 只存储邻接矩阵下三角元素中的非零元素。 稀疏矩阵只是一种存储格式。Matlab 中，普通矩阵使用 sparse 命令变成稀疏矩阵， 稀疏矩阵使用 full 命令变成普通矩阵。 例 1 a=zeros(5); a(1,[2,4])=[3,4]; a(3,[2:4])=[1 3 8]; a(5,[1,5])=[6,7]8 b=sparse(a) %普通矩阵转化成稀疏矩阵 c=full(b) %稀疏矩阵转化成普通矩阵 其他一些特殊矩阵如下表所示： 函数 compan gallery hadamard hankel hilb invhilb 功能 伴随阵 Higham 测试阵 Hardamard 矩阵 Hankel 矩阵 Hilbert 矩阵 反 Hilbert 矩阵 函数 magic rosser toeplitz pascal vander wilkinson 功能 魔方阵 经典对称特征值测试阵 Toeplitz 矩阵 Pascal 矩阵 范德蒙矩阵 Wilkinson’s 特征值测试矩阵4、矩阵的特殊操作 ①变维。 有两种方法，使用冒号（:）和使用函数 reshape 使用“:”表达式对两个矩阵进行变维操作，需要预先定义两个矩阵的维数（例如 A(:,3)=[]） ； reshape 有两种形式，分别为 reshape(X,M,N)和 reshape(X,M,N,P…) ②变向 主要函数如下表所示： 函数 fiplr fipud fipdim Rot90 ③矩阵的抽取 对角线元素抽取函数 diag(X,k)/diag(v,k)，抽取矩阵 X 的第 k 条对角线的元素向量/ 使得向量 v 为所得矩阵的第 k 条对角线元素。 上三角元 素抽取 tril(X,k)和下三角元素抽取 triu(X,k) ④扩展 两种方法：利用对矩阵标示块的赋值命令 X(m1:m2,n1:n2)=a 生成大矩阵， 其中 m2- m1+1 必须等于 a 的行维数， n2- n1+1 必须等于 a 的列维数，生成 m2 × n2 维的 功能 矩阵左右翻转 矩阵上下翻转 矩阵特定维翻转 矩阵反时针 90 翻转 函数 diag tril triu 功能 产生或提取对角阵 产生下三角 产生上三角9 矩阵 X；利用小矩阵组合生成大矩阵，要严格注意矩阵大小的匹配。 （示例： X(1:10,1:10)=zeros(10,10)，LX=[X,X;X,X]） Matlab 中冒号（： ）的使用方法小结： （1）用于生成向量，a:b，一般要求 a&b，否则生成空矩阵。 （2）如果 b-a 不是整数时，则向量的最后一位数为 n+a，其中 n=fix(b-a)(向零取整) （3）a:c:b 表示[a,a+c,…,a+n*c]，其中 n=fix((b-a)/c)，当 c&0 且 a&b 或者 c&0 且 a&b 时出现空矩阵。 （4）A(:)以一列的方式显示 A 中所有的元素 （5）b=A(i,:)表示将 A 中第 i 行存入 b 中 （6）b=A(:,j)表示将 A 中第 j 列存入 b 中 （7）b=A(j:k)表示将 A 中第 j 到第 k 个元素存入 b 中（ Matlab 中矩阵按列存储） （8）b=A(:,c:d)表示将 A 中第 c 列到第 d 列存入 b 中，要求 c，d 不能超过 A 的列 数。 （9）当矩阵很大时，不知道矩阵的维数，可以使用 end 作为矩阵的最后一行或者 一列或者最后一个元素。例如 b=A(1:2,2:end)获取矩阵 A 右上角的元素。2.4 数组及其运算数组与矩阵在形式上完全一致，只是运算与矩阵不同。同型矩阵之间的运算通常称 为数组运算。 （矩阵的数组运算） 1、基本数组运算 ①四则运算。数组的乘除法是指两个同维数组间对应元素之间的乘除法，运算符 为”.*” ,”./”和” .\ ”。数组与常数之间的运算可以加”.”，也可以不加。 ②幂运算。.^,对每个数组元素的幂运算。 ③指数运算 exp,对数运算 log 和开方运算 sqrt。 2、数组函数运算 只要把运算的数组带入到函数中就可以了，通用 形式为 funname(A) 3、数组的逻辑运算和关系运算 指令 & 含义 小于 函数名 lt10 &= & &= == ~= & | ~小于等于 大于 大于等于 等于 不等于 逻辑 与 逻辑 或 逻辑 非le gt ge eq ne and or not指令 xor any all isnan含义 不相同就取 1，否则取 0 只要有非 0 就取 1，否则取 0 全为 1 取 1，否则为 0 为数 NaN 取 1，否则为 0指令 isequal ismember isempty isletter含义 相等取 1，否则取 0 两个矩阵是属于关系取 1， 否则取 0 矩阵为空取 1，否则取 0 是字母取 1， 否则取 0 （可以是字符 串）isinf isfinite ischar find isnumeric为数 inf 取 1，否则为 0 有限大小元素取 1，否则为 0 是字符串取 1，否则为 0 寻找非零元素坐标 判断数值矩阵isstudent isprime isreal isspace islogical学生版取 1 质数取 1，否则取 0 实数取 1，否则取 0 空格位置取 1，否则取 0 判断逻辑数组11 3、Matlab 中的 M 文件及程序调试M 文件分为两种，一种是脚本文件，由一系列 Matlab 的命令组成，可以直接运行； 一种是函数文件，必须由其他 M 文件或者在命令行窗口中调用执行。函数文件具有一 定的通用性，并且可以进行递归调用。3.1、Matlab 中流程控制语句（1）if 语句 有三种形式：if(表达式 ) 语句组 A; end if(表达式) 语句组 A，else 语句组 B，end if(表达式 1) 语句组 A，elseif(表达式 2) 语句组 B, else 语句组 C，end （2）循环语句 while（表达式）语句组 A , end for k=初值：增量：终值 语句组 A，end 循环语句可以结合 break 命令来控制程序的执行顺序。 （3）switch 语句 switch 表达式(标量或者字符串 ) case 值 1 语句组 A case 值 2 语句组 B … otherwise 语句组 N end 用法示例： 示例一：输入数 n，判断其奇偶性 解 Matlab 程序如下： n=input(? n=?); if rem(n,2)==0 %求余数，与 mod 函数的区别，当同符号时两者相同，否则不同 A=?even? ; else A=?odd?,12 end 示例二：求 Matlab 中的最大实数 解：Matlab 程序如下： x=1; while x~=inf, x1=x; x=2*x; %为了获取更接近 realmax 的值，可以改为 x=1.01*x end x1 示例三：求 Matlab 的相对精度 解：matlab 程序如下： y=1; while 1+y&1 y1=y; y=y/2; end y13.2 脚本文件主要特征如下： (1)一般由 clc,clear 命令开始，清除掉屏幕和工作空间中原有的变量和图形，以避免 其他已执行的程序残留数据对本程序的影响。程序中应该添加有注释。 (2)接下来是程序的主体，如果文件中有全局变量，则需要使用 Global 在程序的起 始部分注明。为了提高程序的可读性，注意语句的缩进。 (3)整个程序应按照 Matlab 标识符的要求起文件名，扩展名为 m。 示例：列出求素数的程序 解 Matlab 程序如下： clear,clc N=input('N='),x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end x13 3.3 函数文件函数文件与脚本文件的区别： (1)由 function 开头，后跟的函数名必须与文件名相同； (2)有输入变量和输出变量，可进行变量传递； (3)除非使用 Global 声明，程序中的变量均为局部变量，运行后不保存在工作空间 中。 示例 1：使用函数文件写出求素数的程序 解：Matlab 程序如下： function y=qiuprime(N) %求出 1 到 N 之间的素数，返回素数组成的向量 x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end y=x; 示例 2：写一个递归函数求 n! 解： function f=factor(n) if n&=1 f=1; else f=factor(n-1)*n; end end3.4 函数参数的可调性Matlab 中有两个永久变量 nargin 和 nargout 分别记录调用该函数时的输入实参和输 出实参的个数。只要在函数中包含这两个变量，就可以准确的知道该函数文件被调用时 的输入输出参数个数，从而决定函数如何进行处理。其他类似的变量还有 varargout， varargin。14 示例 1： function fout=examp(a,b,c) if nargin==1 fout=a; elseif nargin==2 fout=a+b; elseif nargin==3 fout=a+b+c; end 示例 2： function [x,y,z]=examp(a,b,c) if nargout==1 x=a; elseif nargout==2 x=a; y=b; elseif nargout==3 x=a; y=b; z=c; end end 示例 3： function vartest(argA, argB, varargin) optargin = size(varargin,2); %可选输入 stdargin = nargin - %标准输入 fprintf('Number of inputs = %d\n', nargin) fprintf(' Inputs from individual arguments(%d):\n', ... stdargin) if stdargin &= 1 fprintf(' %d\n', argA) end if stdargin == 2 fprintf(' %d\n', argB) end fprintf(' Inputs packaged in varargin(%d):\n', optargin) for k= 1 : size(varargin,2) fprintf(' %d\n', varargin{k}) end end15 3.5 程序调试M 文件编辑器中提供了强大的程序调试功能，使用方式如同 VC 程序中的调试器。 在 M 文件编辑器中的工具栏中有如下程序调试图标：这些图标的使用和菜单 Debug 中的一些命令相同。16 4、Matlab 绘图命令4.1 二维绘图命令二维绘图的基本命令有 plot，loglog，semilogx，semilogy 和 polar。它们的使用方法 基本相同，其不同特点是在不同的坐标中绘制图形。plot 命令使用线性坐标空间绘制图 形；loglog 命令在两个对数坐标空间中绘制图形；而 semilogx(或 semilogy)命令使用 x 轴（或 y 轴）为对数刻度，另外一个轴为线性刻度的坐标空间绘制图形； polar 使用极 坐标空间绘制图形。 二维绘图命令 plot 为了适应各种绘图需要，提供了用于控制线色、数据点和线型的 3 组基本参数。它的使用格式如下： plot(x,y,‘color_point_linestyle ’) 该命令是绘制 y 对应 x 的轨迹的命令。 y 与 x 均为向量，且具有相同的元素个数。 用字符串，color_point_linestyle? 完成对上面 3 个参数的设置。 线色：r-red，g-green，b-blue ，w-white ，k-black，i-invisible ，y-yellow，c-cyan （青色） ， m-紫色。 数据点：.（圆点） ，o（小圆圈） ，x（叉号） ，+（加号） ，*（星号） ，s （square 方形） ， h（ hexagram 六角星） ，d（diamond 菱型） ，p（pentagram 五角星） ， v（下三角） ，^ （上 三角） ，&（左三角） ，&（右三角） 。 线型：-（实线） ，--（虚线） ，-.（点画线） ， :（点线） 。 当 plot(x,y)中的 x 和 y 均为 m× n 矩阵时，plot 命令将绘制 n 条曲线。 plot(t,[x1,x2,x3])在同一坐标轴内同时绘制三条曲线。 如果多重曲线 对应不同的 x 轴向量绘制，可使用命令 plot(t1,x1,t2,x2,t3,x3) 式中 x1 对应 t1，x2 对应 t2 等等。在这种情况下，t1，t2 和 t3 可以具有不同的元素 个数，但要求 x1， x2 和 x3 必须分别与 t1，t2 和 t3 具有相同的元素个数。 subplot 命令使得在一个屏幕上可以分开显示 n 个不同坐标系，且可分别在每一个 坐标系中绘制曲线。其命令格式如下 subplot(r,c,p)17 该命令将屏幕分成 r× c 个子窗口，而 p 表示激活第 p 个子窗口。窗口的排号是从 左 到右，自上而下。 下面对几个特殊的坐标系统 进行简要介绍： ①对数坐标曲线，主要有 semilogx,semilogy 和 loglog，前两个分别以 x 坐标和 y 坐标为对数坐标，后一个是双对数坐标。 例如： x=1:0.1*pi:2* y=sin(x); semilogx(x,y,'-*') x = 0:.1:10; semilogy(x,10.^x) x = logspace(-1,2); loglog(x,exp(x),'-s') grid on ②双纵坐标（双 y 轴坐标系）函数 plotyy，调用形式为： plotyy(X1,Y1,X2,Y2) plotyy(X1,Y1,X2,Y2,fun) ③极坐标绘图函数 polar 调用形式为：polar(theta,rho) 例 1：绘制极坐标下的平面曲线 ? ? a cos(b ? n? ) ，并讨论参数 a,b,n 对曲线的影响。 解 Matlab 程序如下： N=100; theta=linspace(0,2*pi,N); for i=1:2 a(i)=input('a='); b(i)=input('b='); n(i)=input('n='); rho(i,:)=a(i)*cos(b(i)+n(i)*theta); subplot(1,2,i),polar(theta,rho(i,:)); end 几个比较漂亮的极坐标系下的图 形： (1)蝴蝶图案 fun 可以是 plot、 semilogx、 semilogy 或 loglog plotyy(X1,Y1,X2,Y2,fun1,fun2) fun1 绘制(X1,Y1)， fun2 绘制(X2,Y2)18 t=0:0.01:36; f=exp(cos(t-pi/2))-2*cos(4*(t-pi/2))+sin((t-pi/2)/12).^5; polar(t,f,'r') (2)枫叶图案 t=-pi/2:0.05:1.5* f=100./(100+(t-pi/2).^8)*2-sin(7*t)-cos(30*t)/2; polar(t,f,'r') 在图形绘制完毕后， 执行如 下命令可以再在图中加入标题、 标号、 说明和分格线等。 这些命令有 title ，xlabel，ylabel，text，gtext，legend 等。它们的命令格式如下 title(“ My Title”)， xlabel(“My X-axis Label”) ，ylabel(“My Y-axis Label”)， xlabel({'first line';'second line'}) ， text(x, y,'Text for annotation')， gtext('Text for annotation')， grid on gtext 命令是使用鼠标器定位的文字注释命令。当输入命令后，可以在屏幕上得到 一个光标， 然后使用鼠标器控制它的位置。 按鼠标器的左键， 即可确定文字设定的位置。 hold on 是图形保持命令，可以把当前图形保持在屏幕上不变，同时在这个坐标系 内绘制另外一个图形。 hold on 命令是一个交替转换命令， 即执行一次， 转变一个状态 （相 当于 hold on、 hold off） 。 还可以设置坐标轴的范围，使用命令 axis ，其格式为： axis([xmin xmax ymin ymax]) axis([xmin xmax ymin ymax zmin zmax cmin cmax]) 示例： x = 0:.01:pi/2; plot(x,tan(x),'-ro') 图形效果如下所示：默认情况下，横纵坐标轴的范围是根据函数自变量和因变量的值自动变化的，有时 效果不好，此时需要设定横纵坐标轴的范围：axis ([0 pi/2 0 5])，效果如上右图所示。19 解 画图的 Matlab 程序如下，画出的图见图 1。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x); plot(x,y1,'.-'); hold on %图形保持命令 plot(x,y2,'*-'); plot(x,y3,'-o'); h=legend('sin($x$)','sin($x+\frac{\pi}{3}$)','cos($x$)') %latex 格式显示 set(h,'Interpreter','latex') %设置 Interpreter 的属性值为 latex ，可以使用数学公式 xlabel('$x$','Interpreter','latex') %latex 格式显示 ylabel('$y$','Interpreter','latex') %latex 格式显示图1解 画图的 Matlab 程序如下，画出的图形结果见图 2。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x);20 subplot(3,1,1), plot(x,y1,'.-'), title('sin($x$)','Interpreter','latex') subplot(312), plot(x,y2,'*-') title('sin($x+\frac{\pi}{3}$)','Interpreter','latex') ylabel('$y$','Interpreter','latex') %latex 格式显示 subplot(313), plot(x,y3,'-o') title('cos($x$)','Interpreter','latex') %latex 格式显示 xlabel('$x$','Interpreter','latex') %latex 格式显示4.2 复数绘图plot(z) ： z 为复数时相当于 plot(real(z),imag(z))；如果是双变量如 plot(t,z)，则 z 中 的虚数部分将丢弃。在复平面中绘出多条曲线，必须使用 hold on 命令，或者把多条曲 线的实部和虚部明确的写出。即 plot(real(z1),imag(z1),real(z2),imag(z2)) Matlab 中专门用来绘制复变量函数图形的相关命令是 cplxmap， cplxgrid， cplxroot， 格式如下： z = cplxgrid(m) %产生(m+1)*(2m+1)的极坐标下的复数数据网格 cplxmap(z,f(z),(optional bound)) %绘制复变函数的图形，以 xy 平面为自变量所在的 复平面，以 z 轴表示复变函数的实部，颜色表示复变函数的虚部。 cplxroot：画复数的 n 次函数曲面21 cplxroot(n) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面 cplxroot(n,m) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面,为（ m+1）* （2m+1）的方阵 例 1：绘制 z=exp((-0.1+i)*t)的复数图形 解 Matlab 程序如下： t=0:0.1:15; z=exp((-0.1+i)*t); subplot(2,2,1) plot(z), % pause 可以暂停程序的执行， 通过按任意键可使程序继续进行， 也 可以在 pause(n)中设置时间使执行结果出现动态效果。 title('复数绘图 plot(z)'); subplot(2,2,2) plot(t,z),pause title('plot(t,z)') subplot(2,2,3) polar(angle(z),abs(z)); title('polar(angle(z),abs(z))') subplot(2,2,4) semilogx(t,z); title('semilogx(t,z)') 注：Matlab 中使用函数 C=complex(A,B)构造复数。 例 2：绘图 y=cos(x+i)的图形 解法一：x=-pi:0.1: fun=@(x)cos(x+i); %如果直接使用 plot(x,y) ，会忽略虚部。 plot(fun(x)) 解法二：plot(cos(x+i)); 例 3：绘制 f=z^4 的图形，其中 z 为复数 解：z=cplxgrid(30); cplxmap(z,z.^4)4.3 显函数，符号函数或隐函数的绘图fplot(fun,lims)绘制由字符串 fun 指定函数名的函数在 x 轴区间为 lims=[xmin, xmax] 的函数图。若 lims=[xmin,xmax,ymin,ymax]，则 y 轴也被限制。22 解 （1）首先用 M 文件 fun1.m 定义函数 f(x)如下 function y=Afun1(x); if x&1 y=x+1; else y=1+1./x; end 在 matlab 命令窗口输入 fplot('Afun1',[-3,3]) 就可画出 函数 f (x ) 的图形。 （2）可以使用匿名函数，编写程序如下 fun2=@(x) (x+1)*(x&1)+(1+1/x)*(x&=1); fplot(fun2,[-3,3]) ezplot(f)绘制符号函数或者隐函数 f(x)的图形，x 轴的近似范围为[2? ,?2??。 ezplot(f,[xmin,xmax])使用输入参数来代替默认横坐标范围 [2? ,?2??。 ezplot 函数的其 他格式有： ? ezplot(fun2)绘制 fun2(x,y)=0 的隐函数曲线，默认 x,y 的范围是[2? ,?2??? ezplot(fun2,[xmin,xmax,ymin,ymax]) 示例： function z = myfun(x,y,k) %建立 M 文件 z = x.^k - y.^k - 1; %在命令窗口中输入：ezplot(@(x,y)myfun(x,y,2));解 ezplot('cot(x)')解 ezplot('x^2+y^2/4=1') 例 6-1：同一坐标系下绘制曲线 x^2+y^2=1 和 x^2-y^2=1 的曲线 解：ezplot('x^2-y^2=1');23 ezplot('x^2+y^2=1'); colormap([0,0,1]);4.4 三维图形在实际工程计算中，最常用的三维绘图是三维曲线图、三维网格图和三维曲面图 3 种基本类型。 与此对应， Matlab 也提供了一些三维基本绘图命令， 如三维曲线命令 plot3， 三维网格图命令 mesh 和三维表面图命令 surf。 1. 三维曲线 plot3(x,y,z)通过描点连线画出曲线， 这里 x,y,z 都是 n 维向量， 分别表示该曲线上点 集的横坐标、纵坐标、竖坐标。 例7 t=0:0.05:20* x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z), title('Line in 3-D Space') xlabel('X'), ylabel('Y'), zlabel('Z'), grid on 2. 三维网格图 命令 mesh(x,y,z)画网格曲面。这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标，命令 mesh(x,y,z)将该数据点在空间中描出，并连成网格。解 x=-3:0.1:3;y=-5:0.1:5; x1=ones(size(y'))*x;y1=y'*ones(size(x)); [x2,y2]=meshgrid(x,y); %meshgrid 函数生成 2D/3D 网格矩阵，用于分割空间，在绘 制 3D 网格图或者表面图时都需要调用该函数对生成绘图时所用的数据。 z1=(sin(x1.*y1)+eps)./(x1.*y1+eps); z2=(sin(x2.*y2)+eps)./(x2.*y2+eps);24 subplot(1,2,1),mesh(x1,y1,z1) subplot(1,2,2),mesh(x2,y2,z2) 例 8-1：绘制两个空间相交的平面 z=x+2y-1 和 z=x- y+2 解 ： x=-3:0.1:3;y=-3:0.1:3; [X,Y]=meshgrid(x,y); mesh(X,Y,X+2*Y-1) hold on mesh(X,Y,X-Y+2) 3. 表面图 命令 surf(x,y,z)画三维表面图，这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标。解 [x,y]=meshgrid([-3:0.2:3]); z=(sin(x.*y)+eps)./(x.*y+eps); surf(x,y,z) 4. 旋转曲面 方法一：使用命令函数 cylinder [X,Y,Z] = cylinder(r) 这里的 r 表示构成旋转曲面的曲线。解 Matlab 程序如下。 x=0:10:600; [X,Y,Z]=cylinder(30*exp(-x/400).*sin((x+25*pi)/100)+130); surf(X,Y,Z) 方法二：将旋转曲面用参数方程表示。25 解 因为这里的函数是隐函数，化成显函数后有两支，必须使用参数方程，旋转面 的参数方程为画图的 Matlab 程序如下： alpha=[0:0.1:2*pi]'; beta=0:0.1:2* x=4*cos(alpha)*ones(size(beta)); y=(5+4*sin(alpha))*cos(beta); z=(5+4*sin(alpha))*sin(beta); surf(x,y,z) 或者利用绘制三维隐函数曲面图形的命令 ezsurf 或者 ezmesh，其命令格式如下： ezmesh(fun) ezmesh(fun,domain) domain 要求是一个向量，指定 x,y 轴的范围 ezmesh(funx,funy,funz) %参数方程绘图 funx(s,t), funy(s,t), and funz(s,t) over the square: - 2π & s & 2 π , -2π & t & 2 π. ezsurf(fun) ezsurf(fun,domain) ezsurf(funx,funy,funz) 画图的 Matlab 程序也可以写成 x=@(alpha,beta) 4*cos(alpha); y=@(alpha,beta) (5+4*sin(alpha))*cos(beta); z=@(alpha,beta) (5+4*sin(alpha))*sin(beta); ezsurf(x,y,z) 对于形如 f(x,y,z)=0 的三维隐函数的图像没有现成的函数可以画，但可以利用 isosurface 函数绘制三角网格图。例如 f=@(x,y,z)x.^2+y.^2+z.^2-10;%定义函数 f=x^2+y^2+z^2-10 [x,y,z]=meshgrid(linspace(-4,4,25));%设定格子大小和范围 val=f(x,y,z); [p,v]=isosurface(x,y,z,val,0);%用 isosurface 得到函数 f=0 图形的点和面26 patch('faces',p,'vertices',v,'facevertexcdata',jet(size(v,1)),'facecolor','w','edgecolor','flat'); %用 patch 绘制三角网格图并设定色彩 view(3); axis equal 5、绘制柱面 柱面平行于某个坐标轴，方程中不出现某个坐标轴的变量，方程表示为 F(x,y)=0 或 者 F(x,z)=0 或者 F(y,z)=0 示例一：画出方程 z ? 2 ? x2 表示的柱面 解：方程中的 x 是自变量矩阵， z 是因变量，则另一个自变量矩阵为 y，自变量平 面是 xoy 面， x 轴是真正的自变量， y 轴是柱面方向。Matlab 程序如下： u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X=u*ones(size(v)); %构成自变量矩阵 X Y=ones(size(u))*v; %构成自变量矩阵 Y Z=2-X.^2; %求因变量 Z mesh(X,Y,Z) 示例二：画出方程 ? x 2 ? y2 ? 1 表示的柱面 4解：方程整理为显示函数： y ? ?2 1 ? x 2 ，分别绘制正负两个分量。 clear,clc u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X1=u*ones(size(v)); %构成自变量矩阵 X1 Y1=2*sqrt(1+X1.^2); %求因变量 Y1 Y2=-2*sqrt(1+X1.^2); Z1=ones(size(u))*v; %求因变量 Z1 mesh(X1,Y1,Z1),hold on mesh(X1,Y2,Z1) 6. 其它二次曲面 对于旋转面，如果母线的方程可以表示成关于旋转轴变量的显式函数，则可以直接 使用 Matlab 工具箱中的命令 cylinde r，否则必须把旋转面化成参数方程，然后使用 ezmesh 或 ezsurf 命令绘图。对于其它的二次曲面，如果可以写成显函数直接使用命令 ezmesh 或 ezsurf，否则必须先化成参数方程。还有一些特殊的绘制函数如27 cylinder,ellipsoid 等命令。解： （1）x=@(s,t) 3*sec(s)*cos(t); y=@(s,t) 3*sec(s)*sin(t); z=@(s,t) 2*tan(s); ezmesh(x,y,z) （2）x=@(s,t) 3*sec(s); y=@(s,t) 2*tan(s)*cos(t); z=@(s,t) 2*tan(s)*sin(t); ezmesh(x,y,z) （3）ezsurf(@(y,z) y.^2,50) %直接调用 ezsurf （4）x=@(s,t) 3*tan(s)*cos(t); y=@(s,t) 2*tan(s)*sin(t); z=@(s,t) tan(s); ezsurf(x,y,z) （5）ellipsoid(0,0,0,3,2,sqrt(6)) %专门绘制椭球面 （6）ezsurf(@(x,y)x*y) （7）x=@(s,t) 3*cos(s); %化为参数方程28 y=@(s,t) 2*sin(s); z=@(s,t) ezmesh(x,y,z) 7、绘制空间两曲面的交线 示例：绘制由水平截面与方程 z ? x 2 ? 2 y 2 构成的马鞍面形成的交线，并讨论等高线 和方向导数（梯度）的意义 解： 水平平面与曲面的交线就是等高线， 在 Matlab 中绘制等高线有两个命令， contour 和 contour3，前者把等高线画在 xoy 平面上，后者把等高线画在一定高度的平面上，使 之成为立体的，与所在曲面对应。 Matlab 程序如下： clc,clear [x,y]=meshgrid(-10:.2:10);%确定计算和绘图的定义域网格 z1=(x.^2-2*y.^2)+ %第一个曲面方程 a=input('a=(-50&a&50)'); z2=a*ones(size(x)); %水平面方程 z2=a， z2 必须与 x， y 具有相同的维数， subplot(1,3,1),mesh(x,y,z1);mesh(x,y,z2);%分别画出两个曲面 v=[-10 10 -10 10 -100 100];axis(v), %确定第一个分图的坐标系 colormap(gray); r0=abs(z1-z2)&=1; %求两曲面 z 坐标只差小于 0.5 的网格 zz=r0.*z2;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标值 subplot(1,3,2),plot3(xx,yy,zz,'x'); %画出这些点 axis(v), %使第二个分图取第一个分图的坐标系 pause,subplot(1,3,3); contour3(x,y,z1,20); %用等高线命令求出 20 条不同高度的交线。 等高线与方向导数和梯度的概念密切相关，函数 z1 在每一点的梯度与该处的等高 线垂直，也就是指向最陡的方向。Matlab 中求梯度用 gradient 函数，quiver 函数画出梯 度向量，这两个函数要求在给定的点阵上求梯度和画梯度向量。Matlab 程序如下： [x,y]=meshgrid(-10:2:10); z1=(x.^2-2*y.^2)+ contour(x,y,z1,20); [px,py]=gradient(z1,2,2); %以步长为 2 求 z1 的梯度的 x， y 分量 quiver(x,y,px,py); %绘制梯度向量 8、离散点的图形绘制相关函数 griddata，TripScatterInterp，scatter（scatter3） griddata 用来对离散数据进行曲面拟合 (1)ZI = griddata(x,y,z,XI,YI) 用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量29 x,y,z。griddata 将返回曲面 z 在点（ XI,YI）处的插值。曲面总是经过这些数据点（ x,y,z） 的。输入参量（XI,YI）通常是规则的格点（像用命令 meshgrid 生成的一样） 。 XI 可以 是一行向量，这时 XI 指定一有常数列向量的矩阵。类似地， YI 可以是一列向量，它 指定一有常数行向量的矩阵。 (2)[XI,YI,ZI] = griddata(x,y,z,xi,yi) 返回的矩阵 ZI 含义同上， 同时， 返回的矩阵 XI,YI 是由行向量 xi 与列向量 yi 用命令 meshgrid 生成的。 (3)[XI,YI,ZI] = griddata(.......,method) 用指定的算法 method 计算： ?linear ?：基于三角形的线性插值（缺省算法） ； ?cubic?： 基于三角形的三次插值； ?n earest? ：最邻近插值法； 示例： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(- x.^2-y.^2); ti = -2:.25:2; [xi,yi] = meshgrid(ti,ti); zi = griddata(x,y,z,xi,yi); mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off 在新版本的 Matlab 中，该函数逐渐被 TriScatteredInterp 函数所代替， TriScatteredInterp 的示例如下： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); F = TriScatteredInterp(x,y,z); ti = -2:.25:2; [qx,qy] = meshgrid(ti,ti); qz = F(qx,qy);%计算指定位置的插值 mesh(qx,qy,qz); plot3(x,y,z,'o'); Matlab 中绘制散点图的命令是 scatter（scatter3） ，命令格式如下： scatter(X,Y,S,C) %X，Y 对应散点值，S 代表标记大小，C 代表颜色值，它都可以是 向量。 scatter3(X,Y,Z,S,C) 示例：30 A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1; 4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1; 1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2; 4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2; 2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3; 2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4; 3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5; 3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6; 3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7]; x=A(:,1);y=A(:,2);z=A(:,3); scatter(x,y,5,z)%散点图4.5 动态可视化图形Matlab 中的动画命令， moviein，getframe 和 movie ，用 getframe 把 Matlab 产生的 图形存储下来，每个图形成一个很大的列向量；再用 N 行这样的列保存 N 幅图，成为 一个大矩阵；movie 命令把他们连接起来重放，产生动画效果；moviein 用来预留存储空 间以加快运行的速度。有些情况也需要借助 pause 命令和循环语句来实现点轨迹的动态 演示。绘制彗星图的命令 comet 或者 comet3 也可以实现点的运动轨迹，不过这种方式 速度较快。 示例 1： axis equal, %把坐标设成相等比例 M=moviein(16);%为变量 M 预留 16 幅图的存储空间 for j=1:16 plot(fft(eye(j+16))); M(:,j)= end movie(M,30) %以每秒 30 帧的速度播放 M 中的图形 示例 2：使用 comet/comet3 绘制彗星图 （1） t = 0:.01:2* %返回当前坐标轴下的一帧图像31 x = cos(2*t).*(cos(t).^2); y = sin(2*t).*(sin(t).^2); comet(x,y); （2） n=10; t=n*pi*(0:0.0005:1); x=sin(t); y=cos(t); plot(x,y,'g'); comet(x,y,0.01); 示例 3：极坐标下一个点沿着图形移动的轨迹 a=0:.01:2* b=3; polar(a,b*(1-cos(a)),'*r'); for a=0:0.05:2*pi cla(findobj(gca,'color','r')); %findobj 函数用来使用特定的属性获取图形句柄， cla 用来清理由函数句柄指定的图形所在的坐标轴。 h2=polar(a,b*(1-cos(a)),'b'); set(h2,'linestyle','.','markersize',30); %设置句柄指定函数图形的属性 pause(0.001) end 示例 4： 普通二维坐标下一个点沿着曲线移动的轨迹， 程序与上面的结果基本相同。 x=-2*pi:0.1:2* y=sin(x); plot(x,y,'b-'); hold on for m=-2*pi:0.05:2*pi cla(findobj(gca,'color','b')); %cla 清除当前坐标轴上面的子对象 h=plot(m,sin(m),'r'); set(h,'LineStyle','.','Marker','*','MarkerSize',20); pause(0.001); end 示例 5：普通三维坐标下一个点沿着曲线移动的轨迹。 解： %三维螺旋线坐标 t1=10*pi*(0:;32 x1=cos(t1); y1=sin(t1); z1=-t1; %绘制曲线下面的水平直线 t2=(0:10)/10; x2=x1(end)*(1-t2); y2=y1(end)*(1-t2); z2=z1(end)*ones(size(x2)); %绘制垂直直线 t3=t2; z3=(1-t3)*z1(end); x3=zeros(size(z3)); y3=x3; %绘制曲线上面的水平直线 t4=t2; x4=t4; y4=zeros(size(x4)); z4=y4; %绘制曲线 x=[x1,x2,x3,x4]; y=[y1,y2,y3,y4]; z=[z1,z2,z3,z4]; plot3(x,y,z,'b','linewidth',3) axis off %绘制移动的点 h=line('Color',[1 0 0],'Marker','.','MarkerSize',40,'EraseMode','xor'); n=length(x); i=1;j=1; while 1 set(h,'xdata',x(i),'ydata',y(i),'zdata',z(i)); drawnow %使 Matlab 暂停目前的任务序列而去刷新屏幕 pause(0.0005); f=getframe(gcf); i=i+1; j=j+1; if i&n & j&450 end end33 5、Matlab 在高等数学中的应用5.0 Matlab 中的符号表示1、建立符号变量 Matlab 提供了 sym 和 syms 两个建立符号对象的函数 1）sym 函数，用来建立单个符号变量，调用格式为： 符号变量名=sym(?符号字符串 ?) 注：符号字符串可以是 常量、变量、函数或者表达式。 2）syms 函数，依次可以定义多个符号变量，调用格式为： syms 符号变量 1 符号变量 2 例如： f=sym(? co s(x)?),f=sym(? sin(x)^2=0? ); f=sin(x)+cos(x), 2、建立符号表达式 ? ? ? 利用单引号来生成表达式 利用 sym 函数建立符号表达式 使用已经定义过的符号变量组成符号表达式 注意：符号表达 …… 符号变量 n 注意 符号变量名之间不能使用任何标点符号，只能用空格隔开。式包括符号函数和符号方程，区别在于是否带有等号。 3、符号表达式的运算 (1)四则运算 四则运算分别使用 symadd, symsub,symmul 和 symdiv 来实现加减乘除运算，使用 sympow 实现幂运算。 （注：6.5 以后的版本中已经没有了这些函数，而直接使用+ - * / 即可） (2)提取分子和分母 如果符号表达式是一个有理分式或可以展开为有理分式，可利用 numden 函数来提 取符号表达式中的分子或者分母，调用格式为：34 [n,d]=numden(s) [n,d] = numden(sym(4/5)) returns n = 4 and d = 5. [n,d] = numden(x/y + y/x) returns (3)因式分解与展开 1）factors (s)分解因式 2）expand(s)展开 3）collect(s)合并同类项 4）collect(s,v)按变量 v 进行合并同类项 (4)符号表达式的化简 1）simplify (S)应用函数规则对 S 进行化简 2）simple (S) 进行综合化简 (5) 符号表达式与数值表达式之间的转换 1）利用函数 sym 可以将数值表达式表示成符号表达式 2）numeric 或者 eval 函数可以将符号表达式转换成数值表达式 3）函数 digits 和 vpa 配合替换函数 subs 进行转换。 digits 函数，digits(D)函数设置有效数字个数为 D 的近似解精度。 vpa 函数 ， R=vpa(S)符号表达式 S 在 digits 函数设置下的精度的数值解。 vpa(S,D) 符号表达式 S 在 digits(D)精度下的数值解。 subs 函数， subs(S, NEW, OLD) 4、符号函数的相关运算 （1）复合函数运算。 compose 函数 （2）反函数运算。 finverse 函数 5、符号代数方程求解 (1)线性方程组的求解，函数 linsolve ，solve ，可以得到方程的精确解 (2)多变量的非线性方程的符号解法，使用函数 fsolve ，调用格式有： ? ? ? ? ? X=fsolve(?fun?, X0) X=fsolve(? f un?,X0,options)options 为选择参数输入向量 X=fsolve(? f un?, X0,o ptio ns,?g radfun ? ),gradfun 为输入函数在 X 处的偏导数 X=fsolve( ? f un?,X0,o ptio ns,?g radfun ? ,P1,P2,… )P1,P2 为问题定性参数 [X,optio n s]=fsolve( ?fun ? , X0,…) 返回使用的优化方法的参数 n = x^2+y^2 , d = y*x35 注意：fun 必须使用@指定 示例：定义一个函数 function F = myfun(x,c) F = [ 2*x(1) - x(2) - exp(c*x(1)) -x(1) + 2*x(2) - exp(c*x(2))]; end 求解该函数： c = -1; %定义输入参数 x = fsolve(@(x) myfun(x,c),[-5;-5]) 6、多项式的表示方法 （1）多项式的表示方法――转化为向量问题 对于多项式P ( x) ? a 0 xn ? a 1 x n ?1 ? ... ? a n ?1 x ? a n用行向量表示： P ? [a , a ,..., a , a ] 0 1 n? 1 n ①系数向量直接输入法，MATLAB 自动将向量元素按降幂顺序分配给各系数值。函数 poly2sym 可以将向量表示的多项式转化为符号多项式表示。 注：由特征多项式生成 的多项式首项系数一定为 1；n 阶矩阵一般产生 n 次多项式。 ③由根创建多项式，由函数 poly 实现。注：若要生成实系数多项式，则根中的复数必 定对应共轭；生成的多项式向量包含很小的虚部时可用 real 命令将其过滤掉。 （2）多项式的运算 ①多项式求值。输入变量值代入多项式计算时以数组为单元的使用函数 polyval （对应 元素计算）;以矩阵(必须为方阵 )为计算单元求多项式的值用函数 ②多项式求根。两种方法，一种是调用函数 roots ，另一种是通过建立多项式的伴随矩 阵再求其特征值的方法得到多项式的所有根。 （使用 compan 和 eig 函数） ③多项式的乘除法运算。乘法使用函数 conv（向量卷积） ，除法使用函数 deconv ④多项式微分。微分函数 polyder ⑤多项式拟合。两种方法，一种是由矩阵的除法求解超定方程来进行，另一种是用拟 合函数 polyfit,调用方式为 polyfit(X,Y,n)和[p,s]=polyfit(X,Y,n)36 5.1 求极限Matlab 求极限的命令为 limit(expr, x, a) limit(expr, a) limit(expr) limit(expr, x, a, 'left') limit(expr, x, a, 'right') 其中 limit(expr, x, a)表示求符号表达式 expr 关于符号变量 x 趋近于 a 时的极限， limit(expr) 求表示缺省变量趋近于 0 时的极限。 例 15 求下列表达式的极限解 （1）clc, clear syms x f=(3*sin(x)+x^2*cos(1/x))/((1+cos(x))*log(1+x)); s=limit(f,x,0) s=simplify(s)（2）clc, clear syms x a s1=limit(((x+2*a)/(x-a))^x,x,inf) s2=limit(((x+2*a)/(x-a))^x,x,-inf)5.2 求导数Matlab 的求导数命令为 diff(expr) diff(expr, v) diff(expr, sym('v')) diff(expr, n) diff(expr, v, n)37 diff(expr, n, v) 其中 diff(expr)表示求表达式 expr 关于默认变量的 1 阶导数，diff(expr, v, n)和 diff(expr, n, v)都表示求表达式，expr 关于符号变量 v 的 n 阶导数。解 （1）clc, clear syms x y=log(x+sqrt(1+x^2)); s=diff(y,x,2) s=simple(s) %对符号函数进行化简 （2）a=[1,0.5,3.5,6]; da=diff(a)5.3 求一元函数极值计算的 Matlab 程序如下 syms x y=x^3+6*x^2+8*x-1; dy=diff(y); dy_zero=solve(dy), dy_zero_num=double(dy_zero) %变成数值类型 ezplot(y) %符号函数画图5.4 求积分1. 求不定积分 Matlab 求符号函数不定积分的命令为 int(expr) int(expr, v)38 解 syms x I=int(1/(1+sqrt(1-x^2))); pretty(I) %写成符号表达式的形式 2. 求定积分 （1）求定积分的符号解 Matlab 求符号函数的定积分命令为 int(expr, a, b) int(expr, v, a, b)解 syms x I=int(cos(x)*cos(2*x),-pi/2,pi/2) （2）求定积分的数值解 例 20 求下列积分的数值解输入 I=quadl(@(t) (t-3*t.^2+2*t.^3).^(-1/3),eps,0.5) 得到 I=1.4396。 注：quad()函数使用 simpson 法，即用二次曲线逼近被积函数， quadl()采用更高阶 的逼近方法。它们的调用方法是：q=quad(fun,a,b,tol)或者 quadl(fun,a,b,tol)， fun 可以是39 符号方程，也可以是匿名函数。 ）输入 I=dblquad(@(x,y)sqrt(1- x.^2-y.^2).*(x.^2+y.^2&=x),0,1,-0.5,0.5) 得到 I= 0.6028。 注： Matlab 中计算二重积分的函数形式为 q= dblquad(fun,xmin,xmax,ymin,ymax,tol)， 一般只用于积分区域为矩形的情况。为了使用该函数，有时需要将非矩形区域化成矩形 区域。?2 2 示例：计算二重积分 ?? ?( x ? y )dxdy 积分区域 ? 由 x=1,y=x 及 y=0 围成。 ?解：首先绘制积分区域，Matlab 程序如下： clear,clc fill([0,1,1,0],[0,0,1,0],'y');hold on%绘制积分区域 fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r') %绘制单元条 gtext('y=x'); gtext('x=1');pause gtext('y=0') 按照矩形区域调用 dblquad 函数，程序如下： I=dblquad(?(x.^2+y.^2).*(y-x&0)?,0,1,0,1) 求得 I=0.3333 （3）输入 fun3=@(x,y,z)z.^2*log(x.^2+y.^2+z.^2+1)./(x.^2+y.^2+z.^2+1).*(z&=0 &... z&=sqrt(1- x.^2-y.^2)); %该语句使用了续行符... I=triplequad(fun3,-1,1,-1,1,0,1) 得到 I= 0.1273。 %三重积分计算函数40 Matlab 符号求解的程序如下 clc, clear syms r z theta I1=int((r^2+z^2)*r, z, r^2,sqrt(2-r^2)); %求最内层积分 I2=int(I1,r,0,1); %求中间层的积分 I=int(I2,theta,0,2*pi) %求最终的积分结果 pretty(I) % 分数线中间显示的格式 Inum=double(I) %把符号解化成数值解 解法二 求数值解的 Matlab 程序如下 clc, clear f=@(x,y,z) (x+y+z).^2.*(z&=x.^2+y.^2 & x.^2+y.^2+z.^2&=2); I=triplequad(f,-sqrt(2),sqrt(2),-sqrt(2),sqrt(2),0,sqrt(2)) 求得积分的值为 2.4486。 注：quad 与 int 的区别，int 可以求具有解析解的不定积分和定积分，quad 只能求解 定积分。5.5 求解微分方程（组）1、微分方程(组)的解析解41 在 MATLAB 中用 D 表示导数， 例如 Dy 表示 y?， D2y 表示 y ?? ， Dy(0)=5 表 示 y ?(0)=5 等，符号常微分方程求解通过函数 dsolve 实现，调用格式为： dsolve( ?e? ,? c? ,? v? )包括微分方程，初始条件和指定变量三个部分，求解常微分方程 e 在初始条件 c 下的特解，v 描述方程中的自变量，省略则默认以 t 为自变量，若没有初 始条件，则获得通解。 例 1：求 du ? 1 ? u 2 的通解 dt解：Matlab 程序为： dsolve('Du=1+u^2','t') ?d2y dy ? 2 ? 4 ? 29 y ? 0 例 2：求微分方程 ??dx 的特解 dx ? ??y (0) ? 0, y?(0) ? 15 解：Matlab 程序为： dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') ? dx ? dt ? 2x ? 3 y ? 3 z ?? ? dy 例 3：求常微分方程组的通解 ? ? 4 x ? 5 y ? 3z ??dt ? dz ? 4x ? 4 y ? 2 z ? ??dt 解：Matlab 程序为： [x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z') 一阶齐次线性常微分方程组和非齐次线性方程组可以直接通过矩阵操作求解。 对于一阶齐次线性微分方程组： X ? ? AX , X ? ? x1 , xn ?, A ? (aij )n ?n 其解的形式为：X (t ) ? eA(t ?t0 ) X 0对于非齐次线性方程组： X ? ? AX ? f (t ), X (t 0 ) ? X 0 解的形式为：X (t ) ? eA(t ?t0 ) X 0 ???? eA (t ? s ) f ( s )ds? 例 4： X ? ? ? 0 2 ?1? X , X (0) ? ? 2 ? ? ? ? ? ? ? ?0 0 3 ? ?1 ?2 1 3? ?1 ? ?42 解：a=[2,1,3;0,2,-1;0,0,2]; x0=[1;2;1];x=expm(a*t)*x0 pretty(x); %化简结果?1 0 0 ? ? 0 ? ?0 ?? ? ? ? ? ? ?? 例 5： X ? ? ? 2 1 ?2 ? X ? ? 0 ? , X (0) ? ?1 ?? t ? ? ? ?3 2 1 ?e cos(t ) ?1解： a=[1 0 0;2 1 -2;3 2 1]; fs=[0;0;exp(s)*cos(2*s)]; x0=[0;1;1]; tx=int(expm(a*(t-s))*fs,s); %先求不定积分，再计算定积分的结果，提高速度 xstar=subs(tx,s,t)-subs(tx,s,0); x=expm(a*t)*x0+ x=simple(s),pretty(x) %simple 直接显示化简结果，simplify 显示化简过程，pretty 修 改显示格式 2、微分方程（组）的数值解 当难以求得微分方程的解析解时，可以求其数值 解。Matlab 中求解微分方程数值解 的方法有以下几个：ode45，ode23，ode113，ode23s ，ode15s 。调用形式如下： [t,x]=solver(‘ f’, ts, x0, options)其中，t 是自变量的值，x 是函数值，solver 表示上述 5 种求解器的一种；‘ f’是由待解方程写成的 m 文件名，ts=[t0, tf]表示自变量的初值和 终值，options 用于设定误差限，通过函数 options=odeset(‘reltol’, rt, ‘abstol’, at)分别设 定相对误差和绝对误差。 不同的函数代表不同的内部算法，ode45 表示 4/5 阶龙格- 库塔- 费尔贝格算法（ RK 算法），是解非刚性常微分方程的首选方法，ode23 采用 2/3 阶 RK 方法，ode113 采用 多步法，效率一般比 ode45 高。ode15s ，ode23s,ode23t,ode23tb 用来解刚性常微分方程。 注 1：在解 n 个未知函数的的方程组时，x0 和 xn 均为 n 为向量， m 文件中的待解 方程组应以 x 分量形式写出； 注 2：使用 Matlab 求数值解时，高阶微分方程必须等价的变换成一阶微分方程组。 ? ? d2 x 2 dx ? ? ?1000(1 ??x ) ? x 0 例 1：求解 van der Pol 刚性微分方程 ??dt 2 的数值解 dt ? x (0) ? 2; x '(0) ? 0 ??43 解：令 y1=x， y2=y1?，则微分方程变为一阶微分方程组： ? y1' ? y 2 ? 2 ??y 2 ' ? 1000(1 ? y1 ) y 2 ? y1 ? y1(0) ? 2; y 2(0) ? 0 ?? 写成 m 文件如下所示： function dy=fun1(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1); 取 [t0 ,tf]=[0,3000] ， 输入命 令： [t,y]=ode15s('fun1',[0,3000],[2,0]); plot(t,y(:,1),'-') ? y1 ' ? y2 y3 ?? ? ? y1 y3 ?y ' 例 2：求解微分方程组 ?? 2 ? y3 ' ? ? 0.51 y1 y2 y1 (0) ? 0, y2 (0) ? 1; y 3 (0) ? 1 解：建立方程组的 m 文件： function dy=fun2(t,y) dy=zeros(3,1); %初始化 dy(1)=y(2)*y(3); dy(2)=y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 取 t0=0，tf=12，Matlab 命令如下： [t,y]=ode45('fun2',[0 12],[0 1 1]) plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+')? ? x ? ? ( y ??x) ?? ? y ? ? x ? y ? xz 在初始条件[x0,y0,z0]=[5,13,17]时的解，其 例 3：求解 Lorez 方程 ? ? ? z ? ?? z ??xy ??中 ? ? 10, ? ? 28, ? ? 8 / 3 。 解：rho=10;beta=29;lamda=8/3;44 f=@(t,Y)[rho*(Y(2)-Y(1));beta*Y(1)-Y(2)-Y(1)*Y(3);- lamda*Y(3)+Y(1)*Y(2)];%方程 组必须使用列向量表示 [t,y]=ode45(f,[0,30],[5,13,17]) subplot(2,2,1) plot(t,y(:,1),'*') %x 曲线 subplot(2,2,2) plot(t,y(:,2),'X') subplot(2,2,3) plot(t,y(:,3),'O') subplot(2,2,4) plot3(y(:,1),y(:,2),y(:,3)) %绘制空间轨迹图 3、边值问题的 Matlab 数值解 Matlab 中使用 bvp4c（bvp5c）和 bvpinit 命令求解常微分方程中的两点边值问题， 其调用格式为： sol = bvp4c(odefun,bcfun,solinit,options) solinit = bvpinit(x,yinit,parameters) %该函数用来给 sol 提供初始猜测解。5.6 级数求和Matlab 级数求和的命令为 r = symsum(expr, v) r = symsum(expr, v, a, b) 其中 expr 为级数的通项表达式， v 是求和变量，a 和 b 分别为求和变量的起始点和 终止点，若没有指明 a 和 b，a 的默认值为 0，b 的默认值为 v-1。 无穷数列的累加称为级数，当取其前面若干有限项时，得到的是部分和，数列 a 累 加形成的新序列可用 s=cumsum(a)实现。 a 的长度为 n， 则 s 的长度也为 n， 即每一个 s(k) 是数组 a 中前 k 项的和。注意 cumsum 与 sum 的区别，sum(a)得到的是一个数，即序列 s 中最后一项。 示例： clc,clear k=1:10; a=1./k.^2;45 s=cumsum(a) ss=sum(a) 例 22 求如下级数的和解 （1）syms n f1=(2*n-1)/2^n; s1=symsum(f1,n,1,inf) （2）syms n f2=1/n^2; s2=symsum(f2,n,1,inf) 补充一点： Matlab 求数组元素乘积的函数 B=prod(A) 如果 A 是一个 m 行一列的 （向 量），则这种用法即返回这 m 个元素的乘积； 如果 A 是一个 m 行 n 列的矩阵，则 A 的每一列都被看做一个 m 行 1 列的向量，分别计算每个向量中元素的乘积，返回给 B， 因此 B 是一个 1 行 n 列的数组。5.7 Matlab 求解优化问题Matlab 求解各种形式的极值问题被集成在优化工具箱中。 下面是优化工具箱 4.0 版本的功能示意图：46 47 变量 f描述 线性规划目标函数 f*X 或二次规划的目标 函数 X’*H*X+F*X 中线性项的系数向量 非线性优化的目标函数. fun 必须为行命令 对象或 M 文件、 嵌入函数、 或 MEX 文件的 名称 二次规划目标函数 X’*H*X+f*X 中的二次 项的系数矩阵 A 矩阵和 b 向量分别为线性不等式约束： AX≤b 中的系数矩阵和右端向量 Aeq 矩阵和 beq 向量分别为线性等式约束 Aeq*X=beq 中的系数矩阵和右端向量 X 的下限和上限向量 迭代初始点坐标 函数最小化的区间 优化选项参数结构调用函数 linprog,quadprog fminbnd,fminsearch,fminunc, fmincon,lsqcurvefit,lsqnonlin, fgoalattain,fminimax quadprog linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax,lsqcurvefit, lsqnonlin 除 fminbnd 外所有函数 fminbnd 所有优化函数funHA,bAeq,beqvlb,vub X0 x1,x2 options变量描述 描述退出条件： exitflag&0 表示目标函数收敛于 x exitflag=0 表示已达到函数评价或迭代的最大次 数 exitflag&0 表示目标函数不收敛 由优化函数求的 x 值 解 x 处的目标函数值调用函数exitflagx所有优化函数 linprog,quadprog,fgoalattain fmincon,fminimax,lsqcurvefit lsqnonlin,fminbndfevaloutput包含优化结果信息的输出结构 所有优化函数 iterations,Algorithm,FuncCount（函数评价次数）进行极值计算时有两种方式，通过命令 optimtool 打开优化工具箱 GUI 求解和使用48 上述提供的命令函数进行求解。 示例一：使用单纯形计算最小值 min f ? ? 4x1 ? x 2 ? ? x1 ? 2x 2 ? 4 ? 2 x ? 3 x ? 12 ?? 2 s.t . ? 1 ? x1 ? x 2 ? 3 x1 , x 2 ? 0 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=[-4;-1]; A=[-1 2;2 3;1 -1]; b=[4;12;3]; [x,fval,exitflag,output]=linprog(f,A,b) 示例二：非线性优化 f ( x ) ? ? x1 x 2 x 3 ? x1 ? 2x 2 ? 2x 3 ? 0 ?? s.t . ? ? x1 ? 2 x 2 ? 2x 3 ? 72 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=@(x)(-x(1)*x(2)*x(3)); A=[-1 -2 -2;1 2 2]; b=[0;72]; [x,fval,exitflag,output]=fmincon(f,[10;10;10],A,b) 其他求解极值的函数主要 有： (1)x=fminbnd(fun,x1,x2)，求函数在区间[x1,x2]上的极小值对应的自 变量值 示例三：求 f(x)=x^4-x^2+x-1 在区间 [-2,1]上的极小值。 解： Matlab 命令如下： [x,feval,exitflag,output]=fminbnd('x^4-x^2+x-1',-2,1) (2) x=fminsearch(fun,x0)无约束非线性极小值求解， x0 为起始搜索点。49 示例四：fminsearch 函数求 f(x)=1/((x1-2)^2+3)-1/(2(x2+1)^2-5),初始点[0,0]， (3) x=fminunc(fun,x0,options)求解函数的一个局部极小值 (4) x=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon 包括非线性约束 ; (5) x=lsqnonlin(fun,x0,lb,ub);非线性最小二乘优化 (6) x=linprog(f,A,b,Aeq,beq,lb,ub,x0)求解线性规划问题 (7) x=bintprog(f,A,b,Aeq,beq,x0)求解 0-1 规划 (8) x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0);求解二次规划5.8 Matlab 求解插值和拟合问题5.8.1 Matlab 插值工具箱 1、一维插值函数 interp1，命令格式为： yi=interp1(x,Y,xi,? method ?)其中 method 指定的方法为 nearest,linear,spline,cubic， pchip （分段三次 Hermite 插值） ， 所有的插值方法都要求 x 是单调的。x 与 Y 是具有相同大小 的向量，求在 xi 点处的插值函数值 yi 一维快速插值函数 interp1q，命令格式为： yi = interp1q(x,Y,xi)，其中 x 是单调递增的列向量，Y 是列向量或者行数为 length(x) 的矩阵，xi 是一个列向量。 yi = interp1(x,Y,xi,method,'extrap') %对于超出 x 范围的 xi 中的分量将执行特殊的外 插值法 extrap。 示例 1：在一 天 24 小时内，从零点开始每间隔 2 小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28 ，27，25，20 ，18，15 ，13，推测中午 12 点（即 13 点）时的温 度． 解： x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') xi=0:1/3600:24; %绘制一天内的温度变化曲线 yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi) 示例 2： x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,'ro',x,y);50 x1=0:0.02:1; %要插值点 y1=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1); y2=interp1(x,y,x1);y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');y5=interp1(x,y,x1,'cubic'); plot(x1,[y2' y3' y4' y5'],':',x,y,'ro',x1,y1); legend('linear','spline','nearest','cubic','样本点','原函数'); %计算各插值的最大误差 [max(abs(y1-y2)),max(abs(y1- y3)),max(abs(y1-y4)),max(abs(y1- y5))]; 2、二维数据内插值 interp2，插值函数为二元函数。 ZI = interp2(X,Y,Z,XI,YI,? method? ) X 和 Y 分别是 m 维和 n 维的向量，表示节点，Z 为 n*m 的矩阵，表示节点值，XI， YI 是为一维数组，表示插值点，XI 与 YI 是方向不同的向量，一个是行向量，一个是列 向量，ZI 是行数为 YI 的维数，列数为 XI 的维数的矩阵，表示得到的插值。若 XI 与 YI 中有在 X 与 Y 范围之外的点，则相应地返回 nan（Not a Number） 。 method 为插值 方法。 如果在某区域测量了若干个节点的高程（节点值） ，为了画出比较精确的等高线图， 需要先插入更多的点（插值点） ，然后计算这些点的高程。 如果是三次样条插值，可以使用命令 pp=csape({x0,y0},z0,conds,valconds);z=fnval(pp,{x,y}) 示例 1： years = 0; service = 10:10:30; wage = [150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; w = interp2(service,years,wage,15,1975) %计算在 1975 年工作了 15 年的员工的工资[x,y]=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); surf(x,y,z);figure [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1); z2=interp2(x,y,z,x1,y1,'cubic'); z3=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1);figure %比较误差 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);51 surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.08])surf(x1,y1,abs(z0-z3)); axis([-3,3,-2,2,0,0.08]) interp2 函数能够较好的进行二维插值运算，但是它只能处理以网格形式给出的数 据，如果数据是随机给出的，此时可以使用 griddata 函数。其调用格式为： 格式:z=griddata(x0,y0,z0,x,y,? method ? ) ? v4 ?:MATLAB4.0 提供的插值算法，公认效果较好； ? linear ?：双线性插值算法（缺省算法） ； ? nearest? ：最临近插值； ?cubic?：双三次插值。clear x=-3+6*rand(199,1);y=-2+4*rand(199,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); plot(x,y,'*');plot3(x,y,z,'*'); [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1); z2=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z2); %误差比较 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); surf(x1,y1,abs(z0-z1));axis([-3,3,-2,2,0,0.15]) surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15]) 3、 三维插值运算函数 interp3,n 维网格插值 interpn， 其调用格式同 interp1 和 interp2， 对应的三维网格生成函数为[x,y,z]=meshgrid(x1,y1,z1)和非网格生成函数 griddata3(),griddatan()，他们同 griddata 调用格式相同。 示例 1：解：[x,y,z]=meshgrid(-1:0.2:1); [x0,y0,z0]=meshgrid(-1:0.05:1); V=exp(x.^2.*z+y.^2.*x+z.^2.*y).*cos(x.^2.*y.*z+z.^2.*y.*x); V0=exp(x0.^2.*z0+y0.^2.*x0+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0); V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:))52 slice(x0,y0,z0,V1,[-0.5,0.3, 0.9],[0.6,-0.1],[-1,-0.5,0.5,1]) title('Slives for Four Dim Figures'); 4、三次样条插值函数 已知平面上 n 个点 xi,yi(i=1...n)，其中 x 1&x 2& … &xn ，S(x)为三次样条函数需要满足 三个条件： (1)S(xi)=yi(i=1...n)，即函数经过样本点； (2)S(x)在每个子区间 [xi,xi+1]上为三次多项式：(3)S(x)在整个区间 [x1,xn]上有连续的一阶及二阶导数。 Matlab 中有如下函数： y=interp1(x,Y,xi,?spline? ); y=spline (x,y,xi); %对于给定的离散的测量数据 x,y（称为断点） ，要寻找一个三项多 项式 y = p(x) ，以逼近每对数据(x,y)点间的曲线。 示例 1：解：x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); sp=csapi({x0,y0},z); fnplt(s p); 示例 2:对离散分布在 y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算 解：x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; y = exp(x).*sin(x); xx = 0:.25:20; yy = spline(x,y,xx); plot(x,y,'o',xx,yy) pp=csape (x,y,conds); %conds 是边界条件选项，可以取以下值： 'complete',给定边界一阶导数. 'not-a-knot',非扭结条件,不用给边界值. 'periodic',周期性边界条件,不用给边界值.53 'second',给定边界二阶导数. 'variational',自然样条 (边界二阶导数为 0) pp=csape (x,y,conds,valconds);y= ppval(pp,x)， 求插值点的函数值必须调用 ppval 函数。 示例 3：给定数据 x=[1 2 4 5]， y=[1 3 4 2]，边界条件 S”(1)=2.5，S”(5)=-3，计算插 值。 解：x=[1 2 4 5];y=[1 3 4 2]; pp=csape(x,y,'second',[2.5,-3]); pp.coefs %显示每个区间上三次多项式的系数 xi=1:0.1:5; yi=ppval(pp,xi); %计算插值点对应的函数值 plot(x,y,?o?,xi,yi) Matlab 还提供了一个生成三次参数样条类的函数 S=csapi(x,y) 其中 x=[x 1,x 2,….,xn], y=[y1,y2,…,yn]为样本点。S 返回样条函数对象的插值结果， 包括子区间点、各区间点三次多项式系数等。 可用 fnplt()绘制出插值结果，其调用格式： fnplt(S) 对给定的向量 xp，可用 fnval()函数计算，其调用格式： yp=fnval(S,xp) 其中得出的 yp 是 xp 上各点的插值结果。 示例 4：解：x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); sp=csapi(x,y); fnplt(sp) c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:),... sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:)] 处理多个自变量的网格数据三次样条插值类的函数： S=csapi({x1,x 2,…,xn},z) 其中 xi 为自变量标识， z 是网格数据的样本点，S 是三次样条曲线的对象。 5.8.2 Matlab 中的拟合方法 1、线性最小二乘法拟合曲线 当用 m 次多项式拟合给定数据时，一般使用函数 polyfit，命令格式为：54 a=polyfit(x,y,m)，输入参数 x，y 为要拟合的数据，m 为拟合多项式的次数，输出参 数 a 为拟合多项式的系数向量。 多项式在 x 处的函数值 y 可以使用函数 polyval，格式为： y=polyval(a,x)，a 为向量表示的多项式 示例：某乡镇企业
年的生产利润如下表所示，试预测 1997 年和 1998 年 的利润。 年份 利润 解： x=[92 95 1996]; 1 122 3 152 5 196 y=[70 122 144 152 174 196 202]; a=polyfit(x,y,1); %根据利润增长的情况，使用一次线性拟合函数 y97=polyval(a,1997); y98=polyval(a,1998); 2、最小二乘优化 常用的求解最小二乘优化的函数 lsqlin,lsqcurvefit,lsqnonlin,lsqnonneg， 这些函数都 位于 Matlab 的优化工具箱中，下面只举例说明它们的用法。 x= lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)求解的模型为：例 1：x=[19 25 31 38 44]'; y=[19 32.3 49 73.3 97.8]'; r=[ones(5,1),x.^2]; %拟合形如 y=a+b*x^2 的多项式，如果拟合的多项式的形式为 y=a0+a1*x+a2*x^2，则 r 应该写成 r=[ones(5,1),x,x.^2] ab=lsqlin(r,y); %计算多项式的系数 x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; plot(x,y, ? o? ,x0, y0, ? r ? ) % 绘制对比曲线查看结果55 lsqcurvefit 函数求解的数学模型为函数的调用格式为： x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)； fun 参数是定义函数 F(x,xdata)的 M 文 件， x0 为初始点， xdata 是输入数据， ydata 是观测数据， 它们都是向量或者矩阵， F(x,xdata) 的维数和 ydata 一致。 例 2：用给定的观测数据拟合函数 y ? e? k1 x1sin(k2 x2 ) ? x23 中的参数 k1 和 k2。解： （1）首先编写 M 文件，定义拟合函数 function y=fun(k,xdata) f=exp(-k(1)*xdata(:,1)).*sin(k(2)*xdata(:,2))+xdata(:,3).^2; end （2）调用函数 lsqcurvefit x0=[23.73 5.49 1.21;…]; %x0 是 25 行 3 列的输入样本 y0=[15.02;….]; %y0 是 25 行 1 列的观测值 k0=rand(2,1); %初始化待求参数的值 lb=zeros(2,1);ub=[20;2]; %给定参数的上下界，没有边界时用 []表示 k=lsqcurvefit(@fun,k0,x0,y0,lb,ub) lsqnonlin 函数用来解决的数学模型是： 命令格式为： x=lsqnonlin(fun,x0,lb,ub,options);其中 fun 是定义向量函数 f(x)的 M 文件或者函数句 柄。 例 3：用最小二乘法拟合 y ? x0,y0 已知。 解：x0=-10:0.01:10; y0=normpdf(x0,0,1); %计算标准正态分布密度函数在 x0 处的取值 F=@(xs)1/sqrt(2*pi)/xs(2)*exp(-(x0-xs(1)).^2)/xs(2)^2/2- y0; %定义函数( x ? ? )2 ? 2 ? 2?1 2? ?e中的位未知参数?, ?，假设数据值56 xs0=rand(2,1);%非线性拟合要有初始值 xs=lsqnonlin（F,xs0）; 例 4： 示例 3:体重约 70kg 的某人在短时间内喝下 2 瓶啤酒后， 隔一定时间测量他的 血液中的酒精含量(mg/100ml)，得到如下数据，使用所给数据用函数 f=atbect 进行拟合， 求出常数 a,b,c。 解： t=[0.25 0.5 0.75 1 1.5 2 2.5 3 h1=log(h); % 进行对数变换 f=inline('a(1)+a(2).*log(t)+a(3).*t','a','t'); [x,r]=lsqcurvefit(f,[1,0.5,-0.5],t,h1); plot(t,h,'*'); ezplot('exp(4.9*log(t)-0.2663*t)',[0.2 16]); 3.5 4 4.5 5:16];h=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4];lsqnonneg 函数求解的数学模型为： 命令格式为： x=lsqnonneg(C,d,options)，该函数的用法同 lsqlin 用法类似。 3、Matlab 曲线拟合的图形用户界面操作 拟合一元函数的命令是 cftool，拟合二元函数的命令 sftool，非线性拟合的命令是 nlintool。 主要拟合步骤是： (1)将 数据导入工作空间； (2)运行 cftool 命令； (3)对数据进行预 处理； (4)选择适当的模型进 行拟合； (5)生成一些相关的统计量，并进行预测操作57 示例 1：利用下表的数据拟合方程 y ? a0 e 1 1 e 2axa x2解： （1）把数据加载到 Matlab 工作空间中， a=[5 2 1.5;5 5 2.3;7 10 3.5;8 4 5.2;8 6 1.8;3 3 2.9;7 8 3.8;7 3 4.5;2 6 4.7]; x1=a(:,1); x2=a(:,2); z=a(:,3); (2)在命令窗口中输入 sftool，打开如下界面（3）输入自变量 X，Y 的数据 x1,x2，因变量 Z 的数据 z，并定义要拟合的函数58 a0*exp(a1*x)*exp(a2*y) （4）用鼠标点击 fit 进行拟合，结果如下59 6、Matlab 在线性代数中的应用6.1 向量组的线性相关性求列向量组 A 的一个最大线性无关组可用命令 rref(A)将 A 化成行最简形，其中单 位向量对应的列向量即为最大线性无关组所含向量，其它列向量的坐标即为其对应向量 用最大线性无关组线性表示的系数。 例 23 求下列矩阵列向量组的一个最大无关组。解 编写 Matlab 程序如下 format rat %数据是有理分数表示 a=[1,-2,-1,0,2;-2,4,2,6,-6;2,-1,0,2,3;3,3,3,3,4]; b=rref(a) format %恢复到短小数的显示格式解 编写 Matlab 程序如下 format rat60 a=[2,2,-1;2,-1,2;-1,2,2];b=[1,4;0,3;-4,2]; c=rref([a,b]) format %恢复到短小数的显示格式6.2 齐次线性方程组在 Matlab 中，函数 null 用来求解零空间，即满足 Ax=0 的解空间，实际上是求出解 空间的一组基（基础解系）。格式如下 : z=null(A) %z 的列向量为方程组的正交规范基，满足。 ZT Z=E z=null(A,? r? ) %z 的列向量是方程 Ax=0 的有理 基。 例 25 求方程组的通解解 编写程序如下 format rat a=[1,2,2,1;2,1,-2,-2;1,-1,-4,-3] b=null(a,'r') %求有理基 syms k1 k2 x=k1*b(:,1)+k2*b(:,2) %写出方程组的通解 format %恢复到短小数的显示格式6.3 非齐次线性方程组Matlab 中解非齐次线性方程组可以使用 “\”。虽然表面上只是一个简简单单的符号， 而它的内部却包含许许多多的自适应算法，如对超定方程（无解）用最小二乘法，对欠 定方程（多解）它将给出范数最小的一个解。61 对于非齐次线性方程组 Ax=b，Matlab 的求解命令为 x=A\b 或者(x=inv(A)*b) 无论 Ax ??b 是有唯一解，无穷多解还是无解时，Matlab 总是给出一个解。 （1）当 Ax ??b 唯一解时，Matlab 给出的就是该唯一解。 （2）当 Ax ??b 有无穷多解时，Matlab 给出的是最小范数解。 （3）当 Ax ??b 无解时，Matlab 给出的是最小二乘解。 另外求解欠定方程组 （多解）可以使用求矩阵 A 的行最简形命令 rref(A)，求出所 有的基础解系。 例 26 求超定方程组解 编写程序如下 a=[2,4;3,-5;1,2;2,1]; b=[11;3;6;7]; solution=a\b 注：上面解超定方程组的“\”可以用伪逆命令 pinv 代替，且 pinv 的使用范围比“ \”更 加广泛，pinv 也给出最小二乘解或最小范数解。 例 27 用最小二乘解法解方程组解 编写程序如下 format rat a=[1,1,0;1,0,1;1,1,1;1,2,1]; b=[1;2;0;-1]; x1=a\b %这里\和 pinv 是等价的 x2=pinv(a)*b62 format %恢复到短小数的显示格式 例 28 求解欠定方程组解 编写程序如下 format rat a=[1,-1,-1,1,0;1,-1,1,-3,1;1,-1,-2,1,-1/2]; b=rref(a) format %恢复到短小数的显示格式 例 29 利用表 1 的数据拟合函数解 拟合的 Matlab 程序如下 clc, clear d=[5 2 1.5 5 5 2.3 7 10 3.5 8 4 5.2 8 6 1.8 3 3 2.9 7 8 3.8 7 3 4.5 2 6 4.7];63 m=size(d,1); %获取 d 的行数 a=[ones(m,1),d(:,1).^3,d(:,1).*d(:,2),d(:,2).^2]; % 求 d(:,3)=xs*a 中的 xs=a\d(:,3) %求的多项式的系数 例 30 利用表 1 的数据拟合函数clc, clear d=[5 2 1.5 5 5 2.3 7 10 3.5 8 4 5.2 8 6 1.8 3 3 2.9 7 8 3.8 7 3 4.5 2 6 4.7]; m=size(d,1); a=[ones(m,1),d(:,1),d(:,2)]; b=log(d(:,3)); xs=a\b; xs(1)=exp(xs(1)) 例 31 使用用户图形界面解法命令 sftool 拟合例 29 和例 30 中的参数。6.4 相似矩阵及二次型解 二次型的矩阵为编写如下程序 A=[0,1,1,-1;1,0,-1,1;1,-1,0,1;-1,1,1,0]; [P,D]=eig(A)64 解 Matlab 程序如下 a=[2 -2 0; -2 4 0; 0 0 5]; b=eig(a); if all(b&0) fprintf('二次型正定 \n'); else fprintf('二次型非正定 \n'); end [c,d]=eig(a) %c 为正交变换的变换矩阵6.5 线性代数中的其它应用计算的 Matlab 程序如下 clc, clear t=[12 10 1; 2 6 1; 4 -10 1]; S=1/2*abs(det(t)) %不考虑顶点的排列，加上绝对值函数 abs。 解法二 利用复数运算的命令计算。 clc, clear a=complex(12,10); b=complex(2,6); c=complex(4,-10); %注意复数的写法。 plot([a,b,c,a]) %画出三角形 theta=angle((b-a)/(c-a)) %利用复数运算求两个向量夹角 S=1/2*abs(b-a)*abs(c-a)*abs(sin(theta))65 解法一 使用 Matlab 的命令 orth，Matlab 程序如下 clc, clear format a=[1 -1 4; 2 3 -1; -1 1 0]; b=orth(a)就是一组规范正交基。 编写程序实现Schimidt正交化过程，程序如下 clc, clearformat a=[1 -1 4; 2 3 -1; -1 1 0]; [m,n]=size(a); %求矩阵 a 的维数 b(:,1)=a(:,1); c(:,1)=b(:,1)/norm(b(:,1)); for i=2:n s=a(:,i); for k=1:i-1 s=s-dot(a(:,i),c(:,k))*c(:,k); %正交化过程 end c(:,i)=s/norm(s); %单位化 end c66 7、Matlab 数据处理7.1 Matlab 中的默认数据文件 mat 文件1、将工作区变量保存到文本文件中 save filename var1 var2 var3… -mat 默认保存的文件扩展名为.mat，如果保存的 数据需要跨平台处理，需要采用如下格式： save filename.dat var1 var2 … -ascii 另外： save filename var1 var2 … -append 附加到已经存在的 Mat 文件中 例 36 把 Matlab 工作空间中的数据矩阵 a，b，c 保存到数据文件 data1.mat 中。 save data1 a b c 注：Matlab 中的默认数据文件 mat 文件可以省略后缀名。 2、将文本文件中的数据导入工作区 load filename 例 37 把例 36 生成的 data1.mat 中的所有数据加载到 Matlab 工作空间中。 load data17.2 纯文本文件可以把 word 文档中整行整列的数据粘贴到纯文本文件，然后调入到 Matlab 工作空 间中。 例 38 把纯文本文件 data2.txt 加载到工作空间。 a=load('data2.txt'); 或者是 a=textread('data2.txt'); 例 39 使用 dlmwrite 命令把矩阵 b 保存到纯文本文件 data3.txt。 dlmwrite('data3.txt',b) 例 40 生成服从标准正态分布随机数的矩阵，然后用 fprintf 命令保存到纯文本文件 data4.txt。67 解 clc, clear fid=fopen('data4.txt','w'); a=normrnd(0,1,100,200); fprintf(fid,'%f\n',a'); fclose(fid); 注：对于高维矩阵，用 dlmwrite 构造的纯文本文件，Lingo 软件不识别；为了 Lingo 软件识别，纯文本文件必须用 fprintf 构造，而且数据之间的分割符为“\ n”。 示例： M = gallery('integerdata', 100, [5 8], 0); dlmwrite('myfile.txt', M, 'delimiter', '\t') dlmread('myfile.txt') dlmread('myfile.txt', '\t', 2, 3) dlmread('myfile.txt', '\t', 'C1..G4')7.3 Excel 文件xlswrite(filename,A,sheet,range) 将阵列 A 写入 Excel 文件 filename 中 sheet 表格的 range 指示的单元格内。 例 41 把一个矩阵写到 Excel 文件 data5.xls 表单 Sheet2 中 B2 开始的域中。 105 a=rand(5,10); xlswrite ('data5.xls',a,'Sheet2','B2') %默认的表格为 Sheet1 [num,txt,raw] = xlsread(filename,sheet,range)从 Excel 文件 filename 中的 sheet 表格 中读取 range 单元格内的数据返回。 例 42 把例 41 生成的 Excel 文件 data5.xls 中表单 Sheet2 的域“C3：F6”中的数据赋 给 b。 b=xlsread('data5.xls','Sheet2','C3:F6') 示例： values = {1, 2, 3 ; 4, 5, 'x' ; 7, 8, 9}; headers = {'First', 'Second', 'Third'}; xlswrite('myExample.xls', [ values]); A = xlsread('myExample.xls') subsetA = xlsread('myExample.xls', 1, 'B2:C3') columnB = xlsread('myExample.xls', 'B:B') [ndata, text, alldata] = xlsread('myExample.xls')7.4 字符串数据例 43 统计下列五行字符串中字符 a、c、g、 t 出现的频数。68 1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggagg 2.cggaggacaaacgggatggcggtattggaggtggcggactgttcgggga 3.gggacggatacggattctggccacggacggaaaggaggacacggcggacataca 4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagctta 5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggc 解 把上述五行复制到一个纯文本数据文件 shuju.txt 中，编写如下程序 clc fid=fopen('shuju.txt','r'); %以只读的形式打开文件 i=1; while (~feof(fid)) %判断是否到文件末尾 data=fgetl(fid); a=length(find(data==97)); b=length(find(data==99)); c=length(find(data==103)); d=length(find(data==116)); e=length(find(data&=97&data&=122)); f(i,:)=[a b c d e a+b+c+d]; i=i+1; end f, he=sum(f) dlmwrite ('pinshu.txt',f); dlmwrite('pinshu.txt',he,'-append'); fclose(fid); %打开的文件使用完之后一定要关闭。 Matlab 中关于文件操作的常用函数： 1) textread 函数可以读取格式化为列数据的 ASCII 码文件， 每一列可以是不同的数 据类型，调用形式为： [a,b,c,...]=textread(filename,format,n) 示例 1：读文件 test_input.dat [first,last,blood,gpa,age,answer]=textread(?t est_input.dat ? ,? %s %s %s %f %d %s?) 通过使用%*s 可以忽略对应的列数据。 示例 2：假设文件 test.dat 中的数据如下： Jam}
产生一个 3x5 随机整数矩阵,每个元素的范围在-5 到5 之间;给出每个元素的符号(+1 -1)表示，用MATLAB做
A=round(rand(3,5)*10-5） ;//产生随机矩阵B=sign(A)//产生符号矩阵
采纳率：71%
来自团队：
为您推荐：
其他类似问题
换一换
回答问题，赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。怎样用MATLAB生成一个4行5列的间的随机矩阵_百度知道
怎样用MATLAB生成一个4行5列的间的随机矩阵
我有更好的答案
rand(4,5)&%&四行五列(0,1)之间的随机矩阵%&如果范围不在0-1之间，可以相应调整ceil(rand(4,5)*10)&%&[1,10]之间随机整数
采纳率：91%
为您推荐：
其他类似问题
随机矩阵的相关知识
换一换
回答问题，赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。当前位置： >>
Matlab 2014软件教程(完美版)
Matlab 软件目录1、Matlab 帮助的使用 ........................................................................ 3 2、Matlab 数据输入与类型 ................................................................ 4 3、Matlab 中的 M 文件及程序调试 ................................................ 12 4、Matlab 绘图命令 .......................................................................... 17 5、Matlab 在高等数学中的应用 ...................................................... 34 6、Matlab 在线性代数中的应用 ...................................................... 60 7、Matlab 数据处理 .......................................................................... 67 9、评价方法 ...................................................................................... 82 10、预测方法 .................................................................................... 97 11、蒙特卡洛方法 .......................................................................... 110 12、智能算法 .................................................................................. 122 13、分形 .......................................................................................... 129 14、Simulink 初步 .......................................................................... 134 15、Matlab 在概率统计中的应用 .................................................. 147 参考文献 .......................................................................................... 1511 作为和 Mathematica 、 Maple 并列的三大数学软件。 其强项就是其强大的矩阵计算以 及仿真能力。要知道 Matlab 的由来就是 Matrix + Laboratory = Matlab，所以这个软件在 国内也被称作《矩阵实验室》 。每次 MathWorks 发布 Matlab 的同时也会发布仿真工具 Simulink。在欧美很多大公司在将产品投入实际使用之前都会进行仿真试验，他们所主 要使用的仿真软件就是 Simulink。 Matlab 提供了自己的编译器： 全面兼容 C++以及 Fortran 两大语言。所以 Matlab 是工程师，科研工作者手上最好的语言，最好的工具和环境。 Matlab 已经成为广大科研人员的最值得信赖的助手和朋友！ 目前 MATLAB 产品族可以用来进行： － 数值分析 － 数值和符号计算 － 工程与科学绘图 － 控制系统的设计与方针 － 数字图像处理 － 数字信号处理 － 通讯系统设计与仿真 － 财务与金融工程... Simulink 是基于 MATLAB 的框图设计环境，可以用来对各种动态系统进行建模、 分析和仿真，它的建模范围广泛，可以针对任何能够用数学来描述的系统进行建模，例 如航空航天动力学系统、卫星控制制导系统、通讯系统、船舶及汽车等等，其中了包括 连续、离散，条件执行，事件驱动，单速率、多速率和混杂系统等等。 Simulink 提供 了利用鼠标拖放的方法建立系统框图模型的图形界面，而且 Simulink 还提供了丰富的 功能块以及不同的专业模块集合，利用 Simulink 几乎可以做到不书写一行代码完成整 个动态系统的建模工作。2 1、Matlab 帮助的使用1.1 helphelp ?????????????%帮助总览 help elfun ???%关于基本函数的帮助信息 help exp ????? %指数函数 exp 的详细信息1.2 lookfor 指令当要查找具有某种功能但又不知道准确名字的指令时， help 的能力就不够了， lookfor 可以根据用户提供的完整或不完整的关键词，去搜索出一组与之相关的指令。 lookfor integral ?? %查找有关积分的指令 lookfor fourier ??%查找能进行傅里叶变换的指令1.3 超文本格式的帮助文件在 Matlab 中，关于一个函数的帮助信息可以用 doc 命令以超文本的方式给出，如?doc ??? doc doc ??? doc eig ????%eig 求矩阵的特征值和特征向量1.4 pdf 帮助文件可从 MathWorks 网站上下载有关的 pdf 帮助文件。 网站地址：http://www.mathworks.com/3 2、Matlab 数据输入与类型2.1 Matlab 中的变量MATLAB 程序中的基本数据单元称为阵列(Array)，是一个分为行与列的数据集合。 变量被看作是只有一行一列的阵列。MATLAB 语言不需要对变量进行事先声明，也不 需要指定变量类型，它会自动根据所赋予变量的值或对变量所进行的操作来确定变量的 类型。其命名规则为： （1）变量名的大小写是敏感的。 （2）变量的第一个字符必须为英文字母，而且不能超过 31 个字符。 （3）变量名可以包含下划线、数字，但不能为空格符、标点。 （4）命名变量时可以取一个容易记忆并且能表达出其含义的名称，如汇率，可以 定义为 exchange_rate。 对于变量作用域，默认情况是局部变量，使用 global 定义全局变量，而且全局变 量常用大写的英文字母表示。 MATLAB 预定义的变量如下表所示： ans eps pi inf NaN i 或 j nargin nargout realmax realmin flops 预设的计算结果的变量名 MATLAB 定义的正的极小值=2.2204e-16 内建的π 值 ∞值，无限大 无法定义一个数目 虚数单位 i=j=√-1 函数输入参数个数 函数输出参数个数 最大的正实数 2 最小的正实数 2 浮点运算次数1023-10224 注 1：在定义变量时要尽量与避免与这些名字相同，以免改变它们的值，如果已经 改变，可以通过 clear 复这些值。 注 2：数字的输入输出格式。所有数据均按 IEEE 浮点标准的长型格式存储。输入 格式沿用了 C 语言的风格和规则；输出格式使用 format 数据格式 命令控制，只影响 在屏幕上的显示结果，不影响内部的存储和运算。 变量名 来恢复它的初始值，也可以通过重新启动 MATLAB 恢2.2 向量及其运算1、向量的生成 ①命令窗口直接输入，使用[ ]，元素之间用空格、逗号或者分号隔开。 ②使用冒号表达式，基本形式为 x=x0:step:xn，其中 xn 为尾元素数值限，而不一定 是尾元素的值。当 step=1 时可省略步长。 ③生成线性等分向量，使用 linspace 函数。Y=linspace(x1,x2,n) ④生成对数等分向量，使用 logspace 函数。Y=logspace(x1,x2,n) 2、向量的基本运算5 数加（减） 、数乘、点积(dot 函数)、叉积 (cross 函数)、混合积 dot(a,cross(b,c))2.3 矩阵及其运算1、简单矩阵的输入 （1）要直接输入矩阵时，矩阵一行中的元素用空格或逗号分隔；矩阵行与行之间 用分号“;”隔离，整个矩阵放在方括号“[ ]”里。 A=[1,2,3;4,5,6;7,8,9] ??? 说明：指令执行后，矩阵 A 被保存在 Matlab 的工作空间中，以备后用。如果用户 不用 clear 指令清除它，或对它进行重新赋值，那么该矩阵会一直保存在工作空间中， 直到本次指令窗关闭为止。 （2）矩阵的分行输入，此时回车键作为分行标志， A=[1,2,3 4,5,6 7,8,9] （3）使用 M 文件创建大矩阵，当矩阵维数非常大时，可以创建 m 文件，在 m 文 件中输入数据或者导出数据文件。 2、矩阵的基本运算 ①矩阵的四则运算。其中乘法运算要注意相乘的双方有相邻公共维，除法分为左除 “\” (A\B=inv(A)*B)和右除”/”（ A/B=A*inv (B)）(需要计算逆矩阵 ) ②矩阵的逆运算。inv 函数。 ③矩阵的幂运算。^ 。 ④矩阵的指数运算。exp(返回每个元素的指数值)，expm([V,D] = EIG(X) and EXPM(X) = V*diag(exp(diag(D)))/V)，expm1 (exp(x)-1) ⑤矩阵的对数运算。logm ⑥矩阵的特征值函数。eig 和 eigs(适合于大型稀疏方阵) ⑦矩阵的奇异值函数。svd（[U,S,V] = SVD(X)， X = U*S*V'）和 svds ⑧矩阵的条件数函数。cond （矩阵 A 的条件数等于 A 的范数与 A 的逆的范数的 乘积，c = cond(A,p)等价于 norm(A,p) * norm(inv(A),p)） ，condest(1 范数的条件数的估计 值)， rcond ⑨特征值的条件数函数。codeig（[V,D,s] = condeig(A) 等价于 [V,D] = eig(A); s6 =condeig(A);） ⑩范数函数。norm （1- 范数：即列范数，矩阵的各列绝对值之和的最大值；2- 范数： 所有元素的平方和开根号(默认)；无穷范数：即行范数，矩阵各行的绝对值之和的最大 值） ， normest (矩阵的 2 范数的估计值) 其他还有秩函数 rank， 迹函数 trace ， 零空间函数 null （又称为核空间， X=null(A)， 则 A*X=0 ，X? *X=I ） ，正 交 空 间 函 数 orth(B = orth(A),则 B'*B = eye(rank(A)))，伪逆函数 pinv 等。 3、特殊向量和特殊矩阵 （1）特殊向量 t=[0:0.1:10] %产生从 0 到 10 的行向量，元素之间间隔为 0.1 t=linspace(n1,n2,n) %产生 n1 和 n2 之间线性均匀分布的 n 个数 (缺省 n 时,产生 100 个数) t=logspace(n1,n2,n) ( 缺省 n 时,产生 50 个数) %在和之间按照对数距离等间距产生 n 个数。 （2）特殊矩阵 i)单位矩阵 eye(m), eye(m,n) 可得到一个可允许的最大单位矩阵而其余处补 0, eye(size(a)) 可以得到与矩阵 a 同样大小的单位矩阵。 ii) 所有元素为 1 的矩阵 ones(n)，ones(size(a))，ones(m， n)。 iii）所有元素为 0 的矩阵 zeros(n)， zeros(m,n)。 iv）空矩阵是 一个特殊矩阵，这在线性代数中是不存在的。例如 q=[ ] 矩阵 q 在工作空间之中，但它的大小为零。通过空矩阵的办法可以删除矩阵的行与 列。例如 a(:,3)=[] 表示删除矩阵 a 的第 3 列。 v）随机数矩阵7 rand(m,n) 产生 m× n 矩阵，其中的元素是服从[0,1]上均匀分布的随机数。 randint(m,n,[min,max]) 产生 m× n 矩阵，其中的元素是[min,max]上的随机整数。 normrnd(mu,sigma,m,n)产生 m× n 矩阵，其中的元素是服从均值为 mu，标准差 为 sigma 的正态分布的随机数。 exprnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的指数分布的随 机 数。 poissrnd(mu,m,n) 产生 m× n 矩阵，其中的元素是服从均值为 mu 的泊松 （Poisson） 分布的随机数。 unifrnd(a,b,m,n) 产生 m× n 矩阵， 其中的元素是服从区间[a,b]上均匀分布的随机数。 r = mvnrnd(MU,SIGMA,cases) 产生 cases 对均值向量为 MU，协方差阵为 SIGMA 的多维正态分布的随机数。 vi）随机置换 randperm(n)产生 1 到 n 的一个随机全排列。 perms([1:n])产生 1 到 n 的所有全排列。 vii) 稀疏矩阵 稀疏矩阵是指矩阵中零元素很多，非零元素很少的矩阵。对于稀疏矩阵，只要存放 非零元素的行标、列标、非零元素的值即可，可以按如下方式存储 （非零元素的行地址，非零元素的列地址），非零元素的值。 在 Matlab 中无向图 和有向图邻接矩阵的使用上有很大差异。 对于有向图，只要写出邻接矩阵，直接 使用 Matlab 的命令 s parse 命令，就可以把 邻接矩阵转化为稀疏矩阵的表示方式。 对于无向图，由于邻接矩阵是对称阵，Matlab 中只需使用邻接矩阵的下三角元素， 即 Matlab 只存储邻接矩阵下三角元素中的非零元素。 稀疏矩阵只是一种存储格式。Matlab 中，普通矩阵使用 sparse 命令变成稀疏矩阵， 稀疏矩阵使用 full 命令变成普通矩阵。 例 1 a=zeros(5); a(1,[2,4])=[3,4]; a(3,[2:4])=[1 3 8]; a(5,[1,5])=[6,7]8 b=sparse(a) %普通矩阵转化成稀疏矩阵 c=full(b) %稀疏矩阵转化成普通矩阵 其他一些特殊矩阵如下表所示： 函数 compan gallery hadamard hankel hilb invhilb 功能 伴随阵 Higham 测试阵 Hardamard 矩阵 Hankel 矩阵 Hilbert 矩阵 反 Hilbert 矩阵 函数 magic rosser toeplitz pascal vander wilkinson 功能 魔方阵 经典对称特征值测试阵 Toeplitz 矩阵 Pascal 矩阵 范德蒙矩阵 Wilkinson’s 特征值测试矩阵4、矩阵的特殊操作 ①变维。 有两种方法，使用冒号（:）和使用函数 reshape 使用“:”表达式对两个矩阵进行变维操作，需要预先定义两个矩阵的维数（例如 A(:,3)=[]） ； reshape 有两种形式，分别为 reshape(X,M,N)和 reshape(X,M,N,P…) ②变向 主要函数如下表所示： 函数 fiplr fipud fipdim Rot90 ③矩阵的抽取 对角线元素抽取函数 diag(X,k)/diag(v,k)，抽取矩阵 X 的第 k 条对角线的元素向量/ 使得向量 v 为所得矩阵的第 k 条对角线元素。 上三角元 素抽取 tril(X,k)和下三角元素抽取 triu(X,k) ④扩展 两种方法：利用对矩阵标示块的赋值命令 X(m1:m2,n1:n2)=a 生成大矩阵， 其中 m2- m1+1 必须等于 a 的行维数， n2- n1+1 必须等于 a 的列维数，生成 m2 × n2 维的 功能 矩阵左右翻转 矩阵上下翻转 矩阵特定维翻转 矩阵反时针 90 翻转 函数 diag tril triu 功能 产生或提取对角阵 产生下三角 产生上三角9 矩阵 X；利用小矩阵组合生成大矩阵，要严格注意矩阵大小的匹配。 （示例： X(1:10,1:10)=zeros(10,10)，LX=[X,X;X,X]） Matlab 中冒号（： ）的使用方法小结： （1）用于生成向量，a:b，一般要求 a&b，否则生成空矩阵。 （2）如果 b-a 不是整数时，则向量的最后一位数为 n+a，其中 n=fix(b-a)(向零取整) （3）a:c:b 表示[a,a+c,…,a+n*c]，其中 n=fix((b-a)/c)，当 c&0 且 a&b 或者 c&0 且 a&b 时出现空矩阵。 （4）A(:)以一列的方式显示 A 中所有的元素 （5）b=A(i,:)表示将 A 中第 i 行存入 b 中 （6）b=A(:,j)表示将 A 中第 j 列存入 b 中 （7）b=A(j:k)表示将 A 中第 j 到第 k 个元素存入 b 中（ Matlab 中矩阵按列存储） （8）b=A(:,c:d)表示将 A 中第 c 列到第 d 列存入 b 中，要求 c，d 不能超过 A 的列 数。 （9）当矩阵很大时，不知道矩阵的维数，可以使用 end 作为矩阵的最后一行或者 一列或者最后一个元素。例如 b=A(1:2,2:end)获取矩阵 A 右上角的元素。2.4 数组及其运算数组与矩阵在形式上完全一致，只是运算与矩阵不同。同型矩阵之间的运算通常称 为数组运算。 （矩阵的数组运算） 1、基本数组运算 ①四则运算。数组的乘除法是指两个同维数组间对应元素之间的乘除法，运算符 为”.*” ,”./”和” .\ ”。数组与常数之间的运算可以加”.”，也可以不加。 ②幂运算。.^,对每个数组元素的幂运算。 ③指数运算 exp,对数运算 log 和开方运算 sqrt。 2、数组函数运算 只要把运算的数组带入到函数中就可以了，通用 形式为 funname(A) 3、数组的逻辑运算和关系运算 指令 & 含义 小于 函数名 lt10 &= & &= == ~= & | ~小于等于 大于 大于等于 等于 不等于 逻辑 与 逻辑 或 逻辑 非le gt ge eq ne and or not指令 xor any all isnan含义 不相同就取 1，否则取 0 只要有非 0 就取 1，否则取 0 全为 1 取 1，否则为 0 为数 NaN 取 1，否则为 0指令 isequal ismember isempty isletter含义 相等取 1，否则取 0 两个矩阵是属于关系取 1， 否则取 0 矩阵为空取 1，否则取 0 是字母取 1， 否则取 0 （可以是字符 串）isinf isfinite ischar find isnumeric为数 inf 取 1，否则为 0 有限大小元素取 1，否则为 0 是字符串取 1，否则为 0 寻找非零元素坐标 判断数值矩阵isstudent isprime isreal isspace islogical学生版取 1 质数取 1，否则取 0 实数取 1，否则取 0 空格位置取 1，否则取 0 判断逻辑数组11 3、Matlab 中的 M 文件及程序调试M 文件分为两种，一种是脚本文件，由一系列 Matlab 的命令组成，可以直接运行； 一种是函数文件，必须由其他 M 文件或者在命令行窗口中调用执行。函数文件具有一 定的通用性，并且可以进行递归调用。3.1、Matlab 中流程控制语句（1）if 语句 有三种形式：if(表达式 ) 语句组 A; end if(表达式) 语句组 A，else 语句组 B，end if(表达式 1) 语句组 A，elseif(表达式 2) 语句组 B, else 语句组 C，end （2）循环语句 while（表达式）语句组 A , end for k=初值：增量：终值 语句组 A，end 循环语句可以结合 break 命令来控制程序的执行顺序。 （3）switch 语句 switch 表达式(标量或者字符串 ) case 值 1 语句组 A case 值 2 语句组 B … otherwise 语句组 N end 用法示例： 示例一：输入数 n，判断其奇偶性 解 Matlab 程序如下： n=input(? n=?); if rem(n,2)==0 %求余数，与 mod 函数的区别，当同符号时两者相同，否则不同 A=?even? ; else A=?odd?,12 end 示例二：求 Matlab 中的最大实数 解：Matlab 程序如下： x=1; while x~=inf, x1=x; x=2*x; %为了获取更接近 realmax 的值，可以改为 x=1.01*x end x1 示例三：求 Matlab 的相对精度 解：matlab 程序如下： y=1; while 1+y&1 y1=y; y=y/2; end y13.2 脚本文件主要特征如下： (1)一般由 clc,clear 命令开始，清除掉屏幕和工作空间中原有的变量和图形，以避免 其他已执行的程序残留数据对本程序的影响。程序中应该添加有注释。 (2)接下来是程序的主体，如果文件中有全局变量，则需要使用 Global 在程序的起 始部分注明。为了提高程序的可读性，注意语句的缩进。 (3)整个程序应按照 Matlab 标识符的要求起文件名，扩展名为 m。 示例：列出求素数的程序 解 Matlab 程序如下： clear,clc N=input('N='),x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end x13 3.3 函数文件函数文件与脚本文件的区别： (1)由 function 开头，后跟的函数名必须与文件名相同； (2)有输入变量和输出变量，可进行变量传递； (3)除非使用 Global 声明，程序中的变量均为局部变量，运行后不保存在工作空间 中。 示例 1：使用函数文件写出求素数的程序 解：Matlab 程序如下： function y=qiuprime(N) %求出 1 到 N 之间的素数，返回素数组成的向量 x=2:N; %列出从 2 到 N 的全部自然数 for u=2:sqrt(N) %依次列出除数 n=find(rem(x,u)==0&x~=u);%找到能被 u 整除而不等于 u 的数序号 x(n)=[]; %剔除该数 end y=x; 示例 2：写一个递归函数求 n! 解： function f=factor(n) if n&=1 f=1; else f=factor(n-1)*n; end end3.4 函数参数的可调性Matlab 中有两个永久变量 nargin 和 nargout 分别记录调用该函数时的输入实参和输 出实参的个数。只要在函数中包含这两个变量，就可以准确的知道该函数文件被调用时 的输入输出参数个数，从而决定函数如何进行处理。其他类似的变量还有 varargout， varargin。14 示例 1： function fout=examp(a,b,c) if nargin==1 fout=a; elseif nargin==2 fout=a+b; elseif nargin==3 fout=a+b+c; end 示例 2： function [x,y,z]=examp(a,b,c) if nargout==1 x=a; elseif nargout==2 x=a; y=b; elseif nargout==3 x=a; y=b; z=c; end end 示例 3： function vartest(argA, argB, varargin) optargin = size(varargin,2); %可选输入 stdargin = nargin - %标准输入 fprintf('Number of inputs = %d\n', nargin) fprintf(' Inputs from individual arguments(%d):\n', ... stdargin) if stdargin &= 1 fprintf(' %d\n', argA) end if stdargin == 2 fprintf(' %d\n', argB) end fprintf(' Inputs packaged in varargin(%d):\n', optargin) for k= 1 : size(varargin,2) fprintf(' %d\n', varargin{k}) end end15 3.5 程序调试M 文件编辑器中提供了强大的程序调试功能，使用方式如同 VC 程序中的调试器。 在 M 文件编辑器中的工具栏中有如下程序调试图标：这些图标的使用和菜单 Debug 中的一些命令相同。16 4、Matlab 绘图命令4.1 二维绘图命令二维绘图的基本命令有 plot，loglog，semilogx，semilogy 和 polar。它们的使用方法 基本相同，其不同特点是在不同的坐标中绘制图形。plot 命令使用线性坐标空间绘制图 形；loglog 命令在两个对数坐标空间中绘制图形；而 semilogx(或 semilogy)命令使用 x 轴（或 y 轴）为对数刻度，另外一个轴为线性刻度的坐标空间绘制图形； polar 使用极 坐标空间绘制图形。 二维绘图命令 plot 为了适应各种绘图需要，提供了用于控制线色、数据点和线型的 3 组基本参数。它的使用格式如下： plot(x,y,‘color_point_linestyle ’) 该命令是绘制 y 对应 x 的轨迹的命令。 y 与 x 均为向量，且具有相同的元素个数。 用字符串，color_point_linestyle? 完成对上面 3 个参数的设置。 线色：r-red，g-green，b-blue ，w-white ，k-black，i-invisible ，y-yellow，c-cyan （青色） ， m-紫色。 数据点：.（圆点） ，o（小圆圈） ，x（叉号） ，+（加号） ，*（星号） ，s （square 方形） ， h（ hexagram 六角星） ，d（diamond 菱型） ，p（pentagram 五角星） ， v（下三角） ，^ （上 三角） ，&（左三角） ，&（右三角） 。 线型：-（实线） ，--（虚线） ，-.（点画线） ， :（点线） 。 当 plot(x,y)中的 x 和 y 均为 m× n 矩阵时，plot 命令将绘制 n 条曲线。 plot(t,[x1,x2,x3])在同一坐标轴内同时绘制三条曲线。 如果多重曲线 对应不同的 x 轴向量绘制，可使用命令 plot(t1,x1,t2,x2,t3,x3) 式中 x1 对应 t1，x2 对应 t2 等等。在这种情况下，t1，t2 和 t3 可以具有不同的元素 个数，但要求 x1， x2 和 x3 必须分别与 t1，t2 和 t3 具有相同的元素个数。 subplot 命令使得在一个屏幕上可以分开显示 n 个不同坐标系，且可分别在每一个 坐标系中绘制曲线。其命令格式如下 subplot(r,c,p)17 该命令将屏幕分成 r× c 个子窗口，而 p 表示激活第 p 个子窗口。窗口的排号是从 左 到右，自上而下。 下面对几个特殊的坐标系统 进行简要介绍： ①对数坐标曲线，主要有 semilogx,semilogy 和 loglog，前两个分别以 x 坐标和 y 坐标为对数坐标，后一个是双对数坐标。 例如： x=1:0.1*pi:2* y=sin(x); semilogx(x,y,'-*') x = 0:.1:10; semilogy(x,10.^x) x = logspace(-1,2); loglog(x,exp(x),'-s') grid on ②双纵坐标（双 y 轴坐标系）函数 plotyy，调用形式为： plotyy(X1,Y1,X2,Y2) plotyy(X1,Y1,X2,Y2,fun) ③极坐标绘图函数 polar 调用形式为：polar(theta,rho) 例 1：绘制极坐标下的平面曲线 ? ? a cos(b ? n? ) ，并讨论参数 a,b,n 对曲线的影响。 解 Matlab 程序如下： N=100; theta=linspace(0,2*pi,N); for i=1:2 a(i)=input('a='); b(i)=input('b='); n(i)=input('n='); rho(i,:)=a(i)*cos(b(i)+n(i)*theta); subplot(1,2,i),polar(theta,rho(i,:)); end 几个比较漂亮的极坐标系下的图 形： (1)蝴蝶图案 fun 可以是 plot、 semilogx、 semilogy 或 loglog plotyy(X1,Y1,X2,Y2,fun1,fun2) fun1 绘制(X1,Y1)， fun2 绘制(X2,Y2)18 t=0:0.01:36; f=exp(cos(t-pi/2))-2*cos(4*(t-pi/2))+sin((t-pi/2)/12).^5; polar(t,f,'r') (2)枫叶图案 t=-pi/2:0.05:1.5* f=100./(100+(t-pi/2).^8)*2-sin(7*t)-cos(30*t)/2; polar(t,f,'r') 在图形绘制完毕后， 执行如 下命令可以再在图中加入标题、 标号、 说明和分格线等。 这些命令有 title ，xlabel，ylabel，text，gtext，legend 等。它们的命令格式如下 title(“ My Title”)， xlabel(“My X-axis Label”) ，ylabel(“My Y-axis Label”)， xlabel({'first line';'second line'}) ， text(x, y,'Text for annotation')， gtext('Text for annotation')， grid on gtext 命令是使用鼠标器定位的文字注释命令。当输入命令后，可以在屏幕上得到 一个光标， 然后使用鼠标器控制它的位置。 按鼠标器的左键， 即可确定文字设定的位置。 hold on 是图形保持命令，可以把当前图形保持在屏幕上不变，同时在这个坐标系 内绘制另外一个图形。 hold on 命令是一个交替转换命令， 即执行一次， 转变一个状态 （相 当于 hold on、 hold off） 。 还可以设置坐标轴的范围，使用命令 axis ，其格式为： axis([xmin xmax ymin ymax]) axis([xmin xmax ymin ymax zmin zmax cmin cmax]) 示例： x = 0:.01:pi/2; plot(x,tan(x),'-ro') 图形效果如下所示：默认情况下，横纵坐标轴的范围是根据函数自变量和因变量的值自动变化的，有时 效果不好，此时需要设定横纵坐标轴的范围：axis ([0 pi/2 0 5])，效果如上右图所示。19 解 画图的 Matlab 程序如下，画出的图见图 1。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x); plot(x,y1,'.-'); hold on %图形保持命令 plot(x,y2,'*-'); plot(x,y3,'-o'); h=legend('sin($x$)','sin($x+\frac{\pi}{3}$)','cos($x$)') %latex 格式显示 set(h,'Interpreter','latex') %设置 Interpreter 的属性值为 latex ，可以使用数学公式 xlabel('$x$','Interpreter','latex') %latex 格式显示 ylabel('$y$','Interpreter','latex') %latex 格式显示图1解 画图的 Matlab 程序如下，画出的图形结果见图 2。 clc, clear x=-2*pi:0.1:2* y1=sin(x); y2=sin(x+pi/3)+2; y3=cos(x);20 subplot(3,1,1), plot(x,y1,'.-'), title('sin($x$)','Interpreter','latex') subplot(312), plot(x,y2,'*-') title('sin($x+\frac{\pi}{3}$)','Interpreter','latex') ylabel('$y$','Interpreter','latex') %latex 格式显示 subplot(313), plot(x,y3,'-o') title('cos($x$)','Interpreter','latex') %latex 格式显示 xlabel('$x$','Interpreter','latex') %latex 格式显示4.2 复数绘图plot(z) ： z 为复数时相当于 plot(real(z),imag(z))；如果是双变量如 plot(t,z)，则 z 中 的虚数部分将丢弃。在复平面中绘出多条曲线，必须使用 hold on 命令，或者把多条曲 线的实部和虚部明确的写出。即 plot(real(z1),imag(z1),real(z2),imag(z2)) Matlab 中专门用来绘制复变量函数图形的相关命令是 cplxmap， cplxgrid， cplxroot， 格式如下： z = cplxgrid(m) %产生(m+1)*(2m+1)的极坐标下的复数数据网格 cplxmap(z,f(z),(optional bound)) %绘制复变函数的图形，以 xy 平面为自变量所在的 复平面，以 z 轴表示复变函数的实部，颜色表示复变函数的虚部。 cplxroot：画复数的 n 次函数曲面21 cplxroot(n) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面 cplxroot(n,m) %画复数 n 次根的函数曲面,复数为最大半径为 1 的圆面,为（ m+1）* （2m+1）的方阵 例 1：绘制 z=exp((-0.1+i)*t)的复数图形 解 Matlab 程序如下： t=0:0.1:15; z=exp((-0.1+i)*t); subplot(2,2,1) plot(z), % pause 可以暂停程序的执行， 通过按任意键可使程序继续进行， 也 可以在 pause(n)中设置时间使执行结果出现动态效果。 title('复数绘图 plot(z)'); subplot(2,2,2) plot(t,z),pause title('plot(t,z)') subplot(2,2,3) polar(angle(z),abs(z)); title('polar(angle(z),abs(z))') subplot(2,2,4) semilogx(t,z); title('semilogx(t,z)') 注：Matlab 中使用函数 C=complex(A,B)构造复数。 例 2：绘图 y=cos(x+i)的图形 解法一：x=-pi:0.1: fun=@(x)cos(x+i); %如果直接使用 plot(x,y) ，会忽略虚部。 plot(fun(x)) 解法二：plot(cos(x+i)); 例 3：绘制 f=z^4 的图形，其中 z 为复数 解：z=cplxgrid(30); cplxmap(z,z.^4)4.3 显函数，符号函数或隐函数的绘图fplot(fun,lims)绘制由字符串 fun 指定函数名的函数在 x 轴区间为 lims=[xmin, xmax] 的函数图。若 lims=[xmin,xmax,ymin,ymax]，则 y 轴也被限制。22 解 （1）首先用 M 文件 fun1.m 定义函数 f(x)如下 function y=Afun1(x); if x&1 y=x+1; else y=1+1./x; end 在 matlab 命令窗口输入 fplot('Afun1',[-3,3]) 就可画出 函数 f (x ) 的图形。 （2）可以使用匿名函数，编写程序如下 fun2=@(x) (x+1)*(x&1)+(1+1/x)*(x&=1); fplot(fun2,[-3,3]) ezplot(f)绘制符号函数或者隐函数 f(x)的图形，x 轴的近似范围为[2? ,?2??。 ezplot(f,[xmin,xmax])使用输入参数来代替默认横坐标范围 [2? ,?2??。 ezplot 函数的其 他格式有： ? ezplot(fun2)绘制 fun2(x,y)=0 的隐函数曲线，默认 x,y 的范围是[2? ,?2??? ezplot(fun2,[xmin,xmax,ymin,ymax]) 示例： function z = myfun(x,y,k) %建立 M 文件 z = x.^k - y.^k - 1; %在命令窗口中输入：ezplot(@(x,y)myfun(x,y,2));解 ezplot('cot(x)')解 ezplot('x^2+y^2/4=1') 例 6-1：同一坐标系下绘制曲线 x^2+y^2=1 和 x^2-y^2=1 的曲线 解：ezplot('x^2-y^2=1');23 ezplot('x^2+y^2=1'); colormap([0,0,1]);4.4 三维图形在实际工程计算中，最常用的三维绘图是三维曲线图、三维网格图和三维曲面图 3 种基本类型。 与此对应， Matlab 也提供了一些三维基本绘图命令， 如三维曲线命令 plot3， 三维网格图命令 mesh 和三维表面图命令 surf。 1. 三维曲线 plot3(x,y,z)通过描点连线画出曲线， 这里 x,y,z 都是 n 维向量， 分别表示该曲线上点 集的横坐标、纵坐标、竖坐标。 例7 t=0:0.05:20* x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z), title('Line in 3-D Space') xlabel('X'), ylabel('Y'), zlabel('Z'), grid on 2. 三维网格图 命令 mesh(x,y,z)画网格曲面。这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标，命令 mesh(x,y,z)将该数据点在空间中描出，并连成网格。解 x=-3:0.1:3;y=-5:0.1:5; x1=ones(size(y'))*x;y1=y'*ones(size(x)); [x2,y2]=meshgrid(x,y); %meshgrid 函数生成 2D/3D 网格矩阵，用于分割空间，在绘 制 3D 网格图或者表面图时都需要调用该函数对生成绘图时所用的数据。 z1=(sin(x1.*y1)+eps)./(x1.*y1+eps); z2=(sin(x2.*y2)+eps)./(x2.*y2+eps);24 subplot(1,2,1),mesh(x1,y1,z1) subplot(1,2,2),mesh(x2,y2,z2) 例 8-1：绘制两个空间相交的平面 z=x+2y-1 和 z=x- y+2 解 ： x=-3:0.1:3;y=-3:0.1:3; [X,Y]=meshgrid(x,y); mesh(X,Y,X+2*Y-1) hold on mesh(X,Y,X-Y+2) 3. 表面图 命令 surf(x,y,z)画三维表面图，这里 x,y,z 是三个同维数的数据矩阵，分别表示数据 点的横坐标、纵坐标、竖坐标。解 [x,y]=meshgrid([-3:0.2:3]); z=(sin(x.*y)+eps)./(x.*y+eps); surf(x,y,z) 4. 旋转曲面 方法一：使用命令函数 cylinder [X,Y,Z] = cylinder(r) 这里的 r 表示构成旋转曲面的曲线。解 Matlab 程序如下。 x=0:10:600; [X,Y,Z]=cylinder(30*exp(-x/400).*sin((x+25*pi)/100)+130); surf(X,Y,Z) 方法二：将旋转曲面用参数方程表示。25 解 因为这里的函数是隐函数，化成显函数后有两支，必须使用参数方程，旋转面 的参数方程为画图的 Matlab 程序如下： alpha=[0:0.1:2*pi]'; beta=0:0.1:2* x=4*cos(alpha)*ones(size(beta)); y=(5+4*sin(alpha))*cos(beta); z=(5+4*sin(alpha))*sin(beta); surf(x,y,z) 或者利用绘制三维隐函数曲面图形的命令 ezsurf 或者 ezmesh，其命令格式如下： ezmesh(fun) ezmesh(fun,domain) domain 要求是一个向量，指定 x,y 轴的范围 ezmesh(funx,funy,funz) %参数方程绘图 funx(s,t), funy(s,t), and funz(s,t) over the square: - 2π & s & 2 π , -2π & t & 2 π. ezsurf(fun) ezsurf(fun,domain) ezsurf(funx,funy,funz) 画图的 Matlab 程序也可以写成 x=@(alpha,beta) 4*cos(alpha); y=@(alpha,beta) (5+4*sin(alpha))*cos(beta); z=@(alpha,beta) (5+4*sin(alpha))*sin(beta); ezsurf(x,y,z) 对于形如 f(x,y,z)=0 的三维隐函数的图像没有现成的函数可以画，但可以利用 isosurface 函数绘制三角网格图。例如 f=@(x,y,z)x.^2+y.^2+z.^2-10;%定义函数 f=x^2+y^2+z^2-10 [x,y,z]=meshgrid(linspace(-4,4,25));%设定格子大小和范围 val=f(x,y,z); [p,v]=isosurface(x,y,z,val,0);%用 isosurface 得到函数 f=0 图形的点和面26 patch('faces',p,'vertices',v,'facevertexcdata',jet(size(v,1)),'facecolor','w','edgecolor','flat'); %用 patch 绘制三角网格图并设定色彩 view(3); axis equal 5、绘制柱面 柱面平行于某个坐标轴，方程中不出现某个坐标轴的变量，方程表示为 F(x,y)=0 或 者 F(x,z)=0 或者 F(y,z)=0 示例一：画出方程 z ? 2 ? x2 表示的柱面 解：方程中的 x 是自变量矩阵， z 是因变量，则另一个自变量矩阵为 y，自变量平 面是 xoy 面， x 轴是真正的自变量， y 轴是柱面方向。Matlab 程序如下： u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X=u*ones(size(v)); %构成自变量矩阵 X Y=ones(size(u))*v; %构成自变量矩阵 Y Z=2-X.^2; %求因变量 Z mesh(X,Y,Z) 示例二：画出方程 ? x 2 ? y2 ? 1 表示的柱面 4解：方程整理为显示函数： y ? ?2 1 ? x 2 ，分别绘制正负两个分量。 clear,clc u=linspace(-5,5,10)';%设定参数列向量 u v=linspace(-5,10,10); %设定参数行向量 v X1=u*ones(size(v)); %构成自变量矩阵 X1 Y1=2*sqrt(1+X1.^2); %求因变量 Y1 Y2=-2*sqrt(1+X1.^2); Z1=ones(size(u))*v; %求因变量 Z1 mesh(X1,Y1,Z1),hold on mesh(X1,Y2,Z1) 6. 其它二次曲面 对于旋转面，如果母线的方程可以表示成关于旋转轴变量的显式函数，则可以直接 使用 Matlab 工具箱中的命令 cylinde r，否则必须把旋转面化成参数方程，然后使用 ezmesh 或 ezsurf 命令绘图。对于其它的二次曲面，如果可以写成显函数直接使用命令 ezmesh 或 ezsurf，否则必须先化成参数方程。还有一些特殊的绘制函数如27 cylinder,ellipsoid 等命令。解： （1）x=@(s,t) 3*sec(s)*cos(t); y=@(s,t) 3*sec(s)*sin(t); z=@(s,t) 2*tan(s); ezmesh(x,y,z) （2）x=@(s,t) 3*sec(s); y=@(s,t) 2*tan(s)*cos(t); z=@(s,t) 2*tan(s)*sin(t); ezmesh(x,y,z) （3）ezsurf(@(y,z) y.^2,50) %直接调用 ezsurf （4）x=@(s,t) 3*tan(s)*cos(t); y=@(s,t) 2*tan(s)*sin(t); z=@(s,t) tan(s); ezsurf(x,y,z) （5）ellipsoid(0,0,0,3,2,sqrt(6)) %专门绘制椭球面 （6）ezsurf(@(x,y)x*y) （7）x=@(s,t) 3*cos(s); %化为参数方程28 y=@(s,t) 2*sin(s); z=@(s,t) ezmesh(x,y,z) 7、绘制空间两曲面的交线 示例：绘制由水平截面与方程 z ? x 2 ? 2 y 2 构成的马鞍面形成的交线，并讨论等高线 和方向导数（梯度）的意义 解： 水平平面与曲面的交线就是等高线， 在 Matlab 中绘制等高线有两个命令， contour 和 contour3，前者把等高线画在 xoy 平面上，后者把等高线画在一定高度的平面上，使 之成为立体的，与所在曲面对应。 Matlab 程序如下： clc,clear [x,y]=meshgrid(-10:.2:10);%确定计算和绘图的定义域网格 z1=(x.^2-2*y.^2)+ %第一个曲面方程 a=input('a=(-50&a&50)'); z2=a*ones(size(x)); %水平面方程 z2=a， z2 必须与 x， y 具有相同的维数， subplot(1,3,1),mesh(x,y,z1);mesh(x,y,z2);%分别画出两个曲面 v=[-10 10 -10 10 -100 100];axis(v), %确定第一个分图的坐标系 colormap(gray); r0=abs(z1-z2)&=1; %求两曲面 z 坐标只差小于 0.5 的网格 zz=r0.*z2;yy=r0.*y;xx=r0.*x; %求这些网格上的坐标值，即交线坐标值 subplot(1,3,2),plot3(xx,yy,zz,'x'); %画出这些点 axis(v), %使第二个分图取第一个分图的坐标系 pause,subplot(1,3,3); contour3(x,y,z1,20); %用等高线命令求出 20 条不同高度的交线。 等高线与方向导数和梯度的概念密切相关，函数 z1 在每一点的梯度与该处的等高 线垂直，也就是指向最陡的方向。Matlab 中求梯度用 gradient 函数，quiver 函数画出梯 度向量，这两个函数要求在给定的点阵上求梯度和画梯度向量。Matlab 程序如下： [x,y]=meshgrid(-10:2:10); z1=(x.^2-2*y.^2)+ contour(x,y,z1,20); [px,py]=gradient(z1,2,2); %以步长为 2 求 z1 的梯度的 x， y 分量 quiver(x,y,px,py); %绘制梯度向量 8、离散点的图形绘制相关函数 griddata，TripScatterInterp，scatter（scatter3） griddata 用来对离散数据进行曲面拟合 (1)ZI = griddata(x,y,z,XI,YI) 用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量29 x,y,z。griddata 将返回曲面 z 在点（ XI,YI）处的插值。曲面总是经过这些数据点（ x,y,z） 的。输入参量（XI,YI）通常是规则的格点（像用命令 meshgrid 生成的一样） 。 XI 可以 是一行向量，这时 XI 指定一有常数列向量的矩阵。类似地， YI 可以是一列向量，它 指定一有常数行向量的矩阵。 (2)[XI,YI,ZI] = griddata(x,y,z,xi,yi) 返回的矩阵 ZI 含义同上， 同时， 返回的矩阵 XI,YI 是由行向量 xi 与列向量 yi 用命令 meshgrid 生成的。 (3)[XI,YI,ZI] = griddata(.......,method) 用指定的算法 method 计算： ?linear ?：基于三角形的线性插值（缺省算法） ； ?cubic?： 基于三角形的三次插值； ?n earest? ：最邻近插值法； 示例： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(- x.^2-y.^2); ti = -2:.25:2; [xi,yi] = meshgrid(ti,ti); zi = griddata(x,y,z,xi,yi); mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off 在新版本的 Matlab 中，该函数逐渐被 TriScatteredInterp 函数所代替， TriScatteredInterp 的示例如下： x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2); F = TriScatteredInterp(x,y,z); ti = -2:.25:2; [qx,qy] = meshgrid(ti,ti); qz = F(qx,qy);%计算指定位置的插值 mesh(qx,qy,qz); plot3(x,y,z,'o'); Matlab 中绘制散点图的命令是 scatter（scatter3） ，命令格式如下： scatter(X,Y,S,C) %X，Y 对应散点值，S 代表标记大小，C 代表颜色值，它都可以是 向量。 scatter3(X,Y,Z,S,C) 示例：30 A=[1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1; 4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1; 1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2; 4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2; 2.757,3.024,0.3;3.439,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3; 2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4; 3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5; 3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6; 3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7]; x=A(:,1);y=A(:,2);z=A(:,3); scatter(x,y,5,z)%散点图4.5 动态可视化图形Matlab 中的动画命令， moviein，getframe 和 movie ，用 getframe 把 Matlab 产生的 图形存储下来，每个图形成一个很大的列向量；再用 N 行这样的列保存 N 幅图，成为 一个大矩阵；movie 命令把他们连接起来重放，产生动画效果；moviein 用来预留存储空 间以加快运行的速度。有些情况也需要借助 pause 命令和循环语句来实现点轨迹的动态 演示。绘制彗星图的命令 comet 或者 comet3 也可以实现点的运动轨迹，不过这种方式 速度较快。 示例 1： axis equal, %把坐标设成相等比例 M=moviein(16);%为变量 M 预留 16 幅图的存储空间 for j=1:16 plot(fft(eye(j+16))); M(:,j)= end movie(M,30) %以每秒 30 帧的速度播放 M 中的图形 示例 2：使用 comet/comet3 绘制彗星图 （1） t = 0:.01:2* %返回当前坐标轴下的一帧图像31 x = cos(2*t).*(cos(t).^2); y = sin(2*t).*(sin(t).^2); comet(x,y); （2） n=10; t=n*pi*(0:0.0005:1); x=sin(t); y=cos(t); plot(x,y,'g'); comet(x,y,0.01); 示例 3：极坐标下一个点沿着图形移动的轨迹 a=0:.01:2* b=3; polar(a,b*(1-cos(a)),'*r'); for a=0:0.05:2*pi cla(findobj(gca,'color','r')); %findobj 函数用来使用特定的属性获取图形句柄， cla 用来清理由函数句柄指定的图形所在的坐标轴。 h2=polar(a,b*(1-cos(a)),'b'); set(h2,'linestyle','.','markersize',30); %设置句柄指定函数图形的属性 pause(0.001) end 示例 4： 普通二维坐标下一个点沿着曲线移动的轨迹， 程序与上面的结果基本相同。 x=-2*pi:0.1:2* y=sin(x); plot(x,y,'b-'); hold on for m=-2*pi:0.05:2*pi cla(findobj(gca,'color','b')); %cla 清除当前坐标轴上面的子对象 h=plot(m,sin(m),'r'); set(h,'LineStyle','.','Marker','*','MarkerSize',20); pause(0.001); end 示例 5：普通三维坐标下一个点沿着曲线移动的轨迹。 解： %三维螺旋线坐标 t1=10*pi*(0:;32 x1=cos(t1); y1=sin(t1); z1=-t1; %绘制曲线下面的水平直线 t2=(0:10)/10; x2=x1(end)*(1-t2); y2=y1(end)*(1-t2); z2=z1(end)*ones(size(x2)); %绘制垂直直线 t3=t2; z3=(1-t3)*z1(end); x3=zeros(size(z3)); y3=x3; %绘制曲线上面的水平直线 t4=t2; x4=t4; y4=zeros(size(x4)); z4=y4; %绘制曲线 x=[x1,x2,x3,x4]; y=[y1,y2,y3,y4]; z=[z1,z2,z3,z4]; plot3(x,y,z,'b','linewidth',3) axis off %绘制移动的点 h=line('Color',[1 0 0],'Marker','.','MarkerSize',40,'EraseMode','xor'); n=length(x); i=1;j=1; while 1 set(h,'xdata',x(i),'ydata',y(i),'zdata',z(i)); drawnow %使 Matlab 暂停目前的任务序列而去刷新屏幕 pause(0.0005); f=getframe(gcf); i=i+1; j=j+1; if i&n & j&450 end end33 5、Matlab 在高等数学中的应用5.0 Matlab 中的符号表示1、建立符号变量 Matlab 提供了 sym 和 syms 两个建立符号对象的函数 1）sym 函数，用来建立单个符号变量，调用格式为： 符号变量名=sym(?符号字符串 ?) 注：符号字符串可以是 常量、变量、函数或者表达式。 2）syms 函数，依次可以定义多个符号变量，调用格式为： syms 符号变量 1 符号变量 2 例如： f=sym(? co s(x)?),f=sym(? sin(x)^2=0? ); f=sin(x)+cos(x), 2、建立符号表达式 ? ? ? 利用单引号来生成表达式 利用 sym 函数建立符号表达式 使用已经定义过的符号变量组成符号表达式 注意：符号表达 …… 符号变量 n 注意 符号变量名之间不能使用任何标点符号，只能用空格隔开。式包括符号函数和符号方程，区别在于是否带有等号。 3、符号表达式的运算 (1)四则运算 四则运算分别使用 symadd, symsub,symmul 和 symdiv 来实现加减乘除运算，使用 sympow 实现幂运算。 （注：6.5 以后的版本中已经没有了这些函数，而直接使用+ - * / 即可） (2)提取分子和分母 如果符号表达式是一个有理分式或可以展开为有理分式，可利用 numden 函数来提 取符号表达式中的分子或者分母，调用格式为：34 [n,d]=numden(s) [n,d] = numden(sym(4/5)) returns n = 4 and d = 5. [n,d] = numden(x/y + y/x) returns (3)因式分解与展开 1）factors (s)分解因式 2）expand(s)展开 3）collect(s)合并同类项 4）collect(s,v)按变量 v 进行合并同类项 (4)符号表达式的化简 1）simplify (S)应用函数规则对 S 进行化简 2）simple (S) 进行综合化简 (5) 符号表达式与数值表达式之间的转换 1）利用函数 sym 可以将数值表达式表示成符号表达式 2）numeric 或者 eval 函数可以将符号表达式转换成数值表达式 3）函数 digits 和 vpa 配合替换函数 subs 进行转换。 digits 函数，digits(D)函数设置有效数字个数为 D 的近似解精度。 vpa 函数 ， R=vpa(S)符号表达式 S 在 digits 函数设置下的精度的数值解。 vpa(S,D) 符号表达式 S 在 digits(D)精度下的数值解。 subs 函数， subs(S, NEW, OLD) 4、符号函数的相关运算 （1）复合函数运算。 compose 函数 （2）反函数运算。 finverse 函数 5、符号代数方程求解 (1)线性方程组的求解，函数 linsolve ，solve ，可以得到方程的精确解 (2)多变量的非线性方程的符号解法，使用函数 fsolve ，调用格式有： ? ? ? ? ? X=fsolve(?fun?, X0) X=fsolve(? f un?,X0,options)options 为选择参数输入向量 X=fsolve(? f un?, X0,o ptio ns,?g radfun ? ),gradfun 为输入函数在 X 处的偏导数 X=fsolve( ? f un?,X0,o ptio ns,?g radfun ? ,P1,P2,… )P1,P2 为问题定性参数 [X,optio n s]=fsolve( ?fun ? , X0,…) 返回使用的优化方法的参数 n = x^2+y^2 , d = y*x35 注意：fun 必须使用@指定 示例：定义一个函数 function F = myfun(x,c) F = [ 2*x(1) - x(2) - exp(c*x(1)) -x(1) + 2*x(2) - exp(c*x(2))]; end 求解该函数： c = -1; %定义输入参数 x = fsolve(@(x) myfun(x,c),[-5;-5]) 6、多项式的表示方法 （1）多项式的表示方法――转化为向量问题 对于多项式P ( x) ? a 0 xn ? a 1 x n ?1 ? ... ? a n ?1 x ? a n用行向量表示： P ? [a , a ,..., a , a ] 0 1 n? 1 n ①系数向量直接输入法，MATLAB 自动将向量元素按降幂顺序分配给各系数值。函数 poly2sym 可以将向量表示的多项式转化为符号多项式表示。 注：由特征多项式生成 的多项式首项系数一定为 1；n 阶矩阵一般产生 n 次多项式。 ③由根创建多项式，由函数 poly 实现。注：若要生成实系数多项式，则根中的复数必 定对应共轭；生成的多项式向量包含很小的虚部时可用 real 命令将其过滤掉。 （2）多项式的运算 ①多项式求值。输入变量值代入多项式计算时以数组为单元的使用函数 polyval （对应 元素计算）;以矩阵(必须为方阵 )为计算单元求多项式的值用函数 ②多项式求根。两种方法，一种是调用函数 roots ，另一种是通过建立多项式的伴随矩 阵再求其特征值的方法得到多项式的所有根。 （使用 compan 和 eig 函数） ③多项式的乘除法运算。乘法使用函数 conv（向量卷积） ，除法使用函数 deconv ④多项式微分。微分函数 polyder ⑤多项式拟合。两种方法，一种是由矩阵的除法求解超定方程来进行，另一种是用拟 合函数 polyfit,调用方式为 polyfit(X,Y,n)和[p,s]=polyfit(X,Y,n)36 5.1 求极限Matlab 求极限的命令为 limit(expr, x, a) limit(expr, a) limit(expr) limit(expr, x, a, 'left') limit(expr, x, a, 'right') 其中 limit(expr, x, a)表示求符号表达式 expr 关于符号变量 x 趋近于 a 时的极限， limit(expr) 求表示缺省变量趋近于 0 时的极限。 例 15 求下列表达式的极限解 （1）clc, clear syms x f=(3*sin(x)+x^2*cos(1/x))/((1+cos(x))*log(1+x)); s=limit(f,x,0) s=simplify(s)（2）clc, clear syms x a s1=limit(((x+2*a)/(x-a))^x,x,inf) s2=limit(((x+2*a)/(x-a))^x,x,-inf)5.2 求导数Matlab 的求导数命令为 diff(expr) diff(expr, v) diff(expr, sym('v')) diff(expr, n) diff(expr, v, n)37 diff(expr, n, v) 其中 diff(expr)表示求表达式 expr 关于默认变量的 1 阶导数，diff(expr, v, n)和 diff(expr, n, v)都表示求表达式，expr 关于符号变量 v 的 n 阶导数。解 （1）clc, clear syms x y=log(x+sqrt(1+x^2)); s=diff(y,x,2) s=simple(s) %对符号函数进行化简 （2）a=[1,0.5,3.5,6]; da=diff(a)5.3 求一元函数极值计算的 Matlab 程序如下 syms x y=x^3+6*x^2+8*x-1; dy=diff(y); dy_zero=solve(dy), dy_zero_num=double(dy_zero) %变成数值类型 ezplot(y) %符号函数画图5.4 求积分1. 求不定积分 Matlab 求符号函数不定积分的命令为 int(expr) int(expr, v)38 解 syms x I=int(1/(1+sqrt(1-x^2))); pretty(I) %写成符号表达式的形式 2. 求定积分 （1）求定积分的符号解 Matlab 求符号函数的定积分命令为 int(expr, a, b) int(expr, v, a, b)解 syms x I=int(cos(x)*cos(2*x),-pi/2,pi/2) （2）求定积分的数值解 例 20 求下列积分的数值解输入 I=quadl(@(t) (t-3*t.^2+2*t.^3).^(-1/3),eps,0.5) 得到 I=1.4396。 注：quad()函数使用 simpson 法，即用二次曲线逼近被积函数， quadl()采用更高阶 的逼近方法。它们的调用方法是：q=quad(fun,a,b,tol)或者 quadl(fun,a,b,tol)， fun 可以是39 符号方程，也可以是匿名函数。 ）输入 I=dblquad(@(x,y)sqrt(1- x.^2-y.^2).*(x.^2+y.^2&=x),0,1,-0.5,0.5) 得到 I= 0.6028。 注： Matlab 中计算二重积分的函数形式为 q= dblquad(fun,xmin,xmax,ymin,ymax,tol)， 一般只用于积分区域为矩形的情况。为了使用该函数，有时需要将非矩形区域化成矩形 区域。?2 2 示例：计算二重积分 ?? ?( x ? y )dxdy 积分区域 ? 由 x=1,y=x 及 y=0 围成。 ?解：首先绘制积分区域，Matlab 程序如下： clear,clc fill([0,1,1,0],[0,0,1,0],'y');hold on%绘制积分区域 fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r') %绘制单元条 gtext('y=x'); gtext('x=1');pause gtext('y=0') 按照矩形区域调用 dblquad 函数，程序如下： I=dblquad(?(x.^2+y.^2).*(y-x&0)?,0,1,0,1) 求得 I=0.3333 （3）输入 fun3=@(x,y,z)z.^2*log(x.^2+y.^2+z.^2+1)./(x.^2+y.^2+z.^2+1).*(z&=0 &... z&=sqrt(1- x.^2-y.^2)); %该语句使用了续行符... I=triplequad(fun3,-1,1,-1,1,0,1) 得到 I= 0.1273。 %三重积分计算函数40 Matlab 符号求解的程序如下 clc, clear syms r z theta I1=int((r^2+z^2)*r, z, r^2,sqrt(2-r^2)); %求最内层积分 I2=int(I1,r,0,1); %求中间层的积分 I=int(I2,theta,0,2*pi) %求最终的积分结果 pretty(I) % 分数线中间显示的格式 Inum=double(I) %把符号解化成数值解 解法二 求数值解的 Matlab 程序如下 clc, clear f=@(x,y,z) (x+y+z).^2.*(z&=x.^2+y.^2 & x.^2+y.^2+z.^2&=2); I=triplequad(f,-sqrt(2),sqrt(2),-sqrt(2),sqrt(2),0,sqrt(2)) 求得积分的值为 2.4486。 注：quad 与 int 的区别，int 可以求具有解析解的不定积分和定积分，quad 只能求解 定积分。5.5 求解微分方程（组）1、微分方程(组)的解析解41 在 MATLAB 中用 D 表示导数， 例如 Dy 表示 y?， D2y 表示 y ?? ， Dy(0)=5 表 示 y ?(0)=5 等，符号常微分方程求解通过函数 dsolve 实现，调用格式为： dsolve( ?e? ,? c? ,? v? )包括微分方程，初始条件和指定变量三个部分，求解常微分方程 e 在初始条件 c 下的特解，v 描述方程中的自变量，省略则默认以 t 为自变量，若没有初 始条件，则获得通解。 例 1：求 du ? 1 ? u 2 的通解 dt解：Matlab 程序为： dsolve('Du=1+u^2','t') ?d2y dy ? 2 ? 4 ? 29 y ? 0 例 2：求微分方程 ??dx 的特解 dx ? ??y (0) ? 0, y?(0) ? 15 解：Matlab 程序为： dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') ? dx ? dt ? 2x ? 3 y ? 3 z ?? ? dy 例 3：求常微分方程组的通解 ? ? 4 x ? 5 y ? 3z ??dt ? dz ? 4x ? 4 y ? 2 z ? ??dt 解：Matlab 程序为： [x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z') 一阶齐次线性常微分方程组和非齐次线性方程组可以直接通过矩阵操作求解。 对于一阶齐次线性微分方程组： X ? ? AX , X ? ? x1 , xn ?, A ? (aij )n ?n 其解的形式为：X (t ) ? eA(t ?t0 ) X 0对于非齐次线性方程组： X ? ? AX ? f (t ), X (t 0 ) ? X 0 解的形式为：X (t ) ? eA(t ?t0 ) X 0 ???? eA (t ? s ) f ( s )ds? 例 4： X ? ? ? 0 2 ?1? X , X (0) ? ? 2 ? ? ? ? ? ? ? ?0 0 3 ? ?1 ?2 1 3? ?1 ? ?42 解：a=[2,1,3;0,2,-1;0,0,2]; x0=[1;2;1];x=expm(a*t)*x0 pretty(x); %化简结果?1 0 0 ? ? 0 ? ?0 ?? ? ? ? ? ? ?? 例 5： X ? ? ? 2 1 ?2 ? X ? ? 0 ? , X (0) ? ?1 ?? t ? ? ? ?3 2 1 ?e cos(t ) ?1解： a=[1 0 0;2 1 -2;3 2 1]; fs=[0;0;exp(s)*cos(2*s)]; x0=[0;1;1]; tx=int(expm(a*(t-s))*fs,s); %先求不定积分，再计算定积分的结果，提高速度 xstar=subs(tx,s,t)-subs(tx,s,0); x=expm(a*t)*x0+ x=simple(s),pretty(x) %simple 直接显示化简结果，simplify 显示化简过程，pretty 修 改显示格式 2、微分方程（组）的数值解 当难以求得微分方程的解析解时，可以求其数值 解。Matlab 中求解微分方程数值解 的方法有以下几个：ode45，ode23，ode113，ode23s ，ode15s 。调用形式如下： [t,x]=solver(‘ f’, ts, x0, options)其中，t 是自变量的值，x 是函数值，solver 表示上述 5 种求解器的一种；‘ f’是由待解方程写成的 m 文件名，ts=[t0, tf]表示自变量的初值和 终值，options 用于设定误差限，通过函数 options=odeset(‘reltol’, rt, ‘abstol’, at)分别设 定相对误差和绝对误差。 不同的函数代表不同的内部算法，ode45 表示 4/5 阶龙格- 库塔- 费尔贝格算法（ RK 算法），是解非刚性常微分方程的首选方法，ode23 采用 2/3 阶 RK 方法，ode113 采用 多步法，效率一般比 ode45 高。ode15s ，ode23s,ode23t,ode23tb 用来解刚性常微分方程。 注 1：在解 n 个未知函数的的方程组时，x0 和 xn 均为 n 为向量， m 文件中的待解 方程组应以 x 分量形式写出； 注 2：使用 Matlab 求数值解时，高阶微分方程必须等价的变换成一阶微分方程组。 ? ? d2 x 2 dx ? ? ?1000(1 ??x ) ? x 0 例 1：求解 van der Pol 刚性微分方程 ??dt 2 的数值解 dt ? x (0) ? 2; x '(0) ? 0 ??43 解：令 y1=x， y2=y1?，则微分方程变为一阶微分方程组： ? y1' ? y 2 ? 2 ??y 2 ' ? 1000(1 ? y1 ) y 2 ? y1 ? y1(0) ? 2; y 2(0) ? 0 ?? 写成 m 文件如下所示： function dy=fun1(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1); 取 [t0 ,tf]=[0,3000] ， 输入命 令： [t,y]=ode15s('fun1',[0,3000],[2,0]); plot(t,y(:,1),'-') ? y1 ' ? y2 y3 ?? ? ? y1 y3 ?y ' 例 2：求解微分方程组 ?? 2 ? y3 ' ? ? 0.51 y1 y2 y1 (0) ? 0, y2 (0) ? 1; y 3 (0) ? 1 解：建立方程组的 m 文件： function dy=fun2(t,y) dy=zeros(3,1); %初始化 dy(1)=y(2)*y(3); dy(2)=y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 取 t0=0，tf=12，Matlab 命令如下： [t,y]=ode45('fun2',[0 12],[0 1 1]) plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+')? ? x ? ? ( y ??x) ?? ? y ? ? x ? y ? xz 在初始条件[x0,y0,z0]=[5,13,17]时的解，其 例 3：求解 Lorez 方程 ? ? ? z ? ?? z ??xy ??中 ? ? 10, ? ? 28, ? ? 8 / 3 。 解：rho=10;beta=29;lamda=8/3;44 f=@(t,Y)[rho*(Y(2)-Y(1));beta*Y(1)-Y(2)-Y(1)*Y(3);- lamda*Y(3)+Y(1)*Y(2)];%方程 组必须使用列向量表示 [t,y]=ode45(f,[0,30],[5,13,17]) subplot(2,2,1) plot(t,y(:,1),'*') %x 曲线 subplot(2,2,2) plot(t,y(:,2),'X') subplot(2,2,3) plot(t,y(:,3),'O') subplot(2,2,4) plot3(y(:,1),y(:,2),y(:,3)) %绘制空间轨迹图 3、边值问题的 Matlab 数值解 Matlab 中使用 bvp4c（bvp5c）和 bvpinit 命令求解常微分方程中的两点边值问题， 其调用格式为： sol = bvp4c(odefun,bcfun,solinit,options) solinit = bvpinit(x,yinit,parameters) %该函数用来给 sol 提供初始猜测解。5.6 级数求和Matlab 级数求和的命令为 r = symsum(expr, v) r = symsum(expr, v, a, b) 其中 expr 为级数的通项表达式， v 是求和变量，a 和 b 分别为求和变量的起始点和 终止点，若没有指明 a 和 b，a 的默认值为 0，b 的默认值为 v-1。 无穷数列的累加称为级数，当取其前面若干有限项时，得到的是部分和，数列 a 累 加形成的新序列可用 s=cumsum(a)实现。 a 的长度为 n， 则 s 的长度也为 n， 即每一个 s(k) 是数组 a 中前 k 项的和。注意 cumsum 与 sum 的区别，sum(a)得到的是一个数，即序列 s 中最后一项。 示例： clc,clear k=1:10; a=1./k.^2;45 s=cumsum(a) ss=sum(a) 例 22 求如下级数的和解 （1）syms n f1=(2*n-1)/2^n; s1=symsum(f1,n,1,inf) （2）syms n f2=1/n^2; s2=symsum(f2,n,1,inf) 补充一点： Matlab 求数组元素乘积的函数 B=prod(A) 如果 A 是一个 m 行一列的 （向 量），则这种用法即返回这 m 个元素的乘积； 如果 A 是一个 m 行 n 列的矩阵，则 A 的每一列都被看做一个 m 行 1 列的向量，分别计算每个向量中元素的乘积，返回给 B， 因此 B 是一个 1 行 n 列的数组。5.7 Matlab 求解优化问题Matlab 求解各种形式的极值问题被集成在优化工具箱中。 下面是优化工具箱 4.0 版本的功能示意图：46 47 变量 f描述 线性规划目标函数 f*X 或二次规划的目标 函数 X’*H*X+F*X 中线性项的系数向量 非线性优化的目标函数. fun 必须为行命令 对象或 M 文件、 嵌入函数、 或 MEX 文件的 名称 二次规划目标函数 X’*H*X+f*X 中的二次 项的系数矩阵 A 矩阵和 b 向量分别为线性不等式约束： AX≤b 中的系数矩阵和右端向量 Aeq 矩阵和 beq 向量分别为线性等式约束 Aeq*X=beq 中的系数矩阵和右端向量 X 的下限和上限向量 迭代初始点坐标 函数最小化的区间 优化选项参数结构调用函数 linprog,quadprog fminbnd,fminsearch,fminunc, fmincon,lsqcurvefit,lsqnonlin, fgoalattain,fminimax quadprog linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax linprog,quadprog,fgoalattain, fmincon,fminimax,lsqcurvefit, lsqnonlin 除 fminbnd 外所有函数 fminbnd 所有优化函数funHA,bAeq,beqvlb,vub X0 x1,x2 options变量描述 描述退出条件： exitflag&0 表示目标函数收敛于 x exitflag=0 表示已达到函数评价或迭代的最大次 数 exitflag&0 表示目标函数不收敛 由优化函数求的 x 值 解 x 处的目标函数值调用函数exitflagx所有优化函数 linprog,quadprog,fgoalattain fmincon,fminimax,lsqcurvefit lsqnonlin,fminbndfevaloutput包含优化结果信息的输出结构 所有优化函数 iterations,Algorithm,FuncCount（函数评价次数）进行极值计算时有两种方式，通过命令 optimtool 打开优化工具箱 GUI 求解和使用48 上述提供的命令函数进行求解。 示例一：使用单纯形计算最小值 min f ? ? 4x1 ? x 2 ? ? x1 ? 2x 2 ? 4 ? 2 x ? 3 x ? 12 ?? 2 s.t . ? 1 ? x1 ? x 2 ? 3 x1 , x 2 ? 0 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=[-4;-1]; A=[-1 2;2 3;1 -1]; b=[4;12;3]; [x,fval,exitflag,output]=linprog(f,A,b) 示例二：非线性优化 f ( x ) ? ? x1 x 2 x 3 ? x1 ? 2x 2 ? 2x 3 ? 0 ?? s.t . ? ? x1 ? 2 x 2 ? 2x 3 ? 72 方法一：通过优化工具箱求解 方法二：输入 Matlab 命令： clc,clear f=@(x)(-x(1)*x(2)*x(3)); A=[-1 -2 -2;1 2 2]; b=[0;72]; [x,fval,exitflag,output]=fmincon(f,[10;10;10],A,b) 其他求解极值的函数主要 有： (1)x=fminbnd(fun,x1,x2)，求函数在区间[x1,x2]上的极小值对应的自 变量值 示例三：求 f(x)=x^4-x^2+x-1 在区间 [-2,1]上的极小值。 解： Matlab 命令如下： [x,feval,exitflag,output]=fminbnd('x^4-x^2+x-1',-2,1) (2) x=fminsearch(fun,x0)无约束非线性极小值求解， x0 为起始搜索点。49 示例四：fminsearch 函数求 f(x)=1/((x1-2)^2+3)-1/(2(x2+1)^2-5),初始点[0,0]， (3) x=fminunc(fun,x0,options)求解函数的一个局部极小值 (4) x=fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon 包括非线性约束 ; (5) x=lsqnonlin(fun,x0,lb,ub);非线性最小二乘优化 (6) x=linprog(f,A,b,Aeq,beq,lb,ub,x0)求解线性规划问题 (7) x=bintprog(f,A,b,Aeq,beq,x0)求解 0-1 规划 (8) x=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0);求解二次规划5.8 Matlab 求解插值和拟合问题5.8.1 Matlab 插值工具箱 1、一维插值函数 interp1，命令格式为： yi=interp1(x,Y,xi,? method ?)其中 method 指定的方法为 nearest,linear,spline,cubic， pchip （分段三次 Hermite 插值） ， 所有的插值方法都要求 x 是单调的。x 与 Y 是具有相同大小 的向量，求在 xi 点处的插值函数值 yi 一维快速插值函数 interp1q，命令格式为： yi = interp1q(x,Y,xi)，其中 x 是单调递增的列向量，Y 是列向量或者行数为 length(x) 的矩阵，xi 是一个列向量。 yi = interp1(x,Y,xi,method,'extrap') %对于超出 x 范围的 xi 中的分量将执行特殊的外 插值法 extrap。 示例 1：在一 天 24 小时内，从零点开始每间隔 2 小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28 ，27，25，20 ，18，15 ，13，推测中午 12 点（即 13 点）时的温 度． 解： x=0:2:24; y=[12 9 9 10 18 24 28 27 25 20 18 15 13]; a=13; y1=interp1(x,y,a,'spline') xi=0:1/3600:24; %绘制一天内的温度变化曲线 yi=interp1(x,y,xi, 'spline'); plot(x,y,'o' ,xi,yi) 示例 2： x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,'ro',x,y);50 x1=0:0.02:1; %要插值点 y1=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1); y2=interp1(x,y,x1);y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');y5=interp1(x,y,x1,'cubic'); plot(x1,[y2' y3' y4' y5'],':',x,y,'ro',x1,y1); legend('linear','spline','nearest','cubic','样本点','原函数'); %计算各插值的最大误差 [max(abs(y1-y2)),max(abs(y1- y3)),max(abs(y1-y4)),max(abs(y1- y5))]; 2、二维数据内插值 interp2，插值函数为二元函数。 ZI = interp2(X,Y,Z,XI,YI,? method? ) X 和 Y 分别是 m 维和 n 维的向量，表示节点，Z 为 n*m 的矩阵，表示节点值，XI， YI 是为一维数组，表示插值点，XI 与 YI 是方向不同的向量，一个是行向量，一个是列 向量，ZI 是行数为 YI 的维数，列数为 XI 的维数的矩阵，表示得到的插值。若 XI 与 YI 中有在 X 与 Y 范围之外的点，则相应地返回 nan（Not a Number） 。 method 为插值 方法。 如果在某区域测量了若干个节点的高程（节点值） ，为了画出比较精确的等高线图， 需要先插入更多的点（插值点） ，然后计算这些点的高程。 如果是三次样条插值，可以使用命令 pp=csape({x0,y0},z0,conds,valconds);z=fnval(pp,{x,y}) 示例 1： years = 0; service = 10:10:30; wage = [150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243]; w = interp2(service,years,wage,15,1975) %计算在 1975 年工作了 15 年的员工的工资[x,y]=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); surf(x,y,z);figure [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1); z2=interp2(x,y,z,x1,y1,'cubic'); z3=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1);figure %比较误差 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);51 surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.08])surf(x1,y1,abs(z0-z3)); axis([-3,3,-2,2,0,0.08]) interp2 函数能够较好的进行二维插值运算，但是它只能处理以网格形式给出的数 据，如果数据是随机给出的，此时可以使用 griddata 函数。其调用格式为： 格式:z=griddata(x0,y0,z0,x,y,? method ? ) ? v4 ?:MATLAB4.0 提供的插值算法，公认效果较好； ? linear ?：双线性插值算法（缺省算法） ； ? nearest? ：最临近插值； ?cubic?：双三次插值。clear x=-3+6*rand(199,1);y=-2+4*rand(199,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); plot(x,y,'*');plot3(x,y,z,'*'); [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1); z2=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z2); %误差比较 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); surf(x1,y1,abs(z0-z1));axis([-3,3,-2,2,0,0.15]) surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15]) 3、 三维插值运算函数 interp3,n 维网格插值 interpn， 其调用格式同 interp1 和 interp2， 对应的三维网格生成函数为[x,y,z]=meshgrid(x1,y1,z1)和非网格生成函数 griddata3(),griddatan()，他们同 griddata 调用格式相同。 示例 1：解：[x,y,z]=meshgrid(-1:0.2:1); [x0,y0,z0]=meshgrid(-1:0.05:1); V=exp(x.^2.*z+y.^2.*x+z.^2.*y).*cos(x.^2.*y.*z+z.^2.*y.*x); V0=exp(x0.^2.*z0+y0.^2.*x0+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0); V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:))52 slice(x0,y0,z0,V1,[-0.5,0.3, 0.9],[0.6,-0.1],[-1,-0.5,0.5,1]) title('Slives for Four Dim Figures'); 4、三次样条插值函数 已知平面上 n 个点 xi,yi(i=1...n)，其中 x 1&x 2& … &xn ，S(x)为三次样条函数需要满足 三个条件： (1)S(xi)=yi(i=1...n)，即函数经过样本点； (2)S(x)在每个子区间 [xi,xi+1]上为三次多项式：(3)S(x)在整个区间 [x1,xn]上有连续的一阶及二阶导数。 Matlab 中有如下函数： y=interp1(x,Y,xi,?spline? ); y=spline (x,y,xi); %对于给定的离散的测量数据 x,y（称为断点） ，要寻找一个三项多 项式 y = p(x) ，以逼近每对数据(x,y)点间的曲线。 示例 1：解：x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); sp=csapi({x0,y0},z); fnplt(s p); 示例 2:对离散分布在 y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算 解：x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; y = exp(x).*sin(x); xx = 0:.25:20; yy = spline(x,y,xx); plot(x,y,'o',xx,yy) pp=csape (x,y,conds); %conds 是边界条件选项，可以取以下值： 'complete',给定边界一阶导数. 'not-a-knot',非扭结条件,不用给边界值. 'periodic',周期性边界条件,不用给边界值.53 'second',给定边界二阶导数. 'variational',自然样条 (边界二阶导数为 0) pp=csape (x,y,conds,valconds);y= ppval(pp,x)， 求插值点的函数值必须调用 ppval 函数。 示例 3：给定数据 x=[1 2 4 5]， y=[1 3 4 2]，边界条件 S”(1)=2.5，S”(5)=-3，计算插 值。 解：x=[1 2 4 5];y=[1 3 4 2]; pp=csape(x,y,'second',[2.5,-3]); pp.coefs %显示每个区间上三次多项式的系数 xi=1:0.1:5; yi=ppval(pp,xi); %计算插值点对应的函数值 plot(x,y,?o?,xi,yi) Matlab 还提供了一个生成三次参数样条类的函数 S=csapi(x,y) 其中 x=[x 1,x 2,….,xn], y=[y1,y2,…,yn]为样本点。S 返回样条函数对象的插值结果， 包括子区间点、各区间点三次多项式系数等。 可用 fnplt()绘制出插值结果，其调用格式： fnplt(S) 对给定的向量 xp，可用 fnval()函数计算，其调用格式： yp=fnval(S,xp) 其中得出的 yp 是 xp 上各点的插值结果。 示例 4：解：x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); sp=csapi(x,y); fnplt(sp) c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:),... sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:)] 处理多个自变量的网格数据三次样条插值类的函数： S=csapi({x1,x 2,…,xn},z) 其中 xi 为自变量标识， z 是网格数据的样本点，S 是三次样条曲线的对象。 5.8.2 Matlab 中的拟合方法 1、线性最小二乘法拟合曲线 当用 m 次多项式拟合给定数据时，一般使用函数 polyfit，命令格式为：54 a=polyfit(x,y,m)，输入参数 x，y 为要拟合的数据，m 为拟合多项式的次数，输出参 数 a 为拟合多项式的系数向量。 多项式在 x 处的函数值 y 可以使用函数 polyval，格式为： y=polyval(a,x)，a 为向量表示的多项式 示例：某乡镇企业
年的生产利润如下表所示，试预测 1997 年和 1998 年 的利润。 年份 利润 解： x=[92 95 1996]; 1 122 3 152 5 196 y=[70 122 144 152 174 196 202]; a=polyfit(x,y,1); %根据利润增长的情况，使用一次线性拟合函数 y97=polyval(a,1997); y98=polyval(a,1998); 2、最小二乘优化 常用的求解最小二乘优化的函数 lsqlin,lsqcurvefit,lsqnonlin,lsqnonneg， 这些函数都 位于 Matlab 的优化工具箱中，下面只举例说明它们的用法。 x= lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)求解的模型为：例 1：x=[19 25 31 38 44]'; y=[19 32.3 49 73.3 97.8]'; r=[ones(5,1),x.^2]; %拟合形如 y=a+b*x^2 的多项式，如果拟合的多项式的形式为 y=a0+a1*x+a2*x^2，则 r 应该写成 r=[ones(5,1),x,x.^2] ab=lsqlin(r,y); %计算多项式的系数 x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; plot(x,y, ? o? ,x0, y0, ? r ? ) % 绘制对比曲线查看结果55 lsqcurvefit 函数求解的数学模型为函数的调用格式为： x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)； fun 参数是定义函数 F(x,xdata)的 M 文 件， x0 为初始点， xdata 是输入数据， ydata 是观测数据， 它们都是向量或者矩阵， F(x,xdata) 的维数和 ydata 一致。 例 2：用给定的观测数据拟合函数 y ? e? k1 x1sin(k2 x2 ) ? x23 中的参数 k1 和 k2。解： （1）首先编写 M 文件，定义拟合函数 function y=fun(k,xdata) f=exp(-k(1)*xdata(:,1)).*sin(k(2)*xdata(:,2))+xdata(:,3).^2; end （2）调用函数 lsqcurvefit x0=[23.73 5.49 1.21;…]; %x0 是 25 行 3 列的输入样本 y0=[15.02;….]; %y0 是 25 行 1 列的观测值 k0=rand(2,1); %初始化待求参数的值 lb=zeros(2,1);ub=[20;2]; %给定参数的上下界，没有边界时用 []表示 k=lsqcurvefit(@fun,k0,x0,y0,lb,ub) lsqnonlin 函数用来解决的数学模型是： 命令格式为： x=lsqnonlin(fun,x0,lb,ub,options);其中 fun 是定义向量函数 f(x)的 M 文件或者函数句 柄。 例 3：用最小二乘法拟合 y ? x0,y0 已知。 解：x0=-10:0.01:10; y0=normpdf(x0,0,1); %计算标准正态分布密度函数在 x0 处的取值 F=@(xs)1/sqrt(2*pi)/xs(2)*exp(-(x0-xs(1)).^2)/xs(2)^2/2- y0; %定义函数( x ? ? )2 ? 2 ? 2?1 2? ?e中的位未知参数?, ?，假设数据值56 xs0=rand(2,1);%非线性拟合要有初始值 xs=lsqnonlin（F,xs0）; 例 4： 示例 3:体重约 70kg 的某人在短时间内喝下 2 瓶啤酒后， 隔一定时间测量他的 血液中的酒精含量(mg/100ml)，得到如下数据，使用所给数据用函数 f=atbect 进行拟合， 求出常数 a,b,c。 解： t=[0.25 0.5 0.75 1 1.5 2 2.5 3 h1=log(h); % 进行对数变换 f=inline('a(1)+a(2).*log(t)+a(3).*t','a','t'); [x,r]=lsqcurvefit(f,[1,0.5,-0.5],t,h1); plot(t,h,'*'); ezplot('exp(4.9*log(t)-0.2663*t)',[0.2 16]); 3.5 4 4.5 5:16];h=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4];lsqnonneg 函数求解的数学模型为： 命令格式为： x=lsqnonneg(C,d,options)，该函数的用法同 lsqlin 用法类似。 3、Matlab 曲线拟合的图形用户界面操作 拟合一元函数的命令是 cftool，拟合二元函数的命令 sftool，非线性拟合的命令是 nlintool。 主要拟合步骤是： (1)将 数据导入工作空间； (2)运行 cftool 命令； (3)对数据进行预 处理； (4)选择适当的模型进 行拟合； (5)生成一些相关的统计量，并进行预测操作57 示例 1：利用下表的数据拟合方程 y ? a0 e 1 1 e 2axa x2解： （1）把数据加载到 Matlab 工作空间中， a=[5 2 1.5;5 5 2.3;7 10 3.5;8 4 5.2;8 6 1.8;3 3 2.9;7 8 3.8;7 3 4.5;2 6 4.7]; x1=a(:,1); x2=a(:,2); z=a(:,3); (2)在命令窗口中输入 sftool，打开如下界面（3）输入自变量 X，Y 的数据 x1,x2，因变量 Z 的数据 z，并定义要拟合的函数58 a0*exp(a1*x)*exp(a2*y) （4）用鼠标点击 fit 进行拟合，结果如下59 6、Matlab 在线性代数中的应用6.1 向量组的线性相关性求列向量组 A 的一个最大线性无关组可用命令 rref(A)将 A 化成行最简形，其中单 位向量对应的列向量即为最大线性无关组所含向量，其它列向量的坐标即为其对应向量 用最大线性无关组线性表示的系数。 例 23 求下列矩阵列向量组的一个最大无关组。解 编写 Matlab 程序如下 format rat %数据是有理分数表示 a=[1,-2,-1,0,2;-2,4,2,6,-6;2,-1,0,2,3;3,3,3,3,4]; b=rref(a) format %恢复到短小数的显示格式解 编写 Matlab 程序如下 format rat60 a=[2,2,-1;2,-1,2;-1,2,2];b=[1,4;0,3;-4,2]; c=rref([a,b]) format %恢复到短小数的显示格式6.2 齐次线性方程组在 Matlab 中，函数 null 用来求解零空间，即满足 Ax=0 的解空间，实际上是求出解 空间的一组基（基础解系）。格式如下 : z=null(A) %z 的列向量为方程组的正交规范基，满足。 ZT Z=E z=null(A,? r? ) %z 的列向量是方程 Ax=0 的有理 基。 例 25 求方程组的通解解 编写程序如下 format rat a=[1,2,2,1;2,1,-2,-2;1,-1,-4,-3] b=null(a,'r') %求有理基 syms k1 k2 x=k1*b(:,1)+k2*b(:,2) %写出方程组的通解 format %恢复到短小数的显示格式6.3 非齐次线性方程组Matlab 中解非齐次线性方程组可以使用 “\”。虽然表面上只是一个简简单单的符号， 而它的内部却包含许许多多的自适应算法，如对超定方程（无解）用最小二乘法，对欠 定方程（多解）它将给出范数最小的一个解。61 对于非齐次线性方程组 Ax=b，Matlab 的求解命令为 x=A\b 或者(x=inv(A)*b) 无论 Ax ??b 是有唯一解，无穷多解还是无解时，Matlab 总是给出一个解。 （1）当 Ax ??b 唯一解时，Matlab 给出的就是该唯一解。 （2）当 Ax ??b 有无穷多解时，Matlab 给出的是最小范数解。 （3）当 Ax ??b 无解时，Matlab 给出的是最小二乘解。 另外求解欠定方程组 （多解）可以使用求矩阵 A 的行最简形命令 rref(A)，求出所 有的基础解系。 例 26 求超定方程组解 编写程序如下 a=[2,4;3,-5;1,2;2,1]; b=[11;3;6;7]; solution=a\b 注：上面解超定方程组的“\”可以用伪逆命令 pinv 代替，且 pinv 的使用范围比“ \”更 加广泛，pinv 也给出最小二乘解或最小范数解。 例 27 用最小二乘解法解方程组解 编写程序如下 format rat a=[1,1,0;1,0,1;1,1,1;1,2,1]; b=[1;2;0;-1]; x1=a\b %这里\和 pinv 是等价的 x2=pinv(a)*b62 format %恢复到短小数的显示格式 例 28 求解欠定方程组解 编写程序如下 format rat a=[1,-1,-1,1,0;1,-1,1,-3,1;1,-1,-2,1,-1/2]; b=rref(a) format %恢复到短小数的显示格式 例 29 利用表 1 的数据拟合函数解 拟合的 Matlab 程序如下 clc, clear d=[5 2 1.5 5 5 2.3 7 10 3.5 8 4 5.2 8 6 1.8 3 3 2.9 7 8 3.8 7 3 4.5 2 6 4.7];63 m=size(d,1); %获取 d 的行数 a=[ones(m,1),d(:,1).^3,d(:,1).*d(:,2),d(:,2).^2]; % 求 d(:,3)=xs*a 中的 xs=a\d(:,3) %求的多项式的系数 例 30 利用表 1 的数据拟合函数clc, clear d=[5 2 1.5 5 5 2.3 7 10 3.5 8 4 5.2 8 6 1.8 3 3 2.9 7 8 3.8 7 3 4.5 2 6 4.7]; m=size(d,1); a=[ones(m,1),d(:,1),d(:,2)]; b=log(d(:,3)); xs=a\b; xs(1)=exp(xs(1)) 例 31 使用用户图形界面解法命令 sftool 拟合例 29 和例 30 中的参数。6.4 相似矩阵及二次型解 二次型的矩阵为编写如下程序 A=[0,1,1,-1;1,0,-1,1;1,-1,0,1;-1,1,1,0]; [P,D]=eig(A)64 解 Matlab 程序如下 a=[2 -2 0; -2 4 0; 0 0 5]; b=eig(a); if all(b&0) fprintf('二次型正定 \n'); else fprintf('二次型非正定 \n'); end [c,d]=eig(a) %c 为正交变换的变换矩阵6.5 线性代数中的其它应用计算的 Matlab 程序如下 clc, clear t=[12 10 1; 2 6 1; 4 -10 1]; S=1/2*abs(det(t)) %不考虑顶点的排列，加上绝对值函数 abs。 解法二 利用复数运算的命令计算。 clc, clear a=complex(12,10); b=complex(2,6); c=complex(4,-10); %注意复数的写法。 plot([a,b,c,a]) %画出三角形 theta=angle((b-a)/(c-a)) %利用复数运算求两个向量夹角 S=1/2*abs(b-a)*abs(c-a)*abs(sin(theta))65 解法一 使用 Matlab 的命令 orth，Matlab 程序如下 clc, clear format a=[1 -1 4; 2 3 -1; -1 1 0]; b=orth(a)就是一组规范正交基。 编写程序实现Schimidt正交化过程，程序如下 clc, clearformat a=[1 -1 4; 2 3 -1; -1 1 0]; [m,n]=size(a); %求矩阵 a 的维数 b(:,1)=a(:,1); c(:,1)=b(:,1)/norm(b(:,1)); for i=2:n s=a(:,i); for k=1:i-1 s=s-dot(a(:,i),c(:,k))*c(:,k); %正交化过程 end c(:,i)=s/norm(s); %单位化 end c66 7、Matlab 数据处理7.1 Matlab 中的默认数据文件 mat 文件1、将工作区变量保存到文本文件中 save filename var1 var2 var3… -mat 默认保存的文件扩展名为.mat，如果保存的 数据需要跨平台处理，需要采用如下格式： save filename.dat var1 var2 … -ascii 另外： save filename var1 var2 … -append 附加到已经存在的 Mat 文件中 例 36 把 Matlab 工作空间中的数据矩阵 a，b，c 保存到数据文件 data1.mat 中。 save data1 a b c 注：Matlab 中的默认数据文件 mat 文件可以省略后缀名。 2、将文本文件中的数据导入工作区 load filename 例 37 把例 36 生成的 data1.mat 中的所有数据加载到 Matlab 工作空间中。 load data17.2 纯文本文件可以把 word 文档中整行整列的数据粘贴到纯文本文件，然后调入到 Matlab 工作空 间中。 例 38 把纯文本文件 data2.txt 加载到工作空间。 a=load('data2.txt'); 或者是 a=textread('data2.txt'); 例 39 使用 dlmwrite 命令把矩阵 b 保存到纯文本文件 data3.txt。 dlmwrite('data3.txt',b) 例 40 生成服从标准正态分布随机数的矩阵，然后用 fprintf 命令保存到纯文本文件 data4.txt。67 解 clc, clear fid=fopen('data4.txt','w'); a=normrnd(0,1,100,200); fprintf(fid,'%f\n',a'); fclose(fid); 注：对于高维矩阵，用 dlmwrite 构造的纯文本文件，Lingo 软件不识别；为了 Lingo 软件识别，纯文本文件必须用 fprintf 构造，而且数据之间的分割符为“\ n”。 示例： M = gallery('integerdata', 100, [5 8], 0); dlmwrite('myfile.txt', M, 'delimiter', '\t') dlmread('myfile.txt') dlmread('myfile.txt', '\t', 2, 3) dlmread('myfile.txt', '\t', 'C1..G4')7.3 Excel 文件xlswrite(filename,A,sheet,range) 将阵列 A 写入 Excel 文件 filename 中 sheet 表格的 range 指示的单元格内。 例 41 把一个矩阵写到 Excel 文件 data5.xls 表单 Sheet2 中 B2 开始的域中。 105 a=rand(5,10); xlswrite ('data5.xls',a,'Sheet2','B2') %默认的表格为 Sheet1 [num,txt,raw] = xlsread(filename,sheet,range)从 Excel 文件 filename 中的 sheet 表格 中读取 range 单元格内的数据返回。 例 42 把例 41 生成的 Excel 文件 data5.xls 中表单 Sheet2 的域“C3：F6”中的数据赋 给 b。 b=xlsread('data5.xls','Sheet2','C3:F6') 示例： values = {1, 2, 3 ; 4, 5, 'x' ; 7, 8, 9}; headers = {'First', 'Second', 'Third'}; xlswrite('myExample.xls', [ values]); A = xlsread('myExample.xls') subsetA = xlsread('myExample.xls', 1, 'B2:C3') columnB = xlsread('myExample.xls', 'B:B') [ndata, text, alldata] = xlsread('myExample.xls')7.4 字符串数据例 43 统计下列五行字符串中字符 a、c、g、 t 出现的频数。68 1.aggcacggaaaaacgggaataacggaggaggacttggcacggcattacacggagg 2.cggaggacaaacgggatggcggtattggaggtggcggactgttcgggga 3.gggacggatacggattctggccacggacggaaaggaggacacggcggacataca 4.atggataacggaaacaaaccagacaaacttcggtagaaatacagaagctta 5.cggctggcggacaacggactggcggattccaaaaacggaggaggcggacggaggc 解 把上述五行复制到一个纯文本数据文件 shuju.txt 中，编写如下程序 clc fid=fopen('shuju.txt','r'); %以只读的形式打开文件 i=1; while (~feof(fid)) %判断是否到文件末尾 data=fgetl(fid); a=length(find(data==97)); b=length(find(data==99)); c=length(find(data==103)); d=length(find(data==116)); e=length(find(data&=97&data&=122)); f(i,:)=[a b c d e a+b+c+d]; i=i+1; end f, he=sum(f) dlmwrite ('pinshu.txt',f); dlmwrite('pinshu.txt',he,'-append'); fclose(fid); %打开的文件使用完之后一定要关闭。 Matlab 中关于文件操作的常用函数： 1) textread 函数可以读取格式化为列数据的 ASCII 码文件， 每一列可以是不同的数 据类型，调用形式为： [a,b,c,...]=textread(filename,format,n) 示例 1：读文件 test_input.dat [first,last,blood,gpa,age,answer]=textread(?t est_input.dat ? ,? %s %s %s %f %d %s?) 通过使用%*s 可以忽略对应的列数据。 示例 2：假设文件 test.dat 中的数据如下： Jam}
我要回帖
更多推荐
- ·发问下美世教育申请初中去新加坡留学利弊咋样？
- ·上海观澜小学是第几梯队同时对口川沙南校的小区里那些？
- ·驾考科目一最牛作弊作弊如何处理
- ·朱志辉是哪个高校的高校党委宣传部是做什么的任职？
- ·研究生考试辅导书系列：马哲、毛概、近现代史纲要及思修毛概马原近代史先学哪个法基高频考点及模拟试题覆盖全面吗？
- ·武汉万达暑假工对三级教育工人年龄要求有什么要求
- ·广州工程职业学院技术职业学院总体如何
- ·cnn互动英语在win10 faster rcnn上安装完后，却不能使用，该怎么办
- ·虫洞存在于dnf异次元裂缝空间的裂缝里，人类的肉眼，看不到，摸不着。相当于量子泡沫，的皱纹。
- ·今天回答了很多个问题，没有一个被采纳，不百度知道提升回答采纳为什么？过去每天都有30-50个采纳的，今天是怎么了呀
- ·家里这几种食物尽量少吃，容易吃出肾结石多吃什么食物的食物，都来
- ·电工kmd电工fr是什么意思思
- ·如何看待插队问题，到底是别人没测试人性的问题，还是你
- ·新浪微博企业认证公函取消认证，在公函中需增加什么取消认证理由？
- ·DeleteFront(La);波胆是什么意思思
- ·五河县公园首府开盘哪里有专业开墙的
- ·每天拖实木地板用什么拖把的时候，水里就会洗出很多这样的虫子，这是什么虫子？
- ·脑子急转弯大全会转弯打一动物
- ·为什么深圳只适合打工？好少人在深圳黄依玲 许晴安居乐业
- ·每天拖实木地板用什么拖把的时候，水里就会洗出很多这样的虫子，这是什么虫子？
- ·怎么用MATLABphp 从数组中随机取值生成一个取值在(5，80)的四阶整数矩阵
- ·哪位大神可以教我个纯物理简单百宝箱鼠标连点器器。（我只是个初一的，最好不要那么麻烦）谢谢！
- ·科幻美剧无垠的太空第一集下载季第一集11分46年左右的bgm是什么
- ·还有mymymy吗？ 戳爷mymymy百度云的。谢谢了
- ·哪位伙伴知道这小咬是什么虫子子？被咬了是什么样子？怎么办
- ·请问这个动作表示人物动作的词语什么意思
- ·我有一个梦想演讲视频想
- ·高中数学辅导，求解
- ·观察从事物的发展变化着眼不断变化是寻找其奥妙的门路
- ·怀孕七个月可以同房吗零几天可以做飞机出国吗
- ·宝宝昨天打小孩生水痘怎么治疗疫苗，今天可以洗澡了吗
- ·得了粉瘤切除手术视频，术后吃韭菜有什么影响吗
- ·疼慢慢发展到一站就疼，咳嗽到心口疼怎么回事事
- ·当罪犯 再犯罪说自己是因为精神疾病而犯罪的，如何鉴定其犯罪
- ·焦虑头疼一个月了是脑瘤会头疼吗吗
- ·颈椎引起手麻手指疼痛能吃脉血康胶囊？
