15.1 What is a chebfun2v?

Chebfun2 objects represent vector-valued functions. We use a lower case letter like $f$ for a chebfun2 and an upper case letter like $F$ for a chebfun2v.

Chebfun2 represents a vector-valued function $F(x,y) = (f(x,y);g(x,y))$ by approximating each component by a low rank approximant, as described in Section 12.8. There are two ways to form a chebfun2v: either by explicitly calling the constructor, or by vertical concatenation of two chebfun2 objects. Here are these two alternatives:

d = [0 1 0 2];
F = chebfun2v(@(x,y) sin(x.*y), @(x,y) cos(y), d);
f = chebfun2(@(x,y) sin(x.*y), d);
g = chebfun2(@(x,y) cos(y), d);
G = [f;g]
G =
   chebfun2v object (Column vector) containing:
   chebfun2 object
       domain                 rank       corner values
[   0,   1] x [   0,   2]        7     [3.6e-23 3.2e-17 1.7e-17 0.91]
vertical scale =   1 
   chebfun2 object
       domain                 rank       corner values
[   0,   1] x [   0,   2]        1     [   1    1 -0.42 -0.42]
vertical scale =   1 

Note that displaying a chebfun2v shows that it is a vector of two chebfun2 objects.

15.2 Algebraic operations

Chebfun2 objects are useful for performing 2D vector calculus. The basic algebraic operations are scalar multiplication, vector addition, dot product and cross product.

Scalar multiplication is the product of a scalar function with a vector function:

f = chebfun2(@(x,y) exp(-(x.*y).^2/20), d);
f*F
ans =
   chebfun2v object (Column vector) containing:
   chebfun2 object
       domain                 rank       corner values
[   0,   1] x [   0,   2]        8     [-9.2e-22 1.7e-16 1.7e-16 0.74]
vertical scale = 0.89 
   chebfun2 object
       domain                 rank       corner values
[   0,   1] x [   0,   2]        7     [   1    1 -0.42 -0.34]
vertical scale =   1 

Vector addition yields another chebfun2v and satisfies the parallelogram law:

plaw = abs((2*norm(F)^2 + 2*norm(G)^2) - (norm(F+G)^2 + norm(F-G)^2));
fprintf('Parallelogram law holds with error = %10.5e\n',plaw)
Parallelogram law holds with error = 6.21725e-15

The dot product combines two vector functions into a scalar function. If the dot product of two chebfun2v objects takes the value zero at some $(x,y)$, then the vector-valued functions are orthogonal there. For example, the following code segment determines a curve along which two vector-valued functions are orthogonal:

LW = 'linewidth';
F = chebfun2v(@(x,y) sin(x.*y), @(x,y) cos(y),d);
G = chebfun2v(@(x,y) cos(4*x.*y), @(x,y) x + x.*y.^2,d);
plot(roots(dot(F,G)),LW,1.6), axis equal, axis(d)

The cross product for 2D vector fields works as follows.

help chebfun2v/cross
 CROSS   Vector cross product.
    CROSS(F, G) returns the cross product of the CHEBFUN2V objects F and G. If F
    and G both have two components, then it returns the CHEBFUN2 representing
        CROSS(F,G) = F(1) * G(2) - F(2) * G(1)
    where F = (F(1); F(2)) and G = (G(1); G(2)). If F and G have three
    components then it returns the CHEBFUN2V representing the 3D cross 
    product.

15.3 Differential operators

Vector calculus also involves various differential operators defined on scalar- or vector-valued functions such as gradient, curl, divergence, and Laplacian.

The gradient of a chebfun2 representing a scalar function $f(x,y)$ represents, geometrically, the direction and magnitude of steepest ascent of $f$. If the gradient of $f$ is $0$ at $(x,y)$, then $f$ has a critical point at $(x,y)$. Here are the critical points of a sum of Gaussian bumps:

f = chebfun2(0);
rng('default')
for k = 1:10
    x0 = 2*rand-1; y0 = 2*rand-1;
    f = f + chebfun2(@(x,y) exp(-10*((x-x0).^2 + (y-y0).^2)));
