91 - Using the 'SymPy' object directly
In caracas: Computer Algebra

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(caracas)

inline_code <- function(x) {
  x
}

if (!has_sympy()) {
  # SymPy not available, so the chunks shall not be evaluated
  knitr::opts_chunk$set(eval = FALSE)

  inline_code <- function(x) {
    deparse(substitute(x))
  }
}

Using `SymPy` directly

First we get the SymPy object:

sympy <- get_sympy()

sympy$diff("2*a*x", "x")
sympy$solve("x**2 - 1", "x")

Elaborate example

How can we minimise the amount of material used to produce a cylindric tin can that contains 1 litre. The cylinder has diameter $d$ and height $h$. The question is therefore: What is $d$ and $h$?

We introduce the variables d (diameter) and h (height):

d <- sympy$symbols('d')
h <- sympy$symbols('h')

The problem is a constrained optimisation problem, and we solve it by a Lagrange multiplier, and therefore we introduce lam (the Lagrange multiplier):

lam <- sympy$symbols('lam')

We now set up the problem:

area_str <- "Pi/2 * d**2 + Pi * h * d"
vol_str <- "Pi/4 * d**2 * h"
lap_str <- paste0("(", area_str, ") - lam*((", vol_str, ") - 1)")
lap <- sympy$parsing$sympy_parser$parse_expr(
  lap_str,
  local_dict = list('d' = d, 'h' = h, 'lam' = lam))

We can now find the gradient:

grad <- sympy$derive_by_array(lap, list(d, h, lam))
grad

And find the critical points:

sol <- sympy$solve(grad, list(d, h, lam), dict = TRUE)
sol

We take the one with the real solution:

sol[[1]]

We now have a short helper function to help getting appropriate R expressions (such a function will be included in later versions of this package):

to_r <- function(x) {
  x <- as.character(x)
  x <- gsub("Pi", "pi", x, fixed = TRUE)
  x <- gsub("**", "^", x, fixed = TRUE)
  x <- parse(text = x)
  return(x)
}

sol_d <- to_r(sol[[1]]$d)
sol_d
eval(sol_d)
sol_h <- to_r(sol[[1]]$h)
sol_h
eval(sol_h)

(It is left as an exercise to the reader to show that the critical point indeed is a minimum.)

Simple example with assumptions

x <- sympy$symbols('x')
x$assumptions0
x <- sympy$symbols('x', positive = TRUE)
x$assumptions0
eq <- sympy$parsing$sympy_parser$parse_expr("x**2 - 1",
                                            local_dict = list('x' = x))
sympy$solve(eq, x, dict = TRUE)

Another example with assumptions

x <- sympy$symbols('x', positive = TRUE)
eq <- sympy$parsing$sympy_parser$parse_expr("x**3/3 - x",
                                            local_dict = list('x' = x))
eq
grad <- sympy$derive_by_array(eq, x)
grad
sympy$solve(grad, x, dict = TRUE)