---
title: "R Notebook"
output: html_notebook
---


Kvadraticke a kubicke rovnice


Velikost výslednice dvou navzajem kolmých sil je 34 N. Jaká je velikost skládaných sil, je-li jedna z nich o 14 N větší než druha?

```{r echo=FALSE, message=FALSE, warning=FALSE}

# 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) {cat("Kvadraticka rovnice ma dva sobe rovne realne koreny (dvojnasobny koren).")}
if (D > 0) {cat("Kvadraticka rovnice ma dva ruzne realne koreny.")}
if (D < 0) {cat("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

# Hornerovo schema
library(sfsmisc)
polyn.eval(dc, x)


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)
# uniroot(fun, c(-1, 1),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9) 
F1$root

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)
F1$root 


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


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

# Vietovy vzorce
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)$root
(-dc[2]/dc[1]) - F1
(dc[3]/dc[1])/F1

# koeficienty kvadraticke rovnice z korenu
library(pracma)
Poly(c(16, -30))


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

## graficke reseni
xx = seq(-50,50,1)
yy = xx^2
zz = -14*xx + 480
plot(xx,yy,type="l")
abline(480,-14,col=2)

plot(xx,yy,type="l",xlim = c(10,20),ylim=c(0,500))
abline(480,-14,col=2)


```


Tlaková láhev s oxidem uhličitým obsahuje 10.0 kg plynu. Jaký objem zaujímá stlačený plyn, když při teplotě 30 °C je tlak v láhvi 13.17e-6 Pa? Výpočet proved'te pomocí van der Waalsovy rovnice.

```{r echo=FALSE, message=FALSE, warning=FALSE}

# (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 idealni plyn
# p*Vm - R*T = 0  
Vm = R*t/p
V = Vm*n; V  # [m3]

# vypocet pro realny 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]
Vm
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)

```


Jake pH má roztok kyseliny mravenčí o koncentraci 8.5 x 10-4 mol/l?

```{r echo=FALSE, message=FALSE, warning=FALSE}

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 

```



Soustavy lineárních rovnic


Ze dvou slitin s 60% a 80% obsahem mědi se ma získat 40 kg slitiny se 75% obsahem mědi. Kolik kg každé slitiny je třeba použít?  [10 a 30 kg]

```{r}

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

```{r}

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)

```





```{r echo=FALSE, message=FALSE, warning=FALSE}

### skalarni soucin (inner product, dot product, scalar product)

u <- rep(3,3)
v <- 1:3
u%*%v # the inner product
library(pracma)
dot(u, v)
dot(c(1, 2, 3), c(4, 5, 6))


### vektorovy soucin (outer product, cross product, vector product)

u <- rep(3,3)
v <- 1:3
u%o%v # The outer product
library(pracma)
cross(u, v)
cross(c(1, 2, 3), c(4, 5, 6))  
A <- matrix(c(1,0,0, 0,1,0), nrow=2, ncol=3, byrow=TRUE)
crossn(A)  


### Kartezsky soucin

expand.grid(u, v)
expand.grid(c(1, 2, 3), c(4, 5, 6))
expand.grid(x=c("a","b","c"),y=c(1,2,3))
library(data.table)
CJ(u, v)
CJ(c(1, 2, 3), c(4, 5, 6))
CJ(x=c("a","b","c"),y=c(1,2,3))
library(tidyr)
x <- data.frame(x=c("a","b","c"))
y <- data.frame(y=c(1,2,3))
crossing(x, y)

```



```{r echo=FALSE, message=FALSE, warning=FALSE}

M <- matrix(data = rnorm(12), ncol = 3)

# transponovana matice
t(M) 

dim(M)
nrow(M) 
ncol(M)
length(M)
library(Matrix)
nnzero(M, na.counted = NA) # pocet nenulovych prvku matice

round(M,2)
library(Rfast)
Round(M,digit=2,na.rm = FALSE)


### prevod symbolicke matice na numerickou
mch = matrix("1",3,3)
mch
as.numeric(mch) # nefunguje
matrix(as.numeric(mch),ncol=ncol(mch))
library(gmp)
asNumeric(mch)


