#### Kombinatorika ####################################################################################

choose(n, k)   # kombinace
choose(3,2)
lchoose(n, k)  # log kombinace
factorial(x)   # faktorial x!
factorial(5)
lfactorial(x)  # faktorial log x!

perm = function(n, x){ factorial(n) / factorial(n-x)}
comb = function(n, x){ factorial(n) / factorial(n-x) / factorial(x)}


## kartezsky soucin mnozin
expand.grid(c(1, 2, 3), c(4, 5, 6))
# https://stackoverflow.com/questions/4309217/cartesian-product-data-frame


my_list <- c('one','two','three','four','five')
combn(my_list, 4, simplify=F)
combn(my_list, 4, simplify=T)
combn(c(1, 2, 3, 4, 5, 6), m=2, FUN = NULL, simplify = TRUE)
combn(5,3)


library(combinat)
permn(3)
length(permn(3))
combn(3, 2)
dim(combn(3,2))[2]


library(gtools)
n=5
r=3
combinations(n, r, v=1:n, set=TRUE, repeats.allowed=FALSE)
permutations(n, r, v=1:n, set=TRUE, repeats.allowed=FALSE)


### Urna obsahuje 3 koule ruzných barev, kolika zpusoby (s vracením) lze vybrat 2 koule z urny ?
x <- c('red', 'blue', 'black')
permutations(n=3,r=2,v=x,repeats.allowed=T)
nrow(permutations(n=3,r=2,v=x,repeats.allowed=T))
library(prob)
nsamp(3, 2, replace = TRUE, ordered = TRUE)
urnsamples(1:3, size = 2, replace = TRUE, ordered = TRUE)
nrow(urnsamples(1:3, size = 2, replace = TRUE, ordered = TRUE))


### Je treba zvolit ctyrmístný PIN. Kolik je mozno vytvorit ruznych kombinací císlic za predpokladu, ze se císlice nebudou opakovat?
n=10
r=4
permutations(n=10,r=4,v=1:n,repeats.allowed=F)
nrow(permutations(n=10,r=4,v=1:n,repeats.allowed=F))


### Kolik je trojciferných císel v nichz se zádná císlice neopakuje ?
n=10
r=3
mpm = permutations(n=10,r=3,v=0:9,repeats.allowed=F)
nrow(mpm)
which(mpm[,1]==0)
mpm = mpm[-which(mpm[,1]==0),]
nrow(mpm)


### Urcete celkový pocet kladných trojciferných císel, která obsahují císlici 6.
n=10
r=3
mpm = permutations(n=10,r=3,v=0:9,repeats.allowed=T)
nrow(mpm)
which(mpm[,1]==0)
mpm = mpm[-which(mpm[,1]==0),]
nrow(mpm)
length(unique(unlist(lapply(c(1:3),function(xx) which(mpm[,xx]==6)))))


### Kolik ruzných tetiv muze být vedeno 20 body po obvodu kruhu?
n=20
r=2
combinations(20, r=2, v=1:n, set=TRUE, repeats.allowed=FALSE)
choose(20, 2) 

### V sérii 12 výrobku jsou práve 3 vadné. Kolika zpusoby z nich lze vybrat 
# 6 libovolných výrobku
n=12
r=6
c1 = combinations(12, r=6, v=1:n, set=TRUE, repeats.allowed=FALSE)
cc1 = length(c1[,1]); cc1
choose(12, 6) 
# 6 výrobku bezvadných
c2 = combinations(9, r=6, v=1:n, set=TRUE, repeats.allowed=FALSE)
cc2 = length(c2[,1])
choose(9, 6) 
# 6 výrobku z nichz práve 1 je vadný
c3a = combinations(3, r=1, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c3a[,1])
c3b = combinations(9, r=5, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c3b[,1])
cc3 = length(c3a[,1])*length(c3b[,1])
choose(3, 1) 
choose(9, 5) 
# 6 výrobku z nichz práve 2 jsou vadné
c4a = combinations(3, r=2, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c4a[,1])
c4b = combinations(9, r=4, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c4b[,1])
cc4 = length(c4a[,1])*length(c4b[,1])
choose(3, 2) 
choose(9, 4) 
# 6 výrobku z nichz práve 3 jsou vadné
c5a = combinations(3, r=3, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c5a[,1])
c5b = combinations(9, r=3, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c5b[,1])
cc5 = length(c5a[,1])*length(c5b[,1])
choose(3, 3) 
choose(9, 3) 
# check - dle kombinatorického pravidla souctu
(cc2 + cc3 + cc4 + cc5) == cc1



### Pravdepodobnost
library(prob)
S <- rolldie(times = 3, makespace = TRUE)
Prob(S, X1+X2 > 9 )
Prob(S, X1+X2 > 9, given = X1+X2+X3 > 7)

#### mnoziny ####################################################################################

(x <- c(sort(sample(1:20, 3))))
(y <- c(sort(sample(3:23, 7))))

# x je prvkem mnoziny y
is.element(x, y) 
x %in% y 
library(prob) # pro cely vektor
isin(x,y, ordered = FALSE)
all(x %in% y)

library(prob)
x <- c(3,3,2,2,3,3,4,4)
isrep(x)  # one pair each of 2s and 4s
countrep(x)
isrep(x, nrep = 4)
countrep(x, nrep = 4)
isrep(x, vals = 4) # one pair of 4s
countrep(x, vals = 4) 

union(x, y) # sjednoceni
intersect(x, y) # prunik
setdiff(x, y) # rozdil mnozin
setdiff(y, x) # rozdil mnozin
setequal(x, y) # shoda mnozin
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/sets.html
subset(x, .)


library(prob)
union(x, y)
setdiff(x, y)
intersect(x, y)
subset(x, ...)

table()



### Vennovy diagramy
library(gplots)
venn( list(A=1:5,B=4:6,C=c(4,8:10)) )
v.table<-venn( list(A=1:5,B=4:6,C=c(4,8:10),D=c(4:12)) ) 
print(v.table)




### numericke metody ###########################################################################################################################

