#' Multivariate Newton method
#'
#' \code{multinewton()} assumes that f is a vector-valued function of vector
#' argument, although both can be one dimensional. It is therefore a
#' generalization of \code{uninewton()}, but has slightly different output.
#'
#' @param f function
#' @param df function; Jacobian matrix of f
#' @param x0 initial value
#' @param tol tolerance, defaults to 10*.Machine$double.eps
#' @param maxit maximum number of iterations
#' @return a list
#' @name multinewton
#' @examples
#'
#' library("kumerical")
#'
#' f <- function(x) x^2 - 2
#' df <- function(x) 2*x
#' x0 <- 2
#' uninewton(f, df, x0)
#' simple_multinewton(f, df, x0)
#' multinewton(f, df, x0)
#' str(multinewton(f, df, x0))
#'
#' # this is easier with mpoly:
#' library("mpoly")
#' (p <- mp("x^2 - 2"))
#' f <- as.function(p)
#' df <- as.function(gradient(p))
#' x0 <- 2
#' simple_multinewton(f, df, x0)
#' multinewton(f, df, x0)
#'
#'
#'
#'
#' jacobian <- function(ps, varorder = vars(ps)) {
#' list_of_mpolyLists <- lapply(ps, deriv, var = varorder)
#' list_of_gradient_functions <- lapply(
#' list_of_mpolyLists, as.function,
#' varorder = varorder, silent = TRUE
#' )
#' J <- function(.) lapply(list_of_gradient_functions, function(f) f(.))
#' function(v) do.call(rbind, J(v))
#' }
#'
#'
#'
#' # intersection of the parabola y = x^2 and circle x^2 + y^2 = 1
#' # algebraically, this is
#' # y + y^2 = 1 => y^2 + y - 1 = 0 =>
#' plus_y <- (-1 + sqrt(1 - 4*(1)*(-1))) / (2*1) # = 0.618034 and
#' minus_y <- (-1 - sqrt(1 - 4*(1)*(-1))) / (2*1) # = -1.618034
#' # so that
#' # x = sqrt( plus_y) = =-0.7861514 and
#' # x = sqrt(minus_y) = +-1.27202i
#' # for solutions (+-0.7861514, 0.618034) and (+-1.27202i, -1.618034)
#' theoretical_solns <- list(
#' c( sqrt(plus_y), plus_y), c( sqrt(-minus_y)*1i, minus_y),
#' c(-sqrt(plus_y), plus_y), c(-sqrt(-minus_y)*1i, minus_y)
#' )
#' ps <- mp(c("y - x^2", "x^2 + y^2 - 1"))
#' f <- as.function(ps, varorder = c("x", "y"))
#' lapply(theoretical_solns, f)
#' df <- jacobian(ps, varorder = c("x", "y"))
#' x0 <- c(2, 2)
#' f(x0)
#' df(x0)
#'
#' simple_multinewton(f, df, x0)
#' out <- multinewton(f, df, x0)
#' str(out, 1)
#' str(out$evals, 1)
#'
#'
#'
#' # intersection of a plane, hyperboloid, and cone
#' # true solutions =
#' # c(-3/sqrt(2), 0, 3/sqrt(2))
#' # c( 3/sqrt(2), 0, -3/sqrt(2))
#' # corresponding to the nonlinear system
#' # x + y + z = 0
#' # x^2 - y^2 + z^2 = 9,
#' # x^2 + y^2 - z^2 = 0
#' ps <- mp(c("x + y + z", "x^2 - y^2 + z^2 - 9", "x^2 + y^2 - z^2"))
#' f <- as.function(ps, varorder = c("x", "y", "z"))
#' df <- jacobian(ps, varorder = c("x", "y", "z"))
#' x0 <- c(2, 2, 2)
#' f(x0)
#' df(x0)
#' out <- multinewton(f, df, x0)
#' str(out, 1)
#' c( 3/sqrt(2), 0, -3/sqrt(2))
#' out$root
#'
#' @rdname multinewton
#' @export
multinewton <- function(f, df, x0, tol = 10*.Machine$double.eps, maxit = 100L) {
# initialize x and fx
px <- length(x0)
fx0 <- f(x0)
pf <- length(fx0)
x <- matrix(NA_real_, nrow = maxit, ncol = px)
fx <- matrix(NA_real_, nrow = maxit, ncol = pf)
dfx <- replicate(maxit, matrix(NA_real_, nrow = pf, ncol = px), simplify = FALSE)
# set norm
norm <- function(v) sum(abs(v))
# check endpoint and return early if root there
x[1,] <- x0;
fx[1,] <- fx0
if(norm(fx[1,]) <= tol) {
return(list(
root = x[1,], f.root = fx[1,],
evals = data_frame(x = x[1,], fx = fx[1,], dfx = dfx[1]),
n_evals = 1
))
}
# loop
for(k in 2:maxit) {
dfx[[k-1]] <- df(x[k-1,])
h <- solve(dfx[[k-1]], -fx[k-1,])
x[k,] <- x[k-1,] + h
fx[k,] <- f(x[k,])
n_evals <- k
if(norm(fx[k,]) <= tol) break
}
# return
list(
root = x[n_evals,], f.root = fx[n_evals,],
evals = list("x" = x[1:n_evals,], "fx" = fx[1:n_evals,], "dfx" = dfx[1:(n_evals-1)]),
n_evals = n_evals
)
}
#' @rdname multinewton
#' @export
simple_multinewton <- function(f, df, x0, tol = 10*.Machine$double.eps, maxit = 100L) {
# initialize x and fx
px <- length(x0)
fx0 <- f(x0)
pf <- length(fx0)
x <- matrix(NA_real_, nrow = maxit, ncol = px)
fx <- matrix(NA_real_, nrow = maxit, ncol = pf)
dfx <- replicate(maxit, matrix(NA_real_, nrow = pf, ncol = px), simplify = FALSE)
# set norm
norm <- function(v) sum(abs(v))
# check endpoint and return early if root there
x[1,] <- x0;
fx[1,] <- fx0
if(norm(fx[1,]) <= tol) return(list(root = x[1,], f.root = fx[1,]))
# loop
for(k in 2:maxit) {
dfx[[k-1]] <- df(x[k-1,])
h <- solve(dfx[[k-1]], -fx[k-1,])
x[k,] <- x[k-1,] + h
fx[k,] <- f(x[k,])
n_evals <- k
if(norm(fx[k,]) <= tol) break
}
# return
list(root = x[n_evals,], f.root = fx[n_evals,])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.