[R] eigenvalues of matrices of partial derivatives with ryacas

[Gabor Grothendieck]

On Thu, Mar 15, 2012 at 7:51 PM, Adam Zeilinger <zeil0006 at umn.edu> wrote:
> Hello,
>
> I am trying to construct two matrices, F and V, composed of partial
> derivatives and then find the eigenvalues of F*Inverse(V).  I have the
> following equations in ryacas notation:
>
>> library(Ryacas)
>> FIh <- Expr("betah*Sh*Iv")
>> FIv <- Expr("betav*Sv*Ih")
>> VIh <- Expr("(muh + gamma)*Ih")
>> VIv <- Expr("muv*Iv")
>
> I successfully found the partial derivatives:
>
>> f11 <- deriv(FIh, "Ih")
>> f12 <- deriv(FIh, "Iv")
>> f21 <- deriv(FIv, "Ih")
>> f22 <- deriv(FIv, "Iv")
>> v11 <- deriv(VIh, "Ih")
>> v12 <- deriv(VIh, "Iv")
>> v21 <- deriv(VIv, "Ih")
>> v22 <- deriv(VIv, "Iv")
>
> Next I would like to put these partial derivatives into two matrices, F and
> V:
>
>> F <- Expr("{{f11, f12}, {f21, f22}}")
>> V <- Expr("{{v11, v12}, {v21, v22}}")
>
> Finally, I would like to find the eigenvalues of F*Inverse(V).  Something
> like:
>
>>  yacas("EigenValues(F*Inverse(V))")
>
> However, this does not work.  I get the following error message:
>
> In function "While" : bad argument number 1 (counting from 1)The offending
> argument $ii49<= $nr49 evaluated to Not Length-1<0CommandLine(1) : Invalid
> argument
>
> According to Mathematica, the correct eigenvalues are:
>
> {-((Sqrt[betah] Sqrt[betav] Sqrt[Sh] Sqrt[Sv])/Sqrt[gamma muv + muh muv]),
>   (Sqrt[betah] Sqrt[betav] Sqrt[Sh] Sqrt[Sv])/Sqrt[gamma muv + muh muv]}
>
> I don't understand the error message.  Any suggestions on how to get the
> correct eigenvalues using R would be greatly appreciated.  I'm using R
> 2.14.0.

Try doing it this way instead:

library(Ryacas)
EigenValues <- function(x) Sym("EigenValues(", x, ")")

Iv <- Sym("Iv"); Ih <- Sym("Ih")
Sv <- Sym("Sv"); Sh <- Sym("Sh")
betav <- Sym("betav"); betah <- Sym("betah")
muv <- Sym("muv"); muh <- Sym("muh")
gamma <- Sym("gamma")

FIh <- betah*Sh*Iv
FIv <- betav*Sv*Ih
VIh <- (muh + gamma)*Ih
VIv <- muv*Iv

f11 <- deriv(FIh, Ih)
f12 <- deriv(FIh, Iv)
f21 <- deriv(FIv, Ih)
f22 <- deriv(FIv, Iv)
v11 <- deriv(VIh, Ih)
v12 <- deriv(VIh, Iv)
v21 <- deriv(VIv, Ih)
v22 <- deriv(VIv, Iv)

F <- List(List(f11, f12), List(f21, f22))
V <- List(List(v11, v12), List(v21, v22))

EigenValues(F * Inverse(V))

For the last line I get:

> EigenValues(F*Inverse(V))
expression(Roots(xx^2 - betav * Sv * muv * (betah * Sh * (muh +
    gamma))/((muh + gamma) * muv)^2))

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com