library(animation) # animation
bisection.method(FUN = function(x) x^2 - 4, rg = c(-1, 10), tol = 0.001, interact = FALSE)
newton.method(FUN = function(x) x^2 - 4, init = 10, rg = c(-1, 10), tol = 0.001, interact = FALSE, col.lp = c("blue", "red", "red"))

library(cmna) # numerical
bisection(f = function(x) x^2 - 4, a=-1, b=10, tol = 0.001, m = 100)
newton(f = function(x) x^2 - 4, fp, x, tol = 0.001, m = 100)



### polynomy ####

library(polynom)
(p <- polynomial(c(1,0,0,3,4)))
str(p)
p + polynomial(1:2)
p*p
f1 <- as.function(p)
f1(pi)
f1(matrix(1:6,2,3))


library(multipol)
(a <- as.multipol(matrix(1:10,nrow=2)))
b <-  as.multipol(matrix(1:10,ncol=2))
a+b
a*b




### kvadraticke a kubicke rovnice ###

# Descartova veta: Pocet kladnych korenu polynomu a0*x^n + a1*x^n-1 + ... + an-1*x + an = 0 je bud roven poctu znamenkovych zmen v posloupnosti a0, a1, ..., an-1, an jeho koeficientu, nebo je o sudy pocet mensi. Pocet zapornych korenu se urci obdobne po dosazeni -x za x.
sum(diff(sign(dc)) != 0)

library(sfsmisc)
nr.sign.chg(dc)

# Horner
library(sfsmisc)
polyn.eval(dc, x)
library(cmna)
horner(x, coefs)
rhorner(x, coefs)
naivepoly(x, coefs)
betterpoly(x, coefs)

### Velikost vyslednice dvou navzajem kolmych sil je 34 N. Jaka je velikost skladanych sil, je-li jedna z nich o 14 N vetsi nez druha?  
# Pythagorova veta: x^2 + (x + 14)^2 = 34^2      
library(Ryacas)
eq <- "x^2 + (x + 14)^2 - 34^2" 
yac_str(paste0("Simplify(", eq, ")"))  # zjednodusení výrazu

# diskriminant
dc = c(1,14,-480)
D = dc[2]^2 - 4*dc[1]*dc[3]
if (D == 0) {print("Kvadraticka rovnice ma dva sobe rovne realne koreny (dvojnasobny koren).")}
if (D > 0) {print("Kvadraticka rovnice ma dva ruzne realne koreny.")}
if (D < 0) {print("Kvadraticka rovnice nema zadny koren v oboru realnych cisel. V oboru komplexnich cisel ma dva imaginarni komplexne sdruzene koreny.")}

# Descartes rule
sum(diff(sign(dc)) != 0)
cat("Pocet kladnych korenu:", sum(diff(sign(dc)) != 0)) # pocet kladnych korenu
library(sfsmisc)
nr.sign.chg(dc)
cat("Pocet kladnych korenu:", nr.sign.chg(dc)) # pocet kladnych korenu

fun <- function (x) x^2 + (x + 14)^2 - 34^2
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9)$root 
F1
F2 = F1 + 14
F2

fun <- function (x) x^2 + 14*x - 480
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9)$root 
F1
F2 = F1 + 14
F2

rr = polyroot(c(-480, 14, 1)) # base
Re(rr)

library(pracma)
rr = roots(c(1, 14, -480))
rr
F1 = rr[rr>=0]
F1
F2 = F1 + 14
F2

# pro kontrolu lze pouzit Vietovy vzorce
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)$root
(-dc[2]/dc[1]) - F1
(dc[3]/dc[1])/F1
# nebo
library(pracma)
Poly(c(16, -30))

## graficke reseni
xx = seq(-50,50,1)
yy = xx^2
zz = -14*xx + 480
plot(xx,yy,type="l")
abline(480,-14,type="l",col=2)
plot(xx,yy,type="l",xlim = c(10,20),ylim=c(0,500))
abline(480,-14,type="l",col=2)

yy = abs(xx^2 + (xx + 14)^2 - 34^2)
plot(xx,yy,type="l")
abline(h=0,type="l",col=2)
xx[which(yy==0)]
xx[yy==0]


library(ssanv)
uniroot.integer(f, interval, lower = min(interval), upper = max(interval),step.power = 6, step.up = TRUE, pos.side = FALSE, print.steps = FALSE,maxiter = 1000)


### Jake pH ma roztok kyseliny mravenci o koncentraci 8.5 x 10-4 mol/l?
pKA = 3.752
KA = 10^-pKA; KA
Kw = 10^-14
cA = 8.5e-4 # [mol/l ] 

# H = sqrt(KA*cA)
H <- sqrt(KA*cA)
pH = -log10(H); pH 

# [H+]^2 + KA*[H+] + KA*cA = 0
CE = c(1,KA,-KA*cA)
fun <- function (x) CE[1]*x^2 + CE[2]*x + CE[3]
uniroot(fun, c(-1e-1, 0),tol = 1e-9)
uniroot(fun, c(0, 1e-1),tol = 1e-9)
H <- uniroot(fun, c(0, 1e-1),tol = 1e-9)$root  # base
pH = -log10(H); pH

library(rootSolve)
rootSolve::uniroot.all(fun, c(-1e-1, 0), tol = 1e-9)
rootSolve::uniroot.all(fun, c(0, 1e-1), tol = 1e-9)
H <- rootSolve::uniroot.all(fun, c(-1e-1, 1e-1), tol = 1e-9)
pH = -log10(H); pH 

library(Rmpfr)
Rmpfr::unirootR(fun, lower=-1e-1, upper=0, tol = 1e-9)
Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)
H <- Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)$root
pH = -log10(H); pH 

# [H+]^3 + KA*[H+]^2 - [H+]*(KA*cA + Kw) - KA*Kw = 0
CE = c(1,KA,-(KA*cA + Kw),-KA*Kw)

library(RConics) 
x0 = cubic(CE); x0
H = Re(x0[which(Im(x0)==0)])
H = Re(x0[which(Re(x0)>=0)])
pH = -log10(H); pH 

