### 12.1 What is a chebfun2?

Chebfun2 is the part of Chebfun that deals with functions of two variables defined on a rectangle $[a,b]\times[c,d]$. Just like Chebfun in 1D, it is an extremely convenient tool for all kinds of computations including algebraic manipulation, optimization, integration, and rootfinding. It also extends to vector-valued functions of two variables, so that one can perform vector calculus.

For example, here is a test function that has been part of MATLAB for many years. MATLAB represents the "peaks" function by a $49\times 49$ matrix:

peaks

z =  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)



The same function is available as a chebfun2 in the Chebfun2 gallery:

f = cheb.gallery2('peaks');
plot(f), axis tight, title('Chebfun2 Peaks')

In Chebfun we can do all sorts of things with functions to a high accuracy, such as evaluate them

f(0.5,0.5)
ans =
0.375375578848315


or compute their maxima,

max2(f)
ans =
8.106213589442337


A chebfun2, with a lower-case "c", is a MATLAB object, the 2D analogue of a chebfun. The syntax for chebfun2 objects is similar to the syntax for matrices in MATLAB, and Chebfun2 objects have many MATLAB commands overloaded. For instance, trace(A) returns the sum of the diagonal entries of a matrix $A$ and trace(f) returns the integral of $f(x,x)$ when $f$ is a chebfun2.

Chebfun2 builds on Chebfun's univariate representations and algorithms. Algorithmic details are given in [Townsend & Trefethen 2013b] and mathematical underpinnings in [Townsend & Trefethen 2014]. For more information, see Section 12.8.

### 12.2 What is a chebfun2v?

Chebfun2 can represent scalar-valued functions, such as $\exp(x+y)$, and vector-valued functions, such as $[\exp(x+y);\cos(x-y)]$. A vector-valued function is called a chebfun2v, and chebfun2v objects are useful for computations of vector calculus. For information about chebfun2v objects and vector calculus, see Chapters 15 and 16 of this guide.

### 12.3 Constructing chebfun2 objects

A chebfun2 can be constructed by supplying the Chebfun2 constructor with a bivariate function handle or string. The default rectangular domain is $[-1,1]\times [-1,1]$. (An example showing how to specify a different domain is given in section 6.) For example, here we construct and plot a chebfun2 representing $\cos(2\pi xy)$ on $[-1,1]\times[-1,1]$.

f = chebfun2(@(x,y) cos(2*pi*x.*y));

We could equally well have constructed chebfun2 objects for the variables $x$ and $y$ first and then computed $f$ from these:

x = chebfun2(@(x,y) x);
y = chebfun2(@(x,y) y);
f = cos(2*pi*x.*y);

There's also a shortcut cheb.xy to constructing these objects x and y, so we could also have executed

cheb.xy
f = cos(2*pi*x.*y);

Here is a plot of $f$:

plot(f), zlim([-2 2])

Along with plot, there are also commands contour and surf for displaying a chebfun2. Here is a contour plot of $f$:

contour(f), axis square

One way to find the rank of the approximant used to represent $f$, discussed in Section 8.8, is like this:

length(f)
ans =
11


f
f =
chebfun2 object
domain                 rank       corner values
[  -1,   1] x [  -1,   1]       11     [   1    1    1    1]
vertical scale =   1


The corner values are the values of the chebfun2 at $(-1,-1)$, $(-1,1)$, $(1,-1)$, and $(1,1)$, in that order. The vertical scale is used by operations to aim for close to machine precision relative to that number.

### 12.4 Basic operations

Once we have a chebfun2, we can compute quantities such as its definite double integral:

sum2(f)
ans =
0.902823333580281


This matches well the exact answer obtained by calculus, which is $(2/\pi)\hbox{Si}(2\pi)$:

exact = 0.9028233335802806267957003779
exact =
0.902823333580281


We can also evaluate a chebfun2 at a point $(x,y)$, or along a line. When evaluating along a line a chebfun is returned because the answer is a function of one variable.

Evaluation at a point:

x = 2*rand - 1; y = 2*rand - 1;
f(x,y)
ans =
-0.997677434419261


Evaluation along the line $y = \pi/6$:

