
#-------------------------------------------------------------------------------------------#

x <- c(1,20,30) # un vector fila x (vector)
sqrt(sum(x * x)) # su norma

library("pracma") # paquete que estudiaremos más a fondo
dot(x,x)^(1/2) # su norma

x_ <- matriz(x) # un vector columna
X_


dot(x_,x_)^(1/2) # su norma

#------------------------------------------------- -------------------------------------------#

x <- 0.0000000002 # valor a probar

p <- 0.000000001 # orden de magnitud de la precisión de la prueba
(- p <= x) && (x<= p) # esta prueba nulidad de x da TRUE

p <- 0.00000000000001 # orden de magnitud de la precisión de la prueba
(- p <= x) && (x<= p) # esta prueba de nulidad de x da FALSE

cat("precisión:", .Machine$doble.eps)
cat("el codificado real más pequeño en 'float':", .Machine$double.xmin)
cat("el códificado real más grande en 'float':", .Machine$double.xmax)
cat("entero demasiado grande:",as.integer(2^31))
cat("entero máximo:",as.integer(2^31-1))
cat("entero mínimo:",as.integer(-(2^31-1)))

e <- 2^39
e

class(e)

e <- 2^40
e
class(e)


#-------------------------------------------------------------------------------------------#

z <- c(1,
       x <- c(1,20,30) # un vector fila x (vector)
       sqrt(sum(x * x)) # su norma
       
       library("pracma") # paquete que estudiaremos más a fondo
       dot(x,x)^(1/2) # su norma
       
       x_ <- matriz(x) # un vector columna
       X_
       
       
       dot(x_,x_)^(1/2) # su norma
       
       #------------------------------------------------- -------------------------------------------#
       
       x <- 0.0000000002 # valor a probar
       1i, # solo "yo"
       -1,
       -1i, # -i
       0i, # si quiere forzar un cero en los complejos
       0)
Real(z) # parte real(z)
Imag(z) # parte imaginaria(z)
Conj(z) # conjugado(z)
abs(z) # módulo(z)

angle(z)/(2*pi) # el argumento
0i ^ (-3:3) # creando infinitos y ceros complejos

#------------------------------------------------- -------------------------------------------#

eye(3) # 3*3 matriz identidad cuadrada
eye(3,2) # bloque de identidad en matriz 3*2
eye(2,3) # bloque de identidad en matriz 3*2
eye(1) # una matriz que contiene un 1


A1 <- ones(2,2)
A2 <- 2*ones(3,2) - eye(3,2)
A3 <- 3*ones(2,2)
B <- diag.blk(A1,A2,A3)
B

C <- matrix(1:12, nrow = 4)
C

fliplr(C) # orden inverso de las columnas
flipud(C) # orden de fila inverso
circshift(C, c(2,0)) # rotar 2 filas
circshift(C, c(0,2)) # rotar 2 columnas
circshift(C, c(4,3)) # no hace nada porque nuestra matriz es un 4 * 3!

#------------------------------------------------- -------------------------------------------#

rand(2,3) # distribución uniforme en ]0,1[
randn(2,3) # distribución normal centrada reducida: N(mu=0,s=1)
randi(50,2,3) # enteros entre 0 y 50 en una matriz de 2*3

set.seed(666)
randsample (10, 7, w = NULL, replacement = FALSE) # 7 de 10, todos diferentes

set.seed(666)
randsample (10, 7, w = NULL, replacement = TRUE) # 7 de 10, con descuento

#------------------------------------------------- -------------------------------------------#

set.seed(0) # iniciar generador de números aleatorios
M <- randi(50,2,3) # una matriz de enteros positivos
M

mod(M,4) # módulo 4 (resto de la división por 4)

idivide(M,4,"floor") # división de enteros (número máximo de veces 4)
idivide(M,4,"round") # división en reales, luego redondea a entero

#-------------------------------------------------------------------------------------------#

# familia de vectores, no ortonormales
v1 <- matrix(c(1,
               1,
               0), ncol = 1)

v2 <- matrix(c(0,
               2,
               0), ncol = 1)

v3 <- matrix(c(0,
               0,
               3), ncol = 1)

M <- cbind(v1,v2,v3) # matriz de 3 vectores columna
Rank(M) # su rango es igual a 3, ¡buenas noticias!

M_ <- orth(M) # crear matriz de vector base
M_

Rank(M_) # la matriz obviamente es de rango 3

# extracción de los 3 vectores de la nueva base
v1_ <- matrix(M_[ ,1])
v2_ <- matrix(M_[ ,2])
v3_ <- matrix(M_[ ,3])