fun <- function (x) CE[1]*x^3 + CE[2]*x^2 + CE[3]*x + CE[4]

uniroot(fun, c(-1e-1,0),tol = 1e-9)
uniroot(fun, c(0, 1e-1),tol = 1e-9)
H <- uniroot(fun, c(0, 1e-1),tol = 1e-9)$root  # base
pH = -log10(H); pH 

library(rootSolve)
rootSolve::uniroot.all(fun, c(-1e-1, 0),tol = 1e-9)
rootSolve::uniroot.all(fun, c(0,1e-1),tol = 1e-9)
H <- rootSolve::uniroot.all(fun, c(-1e-1, 1e-1),tol = 1e-9)
pH = -log10(H); pH 

library(Rmpfr)
Rmpfr::unirootR(fun, lower=-1e-1, upper=0, tol = 1e-9)
Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)
H <- Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)$root
pH = -log10(H); pH 


### Tlakova lahev s oxidem uhlicitym obsahuje 10.0 kg plynu. Jaký objem zaujíma stlaceny plyn, kdyz pri teplote 30 °C je tlak v lahvi 13.17e6 Pa? Vypocet provedte pomoci van der Waalsovy rovnice
# (p + a/Vm^2)(Vm - b) = R*T     [Vm = 0.075 m3/kmol, n = 0.2273 kmol, V = 0.0171 m3]
library(measurements)
p = 13.17e6 # [Pa]
m = conv_unit(10.0, from="kg", to="g")
t = conv_unit(30, from="C", to="K") 
R = 8.3141  # [J / mol K]
R = R*1000  # [J / kmol K]
Mr = 44.01  # [g/mol]
n = m/Mr    # [mol]
n = n/100   # [kmol]

# vypocet pro neidealni plyn
a = 0.365e6 # [m6 Pa / kmol2]
b = 0.0428  # [m3 / kmol]
# Vm**3 - (b + R*t/p)*Vm**2 + (a/p)*Vm - a*b/p 

# realne i komplexni koreny
library(RConics) 
CE = c(1,-(b + R*t/p),a/p,-(a*b/p))
x0 = cubic(CE); x0
Vm = Re(x0[which(Im(x0)==0)])  # [m3 / kmol]
V = Vm*n; V  # [m3]

# jen realne koreny
fun <- function (x) CE[1]*x^3 + CE[2]*x^2 + CE[3]*x + CE[4]
#curve(fun(x),-4,4)
#abline(h = 0, lty = 3)
Vm <- uniroot(fun, c(-4, 4))$root  # base
V = Vm*n; V  # [m3]

library(rootSolve)
Vm <- rootSolve::uniroot.all(fun, c(-4, 4))
V = Vm*n; V  # [m3]

library(Rmpfr)
Vm <- Rmpfr::unirootR(fun, lower=-4, upper=4, tol = 1e-9)$root
V = Vm*n; V  # [m3]

## graficke reseni
# CE = c(1,-(b + R*t/p),a/p,-(a*b/p))
xx <- seq(-4,4, by=0.01)
yy <- CE[1]*xx^3 
zz <- -CE[2]*xx^2 - CE[3]*xx - CE[4]
plot(xx,yy,type="l", col=2)
points(xx,zz,type="l", col=4)
plot(xx, yy, type="l",col=2,xlim=c(0.05,0.1),ylim=c(0,0.001))
points(xx, zz, type="l",col=4)


# vypocet pro idealni plyn
# p*Vm - R*T = 0  
Vm = R*t/p
V = Vm*n; V  # [m3]



### soustavy linearnich rovnic ###

### Ze dvou slitin s 60% a 80% obsahem medi se ma ziskat 40 kg slitiny se 75% obsahem medi. Kolik kg kayde slitiny je treba pouzit?  [10 a 30 kg]

library(rootSolve)
model <- function(x){
  F1 <- 0.6*x[1] + 0.8*x[2] - 30
  F2 <- x[1] + x[2] - 40
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 0.6*m1 + 0.8*m2 = 40*0.75 => m1 = 40*0.75/0.6 - 0.8/0.6 * m2  => m1 = 50 - 1.33 * m2 
# m1 + m2 = 40  => m1 = 40 - m2
xx = seq(0,50,0.1)
yy = 50 - 1.33*xx
zz = 40 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(25,35),ylim=c(0,20))
points(xx, zz, type="l",col=4)


### Do bazenu natece prutokem A za 3 h a prutokem B za 4 h celkem 2150 hl vody. Prutokem A za 4 h a prutokem B za 2 h by nateklo 1700 hl vody. Kolik hl vody natece prutokem A a kolik prutokem B za 1 hodinu?      A 250 hl, B 350 hl

library(rootSolve)
model <- function(x){
  F1 <- 3*x[1] + 4*x[2] - 2150
  F2 <- 4*x[1] + 2*x[2] - 1700
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 3*m1 + 4*m2 = 2150 => m2 = 2150/4 - 3/4 * m1  => m2 = 537.5 - 0.75 * m1
# 4*m1 + 2*m2 = 1700 => m2 = 1700/2 - 4/2 * m1 => m2 = 850 - 2 * m1
xx = seq(0,500,1)
yy = 537.5 - 0.75*xx
zz = 850 - 2*xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(200,300),ylim=c(300,400))
points(xx, zz, type="l",col=4)


### Ze dvou druhu caje v cene 160 Kc a 220 Kc za 1 kg je treba pripravit 20 kg smesi v cene 205 Kc za 1 kg. Kolik kg kazdeho caje je treba smichat?   5 kg levnejsiho a 15 kg drazsiho 

library(rootSolve)   ### opravit zadani
model <- function(x){
  F1 <- 160*x[1] + 220*x[2] - 4100
  F2 <- x[1] + x[2] - 20
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 160*m1 + 220*m2 = 4100 => m1 = 4100/160 - 220/160 * m2  => m1 = 25.625 - 1.375 * m2
# m1 + m2 = 20 => m1 = 20 - m2
xx = seq(0,20,0.1)
yy = 25.625 - 1.375*xx
zz = 20 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(10,20),ylim=c(0,10))
points(xx, zz, type="l",col=4)