end
plot(f), hold on
r = roots(gradient(f));
plot3(r(:,1),r(:,2),f(r(:,1),r(:,2)),'k.','markersize',20)
zlim([0 4]), hold off, colormap(pink)

The curl of 2D vector function is a scalar function defined as follows.

help chebfun2v/curl
 CURL  curl of a CHEBFUN2V
    S = CURL(F) returns the CHEBFUN2 of the curl of F. If F is a CHEBFUN2V with
    two components then it returns the CHEBFUN2 representing
          CURL(F) = F(2)_x - F(1)_y,
    where F = (F(1),F(2)).  If F is a CHEBFUN2V with three components then it
    returns the CHEBFUN2V representing the 3D curl operation.

If the chebfun2v $F$ describes a vector velocity field of fluid flow, for example, then curl(F) is the scalar function equal to twice the angular speed of a particle in the flow at each point. A particle moving in a gradient field has zero angular speed and hence, the curl of the gradient is zero. We can check this numerically:

norm(curl(gradient(f)))
ans =
     0

The divergence of a chebfun2v is defined as follows.

help chebfun2v/divergence
 DIVERGENCE   Divergence of a CHEBFUN2V.
    DIVERGENCE(F) returns the divergence of the CHEBFUN2V i.e.,
        divergence(F) = F_x + F_y.

  See also DIV.

This measures a vector field's distribution of sources or sinks. The Laplacian is closely related and is the divergence of the gradient,

norm(laplacian(f) - divergence(gradient(f)))
ans =
     0

15.4 Line integrals

Given a vector field $F$, we can compute the line integral along a curve with the command integral, defined as follows.

help chebfun2v/integral
 INTEGRAL   Line integration of a CHEBFUN2V.

    INTEGRAL(F, C) computes the line integral of F along the curve C, i.e.,

                            /
        INTEGRAL(F, C) =   |  < F(r), dr > 
                           /
                          C 

    where the curve C is parameterised by the complex curve r(t).  

The gradient theorem says that if $F$ is a gradient field, then the line integral along a smooth curve only depends on the end points of that curve. We can check this numerically:

f = chebfun2(@(x,y) cos(10*x.*y.^2) + exp(-x.^2));  % chebfun2
F = gradient(f);                                    % gradient (chebfun2v)
C = chebfun(@(t) t.*exp(10i*t),[0 1]);              % spiral curve
v = integral(F,C);ends = f(cos(10),sin(10))-f(0,0); % line integral
abs(v-ends)                                         % gradient theorem
ans =
     0

15.5 Phase diagram

A phase diagram is a graphical representation of a system of trajectories for a two-variable autonomous dynamical system. Chebfun2 plots phase diagrams with quiver command, which has been overloaded to plot the vector field. Note that there is a potential terminological ambiguity in that a "phase portrait" can also refer to a portrait of a complex-valued function (see section 12.7).

In addition, Chebfun2 makes it easy to compute and plot individual trajectories of a vector field. If $F$ is a chebfun2v, then ode45(F,tspan,y0) computes a trajectory of the autonomous system $dx/dt=f(x,y)$, $dy/dt=g(x,y)$, where $f$ and $g$ are the first and second components of $F$. Given a prescribed time interval and initial conditions, this command returns a complex-valued chebfun representing the trajectory in the form $x(t) + iy(t)$. For example:

d = 0.04; a = 1; b = -.75;
F = chebfun2v(@(x,y)y, @(x,y)-d*y - b*x - a*x.^3, [-2 2 -2 2]);
[t y] = ode45(F,[0 40],[0,.5]);
plot(y,'r',LW,1.6), hold on,
quiver(F,'b'), axis equal
title('The Duffing oscillator','FontSize',14), hold off

15.6 Representing 2D parametric surfaces in 3D space

So far, we have explored chebfun2v objects with two components, but Chebfun2 can also work with functions with three components, i.e., functions from a rectangle in $R^2$ into $R^3$. For example, we can represent the unit sphere via spherical coordinates as follows:

