SymPy
directlyFirst we get the SymPy
object:
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):
The problem is a constrained optimisation problem, and we solve it by
a Lagrange multiplier, and therefore we introduce lam
(the
Lagrange multiplier):
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.)
x <- sympy$symbols('x')
x$assumptions0
#> $commutative
#> [1] TRUE
x <- sympy$symbols('x', positive = TRUE)
x$assumptions0
#> $positive
#> [1] TRUE
#>
#> $nonpositive
#> [1] FALSE
#>
#> $imaginary
#> [1] FALSE
#>
#> $commutative
#> [1] TRUE
#>
#> $complex
#> [1] TRUE
#>
#> $extended_nonzero
#> [1] TRUE
#>
#> $extended_nonnegative
#> [1] TRUE
#>
#> $zero
#> [1] FALSE
#>
#> $extended_real
#> [1] TRUE
#>
#> $extended_negative
#> [1] FALSE
#>
#> $extended_positive
#> [1] TRUE
#>
#> $negative
#> [1] FALSE
#>
#> $nonzero
#> [1] TRUE
#>
#> $infinite
#> [1] FALSE
#>
#> $real
#> [1] TRUE
#>
#> $hermitian
#> [1] TRUE
#>
#> $finite
#> [1] TRUE
#>
#> $nonnegative
#> [1] TRUE
#>
#> $extended_nonpositive
#> [1] FALSE
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