### Kolik g 60% a kolik g 30% roztoku NaCl je treba smichat pri priprave 100 g 40% roztoku? 20 g 60% a a 80 g 35%

library(rootSolve)
model <- function(x){
  F1 <- 0.6*x[1] + 0.3*x[2] - 40
  F2 <- x[1] + x[2] - 100
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 0.6*m1 + 0.3*m2 = 100*0.4 => m1 = 40/0.6 - 0.3/0.6 * m2  => m1 = 66.67 - 0.5 * m2 
# m1 + m2 = 100  => m1 = 100 - m2
xx = seq(0,100,0.1)
yy = 66.67 - 0.5*xx
zz = 100 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(65,70),ylim=c(30,40))
points(xx, zz, type="l",col=4)

### matice

B = c(0.40*100,100) # [mg]
names(B) = c("60%","30%")
r1 = c(0.60,1) # [mg]
r2 = c(0.30,1) # [mg]
A = cbind(r1,r2)
colnames(A) = c("60%","30%")
rownames(A) = c("r30","r3")
det(A) # matice je regularni n = h
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
#
library(limSolve)
G <-matrix(ncol = 2, byrow = TRUE, data = c(1, 0, 0, 1)) 
H <- c(0, 0)
ldei(A, B, G = G, H = H)$X
#
library(cmna) 
gdls(A, B, alpha = 0.05, tol = 1e-06, m = 1e+05) #  least squares with graident descent
jacobi(A, B, tol = 1e-06, maxiter = 100)  # iterativematrix
gaussseidel(A, B, tol = 1e-06, maxiter = 100) # iterativematrix
# cgmmatrix(A, B, tol = 1e-06, maxiter = 100) # nefunguje
solvematrix(A, B) # refmatrix


### Ze dvou kovu o hustotach 7.4 g/cm3 a 8.2 g/cm3 je treba pripravit 0.5 kg slitiny o hustote 7.6 g/cm3. Kolik g kazdiho z kovu je k tomu potreba? 375 g lehciho a 125 g tezsiho

library(rootSolve)
model <- function(x){
  F1 <- 7.4*x[1] + 8.2*x[2] - 3800
  F2 <- x[1] + x[2] - 500
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root  # [kg]

### graficke reseni
# 7.4*m1 + 8.2*m2 = 0.5*7.6 => m1 = 3800/7.4 - 8.2/7.4 * m2  => m1 = 513.5 - 1.108 * m2 
# m1 + m2 = 500  => m1 = 500 - m2
xx = seq(0,500,1)
yy = 513.5 - 1.108*xx
zz = 500 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(100,150),ylim=c(350,400))
points(xx, zz, type="l",col=4)

### matice
B = c(7.6,1) # [mg]
names(B) = c("kov1","kov2")
r1 = c(7.4,1) # [mg]
r2 = c(8.2,1) # [mg]
A = cbind(r1,r2)
colnames(A) = c("kov1","kov2")
rownames(A) = c("7.4","8.2")
det(A) # matice je regularni n = h
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
rr = matlib::Solve(A, B)
read.table(text = rr[1], fill = TRUE)[[3]]*500 # [g]
read.table(text = rr[2], fill = TRUE)[[3]]*500 # [g]
rr = limSolve::Solve(A, B)
rr*500 # [g]
# 
library(limSolve)
G <-matrix(ncol = 2, byrow = TRUE, data = c(1, 0, 0, 1)) 
H <- c(0, 0)
rr = ldei(A, B, G = G, H = H)$X
rr*500 # [g]



### Vodní nádrz by se naplnila prvním prívodem za 36 minut, druhým za 45 minut. Za jak dlouho se nádrz naplní, pritéká-li voda nejprve 9 minut prvním prívodem a pak obema soucasne?
# x/36 + (x - 9)/45 = 1
vr <- ysym("x/36 + (x-9)/45 - 1") 
solve(vr, "x") 

fun <- function (x) x/36 + (x-9)/45 - 1
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)$root


### V tepelné elektrárne mají zásobu uhlí, která vystací na 24 dní, bude-li v cinnosti pouze první blok, na 30 dní, bude-li v provozu pouze 2. blok a na 20 dní, bude-li v provozu pouze 3. blok. Na jak dloho vystací zásoba, budou-li v provozu vsechny bloky najednou?
# x/24 + x/30 + x/20 = 1
vr <- ysym("x/24 + x/30 + x/20 - 1") 
solve(vr, "x") 

fun <- function (x) x/24 + x/30 + x/20 - 1
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)$root

### Míc byl vyhozen svisle vzhuru rychlostí 25 m/s. Za jak dlouho bude míc 20 m nad zemí? Odpor vzduchu zanedbejte.
# h = v t - 1/2 g t^2
g = 10 # [m/s^2]
h = 20 # [m]
v = 25 # [m/s]
t <- ysym("t")
g <- ysym("g")
h <- ysym("h")
v <- ysym("v")
poly <- ysym("g*t^2 - 2*v*t + 2*h") 
ko = solve(poly, "t")   
ko[1]
ko[2]
eval(as_r((2*v+sqrt(-((-4)*v^2+8*g*h)))/(2*g)), list(g = 10, v = 25, h = 20))
eval(as_r((2*v-sqrt(-((-4)*v^2+8*g*h)))/(2*g)), list(g = 10, v = 25, h = 20))

poly <- ysym("10*t^2 - 50*t + 40") 
ko = solve(poly, "t")   
ko[1]
ko[2]

library(pracma)
rr = roots(c(10, -50, 40))
rr
F1 = rr[rr>=0]
F1

rr = polyroot(c(40, -50, 10)) # base
Re(rr)

g = 10 # [m/s^2]
h = 20 # [m]
v = 25 # [m/s]
fun <- function (t) g*t^2 - 2*v*t + 2*h
#uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0.1, 5),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)$root 
F1 # udava jen jeden koren