### Premena matice na vektor
library(gdata)
unmatrix(M, byrow=FALSE)

### Matice na x, y, z
library(squash)
xyzmat2xyz(M)
expand.grid(c(1, 2, 3), c(4, 5, 6))


library(Rfast2) # Split the matrix in lower,upper triangular and diagonal
lud(M)


# jednotkova matice, identita
diag(3)
MM = diag(x = 1, nrow = 3, ncol = 3, names = TRUE)
diag(MM) <- 2
library(Rfast)
Diag.matrix(3,v=1)
Diag.fill(MM,v=0) 
library(Matrix)
Diagonal(3, x = 1)

upper.tri(matrix(1, 3, 3))
round(upper.tri(matrix(1, 3, 3)))
upper.tri(MM)
MM[upper.tri(MM)] <- 0
lower.tri(matrix(1, 3, 3))
round(lower.tri(matrix(1, 3, 3)))
lower.tri(MM)
MM[lower.tri(MM)] <- 0
library(gdata)
upperTriangle(MM, diag=FALSE, byrow=FALSE)
upperTriangle(MM, diag=TRUE, byrow=FALSE)
upperTriangle(MM, diag=FALSE, byrow=FALSE) <- 1
MM
lowerTriangle(MM, diag=FALSE, byrow=FALSE) 
lowerTriangle(MM, diag=FALSE, byrow=FALSE) <- 3
MM
library(Rfast)
lower_tri(MM, suma = FALSE, diag = FALSE) 
upper_tri(MM, suma = FALSE, diag = FALSE) 
lower_tri.assign(MM, 1, diag = FALSE) 
upper_tri.assign(MM, 3, diag = FALSE)
library(Matrix) # v maticovem tvaru
tril(MM) # lower triangular
triu(MM) # upper triangular

# LU decomposition (decompose a matrix into lower- and upper-triangular matrices)
library(cmna)
lumatrix(MM)

```

```{r}

### nasobeni matic

b <- matrix(nrow = 2, ncol = 2, c(1, 2, 3, 4)); b
a <- matrix(nrow = 2, ncol = 2, c(1, 0, 0, -1)); a
a%*%b
b%*%a
library(Rfast)
mat.mult(a, b)


### stopa (trace) - soucet prvku hlavni diagonaly

library(fBasics)
tr(MM)


### Determinant matice

det(M[,c(1:3)])
library(matlib)
Det(M[,c(1:3)], method = "elimination", verbose = FALSE, fractions = FALSE)
# method = c("elimination", "eigenvalues", "cofactors")


### je matice pozitivne definitni?
library(fBasics)
isPositiveDefinite(MM) 
makePositiveDefinite(MM)


### hodnost matice
# Frobeniova veta: soustava linearnich rovnic ma reseni tehdy a jen tehdy, maji-li matice a rozsirena matice soustavy stejnou hodnost.

qr(M)$rank

library(Matrix)
rankMatrix(M, tol = NULL, method = "qr", sval = svd(M, 0, 0)$d, warn.t = TRUE)[1]
# method = c("tolNorm2", "qr.R", "qrLINPACK", "qr", "useGrad", "maybeGrad")

library(fBasics)
rk(M, method = "qr")
# method = c("qr", "chol")

library(matrixcalc) # JEN PRO CTVERCOVOU MATICI
matrix.rank(M[,c(1:3)], method = "qr")
# method = c("qr", "chol")

library(limSolve)
resolution(M)  # $nsolvable = hodnost matice


### Inverzni matice

## ctvercova matice
M2 <- cbind(c(1,0,1),c(0,1,2),c(0,0,1))
solve(M2)
solve(M2)%*%M2  # check

library(fBasics)
inv(M2)
inv(M2)%*%M2  # check

library(cmna)
invmatrix(M2)
invmatrix(M2)%*%M2  # check

library(matlib)
Inverse(M2, verbose=FALSE, fractions=FALSE)
Inverse(M2)%*%M2 # check

