# Using the ‘SymPy’ object directly

library(caracas)

## Using SymPy directly

First we get the SymPy object:

sympy <- get_sympy()
sympy$diff("2*a*x", "x") #> 2*a sympy$solve("x**2 - 1", "x")
#> [[1]]
#> -1
#>
#> [[2]]
#> 1

## 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 #> [-Pi*d*h*lam/2 + Pi*d + Pi*h, -Pi*d**2*lam/4 + Pi*d, -Pi*d**2*h/4 + 1] And find the critical points: sol <- sympy$solve(grad, list(d, h, lam), dict = TRUE)
sol
#> [[1]]
#> [[1]]$d #> 2**(2/3)/Pi**(1/3) #> #> [[1]]$h
#> 2**(2/3)/Pi**(1/3)
#>
#> [[1]]$lam #> 2*2**(1/3)*Pi**(1/3) #> #> #> [[2]] #> [[2]]$d
#> 2**(2/3)*(-1 + sqrt(3)*I)/(2*Pi**(1/3))
#>
#> [[2]]$h #> 2**(2/3)*(-1 + sqrt(3)*I)/(2*Pi**(1/3)) #> #> [[2]]$lam
#> -2**(1/3)*Pi**(1/3) - 2**(1/3)*sqrt(3)*I*Pi**(1/3)
#>
#>
#> [[3]]
#> [[3]]$d #> -2**(2/3)*(1 + sqrt(3)*I)/(2*Pi**(1/3)) #> #> [[3]]$h
#> -2**(2/3)*(1 + sqrt(3)*I)/(2*Pi**(1/3))
#>
#> [[3]]$lam #> -2**(1/3)*Pi**(1/3) + 2**(1/3)*sqrt(3)*I*Pi**(1/3) We take the one with the real solution: sol[[1]] #>$d
#> 2**(2/3)/Pi**(1/3)
#>
#> $h #> 2**(2/3)/Pi**(1/3) #> #>$lam
#> 2*2**(1/3)*Pi**(1/3)

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 #> expression(2^(2/3)/pi^(1/3)) eval(sol_d) #> [1] 1.083852 sol_h <- to_r(sol[[1]]$h)
sol_h
#> expression(2^(2/3)/pi^(1/3))
eval(sol_h)
#> [1] 1.083852

(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
#> $commutative #> [1] TRUE x <- sympy$symbols('x', positive = TRUE)
x$assumptions0 #>$positive
#> [1] TRUE
#>
#> $nonpositive #> [1] FALSE #> #>$zero
#> [1] FALSE
#>
#> $extended_nonnegative #> [1] TRUE #> #>$imaginary
#> [1] FALSE
#>
#> $negative #> [1] FALSE #> #>$infinite
#> [1] FALSE
#>
#> $extended_nonzero #> [1] TRUE #> #>$hermitian
#> [1] TRUE
#>
#> $extended_negative #> [1] FALSE #> #>$nonzero
#> [1] TRUE
#>
#> $nonnegative #> [1] TRUE #> #>$commutative
#> [1] TRUE
#>
#> $extended_real #> [1] TRUE #> #>$extended_nonpositive
#> [1] FALSE
#>
#> $extended_positive #> [1] TRUE #> #>$complex
#> [1] TRUE
#>
#> $real #> [1] TRUE #> #>$finite
#> [1] TRUE
eq <- sympy$parsing$sympy_parser$parse_expr("x**2 - 1", local_dict = list('x' = x)) sympy$solve(eq, x, dict = TRUE)
#> [[1]]
#> [[1]]$x #> 1 ## 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 #> x**3/3 - x grad <- sympy$derive_by_array(eq, x)
sympy$solve(grad, x, dict = TRUE) #> [[1]] #> [[1]]$x
#> 1