library(rootSolve)
F1 <- rootSolve::uniroot.all(fun, c(0, 10), tol = 1e-9)
F1 

library(Rmpfr)
F1 <- Rmpfr::unirootR(fun, lower=0.1, upper=10, tol = 1e-9)
F1 # udava jen jeden koren

# Descartes rule
cat("Pocet kladnych korenu:", sum(diff(sign(c(10, -50, 40))) != 0)) # pocet kladnych korenu


###### derivace ####################################################################################################

D() 
deriv() 

derivative = deriv(~ x^3, "x") 
x <- 2
eval(derivative )

f = expression(x^3)
dx2x <- D(f,"x")
dx2x 


library(pracma)
f <- sin
xs <- seq(-pi, pi, length.out = 100)
ys <- f(xs)
y1 <- fderiv(f, xs, n = 1, method = "backward")
y2 <- fderiv(f, xs, n = 2, method = "backward")
y3 <- fderiv(f, xs, n = 3, method = "backward")
y4 <- fderiv(f, xs, n = 4, method = "backward")
plot(xs, ys, type = "l", col = "gray", lwd = 2,xlab = "", ylab = "", main = "Sinus and its Derivatives")
lines(xs, y1, col=1, lty=2)
lines(xs, y2, col=2, lty=3)
lines(xs, y3, col=3, lty=4)
lines(xs, y4, col=4, lty=5)
grid()

f <- sin
xs <- seq(-pi, pi, length.out = 100)
numdiff(f, xs, maxiter = 16, h = 1/2)


# derivace v bode
library(pracma)
f <- function(x) sin(x)*sqrt(1+sin(x))
F <- function(x)integrate(f, 0, x, rel.tol = 1e-12)$value
x0 <- 1
(dF0 <- numderiv(F, x0, tol = 6.5e-15))      
f(x0)
x=x0
eval(sin(x)*sqrt(1+sin(x)))

fderiv(F, x0)

library(numDeriv)
numDeriv::grad(F, x0)



### Hmotny bod kona harmonicky pohyb, zavislost vychylky y na case t je dana vztahem y = A * sin(om*t + phi), kde A, om a phi jsou konstanty. Urcete rychlost hmotneho bodu v case t = 0.
library(Deriv)
der = Deriv("A*sin(om*t+phi)", "t")
der
der0 = gsub("t", "0", der, fixed = TRUE)
der0
Simplify(der0)




library(mosaic)

f <- makeFun(m * x + b ~ x, m = 3.5, b = 10)
f(x = 2)

g <- makeFun(A * x * cos(pi * x * y) ~ x + y, A = 3)
# function (x, y, A = 3) A * x * cos(pi * x * y)
g
g(x = 1, y = 2)

plotFun(A * exp(k * t) * sin(2 * pi * t/P) ~ t + k, t.lim = range(0, 10), k.lim = range(-0.3,0), A = 10, P = 4)

library(manipulate)
plotFun(A * exp(k * t) * sin(2 * pi * t/P) ~ t + k, t.lim = range(0, 10),k.lim = range(-0.3,0), A = 10, P = 4, surface = TRUE)



library(mosaic)
library(mosaicCalc)
D(sin(x) ~ x)   # derivace
D(A * x^2 * sin(y) ~ x)  # 2. derivace
D(A * x^2 * sin(y) ~ x + y)  # 3. derivace


findZeros(sin(t) ~ t, t.lim = range(-5, 1))
findZeros(sin(t) ~ t, nearest = 5, near = 10)
findZeros(sin(t) ~ t, near = 0, within = 8)


D = function(f,delta=.000001){function(x){ (f(x+delta) - f(x-delta))/(2*delta)} }
f = function(x){ x^2 + 2*x }
plot(f, 0, 10)



###### integrace ####################################################################################################

# integrate() 

# 1/((x+1)*sqrt(x))
integrand <- function(x) {1/((x+1)*sqrt(x))}
integrate(integrand, lower = 0, upper = Inf)

# 1/sqrt(2*pi)*exp(-x^2/2)
f <- function(x) {1/sqrt(2*pi)*exp(-x^2/2)}
integrate(f, lower = -1.96, upper = 1.96)


### mnohorozmerny integral

library(cubature)
f <- function(x) { 2/3 * (x[1] + x[2] + x[3]) }  
adaptIntegrate(f, lowerLimit = c(0, 0, 0), upperLimit = c(0.5, 0.5, 0.5))



library(sfsmisc)
integrate.xy(x, fx=y, a, b, use.spline=TRUE, xtol=2e-08)


### geometrie ##########################################################################################################################################

library(retistruct)
line.line.intersection(c(0, 0), c(1, 1), c(0, 1), c(1, 0))
line.line.intersection(c(0, 0), c(0, 1), c(1, 0), c(1, 1))


library(PBSmapping)
p1 <- data.frame(PID=rep(1, 4), POS=1:4, X=c(1,1,6,6), Y=c(1,3,3,1))
p2 <- data.frame(PID=rep(2, 5), POS=1:5, X=c(4,4,8,8,6), Y=c(2,4,4,2,1))
p3 <- joinPolys(p1,p2)
x11()
par(mar=c(3,3,1,1))
plot(1,1,ylim=c(0,5),xlim=c(0,9), t="n", xlab="", ylab="")
polygon(p1$X, p1$Y, border=2)
polygon(p2$X, p2$Y)
polygon(p3$X, p3$Y, col=rgb(0,0,1,0.2))


### kuzelosecky 

library(sfsmisc)
ellipsePoints(a, b, alpha = 0, loc = c(0, 0), n = 201)

library(conics)
# 2x2^2 + 2x1x2 + 2x2^2 - 20x1 - 28x2 + 10 = 0 
v <- c(2,2,2,-20,-28,10)
A <- conicMatrix(v)
conicCenter(v)
conicAxes(v)
conicCenter(A)
conicAxes(A)
conicPlot(v, main="conicPlot(v)", xlab="", ylab="")
conicPlot(v, col="red", lty=3)
conicPlot(v, asp=1)
conicPlot(v, center=T, sym.axes=T)

