Matlab使用fsurf根据已知函数绘制曲面图

前面几篇讲述了Matlab如何根据已知数据，来绘制各种曲面图的文章。今天来说一下，在Matlab中如何根据已知函数，绘制相关的曲面图。Matlab中有一个fsurf函数，可以实现这种功能。本文主要讲解fsurf函数的常见用法、语法说明、表达式的三维曲面图、指定曲面图区间并绘制分段表达式、参数化曲面图、添加标题和轴标签以及格式化刻度、指定曲面属性、创建后修改曲面、显示曲面图下的等高线以及控制曲面图的分辨率等相关用法。

下面，我们首先给出Matlab中关于fsurf函数的帮助文档如下：

>> help fsurf
 fsurf   Plot 3-D surface
    fsurf(FUN) creates a surface plot of the function FUN(X,Y). FUN is plotted over
    the axes size, with a default interval of -5 < X < 5, -5 < Y < 5.
 
    fsurf(FUN,INTERVAL) plots FUN over the specified INTERVAL instead of the
    default interval.  INTERVAL can be the vector [XMIN,XMAX,YMIN,YMAX] or the
    vector [A,B] (to plot over A < X < B, A < Y < B).
 
    fsurf(FUNX,FUNY,FUNZ) plots the parametric surface FUNX(U,V),
    FUNY(U,V), and FUNZ(U,V) over the interval -5 < U < 5 and
    -5 < V < 5.
 
    fsurf(FUNX,FUNY,FUNZ,[UMIN UMAX VMIN VMAX]) or
    fsurf(FUNX,FUNY,FUNZ,[A B]) uses the specified interval.
 
    fsurf(AX,...) plots into the axes AX instead of the current axes.
 
    H = fsurf(...) returns a handle to the surface object in H.
 
    Examples:
       fsurf(@(x,y) x.*exp(-x.^2-y.^2))
       fsurf(@(x,y) besselj(1,hypot(x,y)))
       fsurf(@(x,y) besselj(1,hypot(x,y)),[-20,20]) % this can take a moment
       fsurf(@(x,y) sqrt(1-x.^2-y.^2),[-1.1,1.1])
       fsurf(@(x,y) x./y+y./x)
       fsurf(@peaks)
 
       f = @(u) 1./(1+u.^2);
       fsurf(@(u,v) u, @(u,v) f(u).*sin(v), @(u,v) f(u).*cos(v),[-2 2 -pi pi])
 
       A = 2/3;
       B = sqrt(2);
       xfcn = @(u,v) A*(cos(u).*cos(2*v) + B*sin(u).*cos(v)).*cos(u) ./ (B - sin(2*u).*sin(3*v));
       yfcn = @(u,v) A*(cos(u).*sin(2*v) - B*sin(u).*sin(v)).*cos(u) ./ (B - sin(2*u).*sin(3*v));
       zfcn = @(u,v) B*cos(u).^2 ./ (B - sin(2*u).*sin(3*v));
       h = fsurf(xfcn,yfcn,zfcn,[0 pi 0 pi]);
 
    If your function has additional parameters, for example k in myfun:
       %------------------------------%
       function z = myfun(x,y,k1,k2,k3)
       z = x.*(y.^k1)./(x.^k2 + y.^k3);
       %------------------------------%
    then you may use an anonymous function to specify that parameter:
       fsurf(@(x,y)myfun(x,y,2,2,4))

常见用法

fsurf(f)
fsurf(f,xyinterval)
fsurf(funx,funy,funz)
fsurf(funx,funy,funz,uvinterval)
fsurf(___,LineSpec)
fsurf(___,Name,Value)
fsurf(ax,___)
fs = fsurf(___)

语法说明

fsurf(f) 在默认区间 [-5 5]（对于 x 和 y）为函数 z = f(x,y) 创建曲面图。

fsurf(f,xyinterval) 将在指定区间绘图。要对 x 和 y 使用相同的区间，请将 xyinterval 指定为 [min max] 形式的二元素向量。要使用不同的区间，请指定 [xmin xmax ymin ymax] 形式的四元素向量。

fsurf(funx,funy,funz) 在默认区间 [-5 5]（对于 u 和 v）绘制由 x = funx(u,v)、y = funy(u,v)、z = funz(u,v) 定义的参数化曲面。

fsurf(funx,funy,funz,uvinterval) 将在指定区间绘图。要对 u 和 v 使用相同的区间，请将 uvinterval 指定为 [min max] 形式的二元素向量。要使用不同的区间，请指定 [umin umax vmin vmax] 形式的四元素向量。

fsurf(_,LineSpec) 设置线型、标记符号和曲面颜色。例如，’-r’ 指定红色线条。在前面的任何输入参数组合之后使用此选项。

fsurf(_,Name,Value) 使用一个或多个名称-值对组参数指定曲面属性。在前面语法中的任何输入参数组合后使用此选项。

fsurf(ax,_) 将图形绘制到 ax 指定的坐标区中，而不是当前坐标区 (gca) 中。

fs = fsurf(_) 返回 FunctionSurface 对象或 ParameterizedFunctionSurface 对象，具体情况取决于输入。使用 fs 来查询和修改特定曲面的属性。

表达式的三维曲面图

在默认区间 −5<x<5 和 −5<y<5 绘制表达式 sin(x)+cos(y)。

fsurf(@(x,y) sin(x)+cos(y))

指定曲面图区间并绘制分段表达式

绘制分段表达式

• erf(x)+cos(y) −5<x<0
• sin(x)+cos(y) 0<x<5

区间为 −5<y<5.

指定绘图区间为 fsurf 的第二个输入参数。在相同坐标区的不同区间绘制多个曲面时，坐标轴范围会调整以包括所有数据。

