multinewton: Multivariate Newton method

In dkahle/kumerical: Numerical Algorithms in R

Description

multinewton() assumes that f is a vector-valued function of vector argument, although both can be one dimensional. It is therefore a generalization of uninewton(), but has slightly different output.

Usage

multinewton(f, df, x0, tol = 10 * .Machine$double.eps, maxit = 100L)

simple_multinewton(f, df, x0, tol = 10 * .Machine$double.eps, maxit = 100L)

Arguments

f	function
df	function; Jacobian matrix of f
x0	initial value
tol	tolerance, defaults to 10*.Machine$double.eps
maxit	maximum number of iterations

Value

a list

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

dkahle/kumerical documentation built on May 27, 2019, 11:49 p.m.