v <- c(2,2,2,-20,-28,10)
conicPlot(v, center=T, lty=1)
v[6] <- 30 
conicPlot(v, add=T, lty=2)
v[6] <- 50 
conicPlot(v, add=T, lty=4)

v <- c(2,2,-2,-20,20,10) 
conicPlot(v, xlim=c(-10,10), ylim=c(-5,18))
conicPlot(v, asymptotes=T, sym.axes=T,as.col="red", as.lty=2, ax.col="blue", ax.lty=4, xlim=c(-10,10), ylim=c(-5,18))
conicPlot(v, asymptotes=T, sym.axes=T, xlim=c(-10,10), ylim=c(-5,18), asp=1, col="blue", main="Hyperbola", bty="n")

v <- c(-4,0,1,0,0,1) 
cp <- conicPlot(v, sym.axes=TRUE, asymptote=TRUE, asp=1, ax.lty=2, as.col="gray") 
points(cp$foci$x,cp$foci$y,col="red",pch=19) 
text(cp$foci$x,cp$foci$y+0.1,paste("F",2:1,sep="")) 
points(cp$vertices$x,cp$vertices$y,col="blue",pch=19)



#### SYMBOLICKA MATEMATIKA ##################################################################


library(rSymPy)

# výpocet na daný pocet platných císlic
sympy("sqrt(8).evalf()")  # default
sympy("sqrt(8).evalf(50)")
sympy("pi.evalf()")
sympy("pi.evalf(3)")
sympy("pi.evalf(120)")


## Recognizing numbers
sympy("nsimplify(4.242640687119286)")
sympy("nsimplify(cos(pi/6))")
sympy("nsimplify(6.28, [pi], tolerance=0.01)")
sympy("nsimplify(pi, tolerance=1e-5)")
sympy("nsimplify(pi, tolerance=1e-6)")
sympy("nsimplify(29.60881, constants=[pi,E], tolerance=1e-5)")


# 
x <- Var("x")
x+x
x*x/x
y <- Var("x**3")
x/y

z <- sympy("2.5*x**2")
z + y

sympy("one = cos(1)**2 + sin(1)**2")
sympy("(one - 1).evalf()")  # chyba zaokrouhlení
sympy("(one - 1).evalf(chop=True)")  # rouding this type of roundoff errors


# FUNKCE
sympy("sqrt(8)")
sympy("sqrt(9)")  
sympy("sqrt(3.14)")  
sympy("sin(1)")  
# create function 
Exp <- function(x) Sym("exp(", x, ")") 
Exp(-x) * Exp(x) 


# rovnice
sympy("Eq(x**2+2*x+1,(x+1)**2)") # vytvoreni rovnice
# jsou obe strany rovnice ekvivalentni?
sympy("a = x**2+2*x+1")
sympy("b = (x+1)**2")
"0" == sympy("simplify(a-b)")  # if they are equal, the result is zero

# urceni korenu rovnice
sympy("solve(x**2 - 2, x)")  # solve x^2-2=0


# zjednoduseni vyrazu
sympy("simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))")
sympy("factor(x**3 - x**2 + x - 1)")
sympy("(x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)")

x <- Var("x")
y <- Var("y")
z <- Var("z")
sympy("collect(x*y + x - 3 + 2*x**2 - z*x**2 + x**3, x)")  # organize equation around var 'x'
sympy("cancel((x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1))")
sympy("simplify((x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1))")


# roznasobeni zavorek
sympy("(x + 2)*(x - 3)")
sympy("expand((x + 2)*(x - 3))")
sympy("expand_func(gamma(x + 3))") # pro funkci


# limita funkce
sympy("limit(1/x, x, oo)")   # limit eg, x -> Inf
sympy("limit(1/x, x, 0)") 
sympy("limit(sin(2*x)/(2*x), x, 0)") 
sympy("limit((1+(1/(3*x)))**(3*x), x, oo)") 
sympy("limit((exp(x)-E)/(x-1), x, 1)") 
sympy("limit((ln(x)-1)/(x-E), x, E)") 
sympy("limit((exp(x)-exp(-x))/sin(x), x, 0)") 


### Pri redení daného mnozství kyseliny sírové pridáním vody o hmotnosti x vznikne mnozství tepla y = Ax / x+B , kde A, B jsou kladné konstanty. 
# Jaké mnozství tepla se uvolní pro x -> oo ? Kolik vody se musí pridat aby vzniklo 80 % tohoto mnozství tepla ?
x = Var("x")
A = Var("A")
B = Var("B")
sympy("limit(A*x/(x+B), x, oo)")

        
### Okamzitá koncentrace látky C, vzniklé pri nevratné chemické reakci A + B = C, je c(t) = ab(1-e^(b-a)kt) / a-bee^(b-a)kt , kde a, b, k jsou kladné konstanty, t > 0.
# Jaké hodnoty nabývá koncentrace c pro t -> oo v prípadech a > b, b > a. 
a = Var("a")
b = Var("b")
k = Var("k")
t = Var("t")
sympy("limit(a*b*(1-exp((b-a)*k*t))/(a-b*exp((b-a)*k*t)), t, oo)") 


### Molární tepelná kapacita cV dvouatomového ideálního plynu pri konstantním objemu V je dána vztahem cV(T) = 5R/2 + (hv/kT)^2 (R exp(hv/kT))/(exp(hv/kT)-1)^2
R = Var("R")
h = Var("h")
v = Var("v")
k = Var("k")
Tk = Var("Tk")
sympy("limit(R*5/2 + ((h*v/(k*Tk))**2)*(R*exp(h*v/(k*Tk))/(exp(h*v/(k*Tk))-1)**2), Tk, oo)") 
        

### Molární tepelná kapacita cV závisí na teplote T podle vztahu cV(T) = 3R(hv/kT)^2 (exp(hv/kT))/(exp(hv/kT)-1)^2
R = Var("R")
h = Var("h")
v = Var("v")
k = Var("k")
Tk = Var("Tk")
sympy("limit(3*R*((h*v/(k*Tk))**2)*(exp(h*v/(k*Tk))/(exp(h*v/(k*Tk))-1)**2), Tk, oo)") 
        
        
        
        
sympy("y = x*x") # create a variable 'y' in the SymPy persistant state
        
