14.1 Zero contours

Chebfun2 comes with the capability to compute the zero contours of a function of two variables. For example, here we compute a representation of Trott's curve, an example from algebraic geometry [Trott 1997].

x = chebfun2(@(x,y) x); y = chebfun2(@(x,y) y);
trott = 144*(x.^4+y.^4)-225*(x.^2+y.^2) + 350*x.^2.*y.^2+81;
r = roots(trott);
LW = 'linewidth'; MS = 'markersize';
plot(r,LW,1.6), axis([-1 1 -1 1]), axis square

The zero curves are represented as complex-valued chebfuns (see Section 13.4). For example,

r(:,1)
ans =
   chebfun column (1 smooth piece)
       interval       length   endpoint values  
[      -1,       1]      576    complex values 
Epslevel = 1.000000e-05.  Vscale = 9.980190e-01.

The zero contours of a function are computed by Chebfun2 to plotting accuracy: they are typically not accurate to machine precision.

14.2 roots

Chebfun2 also comes with the capability of finding zeros of bivariate systems, i.e., the solutions to $f(x,y) = g(x,y) = 0$. If the roots command is supplied with one chebfun2, it computes the zero contours of that function, as in the last section. However, if it is supplied with two chebfun2 objects, as in roots(f,g), then it computes the roots of the bivariate system. Generically, these are isolated points.

What points on Trott's curve intersect the circle of radius $0.9$?

g = chebfun2(@(x,y) x.^2 + y.^2 - .9^2);
r = roots(trott,g);
plot(roots(trott),'b',LW,1.6), hold on
plot(roots(g),'r',LW,1.6)
plot(r(:,1),r(:,2),'.k',LW,1.6,MS,20)
axis([-1 1 -1 1]), axis square, hold off

The solutions to bivariate polynomial systems and intersections of curves are typically computed to full machine precision.

14.3 Intersections of curves

The problem of determining the intersections of real parameterised complex curves can be expressed as a bivariate rootfinding problem. For instance, here are the intersections between the 'splat' curve [Guettel 2010] and a 'figure-of-eight' curve.

t = chebfun('t',[0,2*pi]);
sp = exp(1i*t) + (1+1i)*sin(6*t).^2;     % splat curve
figof8 = cos(t) + 1i*sin(2*t);           % figure of eight curve
plot(sp,LW,1.6), hold on
plot(figof8,'r',LW,1.6), axis equal

d = [0 2*pi 0 2*pi];
f = chebfun2(@(s,t) sp(t)-figof8(s),d);  % rootfinding
r = roots(real(f),imag(f));              % calculate intersections
spr = sp(r(:,2));
plot(real(spr),imag(spr),'.k',MS,20), ylim([-1.1 2.1])
hold off

Chebfun2 rootfinding is based on an algorithm described in [Nakatsukasa, Noferini & Townsend 2014].

14.4 Global optimisation: max2, min2, and minandmax2

Chebfun2 also provides functionality for global optimisation. Here is an example, where we plot the minimum and maximum as red dots.

f = chebfun2(@(x,y) sin(30*x.*y) + sin(10*y.*x.^2) + exp(-x.^2-(y-.8).^2));
[mn mnloc] = min2(f);
[mx mxloc] = max2(f);
plot(f), hold on
plot3(mnloc(1),mnloc(2),mn,'.r',MS,40)
plot3(mxloc(1),mxloc(2),mx,'.r',MS,30)
zlim([-6 6]), colormap('bone'), hold off

If both the global maximum and minimum are required, it is roughly twice as fast to compute them at the same time by using the minandmax2 command. For instance,

tic; [mn mnloc] = min2(f);  [mx mxloc] = max2(f); t=toc;
fprintf('min2 and max2 separately = %5.3fs\n',t)
tic; [Y X] = minandmax2(f); t=toc;
fprintf('minandmax2 command = %5.3fs\n',t)
min2 and max2 separately = 0.487s
minandmax2 command = 0.219s

Here is a complicated function from the 2002 SIAM 100-Dollar, 100-Digit Challenge. Chebfun2 computes its global minimum in a few seconds:

tic
f = cheb.gallery2('challenge');
[minval,minpos] = min2(f);
minval
toc
minval =
  -3.306868647474791
Elapsed time is 1.147926 seconds.

The result closely matches the correct solution, computed to 10,000 digits Bornemann et al. [2004]:

exact = -3.306868647475237280076113
exact =
  -3.306868647475237

Here is a contour plot of this wiggly function, with the minimum circled in black:

colormap('default'), contour(f), hold on
plot(minpos(1),minpos(2),'ok',MS,20,LW,3), hold off

14.5 Critical points

The critical points of a smooth function of two variables can be located by finding the zeros of $\partial f/ \partial y = \partial f / \partial x = 0$. This is a rootfinding problem. For example,

f = chebfun2(@(x,y) (x.^2-y.^3+1/8).*sin(10*x.*y));
r = roots(gradient(f));                       % critical points
plot(roots(diff(f,1,2)),'b',LW,1.2), hold on  % zero contours of f_x
plot(roots(diff(f)),'r',LW,1.2)               % zero contours of f_y
plot(r(:,1),r(:,2),'k.',MS,24)                % extrema
axis([-1,1,-1,1]), axis square

There is a new command here called gradient that computes the gradient vector and represents it as a chebfun2v object. The roots command then solves for the isolated roots of the bivariate polynomial system represented in the chebfun2v representing the gradient. For more information about gradient, see Chapter 15.

14.6 Infinity norm

The $\infty$-norm of a function is the maximum absolute value in its domain. It can be computed by passing the argument inf to the norm command.

f = chebfun2(@(x,y) sin(30*x.*y));
norm(f,inf)
ans =
   0.999999999999995

14.7 References

[Bornemann et al. 2004] F. Bornemann, D. Laurie, S. Wagon and J. Waldvogel, The SIAM 100-Digit Challenge: A Study in High-Accuracy Nuemrical Computing, SIAM 2004.

[Guettel 2010] S. Guettel, "Area and centroid of a 2D region", http://www.chebfun.org/examples/geom/Area.html.

[Nakatsukasa, Noferini & Townsend 2014] Y. Nakatsukasa, V. Noferini and A. Townsend, "Computing the common zeros of two bivariate functions via Bezout resultants", Numerische Mathematik, to appear.

[Trott 2007] M. Trott, "Applying GroebnerBasis to three problems in geometry", Mathematica in Education and Research, 6 (1997), 15-28.