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,
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.
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:
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. :
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
f = chebfun2(@(x,y) sin(30*x.*y)); norm(f,inf)
ans = 0.999999999999995
[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.