sympy("A = Matrix([[1,x], [y,1]])")
sympy("A**2")
        
sympy("B = A.subs(x,1.1)")  # replace x by 1.1 (btw, SymPy objects are immutable)
sympy("B**2")
        
# more replacement, a subexpression by another:
sympy("expr = sin(2*x) + cos(2*x)")
sympy("expr.subs(sin(2*x), 2*sin(x)*cos(x))")
sympy("expr.subs(x,pi/2)")
        
x <- Var("x")
y <- Var("y")
(x+y)*(x+y)
        

        
### Derivace ###
        
D(expression(sin(2*x)), "x") # also possible in base R
        
sympy("diff(sin(2*x), x, 1)")  # first derivative
sympy("diff(sin(2*x), x, 2)")  # second derivative
sympy("diff(sin(2*x), x, 3)")  # third  derivative
sympy("diff(exp(x*y*z), x, y, y)") # d^3/dxdy^2
sympy("diff(exp(x*y*z), x, z, 3)") # d^4/dxdz^3
        
library(Deriv)
Deriv("sin(2*x)", "x", cache.exp=FALSE) # differentiate only by x 
        
        
### Závislost mnozství m látky získané chemickou reakcí v case t je dána vztahem m = A(1-exp(-kt)), kde A a k jsou konstanty. Jak závisí rychlost reakce na case t ?
D(expression(A*(1-exp(-k*t))), "t") # base R
library(rSymPy)
sympy("diff(A*(1-exp(-k*t)), t, 1)")
library(Deriv)
Deriv("A*(1-exp(-k*t))", "t", cache.exp=FALSE)  
        
        
### Pri chemické reakci A + B = C v níz obe pocátecní koncentrace výchozích látek nabývají stejné hodnoty a, je závislost koncentrace x vznikající látky na case dána vztahem x(t) = (a^2)kt / 1+akt. Jak se mení koncentrace s casem, je-li k > 1 ? kt . 
D(expression((a^2)*k*t/(1+a*k*t)), "t") # base R
library(rSymPy)
sympy("diff((a**2)*k*t/(1+a*k*t), t, 1)")
library(Deriv)
Deriv("(a^2)*k*t/(1+a*k*t)", "t", cache.exp=FALSE)  
       
        

library(Deriv)
        Deriv("sin(x^2) * y", "x") # differentiate only by x 
        Deriv("sin(x^2) * y", "y") # differentiate only by x 
        Deriv("sin(x^2) * y", cache.exp=FALSE) # differentiate by all variables (here by x and y) 
        Deriv("sin(x^2) * y", cache.exp=TRUE) # differentiate by all variables (here by x and y) 
        Deriv("sin(cos(x + y^2))", cache.exp=FALSE) # differentiate by all variables (here by x and y) 
        Deriv("sin(cos(x + y^2))", "x")
        Deriv("sin(cos(x + y^2))", "y")
        # Compound function example (here abs(x) smoothed near 0) 
        fc <- function(x, h=0.1) if (abs(x) < h) 0.5*h*(x/h)**2 else abs(x)-0.5*h 
        Deriv("fc(x)", "x", cache.exp=FALSE) 
        
        Simplify("x^3 - x^2 + x - 1")
        Simplify("x*y + x - 3 + 2*x^2 - z*x^2 + x^3")
        Simplify("(x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)")
        
        
        
### Integral ###
        
integrate(function(x) exp(-x), 0, Inf)  # integration is possible in R
        
sympy("integrate(exp(-x))")    # indefinite integral
sympy("integrate(exp(-x*y),x)")  # indefinite integral  
        sympy("integrate(exp(-x), (x, 0, oo))") # definite integral
        sympy("integrate(x**2 - y, (x, -5, 5), (y, -pi, pi))")  # definite integral
        
        
        ### Mnozství tepla Q (v J) potrebné na zahrátí 10 kg zeleza o 1 K je pro T = <273.15, 473.15> vyjádrena funkcí Q = 4.1686 (105.3T + 0.42T^2). Vypocítejte strední mnozství tepla (v J / kg K) potrebného na zahrátí zeleza z teploty 293.15 K na 373.15 K.  
        qm = integrate(function(Tk) 4.1868*(105.3*Tk + 0.142*Tk^2), 293.15, 373.15)  # integration is possible in R
        unlist(qm)
        0.1*unlist(qm)$value/(373.15-293.15) # [J / kg K]
        
        Tk = Var("Tk")
        smp = sympy("integrate(4.1868*(105.3*Tk + 0.142*Tk**2), (Tk, 293.15, 373.15))")  # definite integral
        0.1*as.numeric(smp)/(373.15-293.15) # [J / kg K]
        
        
        ### Merné teplo c benzínu (pri konstantním tlaku) v závislosti na teplote T je vyjádrené funkcí c = -0.06 + 0.0010228T. Vypocítejte strední merné teplo benzínu pro T = <389.15, 491.15>.
        qm = integrate(function(Tk) -0.06 + 0.0010228*Tk, 389.15, 491.15)  # integration is possible in R
        unlist(qm)
        0.1*unlist(qm)$value/(491.15-389.15) # [J / kg K]
        
        Tk = Var("Tk")
        smp = sympy("integrate(-0.06 + 0.0010228*Tk, (Tk, 389.15, 491.15))") 
        as.numeric(smp)/(491.15-389.15) # [J / K]
        
        
