Решение системы из 50 нелинейных уравнений - ошибка в fn (par, ...): аргумент отсутствует - PullRequest
Новая
       64

Решение системы из 50 нелинейных уравнений - ошибка в fn (par, ...): аргумент отсутствует

0 голосов
/ 19 мая 2019

Я пытаюсь решить систему из 50 нелинейных уравнений следующим образом:

Учитывая вектор, y, который содержит 49 различных значений, я хотел бы перевести эти 49 значений в несколько разные значенияв другом векторе, скажем, x такой, что: 

log(x[1], y[1]) = n
...
log(x[49], y[49]) = n
x[1] + ... + x[49] = 1

Для ясности в приведенных выше уравнениях y[i] подразумевается как основание логарифма.

написал следующий код, который, однако, не работает: 

library(xlsx)
rm(list=ls())
setwd("C:/Users/.../folder")

my_data <- read.xlsx("samplefile.xlsx", 1)
y <- matrix(0:0, nrow=49,ncol=1)

for(i in 1:49) {
  if(my_data[i,1]!=0) {
    y[i,1] = 1/my_data[i,1]
    }
  }

for(i in 1:49) {
  fn <- function(x,n) {
  dummy1 <- log(x[i],y[i])-n
  dummy2 <- sum(x[1:49])-1

  return(c(dummy1,dummy2))
}
}

guess <- matrix(0.5:0.5, nrow = 50, ncol = 1)

nleqslv(guess,fn)

Я бы ожидал, что он решает для x[i] и n.Тем не менее, я получаю следующее сообщение об ошибке:

"Ошибка в fn (par, ...): отсутствует аргумент" n "без значения по умолчанию"

Редактировать: Форматирование

1 Ответ

0 голосов
/ 20 мая 2019

Есть два способа решения вашей проблемы.

Первый использует nleqslv для решения системы из 50 уравнений с 50 переменными: x[1] до x[49] и n. Функция, которую нужно решить с помощью nleqslv, должна возвращать вектор, содержащий 50 элементов; система уравнений должна быть квадратной.

Вторая упрощает вашу задачу до одного уравнения в n.

Первое решение: 

library(nleqslv)
fn <- function(z,y) {

    x <- z[1:49]
    n <- z[50]

    f <- numeric(50)

    for( k in 1:49 ) f[k] <- log(x[k],y[k]) - n
    f[50] <- sum(x) - 1
    f
}

Сгенерируйте некоторые значения для y и начальные значения для x и n и попробуйте решить 

y <- rep(2,49)
guess <- c( rep(.5,49), 1)

res <- nleqslv(guess,fn,y=y)
res

Результат проверки: 

# this should sum to 1
sum(res$x[1:49])
# value of n
res$x[50]

Второе решение:

Используйте тот факт, что log(a,b) = n подразумевает a <- b^n. Так что log(x[k],y[k])=n эквивалентно x[k]=y[k]^n. Таким образом, x[1:49] может быть вычислено немедленно. Нам нужно только определить n из ограничения, что элементы x составляют 1. Это подразумевает простую функцию 

f2 <- function(n,y) {
    x <- y^n  # gives values for x
    sum(x) - 1
}

А теперь используйте uniroot, как это при условии, что значения для n будут лежать между -10 и 10 и что знаки f2 в конечных точках имеют противоположный знак. 

uniroot(f2,c(-10,10), y=y)

Вы можете запустить этот код, как указано, и убедиться, что nleqslv и uniroot дают одинаковый результат для n.

Сработает ли это для ваших данных, зависит от вас.

...