Code covered by the BSD License  

Highlights from
Chebfun V4

image thumbnail

Chebfun V4

by

 

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

Absolute value approximations by rationals

Absolute value approximations by rationals

Nick Trefethen, May 2011

(Chebfun example approx/AbsoluteValue.m)

Peter Lax mentioned to me recently an example that no doubt various people have thought about over the years. Suppose we think of x^2 as a given number and we try to find its square root by solving the equation

r^2 = x^2

for r using Newton's method beginning from the guess r=1. The successive iterates are given by the formula

r := (r^2+x^2)/2r .

After k steps we have a rational function of type (2^k,2^k), and these functions will approach the function abs(x).

Let's see the iteration in action:

x = chebfun('x');
r = chebfun('1');
LW = 'linewidth'; lw = 1.6; FS = 'fontsize'; fs = 12;
for k = 0:5
    subplot(3,2,k+1)
    plot(r,LW,lw), axis([-1 1 -.2 1.2]), grid on
    err = norm(r-abs(x),inf);
    s = sprintf('error=%4.1e   len=%d',err,length(r));
    title(s,FS,fs)
    r = (r.^2+x.^2)./(2*r);
end

The curves look nice, but the exponentially growing chebfun lengths do not. To improve this, we can put a breakpoint at x=0:

x = chebfun('x',[-1 0 1]);
r = chebfun('1',[-1 0 1]);
for k = 0:5
    subplot(3,2,k+1)
    plot(r,LW,lw), axis([-1 1 -.2 1.2]), grid on
    err = norm(r-abs(x),inf);
    s = sprintf('error=%4.1e   length = %d',err,length(r));
    title(s,FS,fs)
    r = (r.^2+x.^2)./(2*r);
end

It's interesting to look at the error. In the outer half of the interval, we've already achieved machine precision, whereas near x=0 the errors remain large.

clf, semilogy(abs(r-abs(x)),LW,lw)
axis([-1 1 1e-18 10]), grid on
xlabel('x',FS,fs)
title('Error',FS,fs)
Warning: Negative data ignored 

Let's take six more steps of the iteration:

for k = 0:5
    subplot(3,2,k+1)
    plot(r,LW,lw), axis([-1 1 -.2 1.2]), grid on
    err = norm(r-abs(x),inf);
    s = sprintf('error=%4.1e   length = %d',err,length(r));
    title(s,FS,fs)
    r = (r.^2+x.^2)./(2*r);
end

Here is the error:

clf, semilogy(abs(r-abs(x)),LW,lw)
axis([-1 1 1e-18 10]), grid on
xlabel('x',FS,fs)
title('Error',FS,fs)
Warning: Negative data ignored 

Evidently we are getting convergence to abs(x), for all x. In the infinity norm, the rate looks pretty disappointing. Donald Newman showed that the optimal type (n,n) rational approximants to abs(x) achieve accuracy O(exp(-C*sqrt(n))) [1,2], whereas here the maximum error is exactly 2^(-k) after k steps, which corresponds to 1/n for the type (n,n) approximation. Away from x=0, however, the accuracy is O(exp(-C*n)), thanks to the quadratic convergence of Newton's method.

Incidentally, note that this last curve is not very close to symmetrical about x=0. I wonder why not?

References:

[1] D. J. Newman, Rational approximation of abs(x), Michigan Mathematical Journal 11 (1964), 11-14.

[2] L. N. Trefethen, Approximation Theory and Approximation Practice, draft book available at http://www.maths.ox.ac.uk/chebfun/ATAP/.

Contact us