Как получить с помощью Ryacas

Я совершенно не использую пакет Ryacas. Я хотел бы получить последовательные производные полиномов, поэтому я попробовал следующие строки:

library(Ryacas)
x <- Sym("x")
P <- 1
for (k in 1:10) {
  P <- (1+k*x)*P + x*(1-x)*deriv(P, x)
  P
}

Я ожидал получить последовательные полиномы:

P = (1 + x)(1) = 1 + x
P = (1 + 2x)(1 +x) + x(1 - x)(1) = 1 + 4x + x^2
P = (1 + 3x)(1 + 4x + x^2) + x(1 - x)(4 + 2x) = 1 + 11x + 11x^2 + x^3

и так далее.

Ошибка не отображается, но вывод не отображается. Так что написание моих строк явно неправильное. Не могли бы вы мне помочь, пожалуйста ?


r derivative yacas
person Andrew    schedule 30.03.2018    source источник
comment
Я бы попробовал P <- Sym(1).   -  person Stéphane Laurent    schedule 30.03.2018
comment
Кроме того, внутри цикла вам нужно print(P), а не только P.   -  person G. Grothendieck    schedule 30.03.2018
comment
С P <- Sym(1) у меня работает, когда я делаю for(k in 1:2), но с for(k in 1:10) расчеты не заканчиваются через 10 минут...   -  person Stéphane Laurent    schedule 30.03.2018
comment
ХОРОШО. Это работает с Sym(1), но вы правы: время расчета - проблема...   -  person Andrew    schedule 30.03.2018
comment
Любая подсказка, чтобы получить более компактные выражения? Например, иметь $1+4x +x^2$ при k=2?   -  person Andrew    schedule 30.03.2018
comment
Попробуйте Simplify(P)   -  person Stéphane Laurent    schedule 30.03.2018
comment
Кажется, это работает. Я пишу ответ, чтобы подытожить.   -  person Stéphane Laurent    schedule 30.03.2018

Ответы (1)


arrow_upward
1
arrow_downward

Используйте P <- Sym(1) в начале. Также используйте Simplify в своем цикле, иначе вычисления будут слишком трудоемкими. Вы также должны использовать print в цикле.

library(Ryacas)
x <- Sym("x")
P <- Sym(1)
for (k in 1:10) {
  P <- Simplify((1+k*x)*P + x*(1-x)*deriv(P, x))
  print(P)
}

Результат:

expression(3 * x^2 - x^3 + 9 * x + 1)
expression(x^4 - 4 * x^3 + 18 * x^2 + 20 * x + 1)
expression(5 * x^4 - x^5 + 2 * x^3 + 94 * x^2 + 43 * x + 1)
expression(x^6 - 6 * x^5 + 27 * x^4 + 196 * x^3 + 411 * x^2 + 
    90 * x + 1)
expression(7 * x^6 - x^7 - 9 * x^5 + 527 * x^4 + 2017 * x^3 + 
    1593 * x^2 + 185 * x + 1)
expression(x^8 - 8 * x^7 + 40 * x^6 + 1000 * x^5 + 8686 * x^4 + 
    14440 * x^3 + 5704 * x^2 + 376 * x + 1)
expression(9 * x^8 - x^9 - 24 * x^7 + 2280 * x^6 + 32058 * x^5 + 
    101190 * x^4 + 86280 * x^3 + 19368 * x^2 + 759 * x + 1)
expression(x^10 - 10 * x^9 + 57 * x^8 + 4368 * x^7 + 112134 * 
    x^6 + 597108 * x^5 + 937350 * x^4 + 461328 * x^3 + 63417 * 
    x^2 + 1526 * x + 1)
expression(11 * x^10 - x^11 - 43 * x^9 + 9249 * x^8 + 371346 * 
    x^7 + 3173370 * x^6 + 8269398 * x^5 + 7454718 * x^4 + 2289231 * 
    x^3 + 202459 * x^2 + 3061 * x + 1)
expression(x^12 - 12 * x^11 + 78 * x^10 + 18068 * x^9 + 1197279 * 
    x^8 + 15664248 * x^7 + 63560580 * x^6 + 94344696 * x^5 + 
    53298207 * x^4 + 10776596 * x^3 + 634926 * x^2 + 6132 * x + 
    1)
person Stéphane Laurent    schedule 30.03.2018
comment
Большое спасибо ; выглядит идеально. - person Andrew; 30.03.2018