第5 章 CHAPTER 5 MATLAB科学绘图 图形绘制与可视化是MATLAB 语言的一大特色。MATLAB 提供了一系列直观、简单 的二维、三维图形绘制命令与函数,可以将实验结果和仿真结果用可视的形式显示出来。 从R2014b 版本开始,MATLAB 提供全新的绘图体系与绘图函数,虽然早期版本的一些 命令仍能使用,但以后的版本中这些命令可能会被逐渐淘汰,所以本书将尽量按照新的体例 介绍图形绘制的方法。 5.1 节介绍简单二维曲线的绘制方法,包括由数据绘制曲线的方法和由数学表达式绘制 曲线的方法,还介绍二维曲线的修饰方法。5.2 节介绍特殊二维曲线的绘制方法,如极坐标曲 线的绘制方法、离散点的表示方法、统计图的绘制方法及动态图形的绘制方法。5.3 节介绍三 维图形的绘制方法,包括三维曲线图与曲面图的绘制方法、视角设置方法与三维动态图形的 处理方法。5.4 节介绍二元与三元隐函数图形的绘制方法。 5.1 简单二维图形绘制 二维图形是科学研究中最常见,也是最实用的图形表示。本节首先介绍将数据用二维曲 线表示出来的方法,然后介绍将数学函数用曲线表示出来的方法。 5.1.1 基于数据的绘图 假设用户已经获得了一些实验数据。例如,已知各个时刻t = t1, t2, · · · , tn,测得这 些时刻的函数值y = y1, y2, · · · , yn,则可以将这些数据输入MATLAB 环境中,构成向量 t = [t1, t2, · · · , tn] 和y = [y1, y2, · · · , yn]。如果用户想用图形表示二者之间的关系,则用 plot(t,y) 即可绘制二维图形。可以看出,该函数的调用是相当直观的。 例5-1 试绘制出显函数y = sin(tan x) ? tan(sin x) 在x ∈ [?π, π] 区间内的曲线。 解解决这种问题的最简捷方法是采用下面的语句直接绘制函数曲线。 >> x=-pi:0.05:pi; % 以0.05 为步距构造自变量向量 y=sin(tan(x))-tan(sin(x)); plot(x,y) % 求出并绘制各个点上的函数值 这些语句可以绘制出该函数的曲线,如图5-1 所示。不过由这里给出的曲线看,得出的曲线似 乎有问题。 值得指出的是,由MATLAB 的plot() 函数绘制出的“曲线”不是真正的曲线,只是给 第5 章MATLAB 科学绘图65 图5-1 给定函数的曲线表示 出各个数值点间的折线。如果给出的数据点足够密,或突变少些,则看起来就是曲线了,故以 后将称之为曲线。 例5-2 试重新生成密集些的数据点,得出例5-1 函数正确的曲线表示。 解仔细观察图5-1 中给出的曲线可以看出,在±π/2 附近图形好像有问题,其他位置还是比 较平滑的。为什么会出现这样的现象呢?观察sin(tan x) 项,由于在±π/2 附近括号内的部分将趋 近于无穷大,所以导致其正弦值变化很不规则,会出现强振荡。 可以考虑全程采用小步距,或在比较粗糙的x ∈ (?1.8,?1.2) 及x ∈ (1.2, 1.8) 两个子区间内 选择小步距,其他区域保持现有的步距,这样可以将上述的语句修改为 >> x=[-pi:0.05:-1.8,-1.799:0.0001:-1.2, -1.2:0.05:1.2,... 1.201:0.0001:1.8, 1.81:0.05:pi]; % 以变步距方式构造自变量向量 y=sin(tan(x))-tan(sin(x)); plot(x,y) % 求出并绘制各个点上的函数值 这样将得出如图5-2 所示的曲线。可见,这样得出的曲线在剧烈变化区域内表现良好。前面解释过, 在±π/2 区域内出现强振荡是正常现象。 如果全程都选择0.0001 这样的小步距,也能得出看起来完全一致的效果。 >> x=-pi:0.0001:pi; % 以0.0001 为步距构造自变量向量 y=sin(tan(x))-tan(sin(x)); plot(x,y) % 求出并绘制各个点上的函数值 图5-2 改变步距后的曲线 从这个例子可以看出,不能过分地依赖于MATLAB 绘制的曲线,需要对曲线的正确性 66 MATLAB/Simulink 实用教程——编程、计算与仿真 做检验。比较有效的方法是选择不同的步距,观察得出的曲线是不是吻合,如果吻合则可以 认为曲线正确,否则需要选择更小的步距绘制曲线后再进行检验,直至得出吻合的结果。 在实际应用中,plot() 函数的调用格式还可以进一步扩展: (1)t 仍为向量,而y 为如下的矩阵,则将在同一坐标系下绘制m条曲线,每一行和t 之 间的关系将绘制出一条曲线。注意,这时要求y 矩阵的列数应该等于t 的长度。 y =26664 y11 y12 · · · y1n y21 y22 · · · y2n ... ... ... ... ym1 ym2 · · · ymn 37775 (2)t 和y 均为矩阵,且假设t 和y 矩阵的行数和列数均相同,则可绘制出t 矩阵每行和y 矩阵对应行之间关系的曲线。 (3)假设有多对这样的向量或矩阵(t1, y1),(t2, y2),· · · ,(tm, ym),则可以用下面的语 句直接绘制出各自对应的曲线: plot(t1,y1,t2,y2,· · · ,tm,ym) (4)曲线的性质,如线型、粗细、颜色等,还可以使用下面的命令进行指定: plot(t1,y1, 选项1,t2,y2, 选项2,· · · ,tm,ym, 选项m) 其中“选项”可以按表5-1 中说明的形式给出,其中的选项可以进行组合。例如,若想绘制红 表5-1 MATLAB 绘图命令的各种选项 曲线线型曲线颜色标记符号 选项意义选项意义选项意义选项意义选项意义 '-' 实线'b' 蓝色'c' 蓝绿色'*' 星号'pentagram' I '--' 虚线'g' 绿色'k' 黑色'.' 点号'o' 圆圈 ':' 点线'm' 红紫色'r' 红色'x' 叉号'square' □ '-.' 点划线'w' 白色'y' 黄色'v' ? 'diamond' ? 'none' 无线'^' △ 'hexagram' A '>' ? '<' ? 色的点划线,且每个转折点上用五角星表示,则选项可以使用组合字符串'r-.pentagram'。 (5)除了表5-1 给出的简洁绘图参数方法之外,plot() 函数还允许以 plot(~, 参数名1, 参数值1, 参数名2, 参数值2,· · · ) 的形式指定绘图参数,其中,~表示上述的正常调用格式。常用的参数名与参数值在表5-2 中 列出。 (6)还可以由h=plot(~) 格式调用plot() 函数,绘图的同时返回曲线的句柄h,以后 可以通过该句柄读取或修改该曲线的属性。如果想在句柄为h的坐标系下绘制曲线,还可以 给出plot(h,t,y,· · · ) 命令。 例5-3 分形树的数学模型:任意选定一个二维平面上的初始点坐标(x0, y0),假设可以生成 一个在[0, 1] 区间上均匀分布的随机数γi,那么根据其取值的大小,可以按下面的公式生成一个新 第5 章MATLAB 科学绘图67 表5-2 常用的绘图控制参数 参数名参数值 LineSpec 曲线线型控制字符串,如'r-.pentagram' 字符串,可以参考表5-1 LineWidth 曲线的线宽,默认的宽度是0.5pt,其中1pt=0.3527mm。 MeshDensity 自动步距计算点的密度,默认值为23,主要用于后面将介绍的fimplicit() 和fplot() 等基于函数表达式的曲线绘制函数,增大该值可以提高曲线精度 Color 曲线颜色,除了表5-1 指定的8 种颜色外,还可以设定为RGB 分量[r,g,b] MarkerEdgeColor 标记的边缘颜色,事实上就是标记自身的颜色 MarkerSize 标记的大小,默认值6pt 的坐标点(x1, y1)[18]。 (x1, y1) ?8>><>>: x1 = 0, y1 = y0/2, γi < 0.05 x1 = 0.42(x0?y0), y1 = 0.2 + 0.42(x0+y0), 0.05 ? γi < 0.45 x1 = 0.42(x0+y0), y1 = 0.2 ? 0.42(x0?y0), 0.45 ? γi < 0.85 x1 = 0.1x0, y1 = 0.2 + 0.1y0, 其他 试用圆点绘制出分形树的结果。 解分形树的数据可以由下面给出的语句直接计算。有了这些语句就可以得出x与y 向量,这 样,由下面的命令可以直接绘制出如图5-3 所示的分形树图。 >> v=rand(10000,1); N=length(v); x=0; y=0; for k=2:N, gam=v(k); % 针对每个γi 值,更新x、y if gam<0.05, x(k)=0; y(k)=0.5*y(k-1); elseif gam<0.45 x(k)=0.42*(x(k-1)-y(k-1)); y(k)=0.2+0.42*(x(k-1)+y(k-1)); elseif gam<0.85 x(k)=0.42*(x(k-1)+y(k-1)); y(k)=0.2-0.42*(x(k-1)-y(k-1)); else, x(k)=0.1*x(k-1); y(k)=0.1*y(k-1)+0.2; end, end plot(x,y,'.','MarkerSize',5) % 注意,不能忽略最后的圆点标记 图5-3 分形树的显示 68 MATLAB/Simulink 实用教程——编程、计算与仿真 5.1.2 基于函数表达式的绘图 如果已知数学函数,还可以利用fplot() 函数绘制函数曲线,其调用格式为fplot(f), 其中,f 可以是匿名函数描述的函数句柄,也可以是描述函数的符号表达式或符号函数,默 认的绘图区间为[?5, 5]。如果想指定绘图区域,还可以给出fplot(f,[xm,xM])。 MATLAB 早期版本提供的ezplot() 函数也可以用于绘制类似的曲线,默认的绘图区 间为[?2π, 2π],不过其绘图函数调用方式与数学函数描述方式与fplot() 函数是不同的,这 里不过多讨论该函数的使用方法。 例5-4 试用fplot() 函数重新绘制例5-1 的函数曲线。 解可以考虑用符号函数描述原来的数学函数,给出下面的语句就可以绘制出函数的二维 曲线,如图5-4 所示。从效果上看,与自选数据点得出的曲线基本一致。该曲线还自动绘制了 x = ?π/2 处的虚线。 >> syms x; f(x)=sin(tan(x))-tan(sin(x)); fplot(f,[-pi,pi]) 图5-4 由fplot() 函数绘制的函数曲线 还可以用匿名函数的形式描述原函数,得出一致的结果。用匿名函数描述已知数学函数时,应 该采用点运算。对本例而言,sin() 与tan() 函数的作用与点运算一致,无须特殊处理。 >> f=@(x)sin(tan(x))-tan(sin(x)); fplot(f,[-pi,pi]) 类似于plot() 函数,fplot() 函数也可以支持不同的调用格式,例如,若想在同一坐标 系内绘制多个函数,则可以给出命令fplot([f1,f2,· · · ,fn]),其中,fi 为第i 个数学函数的 函数句柄或符号表达式,此外,该函数允许带有不同的选项,可以返回图形句柄等。 例5-5 考虑正弦函数f(t) = sin x。由高等数学课程中学习的Taylor 级数展开可以得出函数 的有限项逼近表达式为y(t) = x9/362880 ? x7/5040 + x5/120 ? x3/6 + x。试在同一坐标系下绘 制在x ∈ [?4, 4] 区间内两个函数的曲线,并评价函数逼近的效果。 解可以用符号表达式直接表示原函数与Taylor 级数表达式,然后调用fplot() 函数,同时绘 制两条曲线,如图5-5 所示。为方便起见,这里用手工的方法将Taylor 级数曲线设置为虚线,具体的 设置方法将在后面介绍。 >> syms x; f=sin(x); y=x^9/362880-x^7/5040+x^5/120-x^3/6+x; fplot([f,y],[-4 4]) % 同时绘制两个函数的曲线 第5 章MATLAB 科学绘图69 图5-5 函数的Taylor 逼近对比研究 如果想用匿名函数表示数学函数,应该使用点运算。 >> f=@(x)sin(x); y=@(x)x.^9/362880-x.^7/5040+x.^5/120-x.^3/6+x; 5.1.3 参数方程曲线绘制 如果某函数由参数方程给出 x = x(t), y = y(t), tm ? t ? tM (5-1) 且x(t) 由函数句柄hx 表示,y(t) 由句柄hy 表示,其中,函数句柄可以同时为匿名函数,或同 时为符号表达式,则可以由fplot(hx,hy,[tm,tM]) 命令绘制其轨迹曲线。 例5-6 Lissajous 曲线是有两个不同频率正弦函数构成的参数方程,试绘制出x(t) = sin t, y = sin 1.25t,t ∈ [0, 30] 的函数曲线。 解由符号表达式描述参数方程,并绘制出Lissajous 函数曲线,如图5-6 所示。 >> syms t; x=sin(t); y=sin(1.25*t); fplot(x,y,[0,30]) 图5-6 Lissajous 曲线 默认状态下的fplot() 函数自动选择曲线绘制的参数,例如,给出[tm,tM] 时自动选择时间 步距,有时可能得出完全错误的结果。这时,用户可以自行选择Meshdensity 等控制参数,增大其 值,直至得出正确的结果。 例5-7 重新考虑例5-6 给出的Lissajous 曲线。如果t ∈ [0, 1000] 会可能得出错误的曲线,如何 绘制正确的曲线? 70 MATLAB/Simulink 实用教程——编程、计算与仿真 解直接给出下面的曲线,则可能绘制出如图5-7 所示的错误图形。 >> syms t; x=sin(t); y=sin(1.25*t); fplot(x,y,[0,1000]) 图5-7 错误的Lissajous 曲线 这时,需要人为指定MeshDensity 选项的值,将其设置成300 或更大的值,得出的结果与图 5-6 中给出的完全一致。 >> fplot(x,y,[0,1000],'MeshDensity',300) 和plot() 函数类似,fplot() 函数曲线自动绘制的结果同样也应该检验。设置不同的 MeshDensity 参数是一种检测方法。一般说来,较大的值对应着更精确的结果。选择不同的 MeshDensity 参数值,如果得出一致的曲线,则说明选择合理。 5.1.4 双y 轴曲线 前面介绍的plot() 函数可以在同一坐标系下同时绘制多条曲线,不过在某些特定的应 用中,如果两条曲线的幅值差异过大,则可以为图形设置双y 轴,分别表示不同的曲线。具体 做法为,调用yyaxis left 和yyaxis right 命令设置坐标系,绘制所需的曲线。早期版本 中的plotyy() 在当前版本下仍可以使用,不过这里不推荐该函数。 例5-8 考虑两个函数y1 = sin x 与y2 = 0.01 cos x,试绘制它们的曲线。 解当然可以考虑使用plot() 函数直接绘制曲线,不过由于它们的幅值相差过于悬殊,所以 y2 曲线看起来就像一条直线,分辨度很差。所以,可以考虑给曲线设置两个纵轴,这样得出的曲线 如图5-8 所示。从图中可见,实线的标度在左侧坐标轴给出,虚线的标度在右侧坐标轴给出,用这样 的方法可以很好地显示两条幅值相差悬殊的曲线。 >> x=0:0.01:2*pi; y1=sin(x); y2=0.01*cos(x); yyaxis left; plot(x,y1), yyaxis right; plot(x,y2,'--') 5.1.5 图形修饰与编辑 曲线绘制之后,还可以对绘制的图形进行进一步修饰。表5-3 中列出了常用的图形修饰 命令。用户可以利用这些命令在绘制的图形上做适当的修饰,也可以利用图形窗口提供的工 具对图形进行处理。 第5 章MATLAB 科学绘图71 图5-8 双纵轴的函数曲线 表5-3 常用的图形修饰命令 修饰命令调用格式命令的解释 title() title(str) 给图形加标题,标题的内容由字符串str 表示 xlabel() xlabel(str) 给x 轴加标签,而ylabel(str) 命令给y 轴加标签 text() text(x,y,str) 在图形的(x, y) 坐标处添加文字说明 gtext() gtext(str) 允许用户用鼠标选择添加文字说明的点 legend() legend(s1,s2,· · · ) 与图形对应的线型介绍图例,sk 为第k 条曲线的说明文字 annotation annotation(s,x,y) 曲线加标注,其中,s 为标注类型,可以选择'arrow'(箭头)、'line'(线 段)、'doublearrow'(双向箭头)等,由向量[x1,x2] 和[y1,y2] 表示起 点(x1, y1) 和终点坐标(x2, y2) hold hold on 可以用hold on 或hold off 锁定或释放坐标系。如果坐标系被锁定,再 使用plot() 这类命令将在现有的图形上叠印新的曲线;hold off 命令 解除锁定状态。还可以由key=ishold 命令查询坐标系的锁定状态,返 回的值为0 或1 zoom zoom on 局部放大功能,可以用鼠标选择想放大的区域;zoom off 命令可以取消 局部放大功能;还可以使用zoom xon 和zoom yon 单独放大某坐标轴 MATLAB 图形窗口的“插入”菜单的内容如图5-9 所示。这里为排版方便将整个菜单截 成3 段给出。其中,图5-9(a)还给出了图形窗口的主菜单。可见,表5-3 中的大多数函数或命 令都可以由图形窗口的“插入”菜单直接实现。 (a)菜单系统及“插入”菜单上段(b)中段(c)下段 图5-9 MATLAB 图形窗口的“插入”菜单 MATLAB 图形窗口的“查看”菜单在图5-10(a)中给出,其中,前面3 项控制工具栏的 形式。如果前3 项都选中,则图形窗口的工具栏如图5-10(b)所示。默认状态下,只显示第一 行工具。该工具栏的按钮控制图形的编辑状态。如果选中,则可以用鼠标选择曲线单独编 辑,如果反选,则取消编辑状态。这时,如果把鼠标移入某个坐标系内,则坐标系的右上角出 现如图5-10(c)所示的坐标系工具栏。 72 MATLAB/Simulink 实用教程——编程、计算与仿真 (a)“查看”菜单 (b)完全打开的工具栏 (c)坐标系工具栏 图5-10 MATLAB 图形窗口的“查看”菜单与工具栏 5.1.6 图形数据的提取 MATLAB 的菜单系统允许用户将当前的图形窗口直接存入文件,文件的后缀名为fig。 如果想恢复图形,直接用菜单从.fig 文件读入即可。如果想获得当前图形窗口的数据,或从 打开的.fig 文件获取数据,则需要进入图形编辑状态,用鼠标选中感兴趣的曲线。这时,由 gco(get current object,获得当前对象)命令即可获得选中曲线的句柄,再使用get() 函数 即可提取曲线数据。 x=get(gco,'xData'); y=get(gco,'yData'); MATLAB 提供了一系列类似于gco 的命令,如gcf(获得当前窗口句柄)、gca(获得当 前坐标系句柄)。还可以用clf 命令清空当前的图形窗口。 5.2 特殊二维图形 除了前面介绍的plot() 函数与fplot() 函数之外,还可以采用line() 函数绘制曲线, 其作用于plot() 函数相仿。不同的是,line() 函数在当前的坐标系下叠印曲线。另外, line() 函数只支持两个输入变元。 此外,MATLAB 还支持其他的二维图形绘制命令,常用的特殊曲线绘制命令与调用格 式在表5-4 中给出。本节将介绍这些二维图形的绘制方法。 表5-4 MATLAB 提供的特殊二维曲线绘制函数 函数名意义常用调用格式函数名意义常用调用格式 comet() 彗星轨迹图comet(x,y) bar() 二维条形图bar(x,y) compass() 罗盘图compass(x,y) errorbar() 误差限图形errorbar(x,y,ym,yM) feather() 羽毛状图feather(x,y) fill() 二维填充图fill(x,y,c) hist() 直方图hist(y,n) loglog() 对数图loglog(x,y) quiver() 矢量图quiver(x,y) polarplot() 极坐标图polarplot(x,y) stairs() 阶梯图形stairs(x,y) semilogx() x-半对数图semilogx(x,y) stem() 火柴杆图stem(x,y) semilogy() y-半对数图semilogy(x,y) 5.2.1 极坐标 极坐标是以一个给定的点为原点,以从原点出发的某条线为极轴构造的坐标系。空间上 某一点到原点的距离为ρ,该点与原点的连线和极轴的夹角为θ,该角度以从极轴出发逆时 第5 章MATLAB 科学绘图73 针方向的转角为正方向。这样的有序对(ρ, θ) 称为极坐标。极坐标下的曲线一般可以表示为 显式函数ρ = ρ(θ),称为极坐标方程。传统意义下,极坐标方程中的ρ ? 0。如果将其拓展到 实数空间,再通过直角坐标系变换得出的曲线,该方程可以理解为广义极坐标方程。 MATLAB 提供了polarplot() 函数,其调用格式为polarplot(θ,ρ),其中,θ 和ρ为 给定数据构成的向量。该函数早期版本的polar() 函数调用格式也是一样的,但新版本下不 建议使用。在当前的版本下,即使ρ<0,也可以通过坐标变换绘制极坐标曲线。如果不希望 绘制ρ<0 部分的曲线,则可以将ρ<0 时的函数值设置为NaN,以便绘图时自动排除这些点。 例5-9 试用极坐标绘制函数polarplot() 绘制出ρ = 5 sin(4θ/3) 的极坐标曲线。 解由初中数学课程介绍可以立即得出结论,该函数的周期为3π/2。所以若想绘制极坐标曲 线,则应该先构造一个θ 向量,然后求出ρ 向量,调用polarplot() 函数就可以绘制出所需的极坐 标曲线,如图5-11(a)所示。 >> theta=0:0.01:3*pi/2; rho=5*sin(4*theta/3); % 生成极坐标向量 polarplot(theta,rho) % 极坐标图曲线 (a)θ ∈ (0, 3π/2) (b)ρ > 0 的新极坐标曲线 图5-11 极坐标曲线 如果想排除ρ < 0 的部分,则可以将满足ρ < 0 的点强行置为NaN,这样,绘图时自动排除这些 点,绘制的新极坐标曲线如图5-11(b)所示。 >> rho(rho<0)=NaN; polarplot(theta,rho) % 屏蔽掉NaN 点 观察得出的曲线,绘制的极坐标图好像不完全。如何绘制完整的极坐标曲线呢?很多极坐标函 数都是周期函数,如果能确定函数的周期则可绘制完整的曲线,这样就导致了新的问题——如何 确定周期?对MATLAB 这样的工具而言,甚至不必考虑周期问题。若想绘制完整的极坐标曲线,选 一个较大的θ 范围,如0?θ?20π 或更大范围,其完整的极坐标曲线如图5-12(a)所示。通过试凑方 法可见,该函数的实际周期是6π。如果只考虑ρ ? 0 部分,则得出的结果如图5-12(b)所示。 >> theta=0:0.01:20*pi; rho=5*sin(4*theta/3); % 重新生成数据 polarplot(theta,rho) % 极坐标图曲线 figure; rho(rho<0)=NaN; polarplot(theta,rho) 例5-10 试绘制非周期极坐标函数ρ = e?0.1θ sin 3θ 的函数曲线。 解选择自变量θ 的变化范围θ ∈ (0, 10π),则可以生成函数的极坐标数据,并绘制出极坐标图 形,如图5-13(a)和图5-13(b)所示。 74 MATLAB/Simulink 实用教程——编程、计算与仿真 (a)θ ∈ (0, 20π) (b)ρ ? 0 部分 图5-12 更大范围的极坐标曲线 >> theta=0:0.001:10*pi; rho=exp(-0.1*theta).*sin(3*theta); polarplot(theta,rho) figure; rho(rho<0)=NaN; polarplot(theta,rho) (a)广义极坐标曲线(b)ρ ? 0 部分 图5-13 非周期函数的极坐标曲线 很多极坐标函数是周期性函数,选择一个周期内的θ 值就可以把极坐标曲线绘制出来,这里的 函数不是周期性的,所以不论取多大范围都不能绘制完整的极坐标曲线。 5.2.2 离散数据的图形表示 本节将首先给出离散信号的定义,然后介绍基于MATLAB 的离散信号表示方法,还将 介绍经过零阶保持器后的输出信号。 一般的离散信号可以表示为时间序列y1,y2,· · · ,yn。离散信号当然可以用plot() 函数直接绘制,不过更恰当的是用stem() 绘图语句绘制出来的火柴杆图,其调用格式为 stem(t,y),其中,t 为时间点构成的向量。如果离散信号后面跟一个零阶保持器(zero-order hold,ZOH),则该信号将变成连续信号且在每一个采样周期内都保持常值,那么该函数可以 用stairs() 函数绘制出来的阶梯信号表示,其调用格式为stairs(t,y)。 例5-11 假设已知离散信号的数学表示为f(t) = sin t sin 7t,且t = kT,k = 0, 1, 2, · · · , 31, 第5 章MATLAB 科学绘图75 T = 0.1 s 称为采样周期,表示每间隔0.1 s 采集一次函数信号。试用图形方式表示该序列信号。 解根据给出的采样周期,可以由下面命令生成时间序向量t,然后计算出离散信号的函数值, 并直接绘制出其示意图,如图5-14(a)所示。 >> T=0.1; t=(0:31)*T; % 生成时间采样点向量 f=sin(t).*sin(7*t); stem(t,f) % 计算函数值并绘图 (a)火柴杆图(b)阶梯信号图 图5-14 离散信号表示 如果用stairs() 函数取代stem(),则可以绘制出如图5-14(b)所示的阶梯图。 >> stairs(t,f) % 绘制阶梯信号图 5.2.3 统计图形绘制 直方图与饼图是统计学领域经常使用的绘图工具。本节将先给出直方图与频度的定义, 然后通过例子介绍直方图与饼图的绘制方法与技巧。 假设已知一组离散的检测数据x1,x2,· · · ,xn,并且这组数据都位于(a, b) 区间内,则可 以将这个区间分成等间距的m个子区间,使得b1 = a,bm+1 = b。将每个随机量xi 依其大小 投入相应的子区间,并记子区间(bj , bj+1) 落入的数据个数为kj,j = 1, 2, · · · ,m,则可以得 出fj = kj/n,称为频度。 还可以利用histogram() 函数求取各个子区间的频度,该函数的调用格式如下: k=histogram(x,b); % 该函数返回的k 是结构体型数据 f=k.Values/n; bar(b(1:end-1)+δ/2,f/δ); 其中δ = x2 ? x1 为等间距子区间宽度。选择向量b 和f,可以绘制出频度的直方图。注意,直 方图得出的向量长度比b 向量长度短1,另外,调用bar() 函数之前应该将b 向量后移半个 子区间宽度。计算落入每个子区间数据点个数时建议使用新的histogram() 函数,不建议使 用早期版本的hist() 函数。下面将通过例子演示直方图的表示方法。 例5-12 生成满足参数为b = 1 的Rayleigh 分布的30000×1 伪随机数向量,并用直方图验证 生成的数据是否满足期望的分布。 解可以由raylrnd() 函数生成30000 × 1 的伪随机数向量,选择向量x,这样可以通过 histogram() 计算每个子区间落入的数据个数,其结果的Values 属性为落入每个子区间的点数。 可以用函数bar() 近似概率密度函数,如图5-15 所示。该图还叠印了Rayleigh 分布的概率密度函 数理论值,可以看出,二者的吻合度比较高。 76 MATLAB/Simulink 实用教程——编程、计算与仿真 >> b=1; p=raylrnd(1,30000,1); x=0:0.1:4; % 子区间划分 y=histogram(p,x); yy=y.Values/(30000*0.1); % 直方图数据 x0=x(1:end-1)+0.05; bar(x0,yy), y=raylpdf(x,1); line(x,y) % 理论值 图5-15 Rayleigh 分布的概率密度函数及其近似 由前面介绍的频度向量f 还可以绘制出饼图,调用格式为pie(f)。用饼图可以大致显 示落入每个子区间点数的占比。 例5-13 仍考虑例5-12 中的数据。绘制饼图不能像前面例子中分那么多子区间,否则绘制的 饼图意义不大,可以将子区间分为[0, 0.5],(0, 5, 1],· · · , 图5-16 Rayleigh 分布的饼图表示 (3.5, 4],试绘制饼图表示数据的期间分布。 解可以仿照上述方法先将频度向量计算出来,然后根 据频度向量绘制出饼图,如图5-16 所示。 >> b=1; p=raylrnd(1,30000,1); x=0:0.5:4; y=histogram(p,x); f=y.Values/30000; pie(f) f1=f*100 饼图显示虽然直观,但如果不对照频度数据很难看出哪 个饼图分区对应于哪个子区间,所以同时还应该显示占比向 量f1 = [11.5, 28.2, 27.8, 19.1, 8.9, 3.4, 0.8, 0.2]%。 5.2.4 填充图 如果有一组坐标点A1(x1, y1),A2(x2, y2),· · · ,An(xn, yn),由A1~An 做折线,再由An~ A1 做折线,构成一个封闭的形状,MATLAB 提供的fill() 函数可以对封闭形状的内部进 行填充,得出填充图。该函数的调用格式为fill(x,y,c),其中c 是颜色标识,例如,可以 用'g' 表示绿色,参考表5-1,c 也可以是[1, 0, 0] 这类的三原色表示,对应红色。 如果想让Ai 这些坐标点与x 轴围成封闭图形,且x向量是从小到大或从大到小排列的 向量,则可以在左右各补充一个点,分别为(x1, 0) 和(xn, 0),使得x = [x1, x, xn],y = [0, y, 0],这样就可以由fill() 函数获得填充图形了。 例5-14 考虑例5-12 中的Rayleigh 分布,试用填充颜色的方法表示出面积达到95%的概率密 度函数曲线。 第5 章MATLAB 科学绘图77 解生成一个x ∈ (0, 4) 的行向量,得出相应的Rayleigh 概率密度函数值。本例一个关键的步 骤是得出95%面积的关键点,该点可以由逆概率密度函数raylinv() 直接求出,记为x0,由下面的 语句可以得出x0 = 2.4477。由于左端x = 0 时概率密度的值为0,所以左侧不必补充点,现在看填 充向量x1 的右侧。先从x向量中提取出x < x0 的点,横坐标再补上两次x0 的值,在这两个x0 处 相应的y 值一个是概率密度函数计算出来的函数值;另一个是0,确保围成区域是期望的区域。这 样得出的曲线如图5-17 所示。 >> x=0:0.1:4; b=1; y=raylpdf(x,b); x0=raylinv(0.95,b) ii=x<=x0; x1=[x(ii) x0, x0]; y1=[y(ii),raylpdf(x0,b),0]; plot(x,y), hold on; fill(x1,y1,'g'), hold off x0 = 2.4477 图5-17 Rayleigh 分布的概率密度与95%区域 5.2.5 对数图绘制 在一些特定的领域内,如数字信号处理与自动控制等领域,经常需要对信号与系统做频 域分析,而Bode 图分析方法就是一种常用的频域分析方法。 Bode 图是对系统G(s) 在s = jω1, jω2, · · · , jωm 处的增益的描述,其中ωk 称为频率点。 增益G(jω) 为复数向量,而复数采用的是幅值|G(jω)|与相位∠G(jω) 形式描述的。正常情 况下,Bode 图采用上下两幅图表示,分别表示幅值与频率的关系(幅频特性)、相位与频率的 关系(相频特性)。频率坐标轴采用对数型坐标轴,幅值通过20 lg |G(jω)|变换,单位为分贝 (dB),相位采用角度为单位。 如果绘制横坐标为对数坐标、纵坐标为线性坐标,则可以使用semilogx() 函数直接绘 制,而绘制纵坐标为对数坐标、横坐标为线性坐标的函数为semilogy(),两个坐标轴都是对 数坐标的图形可以利用loglog() 函数绘制。 例5-15 假设系统的传递函数如下,选择频率范围ω ∈ (0.01, 1000),试绘制幅值与频率之间 的Bode 图。 G(s) = 2(s0.4 ? 2)0.3 √s(s0.3 + 3)0.8(s0.4 ? 1)0.5 解一般情况下,频率点按照对数等间距的方式直接生成。这样,由给定的传递函数模型可以 直接计算出以分贝为单位的幅值数据,半对数幅频特性曲线如图5-18 所示。 >> G=@(s)2*(s.^0.4-2).^0.3./sqrt(s)./(s.^0.3+3).^0.8./(s.^0.4-1).^0.5; w=logspace(-2,3,100); M=20*log10(abs(G(1i*w))); % 变换为分贝 78 MATLAB/Simulink 实用教程——编程、计算与仿真 semilogx(w,M) % 绘制半对数图,横轴为对数坐标,纵轴为线性坐标 图5-18 Bode 幅频特性曲线 5.2.6 动态轨迹绘制与动画制作 前面所介绍的都是静态曲线的绘制方法。如果将一条曲线看作一个粒子的运动轨迹,则 用前面介绍的方法只能显示运动的最终结果,并不能看出粒子是如何运动的。如果将普通的 曲线绘制函数plot() 替换成comet(),则可以动态地显示粒子的运动轨迹。 例5-16 试动态显示例5-1 中粒子的运动轨迹。 解选择步距为0.001,则可以由下面语句直接动态地显示粒子的运动轨迹。 >> x=-pi:0.001:pi; y=sin(tan(x))-tan(sin(x)); comet(x,y) 前面所给出的绘图命令似乎可以直接绘制出所需的图形。不过可以考虑这样一种场景: 如果一个绘图命令后面紧接一组耗时计算命令,则由MATLAB 现有的执行机制看,图形绘 制往往不能立即进行,需要在计算完成后才能将图形绘制出来,这样的执行机制不利于动画 的处理。MATLAB 提供了drawnow 命令,强行暂缓后面的计算命令,直接完成图形绘制任务 后再开始后续命令,利用这样的方法可以实现动画的处理。 动画处理的另一个关键是如何更新图形的数据点位置,让其动起来。如果调用plot() 函数返回句柄,则曲线对象的数据存储在XData 和YData 属性中,可以更新图形数据点的位 置信息,实现动画的效果。下面通过例子演示动画处理方法。 例5-17 考虑Brown 运动的一群粒子,粒子个数n = 30,观察的区域[?30, 30],每个粒子的位 置可以由xi+1,k = xi,k + σΔxi,k,yi+1,k = yi,k + σΔyi,k,k = 1, 2, · · · , n,其中,σ 为比例因子, 增量Δxi,k 和Δyi,k 满足标准正态分布。试用动画的方法模拟粒子的Brown 运动。 解标准正态分布的随机数可以由randn() 函数直接生成,取比例因子σ = 0.3,则可以在死 循环下由向量化方法模拟。由于这里用到了死循环,所以可按Ctrl+C 组合键强行终止程序的运行。 >> n=30; x=randn(1,n); y=randn(1,n); s=0.3; % 生成伪随机数 figure(gcf), hold off; % 当前的窗口提前,若没有当前窗口则打开新窗口 h=plot(x,y,'o'); axis([-30,30,-30,30]) % 固定坐标系范围 while (1) % 死循环结构模拟动画 x=x+s*randn(1,n); y=y+s*randn(1,n); % 计算粒子新位置 第5 章MATLAB 科学绘图79 h.XData=x; h.YData=y; drawnow % 改写粒子位置并立即更新 end 5.2.7 图形窗口的分割 在实际应用中,还可以根据需要将MATLAB 的图形窗口划分为若干个区域,在每个区 域内绘制出不同的图形。本节将介绍规范的分区方法与不规则的分区方法,并通过例子演示 这样的方法及其应用。 规范分区就是将整个图形窗口分割为m × n个均匀的分区,以便在每个分区绘制出不 同的图形。规范分割方法在实际应用中是很常用的。MATLAB 提供的subplot() 函数可以 直接用于图形窗口的分割,其调用格式为subplot(m,n,k),其中k 是需要绘图的分区编号, 该编号是按行计算的。该函数还可以带一个返回变量h=subplot(m,n,k),其中h为该分 区坐标系的句柄。 例5-18 试绘制例5-15 中分数阶系统G(s) 的Bode 图。 解Bode 图将图形窗口分为上下两个部分,所以比较适合使用subplot() 函数分割分割 成2 × 1 的区域,上面的区域是1,下面的是2。分割完成之后就可以在两个部分分别绘制幅频特性 与相频特性曲线了。幅频特性可以完全使用例5-15的代码,相频特性需要重新计算,得出的结果如 图5-19 所示。 >> G=@(s)2*(s.^0.4-2).^0.3./sqrt(s)./(s.^0.3+3).^0.8./(s.^0.4-1).^0.5; w=logspace(-2,3,100); subplot(211) % 或subplot(2,1,1) G0=G(1i*w); M=20*log10(abs(G0)); semilogx(w,M) % 幅频特性 subplot(212), P=angle(G0)*180/pi; semilogx(w,P) % 相频特性 图5-19 分数阶系统的Bode 图 不规则分区是指利用图形窗口的“插入→坐标区”菜单项,允许用户用鼠标在图形窗 口中画出任意的坐标系,这样,用户就可以在新坐标系绘制图形。 5.3 MATLAB三维绘图 有些数学函数和数据的三维图可以由三维坐标系下的曲线表示,有些需要用三维曲面 表示,这完全取决于数学函数与数据的形式和意义。本节将侧重介绍三维曲线的表示方法。 80 MATLAB/Simulink 实用教程——编程、计算与仿真 5.3.1 三维曲线绘制 考虑一个在三维空间运动的质点,如果这个质点在t 时刻的空间位置由参数方程x(t)、 y(t)、z(t) 表示,则这个质点的轨迹就可以看成一条三维曲线。 MATLAB 中的二维曲线绘制函数plot() 可以扩展到三维曲线的绘制中。这时可以用 plot3() 函数绘制三维曲线。该函数的调用格式为 plot3(x,y,z) plot3(x1,y1,z1, 选项1,x2,y2,z2, 选项2,· · · ,xm,ym,zm, 选项m) 其中“选项”和二维曲线绘制的完全一致,如表5-1 所示。x、y、z 为时刻t 的空间质点的坐标 构成的向量。 相应地,类似于二维曲线绘制函数,MATLAB 还提供了其他的三维曲线绘制函数,如 stem3() 可以绘制三维火柴杆型曲线,fill3() 可以绘制三维的填充图形,bar3() 可以绘制 三维的直方图等。如果采用comet3() 函数将得出动态的轨迹显示。这些函数的调用格式可 以参见其二维曲线绘制函数原型。 例5-19 试绘制参数方程x(t)=t3e?t sin 3t,y(t)=t3e?t cos 3t,z=t2 的三维曲线。 解若想绘制该参数方程的曲线,可以先定义一个时间向量t,由其计算出x、y、z 向量,并用函 数plot3() 绘制出三维曲线,如图5-20 所示。注意,这里应该采用点运算。 >> t=0:0.01:2*pi; % 构造时间向量,注意下面的点运算 x=t.^3.*exp(-t).*sin(3*t); y=t.^3.*exp(-t).*cos(3*t); z=t.^2; plot3(x,y,z), grid % 三维曲线绘制,并绘制坐标系网格 图5-20 三维曲线的绘制 5.3.2 三维参数方程的曲线绘制 如果已知三维函数的参数方程x(t)、y(t)、z(t),还可以使用fplot3() 函数直接绘制三 维函数的曲线,该函数的调用格式为 fplot3(fx,fy,fz), fplot3(fx,fy,fz,[tm,tM]) 其中fx、fy 和fz 为参数方程的数学表示形式,可以是符号表达式,也可以是匿名函数表达 式。参数t 的默认区间为[0, 5]。 第5 章MATLAB 科学绘图81 例5-20 重新考虑例5-19 的空间质点函数,试由数学公式直接绘制三维曲线。 解将参数方程用符号表达式表示,则可以给出如下的命令,这样,由fplot3() 函数可以得出 与例5-19 完全一致的结果。 >> syms t; x=t^3*exp(-t)*sin(3*t); y=t^3*exp(-t)*cos(3*t); z=t^2; fplot3(x,y,z,[0,2*pi]) 参数方程还可以由匿名函数表示,由下面的语句可以绘制出同样的三维曲线。 >> x=@(t)t.^3.*exp(-t).*sin(3*t); y=@(t)t.^3.*exp(-t).*cos(3*t); z=@(t)t.^2; fplot3(x,y,z,[0,2*pi]); 5.3.3 三维曲面绘制 如果已知二元函数z = f(x, y),则可以考虑先在xy 平面生成一些网格点,然后求出每 个点处的函数值z,这样就可以由这些信息绘制三维图了。 例3-13 演示了MATLAB 提供的meshgrid() 函数生成网格的方式与矩阵的生成格式, 由该函数可以生成两个矩阵x与y,将两个矩阵重叠在一起正好形成了每个网格点的x 与y 坐标值,这时如果函数z = f(x, y) 已知,则可以直接用点运算的方式计算出每个网格点的函 数值矩阵z。有了这3 个矩阵,就可以调用MATLAB 提供的mesh() 与surf() 函数直接绘制 三维数据的网格图与表面图了。这两个函数的调用格式为mesh(x,y,z) 或surf(x,y,z)。 surf() 函数还可以返回曲面的句柄,这样就可以对得出的曲面进行进一步的操作处理。 例5-21 给出二元函数z = f(x, y) = (x2 ? 2x)e?x2?y2?xy,其中,?3 ? x ? 2,?2 ? y ? 2, 试绘制该函数的三维表面图形。 解选择步距为0.1,可以调用meshgrid() 函数生成xy 平面的网格矩阵x、y。然后由给出的 公式计算出曲面的z 矩阵。最后调用surf() 函数绘制曲面的三维表面图,如图5-21 所示。 >> [x,y]=meshgrid(-3:0.1:2,-2:0.1:2); % 生成xy 平面的网格矩阵x、y z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 计算高度矩阵z surf(x,y,z) % 绘制三维表面图 图5-21 三维函数的表面图 如果二元函数的解析表达式已知,还可以采用fsurf() 直接绘制函数的曲面,该函数的 调用格式为fsurf(f) 或fsurf(f,[xm,xM,ym,yM]),其中,f 为二元函数的符号表达式或 匿名函数句柄,函数绘制的默认区域为?5 ? x, y ? 5。 82 MATLAB/Simulink 实用教程——编程、计算与仿真 例5-22 试用fsurf() 函数绘制例5-21 函数的三维曲面。 解可以用符号表达式描述给定的函数,然后调用fsurf() 函数绘制函数的表面图,其结果与 图5-21 中给出的结果完全一致。 >> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y); % 给出符号表达式 fsurf(f,[-3,2,-2,2]) % 在指定区域内绘制曲面图 当然,还可以由下面的匿名函数表示给定的函数,替代上面语句的f。由该函数调用上述的 fsurf() 函数,也能得出完全一致的结果。 >> f=@(x,y)(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); fsurf(f,[-3,2,-2,2]) 5.3.4 视角设置 单击MATLAB 图形窗口的坐标系工具栏(图5-10(c))中的图标,则可以用拖动鼠标 的方法直接修改三维图的视角,用可视的方法直接调整为期望的视角。也可以调用view() 函数进行设置。 MATLAB 三维图形视角的定义如图5-22 所示,视 α y x z β 视点 图5-22 视角定义示意图 角是由两个角度唯一描述的,这两个角度分别为方位 角与仰角,其中,方位角α定义为视点与原点连线在 xy 平面投影线与y 轴负方向之间的夹角,默认值为α= ?37.5?,仰角β 定义为视点与原点连线和xy 平面的夹 角,默认值为β = 30?。MATLAB 中的当前的视角可以 由[α,β]=view(3) 语句读出,如果想改变视角观察曲 面,则可以给出view(α, β) 命令。 在工程制图领域经常需要绘制物体的三视图,其 中,俯视图是从上往下看,显然,这时的仰角为90?,方位角为0?,所以,可以由命令 view(0,90) 直接设置视角;类似地,还可以分别由view(0,0) 和view(90,0) 命令定义主 视图和侧左视图。这样由surf() 等命令绘制函数曲面后,给出修改视角的命令,就可以得出 其三视图。下面通过例子演示三视图的绘制方法。 例5-23 仍考虑例5-21 中的函数表达式,试绘制其曲面的三视图。 解可以将整个图形窗口分割成2 × 2 的四个分区,这样就可以在不同的命令在相应的分区绘 制三视图,如图5-23 所示。 >> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y); subplot(221), fsurf(f,[-3,2,-2,2]), view(0,90) % 俯视图 subplot(222), fsurf(f,[-3,2,-2,2]), view(-90,0) % 右侧视图 subplot(223), fsurf(f,[-3,2,-2,2]), view(0,0) % 主视图 subplot(224), fsurf(f,[-3,2,-2,2]) % 三维表面图 第5 章MATLAB 科学绘图83 (a)俯视图(b)右侧视图 (c)主视图(d)三维表面图 图5-23 三维曲面的三视图 5.3.5 二元参数方程的曲面绘制 前面介绍了带有单个自变量的参数方程,该参数方程对应的是三维曲线。如果参数方程 带有两个自变量u、v,三维参数方程的数学形式为 x = fx(u, v), y = fy(u, v), z = fz(u, v) (5-2) 且um ?u?uM,vm ?v?vM,则由fsurf(fx,fy,fz,[um,uM,vm,vM]) 函数可以直接绘制 三维表面图,其中,u、v 变量的默认区间为(?5, 5)。使用早期版本的MATLAB 还可以尝试 ezsurf() 函数绘制表面图,不过不建议使用该函数。 例5-24 著名的M?bius 带可以由数学模型x=cos u+v cos u cos u/2,y=sin u+v sin u cos u/2, z = v sin u/2 描述。如果0 ? u ? 2π,?0.5 ? v ? 0.5,试绘制M?bius 带的三维表面图。 解首先需要声明两个符号变量u、v,并将参数方程输入MATLAB 环境中,这样就可以由下 面的语句直接绘制M?bius 带,得出如图5-24 所示的表面图。 >> syms u v; x=cos(u)+v*cos(u)*cos(u/2); y=sin(u)+v*sin(u)*cos(u/2); z=v*sin(u/2); fsurf(x,y,z,[0,2*pi,-0.5,0.5]) % M?bius 带的绘制 图5-24 M?bius 带的表面图(图形经过了旋转处理) 84 MATLAB/Simulink 实用教程——编程、计算与仿真 5.3.6 三维动画的制作与播放 若某个三维曲面图是时间的函数,则由getframe 命令提取每个时间样本点绘制的三维 曲面句柄,这样就可以提取一系列句柄。有了句柄就可以调用movie() 函数制作三维动画的 视频。本节将通过例子演示三维动画的制作与播放方法。 例5-25 考虑一个时变的函数z(x, y, t) = sin(x2t + y2),其中,0 ? t ? 1,?2 ? x, y ? 2,试用 动画的方式表示函数表面图随时间t 的变化效果。 解三维动画的处理分为两个部分:第一部分是动画的制作,需要计算出各个时刻的曲面图数 据,由每个时刻的数据绘制三维表面图,然后用getframe() 函数提取出一帧图像的句柄,通过这 样的方法可以获得一系列句柄。为使得动画变化平稳,可以考虑用axis() 函数将每帧动画固定在 相同的坐标系范围内。第二部分是动画播放,有了动画的一系列句柄,则可以调用movie() 函数播 放动画。前面的叙述可以由下面的语句直接实现,获得期望的三维动画的结果。 >> t=linspace(0,1); [x,y]=meshgrid(-2:0.1:2); for i=1:length(t) % 用循环的方式对每个时刻单独处理 z=sin(x.^2*t(i)+y.^2); surf(x,y,z); % 绘制三维表面图 axis([-2,2,-2,2,-1,1]); h(i)=getframe; % 提取一帧图像 end figure, movie(h) % 三维动画的直接播放 值得指出的是,getframe() 函数不仅可以提取三维的帧,如果绘制二维曲线图也可 以用该函数提取一帧,所以可以用该函数将二维图形制作成动画的形式显示出来。此外, VideoWriter() 函数可以打开一个视频文件,writeVideo() 函数可以往视频文件中写入一 帧视频,下面将通过例子演示动画视频的制作过程。 例5-26 试将例5-17 中的Brown 二维运动动画制作成视频文件。 解假设运动200 步,则由下面的命令生成动画数据,并制作动画的视频文件,这些语句调用 结束后将在当前文件夹下生成brown.avi 文件,可以用任意媒体播放器播放。 >> n=30; x=randn(1,n); y=randn(1,n); s=0.3; % 生成初始位置 figure(gcf), hold off; % 当前的窗口提前,若没有当前窗口则打开新窗口 h=plot(x,y,'o'); axis([-30,30,-30,30]) vid=VideoWriter('brown.avi'); open(vid); % 打开空白的视频文件 for k=1:200 % 移动200 步的模拟动画 x=x+s*randn(1,n); y=y+s*randn(1,n); % 更新粒子位置 h.XData=x; h.YData=y; drawnow % 生成并绘制一帧图片 hVid=getframe; writeVideo(vid,hVid); % 获得并写入一帧视频 end, close(vid) % 关闭完成视频文件 5.4 隐函数绘制 前面介绍的都是显函数的曲线与曲面绘制,隐函数用前面介绍的方法是不能直接绘制 的,必须使用专门的隐函数图形绘制方法。本节将介绍二元与三元隐函数的图形绘制方法。 第5 章MATLAB 科学绘图85 5.4.1 二维隐函数曲线绘制 隐函数即满足f(x, y) = 0 方程的x 和y 之间的关系式。用前面介绍的曲线绘制方法显 然会有问题。例如,很多隐函数无法求出x 和y 之间的显式关系,所以无法先定义一个x向 量再求出相应的y 向量,从而不能采用plot() 函数绘制曲线。另外,即使能求出x、y 之间 的显式关系,但由于不是单值函数,因此绘制也是很麻烦的。前面介绍的显函数绘制命令 fplot() 也不能绘制隐函数曲线。 MATLAB 提供的fimplicit() 函数可以直接绘制隐函数曲线,该函数的调用格式为 fimplicit(隐函数表达式),其中,“隐函数表达式”既可以是符号表达式,也可以是点运算 描述的匿名函数。用户还可以指定绘图范围fimplicit(隐函数表达式,[xm,xM]),得出可 读性更好的曲线。坐标轴范围的默认区间为[?5, 5]。 早期版本的MATLAB 还提供了绘制二元隐函数曲线的实用函数ezplot(),其调用格 式与fimplicit() 函数接近,还可以用字符串描述隐函数方程。ezplot() 函数不能处理由 piecewise() 语句描述的分段函数模型。下面将通过例子演示该函数的使用方法。 例5-27 试绘制隐函数y2 cos(x + y2) + x2ex+y = 0 在(?2π, 2π) 的曲线。 解从给出的函数可见,无法用解析的方法写出该函数,所以不能用前面给出的plot() 函数 绘制出该函数的曲线。可以给出如下的MATLAB 命令,绘制出如图5-25(a)所示的隐函数曲线。可 见,隐函数绘制是很简单的,只须将隐函数原原本本地表示出来,就能直接得出相应的曲线。 >> syms x y; f=y^2*cos(x+y^2)+x^2*exp(x+y); % 符号表达式描述 fimplicit(f,[-2*pi,2*pi]) % 在指定的范围内绘制隐函数曲线 (a)默认参数直接绘制(b)光滑的隐函数曲线 图5-25 隐函数曲线绘制 从得出的曲线看,默认设置下得到的曲线不平滑,有些地方有毛刺,故可以修改MeshDensity 参数,将其设置为1000,绘制出如图5-25(b)所示的光滑隐函数曲线。 >> fimplicit(f,[-2*pi,2*pi],'MeshDensity',1000) % 绘制光滑的隐函数曲线 还可以使用下面的匿名函数形式描述隐函数,注意,应该使用点运算描述向量运算。描述了函 数之后,绘制命令fimplicit() 与前面介绍的是完全一致的。 >> f=@(x,y)y.^2.*cos(x+y.^2)+x.^2.*exp(x+y); % 匿名函数,点运算 86 MATLAB/Simulink 实用教程——编程、计算与仿真 例5-28 试用隐函数绘制的方法求解联立方程组在?2π?x, y?2π 范围内的解。 (x2e?xy2/2 + e?x/2 sin(xy) = 0 y2 cos(y + x2) + x2ex+y = 0 解上面的每个方程都可以看作一个隐函数,这样就可以用fimplicit() 函数同时绘制这两 个隐函数,得出如图5-26 所示的隐函数曲线。这样,两组曲线的每个交点都是联立方程的解(这是 习题1.8 的解答。) >> syms x y; f1=x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y); f2=y^2*cos(y+x^2)+x^2*exp(x+y); fimplicit([f1,f2],[-2*pi,2*pi],'MeshDensity',1000) 图5-26 联立方程图解法示意图 如果对某个交点处的解感兴趣,可以考虑采用局部放大的方法读出解的x,y 值,这种方 法称为图解法。不过,若想利用图解法逐一求出图中所有交点处的解是极其麻烦的,第9 章 将介绍多解方程的数值求解方法,一次性得出方程在指定区域内所有的解。 5.4.2 三维隐函数曲面绘制 如果某三维曲面由隐函数g(x, y, z) = 0 表示,则可以利用MATLAB 的fimplicit3() 绘制其曲面图形,该函数的调用格式为fimplicit3(fun,[xm,xM,ym,yM,zm,zM]),其中 fun 可以为匿名函数、也可以是符号表达式,坐标轴范围向量xm、xM、ym、yM、zm、zM 的默认 值为±5。如果只给出一对上下限xm、xM,则表示3 个坐标轴均同样设置。该函数的核心部分 是等高面绘制函数。 例5-29 假设某三维隐函数的数学表达式为 x(x, y, z) = x sin ??y + z2+ y2 cos ??x + z+ zx cos ??z + y2= 0 且感兴趣的区域为x, y, z ∈ (?1, 1),试绘制其三维曲面。 解用符号表达式或匿名函数的方式都可以描述原始的隐函数,二者作用相同。用下面语句就 可以直接绘制出该隐函数的三维曲面图,如图5-27(a)所示。 >> syms x y z; f=x*sin(y+z^2)+y^2*cos(x+z)+z*x*cos(z+y^2); fimplicit3(f,[-1 1]) % 三维隐函数曲面绘制 其实,三元隐函数还可以用匿名函数描述,得出的结果是一致的。 第5 章MATLAB 科学绘图87 (a)直接绘制的三维曲面(b)叠印球面的三维曲面 图5-27 隐函数三维曲面绘制 >> f=@(x,y,z)x.*sin(y+z.^2)+y.^2.*cos(x+z)+z.*x.*cos(z+y.^2); fimplicit3(f,[-1,1]) 还可以在原曲面上叠印单位球面x2 + y2 + z2 = 1,如图5-27(b)所示。 >> f1=x^2+y^2+z^2-1; fimplicit3([f f1],[-1 1]); % 叠印两个曲面 5.5 习 题 5.1 试绘制函数曲线y(x) = sin πx/(πx),其中x ∈ (?4, 4)。 5.2 选择合适的步距绘制出图形sin 1/t,其中t ∈ (?1, 1)。 5.3 选择合适的步距,绘制tan t 曲线,t ∈ (?π, π),并观察不连续点的处理方法。 5.4 试绘制下面的函数曲线: (1)f(x) = x sin x,x ∈ (?50, 50),(2)f(x) = x sin 1/x,x ∈ (?1, 1) 5.5 试选择合适的t 范围,绘制x = sin t,y = sin 2t 的曲线,如果有一个质点在该曲线上运动,试 绘制其运动的动态显示。 5.6 试在t ∈ (0, 2π) 区间内在同一坐标系下绘制3 条曲线sin x、sin 2x、sin 3x。 5.7 用MATLAB 语言的基本语句显然可以立即绘制一个正三角形。试结合循环结构,编写一个 小程序,在同一个坐标系下绘制出该正三角形绕其中心旋转后得出的一系列三角形,还可以 调整旋转步距观察效果。 5.8 试在区间?50 ? x, y ? 50 内绘制x sin x + y sin y = 0 的曲线。 5.9 试绘制分段函数的曲线 y(t) = sin t + cos t, t ? 0 tan t, t > 0 5.10 已知正态分布的概率密度函数为p(x) = 1 √2πσ e?(x?μ)2/(2σ2),其中μ为均值,σ 为方差,试 绘制不同μ、σ 参数下的概率密度函数曲线。 5.11 试按照下面的规则构造某序列的前40 项,再用stem() 函数显示序列变化趋势。 xk = 1 + 1 2 + 1 3 + · · · + 1 k ? ln k 5.12 已知迭代模型如下,且迭代初值为x0 = 0,后续各点可以递推求出y0 = 0。 ( xk+1 = 1 + yk ? 1.4x2 k yk+1 = 0.3xk 88 MATLAB/Simulink 实用教程——编程、计算与仿真 如果取迭代初值为x0 = 0,y0 = 0,那么进行30000 次迭代求出一组x和y 向量,然后在所有 的xk 和yk 坐标处点亮一个点(注意不要连线),最后绘制出所需的图形。(提示:这样绘制出 的图形又称为Hénon 引力线图,它将迭代出来的随机点吸引到一起,最后得出貌似连贯的引 力线图。) 5.13 假设某幂级数展开表达式为 f(x) = lim N→∞ NXn=1 (?1)n x2n (2n)! 若N 足够大,则幂级数f(x) 收敛为某个函数? f(x)。试写出一个MATLAB 程序,绘制出 x ∈ (0, π) 区间的? f(x) 的函数曲线,观察并验证? f(x) 是什么样的函数。 5.14 试在t ∈ (0, π) 区间内绘制sin t2 的函数曲线,若某个区域曲线关系不清晰,则可以考虑采用 局部放大的方法显示细节。 5.15 分别选取合适的θ 范围,绘制出下列极坐标图形: (1)ρ = 1.0013θ2,(2)ρ = cos 7θ/2,(3)ρ = sin θ/θ,(4)ρ = 1 ? cos3 7θ 5.16 试绘制参数方程曲线x = (1 + sin 5t/5) cos t,y = (1 + sin 5t/5) sin t,t ∈ (0, 2π)。如果将参 数方程中的5 替换成其他数值会得出什么结果。 5.17 用图解的方式求解下面联立方程的近似解: (1)(x2 + y2 = 3xy2 x3 ? x2 = y2 ? y (2)( e?(x+y)2+π/2 sin(5x + 2y) = 0 (x2 ? y2 + xy)e?x2?y2?xy = 0 5.18 已知正弦函数y = sin(ωt + 20?),t ∈ (0, 2π),ω ∈ (0.01, 10),试绘制当ω 变化时正弦函数曲 线的动画。 5.19 已知某空间质点的运动方程为x(t) = cos t + t sin t,y(t) = sin t?t cos t,z(t) = t2,且 t ∈ (0, 2π),试绘制该空间质点的运动轨迹,并绘制出该质点随时间变化的动态轨迹。 5.20 试分别绘制出xy、sin xy 和e2x/(x2+y2) 的三维表面图。 5.21 试绘制函数的三维表面图f(x, y) = sinpx2 + y2/px2 + y2,?8 ? x, y ? 8。 5.22 试绘制下列参数方程的三维表面图[19]。 (1)x = 2 sin2 u cos2 v,y = 2 sin u sin2 v,z = 2 cos u sin2 v,?π/2 ? u, v ? π/2 (2)x = u ? u3/3 + uv2,y = v ? v3/3 + vu2,z = u2 ? v2,?2 ? u, v ? 2 5.23 试绘制下面函数的表面图,还可以使用waterfall()、surfc()、surfl() 等函数并观察效果。 (1)z = xy,(2)z = sin x2y3,(3)z = (x ? 1)2y2 (x ? 1)2 + y2 ,(4)z = ?xy e?2(x2+y2) 5.24 已知x = u sin t,y = u cos t,z = t/3,t ∈ (0, 15),u ∈ (?1, 1),试绘制参数方程的曲面。 5.25 在图形绘制语句中,若函数值为不定式NaN,则相应的部分不绘制出来。试利用该规律绘制 z = sin xy 的表面图,并剪切x2 + y2 ? 0.52 的部分。 5.26 试绘制x(z, y) = (z2?2z)e?z2?y2?zy 的三维曲面。 5.27 若已知x = cos t(3 + cos u),sin t(3 + cos u),z = sin u,且t ∈ (0, 2π),u ∈ (0, 2π),试绘制表 面图。 第5 章MATLAB 科学绘图89 5.28 已知二元函数f(x, y) = x sin(1/y) + y sin(1/x),试用图形方法研究(x, y) 在(0, 0) 点附近的 行为。 5.29 试绘制复杂隐函数曲线(r?3)√r+0.75+sin 8√r cos 6θ?0.75 sin 5θ = 0,其中r = x2+y2, θ = arctan(y/|x|)。 5.30 试绘制出三维隐函数(x2 + xy + xz)e?z + z2yx + sin(x + y + z2) = 0 的曲面。 5.31 试绘制两个曲面x2 + y2 + z2 = 64,y + z = 0 并观察其交线。 5.32 MATLAB 提供的treeplot() 函数可以绘制树图。例如,图4-7 给出的二叉树图可以由下面 语句直接绘制。 >> nodes=[0,1,1,2,2,3,3,4,4,5,5,6,6,8,8]; treeplot(nodes) 试理解上面语句的实际含义,并扩展该图形,绘制Fibonacci 序列的7 级二叉树结构,并编写 函数绘制k 级二叉树图。