f(:,pi/6)
ans =
chebfun row (1 smooth piece)
interval       length     endpoint values
[      -1,       1]       23     -0.99    -0.99
vertical scale =   1


There are plenty of other questions that may be of interest. For instance, what are the zero contours of $f(x,y) - .95$?

r = roots(f-.95);
plot(r), axis([-1 1 -1 1])
axis square, title('Zero contours of f-.95')

What is the partial derivative $\partial f/\partial y$?

fy = diff(f,1,1);
plot(fy)

The syntax for the diff command can cause confusion because we are following the matrix syntax in MATLAB. Chebfun2 also offers the more easily remembered diffx(f,k) and diffy(f,k), which differentiate $f(x,y)$ $k$ times with respect to the first and second variable, respectively.

What is the mean value of $f$ on $[-1,1]\times[-1,1]$?

mean2(f)
ans =
0.225705833395070


### 12.5 Chebfun2 methods

There are over 100 methods that can be applied to chebfun2 objects. For a complete list type:

methods chebfun2
Methods for class chebfun2:

abs                   horzcat               quiver3
biharm                imag                  rank
biharmonic            integral              rdivide
cdr                   integral2             real
chebcoeffs2           isPeriodicTech        restrict
chebfun2              isempty               roots
chebpolyplot          isequal               roots0
chebpolyplot2         isreal                rootsFeb19
chebpolyval2          iszero                rootsclean
chol                  jacobian              rootsfine
coeffs2               lap                   rootsnew
complex               laplacian             sample
compose               ldivide               sampleTest
conj                  length                simplify
constructor           log                   sin
contour               lu                    singleSignTest
contourf              max                   sinh
cos                   max2                  size
cosh                  mean                  sph2cart
ctranspose            mean2                 sphere
cumprod               median                sqrt
cumsum                mesh                  squeeze
cumsum2               min                   std
diag                  minandmax2            subsref
diff                  minandmax2est         sum
diffx                 minus                 sum2
diffy                 mldivide              sumdisk
discriminant          mrdivide              surf
disp                  mtimes                surface
display               norm                  surfacearea
domainCheck           normalizePivots       svd
domainarea            normalizeRowsAndCols  tan
eig                   pivotplot             tand
ellipsoid             pivots                tanh
end                   plot                  times
exp                   plotcoeffs            trace
explain               plotcoeffs2           transpose
feval                 plus                  uminus
fevalm                pol2cart              uplus
flipdim               poldec                vertcat
fliplr                potential             volt
flipud                power                 vscale
fred                  prod                  waterfall
get                   qr

Static methods:

coeffs2vals           poisson
outerProduct          vals2coeffs



Most of these commands have been overloaded from MATLAB. More information about a Chebfun2 command can be found with help:

help chebfun2/max2
 MAX2   Global maximum of a CHEBFUN2.
Y = MAX2(F) returns the global maximum of F over its domain.

[Y, X] = MAX2(F) returns the global maximum in Y and its location X.

This command may be faster if the OPTIMIZATION TOOLBOX is installed.



### 12.6 Composition of chebfun2 objects

New chebfun2 objects can be constructed from existing ones by composing them with operations such as +, -, .*, and .^. For example,

x = chebfun2(@(x,y) x, [-4 4 -2 2]);
y = chebfun2(@(x,y) y, [-4 4 -2 2]);
f = 1./( 2 + cos(.25 + x.^2.*y + y.^2) );
contour(f), axis equal

### 12.7 Analytic functions

An analytic function $f(z)$ can be thought of as a complex-valued function of two real variables, $f(x,y) = f(x+iy)$. If the Chebfun2 constructor is given an anonymous function with one argument, it assumes that argument is a complex variable. For instance,

f = chebfun2(@(z) sin(z));
f(1+1i), sin(1+1i)
ans =
1.298457581415977 + 0.634963914784736i
ans =
1.298457581415977 + 0.634963914784736i


These functions can be visualised by using a technique known as phase portrait plots. Given a complex number $z = re^{i\theta}$, the phase $e^{i\theta}$ can be represented by a colour. We follow Wegert's colour recommendations [Wegert 2012], using red for a phase $i$, then yellow, green, blue, and violet as the phase moves clockwise around the unit circle. For example,

f = chebfun2(@(z) sin(z)-sinh(z),2*pi*[-1 1 -1 1]);
plot(f)

