CODEWITHSUNDEEP
pythonrmathsympy
Felipe Queiroz
Felipe Queiroz

Reputation: 123

Is it actually possible to use any Computer Algebra System (CAS) in R?

I was searching for some way to use Sympy in R, but I didn't find anything that could solve the equation below without errors.

Equation

Is is actually possible to use something like Sympy in R?

Upvotes: 1

Views: 479

Answers (2)

Ben Bolker
Ben Bolker

Reputation: 226662

It will definitely be easiest to do this numerically, if you can live with that:

uniroot(function(x) 1000/(1+x)^(25/252) - 985, c(0,10))

However: I tried this in Wolfram Alpha, which gave me the exact result

x = (102400000000000000000000 2^(6/25) 5^(4/25) 197^(23/25) - 
     17343170265605241347130653)/17343170265605241347130653

and then moments later reverted to giving me the numerical approximation (x ≈ 0.164562) -- I had to quickly copy the text before it disappeared.

When I tried this with sympy, it quickly gave me a numerical answer:

from sympy import *
x = Symbol("x")
solve(1000/((1+x)**(25/252)) - 985, x)
[0.164562487329317]

... but trying to get an exact solution by setting numerical = False was no good: it seemed to have to work fairly hard (pegged CPU at 100% for several minutes before I gave up). This is consistent with my results from Ryacas,

library(Ryacas)
yac("Solve(1000/((1+x)^(25/252)) == 985, x)")

Error in yac_core(x) : Yacas returned this error: CommandLine(1) : Max evaluation stack depth reached. Please use MaxEvalDepth to increase the stack size as needed.

However, there is a surprising shortcut. The hard part of this problem is actually doing the deep recursion needed to unravel the (25/252) power exactly. If you're happy with a more general solution, you can get it quickly:

In sympy:

a = Symbol("a")
solve(1000/((1+x)**(a)) - 985, x, numerical = False)
[-1 + 200**(1/a)/197**(1/a)]

For what it's worth, this looks simpler (to a human) as -1 + (200/197)**(1/a) (there's probably some simplification setting that would convince sympy to do this but ...)

This works in caracas (an R wrapper for sympy) too, once everything is set up:

library(caracas)
x <- symbol('x')
a <- symbol('a')
solve_sys(1000/((1+x)^a), 985, list(x))
## Solution 1:
##   x =          -1         
##               ───        
##                a  a _____
##       -1 + 197   ⋅╲╱ 200

Unfortunately, the same trick doesn't seem to work for Ryacas.

Upvotes: 1

Dirk is no longer here
Dirk is no longer here

Reputation: 368449

There are several packages at CRAN that may help:

They may all have different levels of difficulty in terms installing the underlying 'engine' but this should get you started.

Upvotes: 3

Related Questions