## zobecnena inverze
library(MASS)  #  Moore-Penrose generalized inverse 
ginv(M2)
round(ginv(M2),0)
ginv(M2)%*%M2   # check
round(ginv(M2)%*%M2,0)

library(matlib) 
Ginv(M2, fractions=FALSE) 
Ginv(M2) %*% M2 # check

library(limSolve)
limSolve::Solve(M2)
Solve(M2) %*% M2 # check


### Eigenvalues, eigenvectors
eigen(M2)

library(matlib)
Eigen(M2,tol = sqrt(.Machine$double.eps), max.iter = 100, retain.zeroes = TRUE)
C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric C 
EC <- Eigen(C) # eigenanalysis of C 
EC$vectors %*% diag(EC$values) %*% t(EC$vectors) # check

```


Sestavte matici konstitučních koeficientů

```{r echo=FALSE, message=FALSE, warning=FALSE}

form =c("CH4","H2O","CO2","CO")
form =c("CO","COS","CH4","S2","H2O")
form = c("KMnO4", "MnSO4", "H2O", "MnO2", "K2SO4", "H2SO4")
form = c("Fe2O3", "Al", "Al2O3", "Fe")

library(CHNOSZ)
chars = lapply(form,function(x){names(count.elements(x))})
nums = lapply(form,function(x){as.vector(count.elements(x))})
ulch = sort(unique(unlist(chars)))
Mrl = as.data.frame(matrix(0,length(chars),length(ulch)))
colnames(Mrl)=ulch
for(X in c(1:length(chars))){  
  for(i in c(1:length(chars[[X]]))){ 
    cn1 = which(ulch == chars[[X]][i])
    cn2 = which(chars[[X]] == chars[[X]][i])
    Mrl[X,cn1] = as.numeric(nums[[X]][cn2])
  }
}
rownames(Mrl)=form
Mrl

```


Určete počet lineárně nezávislých reakcí v soustavě obsahující dané složky.
  Gibbsovo stechiometricke pravidlo: maximální počet lineárně nezávislých reakcí r je menší nebo roven rozdílu počtu složek v systému (n) a hodnosti matice konstitučních koeficientů (h); v n-složkovém systému existuje n - h stechiometrických vztahů.

```{r echo=FALSE, message=FALSE, warning=FALSE}

form = c("CH4","H2O","CO","CO2","H2")
form = c("NH3","N2","NO","NO2","H2O","O2")

library(CHNOSZ)
chars = lapply(form,function(x){names(count.elements(x))})
nums = lapply(form,function(x){as.vector(count.elements(x))})
ulch = sort(unique(unlist(chars)))
Mrl = as.data.frame(matrix(0,length(chars),length(ulch)))
colnames(Mrl)=ulch
for(X in c(1:length(chars))){  
  for(i in c(1:length(chars[[X]]))){ 
    cn1 = which(ulch == chars[[X]][i])
    cn2 = which(chars[[X]] == chars[[X]][i])
    Mrl[X,cn1] = as.numeric(nums[[X]][cn2])
  }
}
rownames(Mrl)=form
Mrl
M = Mrl

# hodnost matice
qr(M)$rank
library(Matrix)
rankMatrix(M)[1]

# pocet linearne nezavisl}ch reakcm
(r = length(form) - qr(M)$rank)

```


Vyčíslete rovnici: Fe2O3 + Al = Al2O3 + Fe

```{r echo=FALSE, message=FALSE, warning=FALSE}

Fe = c(2, 0, 0, -1)
O = c(3, 0, -3, 0)
Al = c(0, 1, -2, 0)
A = as.matrix(rbind(Fe, O, Al))
colnames(A) = c("Fe2O3", "Al", "Al2O3", "Fe")
B = c(rep(0,length(A[,1])))
C = cbind(A,B)
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)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = TRUE)
as.character(GE1)
# cim nasobit koeficienty
nnv1 = 1/min(abs(as.numeric(GE1[,length(A[1,])])))
nnv = as.numeric(GE1[,length(A[1,])])*nnv1
NN = -1*(as.matrix(GE1)*nnv)[,c(1:(length(A[1,])-1))]
rownames(NN) = rownames(A)
colnames(NN) = colnames(A)[-length(colnames(A))]
apply(NN,2,sum)