Many properties of analytic functions can be visualised by these types of plots, such as the location of zeros and their multiplicities. Can you work out the multiplicity of the root at $z=0$ from this plot? For another example, try cheb.gallery2('airycomplex').

Since Chebfun2 only represents smooth functions, a trick is required to draw pictures like these for functions with poles [Trefethen 2013]. For functions with branch points or essential singularities, it is currently not possible in Chebfun2 to draw phase plots.

### 12.8 Chebfun2 low rank approximations

Chebfun2 exploits the observation that many functions of two variables can be well approximated by low rank approximants. A rank $1$ function, also known as separable, is of the form $u(y)v(x)$, and a rank $k$ function is one that can be written as the sum of $k$ rank $1$ functions. Smooth functions tend to be well approximated by functions of low rank. Chebfun2 determines low rank function approximations automatically by means of an algorithm that can be viewed as an iterative application of Gaussian elimination with complete pivoting [Townsend & Trefethen 2013]. The underlying function representations are related to work by Carvajal, Chapman and Geddes [Carvajal, Chapman, & Geddes 2008] and others including Bebendorf [Bebendorf 2008], Hackbusch, Khoromskij, Oseledets, and Tyrtyshnikov. For further aspects of low-rank representations see [Trefethen 2017] and [Beckermann and Townsend 2019].

Here is an example adapted from [Townsend & Trefethen 2013] and cheb.gallery2('smokering'). The function $$f(x,y) = \exp( -40(x^2-xy+2y^2 - 1/2)^2)$$ has the shape of an elliptical ring in the unit square, and Chebfun2 represents it by an approximation of reasonably high rank:

ff = @(x,y) exp(-40*(x.^2 - x.*y + 2*y.^2 - 1/2).^2);
f = chebfun2(ff);
levels = 0.1:0.1:0.9;
contour(f,levels), axis([-1 1 -1 1]), axis square
title(['rank ' int2str(length(f))],'fontsize',12)

To illustrate the nature of low-rank approximations, rather than letting Chebfun2 determine the rank adaptively, we can force it to take ranks $1,2,\dots ,9$. Here are the results, plotted with black level curves at heights $0.2,0.4,0.6,0.8$:

levels = 0.2:0.2:0.8;
clf
for k = 1:9
axes('position',[.03+.33*mod(k-1,3) .67-.3*floor((k-1)/3) .28 .28])
contour(chebfun2(ff,k),levels,'k')
xlim([-1 1]), axis equal, axis off
end

For this function, "plotting accuracy" is achieved approximately at rank 16; the remaining terms are then required to get from 2-3 digits to 15.

### 12.9 References

[Bebendorf 2008] M. Bebendorf, Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems, Springer, 2008.

[Beckermann & Townsend 2019] "Bounds on the singular values of matrices with displacement structure", SIAM Rev. 61 (2019), 319--344.

[Carvajal, Chapman, & Geddes 2008] O. A. Carvajal, F. W. Chapman and K. O. Geddes, "Hybrid symbolic-numeric integration in multiple dimensions via tensor-product series", Proceedings of ISSAC'05, M. Kauers, ed., ACM Press, 2005, 84-91.

[Townsend & Trefethen 2013] A. Townsend and L. N. Trefethen, "Gaussian elimination as an iterative algorithm", SIAM News, March 2013.

[Townsend & Trefethen 2013b] A. Townsend and L. N. Trefethen, "An extension of Chebfun to two dimensions", SIAM Journal on Scientific Computing, 35 (2013), C495-C518.

[Townsend & Trefethen 2014] A. Townsend and L. N. Trefethen, "Continuous analogues of matrix factorizations", Proceedings of the Royal Society A 471 (2014) 20140585.

[Trefethen 2013] L. N. Trefethen, "Phase Portraits for functions with poles", http://www.chebfun.org/examples/complex/PortraitsWithPoles.html.

[Trefethen 2017] L. N. Trefethen, "Cubature, approximation, and isotropy in the hypercube", SIAM Rev. 59 (2017), 469--491.

[Wegert 2012] E. Wegert, Visual Complex Functions: An Introduction with Phase Portraits, Birkhauser/Springer, 2012.