f1 = @(x,y) erf(x)+cos(y);
fsurf(f1,[-5 0 -5 5])
hold on
f2 = @(x,y) sin(x)+cos(y);
fsurf(f2,[0 5 -5 5])
hold off

参数化曲面图

绘制参数化曲面

• x=rcos(u)sin(v)
• y=rsin(u)sin(v)
• z=rcos(v)
• where r=2+sin(7u+5v)

对于 0<u<2π，0<v<π。使用 camlight 为曲面添加光。

r = @(u,v) 2 + sin(7.*u + 5.*v);
funx = @(u,v) r(u,v).*cos(u).*sin(v);
funy = @(u,v) r(u,v).*sin(u).*sin(v);
funz = @(u,v) r(u,v).*cos(v);
fsurf(funx,funy,funz,[0 2*pi 0 pi]) 
camlight

添加标题和轴标签以及格式化刻度

在从 −2π 到 2π 的区间为 x 和 y 绘制三维曲面 ysin(x)−xcos(y)。添加标题和轴标签，并显示轴轮廓。

fsurf(@(x,y) y.*sin(x)-x.*cos(y),[-2*pi 2*pi])
title('ysin(x) - xcos(y) for x and y in [-2\pi,2\pi]')
xlabel('x');
ylabel('y');
zlabel('z');
box on

使用坐标区对象的 XTickLabel 和 XTick 属性设置 x 轴刻度值和关联的标签。使用 gca 访问坐标区对象。按照同样的方式设置 y 轴刻度值和关联的标签。

ax = gca;
ax.XTick = -2*pi:pi/2:2*pi;
ax.XTickLabel = {'-2\pi','-3\pi/2','-\pi','-\pi/2','0','\pi/2','\pi','3\pi/2','2\pi'};

ax.YTick = -2*pi:pi/2:2*pi;
ax.YTickLabel = {'-2\pi','-3\pi/2','-\pi','-\pi/2','0','\pi/2','\pi','3\pi/2','2\pi'};

指定曲面属性

使用不同的线型为不同的 v 值绘制参数化曲面 x=usin(v)、y=−ucos(v) 和 z=v。对于 −5<v<−2，使用绿色虚线绘制曲面边。对于 −2<v<2，通过将 EdgeColor 属性设置为 ‘none’ 来关闭边。

funx = @(u,v) u.*sin(v);
funy = @(u,v) -u.*cos(v);
funz = @(u,v) v;

fsurf(funx,funy,funz,[-5 5 -5 -2],'--','EdgeColor','g')
hold on
fsurf(funx,funy,funz,[-5 5 -2 2],'EdgeColor','none')
hold off

创建后修改曲面

绘制参数化曲面

• x=e^−|u|/10sin(5|v|)
• y=e^−|u|/10cos(5|v|)
• z=u.

将参数化函数曲面对象指定给变量。

x = @(u,v) exp(-abs(u)/10).*sin(5*abs(v));
y = @(u,v) exp(-abs(u)/10).*cos(5*abs(v));
z = @(u,v) u;
fs = fsurf(x,y,z)

输出结果如下：

fs = 

  ParameterizedFunctionSurface (具有属性):

    XFunction: @(u,v)exp(-abs(u)/10).*sin(5*abs(v))
    YFunction: @(u,v)exp(-abs(u)/10).*cos(5*abs(v))
    ZFunction: @(u,v)u
    EdgeColor: [0 0 0]
    LineStyle: '-'
    FaceColor: 'interp'

通过设置对象的 URange 属性，将 u 的绘图区间更改为 [-30 30]。通过将 FaceAlpha 属性设置为 0（透明）到 1（不透明）之间的值，为曲面添加透明度。

fs.URange = [-30 30];

fs.FaceAlpha = .5;

显示曲面图下的等高线

通过将 ‘ShowContours’ 选项设置为 ‘on’，显示曲面图下的等高线。

f = @(x,y) 3*(1-x).^2.*exp(-(x.^2)-(y+1).^2)...
    - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2)...
    - 1/3*exp(-(x+1).^2 - y.^2);
fsurf(f,[-3 3],'ShowContours','on')

控制曲面图的分辨率

使用 ‘MeshDensity’ 选项控制曲面图的分辨率。增大 ‘MeshDensity’ 可以使绘图更平滑、更准确，而减小该值可以提高绘图速度。

在一个分块图布局中创建两个绘图。在第一个绘图中，显示参数化曲面 x=sin(s)、y=cos(s) 和 z=(t/10)sin(1/s)。曲面存在很大间隙。要解决此问题，可在第二个绘图中将 ‘MeshDensity’ 增大到 40。fsurf 将填补间隙，表明增大 ‘MeshDensity’ 可以提高分辨率。

tiledlayout(2,1)
nexttile
fsurf(@(s,t) sin(s), @(s,t) cos(s), @(s,t) t/10.*sin(1./s))
view(-172,25)
title('Default MeshDensity = 35')

nexttile
fsurf(@(s,t) sin(s), @(s,t) cos(s), @(s,t) t/10.*sin(1./s),'MeshDensity',40)
view(-172,25)
title('Increased MeshDensity = 40')

由于我用的是Matlab2016，还没有tiledlayout以及nexttile函数（Matlab2019才有这函数），因此我用subplot来代替，结果是一致的。代码如下：

subplot(2,1,1)
fsurf(@(s,t) sin(s), @(s,t) cos(s), @(s,t) t/10.*sin(1./s))
view(-172,25)
title('Default MeshDensity = 35')

subplot(2,1,2)
fsurf(@(s,t) sin(s), @(s,t) cos(s), @(s,t) t/10.*sin(1./s),'MeshDensity',40)
view(-172,25)
title('Increased MeshDensity = 40')