library(MASS)
nn = length(A[1,])
a = A[,-nn]
b = -A[,nn]
ge1 <- MASS::ginv(a) %*% b
ge1*2

library(limSolve)
Ls = limSolve::Solve(A, B = diag(nrow = nrow(A)), tol = sqrt(.Machine$double.eps))
library(MASS)
fractions(Ls)
fractions(Ls)*30

```


Vyčíslete rovnici: KMnO4 + MnSO4 + H2O = MnO2 +  K2SO4 + H2SO4

```{r}

K = c(1, 0, 0, 0, -2, 0)
Mn = c(0, 1, 0, -1, 0, 0)
O = c(4, 4, 1, -2, -4, -4)
S = c(0, 1, 0, 0, -1, -1)
H = c(0, 0, 2, 0, 0, -2)

A = as.matrix(rbind(K, Mn, O, S, H))
colnames(A) = c("KMnO4", "MnSO4", "H2O", "MnO2", "K2SO4", "H2SO4")
B = c(rep(0,length(A[,1])))
C = cbind(A,B)
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)
GE2 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = TRUE)
as.character(GE2)
# cim nasobit koeficienty
nnv1 = 1/min(abs(as.numeric(GE2[,length(A[1,])])))
nnv = as.numeric(GE2[,length(A[1,])])*nnv1
NN = -1*(as.matrix(GE2)*nnv)[,c(1:(length(A[1,])-1))]
rownames(NN) = rownames(A)
colnames(NN) = colnames(A)[-length(colnames(A))]
apply(NN,2,sum)

```


Kolik 96% kyseliny sírové a kolik vody je potřeba na přípravu 1 l 78% kyseliny sírové?

```{r echo=FALSE, message=FALSE, warning=FALSE}

r1 = c(0, 1)
r2 = c(0.96, 0.04)
A = cbind(r1, r2); rownames(A) = c("H2SO4","H2O")
B = c(0.78, 0.22); names(B) = c("H2SO4","H2O")
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks, viz Frobeniova veta
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = FALSE)
NN = GE1[,3] # [l]
names(NN) = rownames(A)
NN
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x
library(cmna)
cgmmatrix(A, B, tol = 1e-06, maxiter = 100) # iterativematrix
solvematrix(A, B) # refmatrix

```


Slitina A obsahuje 1.5 % Si, 1.4 % Mn, 0.4 % P a 0.3 % S. Slitina B obsahuje 0.5 % Si, 1.6 % Mn, 0.2 % P a 0.2 % S. Slitina C obsahuje 3 % Si, 0.5 % Mn, 0.5 % P a 0.05 % S. Kolik každé slitiny A, B, C je třeba na výrobu 100 kg slitiny obsahující 2 % Si, 1 % Mn, 0.4 % P a 0.15 % S?

```{r echo=FALSE, message=FALSE, warning=FALSE}

sl1 = c(0.015, 0.014, 0.004, 0.003)
sl2 = c(0.005,0.016, 0.002, 0.002)
sl3 = c(0.03, 0.005, 0.005, 0.0005)
A = cbind(sl1, sl2, sl3); rownames(A) = c("Si","Mn","P","S")
B = c(0.02, 0.01, 0.004, 0.0015); names(B) = c("Si","Mn","P","S")
# B = c(2, 1, 0.4, 0.15); names(B) = c("Si","Mn","P","S")
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks, viz Frobeniova veta
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = FALSE)
GE1[,4]*100 # [kg]
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x

```


Kolik g 60% a kolik g 30% roztoku NaCl je třeba smíchat při přípravě 100 g 40% roztoku? [20 g 60% a 80 g 35%]

```{r echo=FALSE, message=FALSE, warning=FALSE}

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
solvematrix(A, B) # refmatrix

```


Ze dvou kovu o hustotách 7.4 g/cm3 a 8.2 g/cm3 je třeba připravit 0.5 kg slitiny o hustotě 7.6 g/cm3. Kolik g každého z kovů je k tomu potřeba? [375 g lehčího a 125 g těžšího]

```{r echo=FALSE, message=FALSE, warning=FALSE}

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]