# comprobar que la base es ortonormal a través 
# del cálculo de los productos escalares de dos en dos
l <- list(v1_,v2_,v3_)
d <- NULL
for (v in l) {for (w in l) {d <-c(d, dot(v,w))}}
matrix(round(d,10), nrow =3) # matrix de los productos escalares de vi_

#-------------------------------------------------------------------------------------------#

# base del subespacio lineal
# es una linea así que solo un vector
vh1 <- matrix(c(1,
                1,
                0), ncol = 1)

# matriz de vectores-columna correspondientes
# aquí solo un vector
M_hYperplan <- cbind(vh1)

# Matriz de columnas-vectores a proyectar
# varios casos didácticos
V <- cbind(vh1,v2,v3,v1_,v2_,v3_,10*v3_,v1_+ 100*v2_)
V

# proyección
l <- linearproj(M_hyperplan, V)
Vp <- l$Q # extracción - matriz de vectores proyectados
Vp

#------------------------------------------------- -------------------------------------------#

# definición de las dos funciones
f <- function(x){exp(-x)*sen(x)}
g <- function(x){exp(x)*sen(x)+1000}

curve(f,-10,10, col = "red", ylab = "f y g")
curve(g,-10,10, col = "blue", sum = TRUE)
grid()

x_ <- linspace(-10,10,100 + 1) # discretización de x a través de pracma::linspace
# -10 -9.8 ... 0 ... 9.8 10

plotyy(x_, f(x_), # f y g no tienen la misma escala
       x_, g(x_))

library (ggplot2)
f <- function(x){exp(-x)*sin(x)} # define su función
x_ <- linspace(-10,10,100 + 1) # discretización de x a través de pracma::linspace
# -10 -9.8 ... 0 ... 9.8 10
p1 <- ggplot(data.frame(x = x_), aes(x = x)) # observe el uso de data.frame
p1 <- p1 + stat_function(fun = f) # así que este será un gráfico de función


g <- function(x){exp(x)*sin(x)+1000} # define su función
p2 <- ggplot(data.frame(x = x_), aes(x = x)) # observe el uso de data.frame
p2 <- p2 + stat_function(fun = g) # así que este será un gráfico de función

library (patchwork)

p1 + p2

p1 + stat_function(fun = g) # agregando g a la gráfica de f


x <- logspace(-2, 2) # de 1/100 a 100
loglog(x, 1/tanh(x), # escala logarítmica en ambos ejes
       type='l',
       ylab="1/tangente_hiperbolica")

r <- function(t) 1-sqrt(2*t)
ezpolar(r)

#-------------------------------------------------------------------------------------------#

h <- function(x) {x * cos(pi*exp(-x)) * sin(pi*exp(x))}
ezplot(h, 0, 1.705, 
       n = 100, 
       fill = TRUE,
       ylab = "oscilador libre",
       main = "")

paso <- 0.001 # el incremento entre cada iteración
test_0 <- 0.000001 # precisión relativa a cero pruebas
mini <- 0 # límites mínimos y máximos de búsqueda raíz
máx <- 1.705
x <- mini # inicializaciones del algoritmo de búsqueda
raices <- NULL

while (x <= maxi) {
  xmin    <- x
  x       <- x + pas
  c_una_raiz      <- abs(h(xmin))  <= test_0     # raíz evidente
  signos_diferentes <- sign(h(xmin)) != sign(h(x)) # h crece en el eje de las x
  if (c_una_raiz) {raices  <- c(raices,xmin)}
  if (!c_una_raiz && signos_diferentes) {
    raices  <- c(raices,ridders(h,xmin,x)$root)  # búsqueda de raíz    
  }
}

raices

library(rootSolve)
raices=uniroot.all(h,c(0,1.705))
raices

library(mosaic)
findZeros(h(x) ~ x , xlim = c(0,1.705) ) 

#-------------------------------------------------------------------------------------------#

library(tidyverse)
library(plotly)

f <- function(x,y){sin(x) * sinh(y)}         # una función de dos variables

x_ <- y_ <-  pracma::linspace(-5,5,100 + 1)  # discretizar x y y
z_ <- outer(x_,y_,f)                         # cálculo de las f(x_,y_)
plot_ly() %>% add_surface(x = ~x_, y = ~y_, z = ~z_) # plot interactivo

# una función frenética
f <- function(x){erf(x[1]+x[2])*exp(sin(x[1]))*(x[2]+1)} 
m1 <- fminunc(c( 5,-5), f)                           # búsqueda de un mínimo
m2 <- fminunc(c(-5, 5), f)                           # búsqueda de otro mínimo