th = chebfun2(@(th,phi) th, [0 pi 0 2*pi]);
phi = chebfun2(@(th,phi) phi, [0 pi 0 2*pi]);

x = sin(th).*cos(phi);
y = sin(th).*sin(phi);
z = cos(th);

F = [x;y;z];
surf(F), colormap('default'), camlight, axis equal

Above, we have formed a chebfun2v with three components by vertical concatenation of chebfun2 objects. However, for the familiar surfaces cylinders, spheres, and ellipsoids, Chebfun2 has overloads of the commands cylinder, sphere, and ellipsoid to make things simpler. For example, a cylinder of radius $1$ and height $5$ can be constructed like this:

h = 5;
r = chebfun(@(th) 1+0*th,[0 h]);
F = cylinder(r);
surf(F), camlight

An important class of parametric surfaces are surfaces of revolution, which are formed by revolving a curve in the left half plane about the $z$-axis. The cylinder command can be used to generate surfaces of revolution. For example:

f = chebfun(@(t) (sin(pi*t)+1.1).*t.*(t-10),[0 5]);
F = cylinder(f);
surf(F), axis([-70 70 -70 70 -2 6]), camlight

Here as another example is a torus with a gap in it.

x = chebfun2(@(x,y) x); y = chebfun2(@(x,y) y);
theta = 0.9*pi*x; phi = pi*y;
F = [-(1+.3*cos(phi)).*sin(theta);
      (1+.3*cos(phi)).*cos(theta);
      .3*sin(phi)];
surf(F), axis equal, camlight

15.7 Surface normals and the divergence theorem

Given a chebfun2v representing a surface, the normal can be computed by the Chebfun2 normal command. Here are the normal vectors of another torus:

r1 = 1; r2 = 1/3;       % inner and outer radius
d = [0 2*pi 0 2*pi];
u = chebfun2(@(u,v) u,d);
v = chebfun2(@(u,v) v,d);
F = [-(r1+r2*cos(v)).*sin(u); (r1+r2*cos(v)).*cos(u); r2*sin(v)];  % torus

surf(F), camlight, hold on
quiver3(F(1),F(2),F(3),normal(F,'unit'),'numpts',10)
axis equal, hold off

Once we have the surface normal vectors we can compute, for instance, the volume of the torus by applying the divergence theorem:

$$ \int\int_V\int \hbox{div}(G) dV = \int_S\int G\cdot d\mathbf{S}, $$

where $\hbox{div}(G)=1$. Instead of integrating over the 3D volume, which is currently not possible in Chebfun2, we integrate over the 2D surface:

G = F./3;  % full 3D divergence of G is 1 because F = [x;y;z].
integral2(dot(G,normal(F)))
exact = 2*pi^2*r1*r2.^2
ans =
   2.193245422464303
exact =
   2.193245422464302

Chebfun2v objects with three components come with a warning. Chebfun2 works with functions of two real variables and therefore, operations such as curl and divergence (in 2D) have little physical meaning to the represented 3D surface. The reason we can compute the volume of the torus (above) is because we are using the divergence theorem and circumventing the 3D divergence.

To finish this section we represent the Klein Bagel. The solid black line shows the parameterisation seam and is displayed with the syntax surf(F,'-'). See [Platte 2013] for more on parameterised surfaces.

u = chebfun2(@(u,v) u, [0 2*pi 0 2*pi]);
v = chebfun2(@(u,v) v, [0 2*pi 0 2*pi]);
x=(3+cos(u/2).*sin(v)-sin(u/2).*sin(2*v)).*cos(u);
y=(3+cos(u/2).*sin(v)-sin(u/2).*sin(2*v)).*sin(u);
z=sin(u/2).*sin(v)+cos(u/2).*sin(2*v);
surf([x;y;z],'-k','FaceAlpha',.6), camlight left, colormap(hot)
axis tight equal off

15.8 References

[Platte 2013] R. Platte, "Parameterizable surfaces," http://www.chebfun.org/examples/geom/ParametricSurfaces.html.