```


Derivace

```{r echo=FALSE, message=FALSE, warning=FALSE}
 
# deriv() 
derivative = deriv(~ x^3, "x") 
x <- 2
eval(derivative)

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


library(Deriv)
Deriv("sin(2*x)", "x", cache.exp=FALSE) # differentiate only by x 

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)")


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()


# 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)


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


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)


```


Hmotný bod koná harmonicky pohyb, závislost výchylky y na čase t je dána vztahem y = A * sin(om*t + phi), kde A, om a phi jsou konstanty. Určete rychlost hmotného bodu v čase t = 0.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(Deriv)
der = Deriv("A*sin(om*t+phi)", "t")
der
der0 = gsub("t", "0", der, fixed = TRUE)
der0
Simplify(der0)

```


Při chemické reakci A + B = C v níž obě počáteční koncentrace výchozích látek nabývají stejné hodnoty a, je závislost koncentrace x vznikající látky na čase dána vztahem x(t) = (a^2)*k*t / 1+a*k*t. Jak se mění koncentrace s časem, je-li k > 1 ? 

```{r echo=FALSE, message=FALSE, warning=FALSE}

D(expression((a^2)*k*t/(1+a*k*t)), "t") # base R

library(Deriv)
Deriv("(a^2)*k*t/(1+a*k*t)", "t", cache.exp=FALSE)  

```


Integrace

```{r echo=FALSE, message=FALSE, warning=FALSE}

# 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)

# exp(-x)
integrate(function(x) exp(-x), 0, Inf)  # integration is possible in R

library(sfsmisc)
x <- seq(1,4,0.1)
integrate.xy(x, exp(x))
integrate.xy(x, fx=exp(x), a=2, b=3, use.spline=TRUE, xtol=2e-08)
integrate.xy(x, fx=exp(x), a=1, b=4, use.spline=TRUE, xtol=2e-08)


### 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))

```


Měrné teplo c benzínu (při konstantním tlaku) v závislosti na teplotě T je vyjádřené funkcí c = -0.06 + 0.0010228T. Vypočítejte střední měrné teplo benzínu pro T = <389.15, 491.15>.

```{r echo=FALSE, message=FALSE, warning=FALSE}

qm = integrate(function(Tk){-0.06 + 0.0010228*Tk}, 389.15, 491.15)
unlist(qm)

0.1*unlist(qm)$value/(491.15-389.15) # [J / kg K]

```


Množství tepla Q (v J) potřebné na zahřátí 10 kg železa o 1 K je pro T = <273.15, 473.15> vyjádřena funkcí Q = 4.1686 (105.3T + 0.42T^2). Vypočítejte střední množství tepla (v J / kg K) potřebného na zahřátí železa z teploty 293.15 K na 373.15 K.  

```{r}

qm = integrate(function(Tk) 4.1868*(105.3*Tk + 0.142*Tk^2), 293.15, 373.15)
unlist(qm)

0.1*unlist(qm)$value/(373.15-293.15) # [J / kg K]

```


SYMBOLICKÁ MATEMATIKA 

