/
AbsoluteValue.m
98 lines (87 loc) · 2.88 KB
/
AbsoluteValue.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
%% Absolute value approximations by rationals
% Nick Trefethen, May 2011
%%
% (Chebfun example approx/AbsoluteValue.m)
% [Tags: #rational, #Newton, #ABS]
%%
% 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 $|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)
%%
% 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)
%%
% Evidently we are getting convergence to $|x|$, for all $x$. In the
% $\infty$-norm, the rate looks pretty disappointing. Donald Newman showed
% that the optimal type $(n,n)$ rational approximants to $|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(-Cn))$, 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 $|x|$, _Michigan Mathematical
% Journal_, 11 (1964), 11-14.
%
% 2. L. N. Trefethen, _Approximation Theory and Approximation Practice_,
% SIAM, 2013.