The function \[ f(x) = \frac{1}{2} \log(x^2+1) + x \arctan(x) \]
is convex with a minimum at \(x = 0\):
f <- function(x) 1/2 * log(1 + x^2) + x * atan(x)
x <- seq(-5, 5, len = 101)
plot(x, f(x), type = "l")
The first derivative is \[ f'(x) = d_1(x) = \arctan(x) \]
d1 <- atan
plot(x, d1(x), type = "l")
The second derivative is \[ f''(x) = d_2(x) = \frac{1}{1 + x^2} \]
d2 <- function(x) 1 / (1 + x^2)
plot(x, d2(x), type = "l")
A single Newton step is taken by
ns <- function(x) x - d1(x) / d2(x)
A set of n
steps is returned by
run <- function(x, step, n = 10) {
v <- matrix(0, nrow = n, ncol = length(x))
for (i in seq_len(n)) {
x <- step(x)
v[i, ] <- x
}
v
}
For a starting point close enough to zero convergence is very fast:
run(1.3, ns)
## [,1]
## [1,] -1.161621e+00
## [2,] 8.588964e-01
## [3,] -3.742407e-01
## [4,] 3.401887e-02
## [5,] -2.624025e-05
## [6,] 1.204517e-14
## [7,] 0.000000e+00
## [8,] 0.000000e+00
## [9,] 0.000000e+00
## [10,] 0.000000e+00
Start a little too far away:
run(1.4, ns)
## [,1]
## [1,] -1.413619e+00
## [2,] 1.450129e+00
## [3,] -1.550626e+00
## [4,] 1.847054e+00
## [5,] -2.893562e+00
## [6,] 8.710326e+00
## [7,] -1.032498e+02
## [8,] 1.654056e+04
## [9,] -4.297215e+08
## [10,] 2.900641e+17
This plot shows the range of convergence:
plot(x, ns(x), type = "l", xlim = c(-2, 2), ylim = c(-5, 5))
lines(x, -x, lty = 2)