Если вы заинтересованы в численной оптимизации, у вас есть несколько возможных подходов. Самый прямой - использовать optim. По умолчанию это симплекс-метод Nelder-Mead, но другие реализованы.
Вам нужно будет начать с разных начальных значений, чтобы сходиться к разным конечным точкам. Я могу предложить вам следующее:
func <- function(a){
x <- a[1]
y <- a[2]
return(
0.5*(x^4 - 16*x^2 + 5*x + y^4 - 16*y^2 + 5*y)
)
}
t0 <- rnorm(100, sd = 20)
t1 <- rnorm(100, sd = 20)
points <- do.call(rbind, lapply(1:100, function(i) optim(par = c(t0[i],t1[i]), fn = func)$par))
А если вы хотите графически увидеть ваши решения:
library(ggplot2)
ggplot(data.frame(points)) + geom_point(aes(x = X1, y = X2))
У вас есть четыре локальных минимума в этом выводе