m  <- rbind(c(m1$par,m1$value),c(m2$par,m2$value))
m

for (i in 1:nrow(m)) { print(abs(f(c(m[i,1],m[i,2])) - m[i,3]) <= 0.00000001) }

library(tidyverse)
library(plotly)

f_ <- Vectorize(function(x,y) f(c(x,y)))       # hace f compatible con  plotly
x_ <- y_ <-  pracma::linspace(-20,20,300 + 1)  # discretizar x e y
z_ <- outer(x_,y_,f_)                          # cálculos de las f_(x_,y_)
plot_ly() %>% add_surface(x = ~x_, y = ~y_, z = ~z_) # plot inteactivo



#-------------------------------------------------------------------------------------------#

A <- matrix(c(73,85,9,87), ncol = 2)   # una matriz
A

b <- matrix(c(1,2), ncol =1)           # un vector
b

det(A)                                 # el determinante

A_1 <- solve(A)                        # el determinante no era nulo
A_1

round(A_1 %*% A,2)                     # verificación opcional 

x_solution <- A_1 %*% b                # cálculo de la solución
x_solution


#-------------------------------------------------------------------------------------------#

A <- matrix(c(7,-4,2,3,4,2,2,10:17), ncol=3, byrow=TRUE)
A

b <- matrix(c(-1,2,-3,4,-5), ncol=1)
b

x_solution <- qrSolve(A, b)
x_solution

#-------------------------------------------------------------------------------------------#


# f es una función de un vector x
# que devuelve un vector f(x) de la misma dimensión
f <- function(x) { 
  x1 <- x[1]
  x2 <- x[2]
  x3 <- x[3]
  
  y1 <- x1 - 2 * x2
  y2 <- x2 + 1
  y3 <- x3 * x1 + 1
  
  c(y1,y2,y3)
}

x0 <- c(1, 1, 1)

broyden(f, x0)$zero
f(c(-2.0,-1.0 ,0.5))

#-------------------------------------------------------------------------------------------#

e <- expression(x^2+3*cos(x)- 1/x)
d <- D(e,'x')
d

x <- pi
eval(d)

Df <- function(x){eval(bquote(.(d)))}
Df(pi)
Df(x)
Df(2*pi)

Df<- Vectorize(Df)  # inútil aquí pero importante para los casos complicados
mes_data <- pi*1:5
Df(mes_data)

# Mi función de derivación
Df_  <- function (f, dx){
  d <- function(x)(f(x + dx)-f(x))/dx
}


f <- function(x){x^2+3*cos(x)- 1/x}

Df_(f,0.000000001)(x)


Df <- function(x){Df_(f,0.000000001)(x)}
Df<- Vectorize(Df)  # inútil aquí pero importante para los casos complicados
mes_data <- pi*1:5
Df(mes_data)

#-------------------------------------------------------------------------------------------#

library(numDeriv)
grad(f,mes_data)

library(Deriv)

Df <- Deriv(f,"x")
Df

Df(mes_data)

Df6 <- Deriv(f,"x",nderiv = 6)
Df6

#-------------------------------------------------------------------------------------------#

Deriv(gamma,"x")

Deriv(gamma,"x", nderiv = 2)

pracma::zeta(-3)-1/120
pracma::zeta(3)

fderiv(zeta,-7)

#-------------------------------------------------------------------------------------------#

f <- function(x,y,z) x + y + 10*z    # la función a integrar

# definición del dominio D
z1 <- 0
z2 <- 1

y1 <- function(z) z
y1 <- Vectorize(y1) # vectoriza y1, inútil en este caso, porque y1 es muy
# simple (hecho aquí para llamar su atención)
y2 <- 1

x1 <- function(y,z) y
x1 <- Vectorize(x1)     # vectoriza x1

x2 <- function(y,z) z
x2 <- Vectorize(x2)     # vectoriza x2

integral3(f, 
          z1, z2,
          y1, y2, 
          x1, x2)

#-------------------------------------------------------------------------------------------#


library(Bessel)

x <- linspace(-6.5,3,100) 
plotyy(x, AiryA(x),               
       x, AiryB(x), 
       xlab = "pequeños valores reales",
       ylab = "Ai (A la izquierda/azul), Bi (A la derecha/rojo)")

#-------------------------------------------------------------------------------------------#

erf(sqrt(2))     # un valor conocido