##################################################################################################################
        
        library(Ryacas)
        # https://www.andrewheiss.com/blog/2019/02/16/algebra-calculus-r-yacas/
        
        yac_str("N(Pi, 50)")
        
        yac_str("x+x+x+x")
        yac_expr("x+x+x+x")
        eval(yac_expr("x+x+x+x"), list(x = 4))
        
        yac_str("Factor(x^2+x-6)")
        yac_expr("Factor(x^2+x-6)")
        
        eq <- "x^2 + 4 + 2*x + 2*x"
        yac_str(eq) # Evaluate yacas command x (a string) and get result as string/character. No task was given to yacas, so we simply get the same returned
        
        # zjednodusení výrazu
        yac_str(paste0("Simplify(", eq, ")"))  # Simplify(x): Simplify an expression
        
        eq <- "x^2 + (x + 14)^2 - 34^2" 
        rr = yac_str(paste0("Simplify(", eq, ")"))  # zjednodusení výrazu
        x = 2
        eval(parse(text="x^2+14*x-480"))
        
        yac_str(paste0("Factor(", eq, ")"))    # Factor(x): Factorise an expression
        yac_str(y_fn(eq, "Factor"))
        
        y_fn(eq, "Factor")
        yac_expr(paste0("Factor(", eq, ")"))
        
        zz = y_fn(eq, "Factor")
        yac_expr(paste0("Expand(", zz, ")"))
        
        #
        expr <- yac_expr(paste0("Factor(", eq, ")"))
        expr
        eval(expr, list(x = 2))
        
        
        x <- ysym("x")
        2*x^2 - 5
        c(-2, 5)*x
        c(2*x, -x^3)
        as_r(c(-2, 5)*x) # or yac_expr(c(-2, 5)*x)
        
        
        xs <- ysym("x")
        poly <- xs^2 - xs - 6
        poly
        zeroes <- solve(poly, "x") # Solve(x^2 - x - 6 == 0, x)
        zeroes
        
        solve(poly, 3, "x") # Solve(x^2 - x - 6 == 3, x)
        
        
        x <- ysym("x")
        y <- ysym("y")
        lhs <- c(3*x*y - y, x)
        rhs <- c(-5*x, y+4)
        sol <- solve(lhs, rhs, c("x", "y"))
        sol
        sol_vals <- lapply(seq_len(nrow(sol)), function(sol_no){y_rmvars(sol[sol_no, ])})
        sol_vals
        sol_envir <- lapply(sol_vals, function(l){list(x = as_r(l[1]), y = as_r(l[2]))})
        sol_envir
        do.call(rbind, lapply(seq_along(sol_envir), function(sol_no) {sol_val <- sol_envir[[sol_no]]
        data.frame(sol_no = sol_no,
                   eq_no = seq_along(sol_val),
                   lhs = eval(as_r(lhs), sol_val),
                   rhs = eval(as_r(rhs), sol_val))}))
        
        
        # Rosenbrock function:
        x <- ysym("x")
        y <- ysym("y")
        f <- (1 - x)^2 + 100*(y - x^2)^2
        f
        
        N <- 30
        x <- seq(-1, 2, length=N)
        y <- seq(-1, 2, length=N)
        f_r <- as_r(f)
        f_r
        z <- outer(x, y, function(x, y) eval(f_r, list(x = x, y = y)))
        levels <- c(0.001, .1, .3, 1:5, 10, 20, 30, 40, 50, 60, 80, 100, 500, 1000)
        cols <- rainbow(length(levels))
        contour(x, y, z, levels = levels, col = cols)
        g <- deriv(f, c("x", "y")) # minimum
        g
        crit_sol_all <- solve(g, c("x", "y"))
        crit_sol_all
        
        
        Solve(expr, var) solve an equation (refer to the Ryacas function y_rmvars())
        Variables(): List yacas variables
        # yacas( "OldSolve({a*x+y==0,x+z==0},{x,y})" )
        #  equation $a x = b$ for $x$.
        solve(a,b)   # ???nd roots of a for variable b, i.e. yacas Solve(a == 0,b) 
        solve(a,b,v) #  ???nd solutions to a == b for variable v, i.e. yacas Solve(a == b,v) 
        #  also works for a system of equations (when a is a vector)
        A <- outer(0:3, 1:4, "-") + diag(2:5) 
        a <- 1:4 
        B <- ysym(A) 
        b <- ysym(a) solve(A) 
        solve(B) 
        solve(A, a) 
        solve(B, b)
        poly <- ysym("x^2 - x - 6") 
        solve(poly, "x")   # Solve(poly == 0, x) 
        solve(poly, 3, "x") # Solve(poly == 3, x)
        
        
        sum(expr, var, lower, upper, ..., na.rm = FALSE)
        
        ## D(x) expr: Take the derivative of expr with respect to x
        
        yac_str("D(x) x^2+x-6")
        yac_expr("D(x) x^2+x-6")
        
        L <- ysym("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
        L
        deriv(L, "x")
        
        my_log <- function(x){return(sin(log(2 + x)))}
        deriv(my_log(x), x)
        integrate(my_log(x), x)
        integrate(f, var, lower, upper )
        
        # Hessian(function, var) # Create the Hessian matrix
        Hessian(L, "x")
        deriv(L, c("x", "y", "a"))
        H <- Hessian(L, c("x", "y", "a")); H
        as_r(H)
        eval(as_r(H), list(x = 2, y = 2, a = 2))
        # Jacobian(function, var) # Create the Jacobian matrix
        L2 <- ysym(c("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)", "x^3 + 4*a^2")) # just some function
        L2
        Jacobian(L2, "x")
        Jacobian(L2, c("x", "y", "a"))
        
        
        # lim(f, var, val, from_left, from_right) # Limit of f(n) for n going towards a (e.g. Infinity or 0)
        cmd <- "Limit(n, Infinity) (1+(1/n))^n"
        yac_str(cmd)
        yac_expr(cmd)
        
        yac_str("Limit(h, 0) (Sin(x+h)-Sin(x))/h") 
        
        
        # sum(k, a, b, f(k)) # Sum of f(k) for k from a to b # urcity integral
        yac_str("Sum(k, 0, n, a^k)")
        yac_expr("Sum(k, 0, n, a^k)")
        
        
        yac_str("poly := (x-3)*(x+2)")
        yac_silent("poly := (x-3)*(x+2)")
        yac_str("Expand(poly)")
        
        
      

        
        