14.1 Zero contours of a bivariate function: roots
The Chebfun2 roots command can compute the zero contours of a function of two variables. For example, here we compute the zero contours of Trott's curve, an example from algebraic geometry [Trott 1997].
cheb.xy; trott = 144*(x.^4+y.^4) - 225*(x.^2+y.^2) + 350*x.^2.*y.^2 + 81; r = roots(trott); plot(r), axis([-1 1 -1 1]), axis square
The curves are represented as complex-valued chebfuns (see Section 13.4). For example, here is one of the four components:
r(:,1)
ans =
chebfun column (1 smooth piece)
interval length endpoint values
[ -1, 1] 576 complex values
vertical scale = 1
Although roots can be outwitted, it is often remarkably accurate. For example, here is the perimeter of a circle of radius 1/2, which we measure by using the fact that the arc length is equal to the 1-norm of the derivative.
f = x.^2 + y.^2 - 1/4; perimeter = norm(diff(roots(f)),1) exact_perimeter = pi
perimeter = 3.141592653589794 exact_perimeter = 3.141592653589793
For some more exotic examples of zero curves computed by roots, see the 2D Approximation section of the Chebfun examples collection at www.chebfun.org.
14.2 Zeros of a pair of bivariate functions: roots again
Chebfun2 roots can also find zeros of bivariate systems, i.e., solutions to $f(x,y) = g(x,y) = 0$. Generically, these are isolated points.
For example, which points on Trott's curve intersect the circle of radius $0.9$?
g = x.^2 + y.^2 - .9^2; r = roots(trott,g) plot(roots(trott),'b'), hold on plot(roots(g),'r') MS = 'markersize'; plot(r(:,1),r(:,2),'.k',MS,20) axis([-1 1 -1 1]), axis square, hold off
r = -0.799441089368585 -0.413393208252350 -0.799441089368583 0.413393208252352 -0.413393208252347 -0.799441089368586 -0.413393208252346 0.799441089368587 0.413393208252345 -0.799441089368587 0.413393208252344 0.799441089368588 0.799441089368588 -0.413393208252343 0.799441089368587 0.413393208252346
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 [Güttel 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), hold on
plot(figof8,'r'), 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 = 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.306s minandmax2 command = 0.093s
Here is a complicated function from the 2002 SIAM 100-Dollar, 100-Digit Challenge [Bornemann et al. 2004]. Chebfun2 computes its global minimum in a fraction of a second:
tic
f = cheb.gallery2('challenge');
[minval,minpos] = min2(f);
minval
toc
minval = -3.306868647474791 Elapsed time is 0.182211 seconds.
The result closely matches the correct solution, computed to 10,000 digits by 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), 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 = (x.^2-y.^3+1/8).*sin(10*x.*y); r = roots(gradient(f)); % critical points plot(roots(diff(f,1,2)),'b'), hold on % zero contours of f_x plot(roots(diff(f)),'r') % 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 = sin(30*x.*y); norm(f,inf)
ans = 1.000000000000003
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 Numerical Computing, SIAM 2004.
[Güttel 2010] S. Güttel, "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 129 (2015), 181--209.
[Trott 2007] M. Trott, "Applying GroebnerBasis to three problems in geometry", Mathematica in Education and Research, 6 (1997), 15-28.