# la inversa  de erf funciona en su dominio
# de definición: -1 Ã  1
erf(erfinv(c(-1.1,-1,0,1,1.1)))  

erfz(1i)         # erf de un valor complejo


#-------------------------------------------------------------------------------------------#

f  <- function(t,y) {(t-sin(y))/(1 + sin(y))} # función de t y y
y0 <- 0                                       # valor inicial de y a t0
t0 <- 0                                       # t0
tf <- 10                                      # t final
df <- data.frame(ode45(f,t0,tf,y0))           # resolución y almacenamiento de la solución
ggplot2::qplot(t,y,data = df, geom = "line", col = y)

#-------------------------------------------------------------------------------------------#

f <- Vectorize(function(x) 2*x^2 / (1 + x^3))
g <- Vectorize(function(x)  -2 / (1 + x^2 + 10 * sin(x)))
h <- Vectorize(function(x) {1}) 

a <- -3.5
b <-  3.5
alpha <-  3.0
beta  <- -3.0

# resolución y almacenamiento de la solución
df <- data.frame(
  bvp(f, g, h, c(a,b), c(alpha,beta))
)

ggplot2::qplot(xs,ys,data = df, geom = "line", col = ys)

#-------------------------------------------------------------------------------------------#



f <- function(t, y1, y2) -y2^2 / y1 + t
h <- function(u, v) u - 30
t0 <- 0 
tf <- 1
y0 <- 2

y <-shooting(f, t0, tf, y0, h, 0, 1)
y1 <- y$y[ ,1]
y2 <- y$y[ ,2]
t  <- y$t
df <- data.frame(y1,y2,t)

dplyr::glimpse(df, width = 70)
qplot(t,y1,geom ="line")
qplot(y1,y2,geom ="line")

library(lattice)
cloud(t~y1+y2,df, type = c("h","l"),drape = TRUE)


library(plotly)
plot_ly(df, 
        x = ~y1, y = ~y2, z = ~t, 
        type = 'scatter3d', mode = 'lines',
        line = list(width = 6))


#-------------------------------------------------------------------------------------------#

f <- Vectorize(function(x){exp(-x)*sin(x)+ x^3})
complexstep(f,1)          # derivada robusta
complexstep(f,c(1,2,3))   # derivada de cada componente del vector


f1 <- function(x) abs(x^2-4)
curve(f1, -4, +4); grid()
fderiv(f1,2)       # resultado engañoso
fderiv(f1, 2, method = "backward") # A izquierda
fderiv(f1, 2, method = "forward")  # A derecha

library(mosaicCalc)
# otra forma de crear una derivada
# reutilizable con el paquete mosaicCalc
f_prime <- D(f(x) ~ x)
f_prime(c(1,2,3)) # tenemos una función derivada vectorizada

F <- antiD(f(x) ~ x) # primitiva, f(x).dx
I <- F(2)-F(1) # integral de 1 a 2 (entonces F(2)-F(1))

fderiv(f, c(1,2,3),2) # segunda derivada vectorizada  


#-------------------------------------------------------------------------------------------#

# h(x) fonction d'un vecteur de dim n=2
h <- function(x) (f(x[1]) + 2* f(x[2]))/3

# derivadas parciales de h(x,y)
f_deriv_partiel_x <- D(h(c(x,y))~x)    
args(f_deriv_partiel_x)                # función de dos variables
f_deriv_partiel_x(0,0)     

f_deriv_partiel_y <- D(h(c(x,y))~y)
f_deriv_partiel_y(0,0)   

D(h(c(x,y))~ x + x )(0,0)
D(h(c(x,y))~ x + y )(0,0)
D(h(c(x,y))~ y + x )(0,0)
D(h(c(x,y))~ y + y )(0,0)


#-------------------------------------------------------------------------------------------#

grad_csd(h, c(0,0))
hessian_csd(h, c(0,0))
laplacian_csd(h, c(0,0))

sum(Diag(hessian_csd(h, c(0,0))))


#-------------------------------------------------------------------------------------------#

g <- function(x) {y1 <- x[1] + x[2]^2
y2 <- sin((x[3]))
y3 <- x[2] - x[1]^2 + 1
y4 <- cos(x[2])
c(y1,y2,y3,y4)
}

x0 <- c(0,0,0)
g(x)


jacobian_csd(g,c(-1,0,0))

#-------------------------------------------------------------------------------------------#

library(rSymPy)  

x <- Var("x")                    # No olvidar las ""
x+3*x

x/x

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

sympy("limit(x/exp(x), x, oo)")

#-------------------------------------------------------------------------------------------#