```{r}

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)")
eval(yac_expr("Factor(x^2+x-6)"), list(x = 4))

eq <- "x^2 + 4 + 2*x + 2*x"
yac_str(eq) # Evaluate yacas command x (a string)
# zjednodusení výrazu
yac_str(paste0("Simplify(", eq, ")"))

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, ")"))  # 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))


yac_str("poly := (x-3)*(x+2)")
yac_str("Expand(poly)")

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)


# Linear algebra 

A <- outer(0:3, 1:4, "-") + diag(2:5)  # numericka matice
a <- 1:4
B <- ysym(A)                           # yacas matice
B
b <- ysym(a)
b

as_r(B) # objekt yacas na objekt R 
as_y(A)


dim(A)
dim(B)

cbind(b,b)
rbind(b,b)

A[, 2:3]
B[, 2:3]

diag(A)
diag(B)
lower.tri(A, diag = FALSE)
lower.tri(B, diag = FALSE)
as_y(lower.tri(B, diag = FALSE))
upper.tri(A, diag = FALSE)
upper.tri(B, diag = FALSE)
as_y(upper.tri(B, diag = FALSE))
A[upper.tri(A)] <- 1
B[upper.tri(B)] <- 1

C <- B
C[lower.tri(C)] <- "x"
C[upper.tri(C)] <- "1"
diag(C) <- 2
C[2, 3] <- "ou"
C

exp(A)
exp(B)

2*A - A
2*B - B

A %*% a
B %*% b
B * B 
B %*% B

# transpose of a matrix
t(A)
y_fn(B, "Transpose")
t(B)

# inverse of a matrix
library(cmna)
invmatrix(A)

y_fn(B, "Inverse")
solve(B)
A %*% solve(A)
round(A %*% solve(A),0)
B %*% solve(B)
solve(A %*% t(A))
solve(B %*% t(B))
solve(A, a)
solve(B, b)

# trace of a matrix, soucet hodnot na hlavni dagonale
y_fn(B, "Trace")


## Limita funkce
# lim(f, var, val, from_left, from_right) 
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") 


## Derivace
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")


## Integrace

integrate(L, "x")
integrate(L, "x", lower=1, upper=4)
yac_str(paste0("Simplify(", integrate(L, "x", lower=1, upper=4), ")"))

# 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)")

```


V tepelné elektrárne mají zásobu uhlí, která vystačí na 24 dní, bude-li v činnosti 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 vystačí zásoba uhlí, budou-li v provozu všechny 3 bloky najednou?

```{r echo=FALSE, message=FALSE, warning=FALSE}

# x/24 + x/30 + x/20 = 1

library(Ryacas)
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

```


Kolik g 60% a kolik g 30% roztoku NaCl je třeba smíchat při přípravě 100 g 40% roztoku? [20 g 60% a a 80 g 35%]

```{r}

r1 = c(0.60, 1)
r2 = c(0.30, 1)
A = rbind(r1, r2); rownames(A) = c("NaCl","H2O")
A = ysym(A)
B = c(0.40,100); names(B) = c("NaCl","H2O")
B = ysym(B)
solve(A,B)

r1 = c(0.3, 1)
r2 = c(0.6, 1)
A = ysym(cbind(r1, r2))
B = ysym(c(0.4, 1))
solve(A,B)
as_r(solve(A,B))

```



Míč byl vyhozen svisle vzhuru rychlostí 25 m/s. Za jak dlouho bude míč 20 m nad zemí? Odpor vzduchu zanedbejte.

```{r echo=FALSE, message=FALSE, warning=FALSE}

# h = v t - 1/2 g t^2
g = 10 # [m/s^2]
h = 20 # [m]
v = 25 # [m/s]

library(Ryacas)
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]
eval(as_r((2*v+sqrt(-((-4)*v^2+8*g*h)))/(2*g)), list(g = 10, v = 25, h = 20))
ko[2]
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]

# Reseni kvadraticke rovnice

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

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

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 

```


Máme směs kyslíku a oxidu dusnatého. Při jaké koncentraci kysliku (v %) ve směsi se oxid dusnatý oxiduje nejrychleji? 

```{r echo=FALSE, message=FALSE, warning=FALSE}

# 2 NO + O2 = 2 NO2
# v = k*(x^2)*y = k*(x^2)*(100-x)

library(Ryacas)
L <- ysym("(x^2)*(100-x)"); L
der <- deriv(L, "x"); der
as_r(der)
solve(der, "x")
round(200/3,1) # [%]

```


Určete velikost hydrostatické tlakové síly, která pusobí na plášť válce výšky v a poloměru r, zcela naplněněho kapalinou o hustotě rho. 

```{r}

# dF = p*dS = p*2*pi*r*dh = rho*g*h*2*pi*r*dh

library(Ryacas)
L <- ysym("2*pi*r*rho*g*h*dh"); L
int = integrate(L, "h"); int
as_r(int)

```