knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# Cargar librerías necesarias
library(Deriv)      # Para derivadas simbólicas
library(Ryacas)     # Para cálculo simbólico general
library(pracma)     # Para funciones matemáticas adicionales
library(ggplot2)    # Para visualizaciones
library(gridExtra)  # Para múltiples gráficos

cat("EJEMPLOS DE DERIVADAS BÁSICAS\n")
cat("==============================\n\n")

# Función polinomial
f1 <- function(x) x^3 + 2*x^2 - 5*x + 7
f1_prime <- Deriv(f1, "x")
cat("1. Función polinomial:\n")
cat("   f(x) = x³ + 2x² - 5x + 7\n")
cat("   f'(x) =", capture.output(print(f1_prime))[1], "\n\n")

# Función trigonométrica
f2 <- function(x) sin(x) * cos(x)
f2_prime <- Deriv(f2, "x")
cat("2. Función trigonométrica:\n")
cat("   f(x) = sin(x) * cos(x)\n")
cat("   f'(x) =", capture.output(print(f2_prime))[1], "\n\n")

# Función exponencial y logarítmica
f3 <- function(x) exp(x) * log(x)
f3_prime <- Deriv(f3, "x")
cat("3. Función exponencial-logarítmica:\n")
cat("   f(x) = e^x * ln(x)\n")
cat("   f'(x) =", capture.output(print(f3_prime))[1], "\n\n")

# Función compuesta
f4 <- function(x) sin(exp(x^2))
f4_prime <- Deriv(f4, "x")
cat("4. Función compuesta:\n")
cat("   f(x) = sin(e^(x²))\n")
cat("   f'(x) =", capture.output(print(f4_prime))[1], "\n")

cat("\nDERIVADAS PARCIALES\n")
cat("===================\n\n")

# Función de dos variables
f_xy <- function(x, y) x^2 * y + sin(x*y) + exp(x+y)

# Derivada parcial respecto a x
df_dx <- Deriv(f_xy, "x")
cat("Función: f(x,y) = x²y + sin(xy) + e^(x+y)\n\n")
cat("∂f/∂x =", capture.output(print(df_dx))[1], "\n")

# Derivada parcial respecto a y
df_dy <- Deriv(f_xy, "y")
cat("∂f/∂y =", capture.output(print(df_dy))[1], "\n\n")

# Derivada parcial de segundo orden
d2f_dx2 <- Deriv(df_dx, "x")
cat("∂²f/∂x² =", capture.output(print(d2f_dx2))[1], "\n")

# Derivada mixta
d2f_dxdy <- Deriv(df_dx, "y")
cat("∂²f/∂x∂y =", capture.output(print(d2f_dxdy))[1], "\n")

cat("\nOPTIMIZACIÓN CON DERIVADAS SIMBÓLICAS\n")
cat("======================================\n\n")

# Función a optimizar
f_opt <- function(x) x^3 - 6*x^2 + 9*x + 2

# Calcular la derivada
f_opt_prime <- Deriv(f_opt, "x")
cat("Función: f(x) = x³ - 6x² + 9x + 2\n")
cat("Derivada: f'(x) =", capture.output(print(f_opt_prime))[1], "\n\n")

# Resolver f'(x) = 0 para encontrar puntos críticos
# La derivada es 3x² - 12x + 9 = 0
# Simplificando: x² - 4x + 3 = 0
# Soluciones: x = 1 y x = 3

# Calcular la segunda derivada para clasificar
f_opt_double_prime <- Deriv(f_opt_prime, "x")
cat("Segunda derivada: f''(x) =", capture.output(print(f_opt_double_prime))[1], "\n\n")

# Evaluar en los puntos críticos
x_criticos <- c(1, 3)
for (x_c in x_criticos) {
  f_val <- f_opt(x_c)
  f_double_val <- f_opt_double_prime(x_c)
  tipo <- ifelse(f_double_val > 0, "Mínimo local", 
                 ifelse(f_double_val < 0, "Máximo local", "Punto de inflexión"))
  cat(sprintf("En x = %d: f(%d) = %.1f, f''(%d) = %.1f → %s\n", 
              x_c, x_c, f_val, x_c, f_double_val, tipo))
}

# Visualización
x_vals <- seq(-1, 5, length.out = 200)
y_vals <- sapply(x_vals, f_opt)
df_plot <- data.frame(x = x_vals, y = y_vals)

ggplot(df_plot, aes(x = x, y = y)) +
  geom_line(size = 1.2, color = "blue") +
  geom_point(data = data.frame(x = x_criticos, y = sapply(x_criticos, f_opt)),
             aes(x = x, y = y), color = "red", size = 4) +
  labs(title = "Optimización usando derivadas simbólicas",
       x = "x", y = "f(x)") +
  theme_minimal() +
  annotate("text", x = 1, y = f_opt(1) + 2, label = "Máximo local", color = "red") +
  annotate("text", x = 3, y = f_opt(3) - 2, label = "Mínimo local", color = "red")

cat("\nINTEGRACIÓN SIMBÓLICA\n")
cat("=====================\n\n")

# Inicializar Yacas
yacas("Clear(x)")

# Ejemplo 1: Integral de un polinomio
cat("1. Integral de polinomio:\n")
integral1 <- yacas("Integrate(x) x^3 + 2*x - 1")
cat("   ∫(x³ + 2x - 1)dx =", as.character(integral1), "\n\n")

# Ejemplo 2: Integral trigonométrica
cat("2. Integral trigonométrica:\n")
integral2 <- yacas("Integrate(x) Sin(x)*Cos(x)")
cat("   ∫sin(x)cos(x)dx =", as.character(integral2), "\n\n")

# Ejemplo 3: Integral definida
cat("3. Integral definida:\n")
integral3 <- yacas("Integrate(x, 0, 1) x^2")
cat("   ∫₀¹ x² dx =", as.character(integral3), "\n\n")

# Verificación numérica de la integral definida
verificacion <- integrate(function(x) x^2, 0, 1)
cat("   Verificación numérica:", verificacion$value, "\n")

cat("\nSERIES DE TAYLOR\n")
cat("=================\n\n")

# Función para calcular coeficientes de Taylor
taylor_coefs <- function(f, x0, n) {
  coefs <- numeric(n + 1)
  f_deriv <- f
  
  for (i in 0:n) {
    if (i == 0) {
      coefs[i + 1] <- f(x0)
    } else {
      f_deriv <- Deriv(f_deriv, "x")
      coefs[i + 1] <- f_deriv(x0) / factorial(i)
    }
  }
  return(coefs)
}

# Ejemplo: Serie de Taylor de e^x alrededor de x=0
f_exp <- function(x) exp(x)
n_terms <- 5
x0 <- 0
coefs_exp <- taylor_coefs(f_exp, x0, n_terms)

cat("Serie de Taylor de e^x alrededor de x=0:\n")
cat("e^x ≈")
for (i in 0:n_terms) {
  if (i == 0) {
    cat(sprintf(" %.0f", coefs_exp[i + 1]))
  } else {
    cat(sprintf(" + %.4f*x^%d", coefs_exp[i + 1], i))
  }
}
cat("\n\n")

# Visualización de la aproximación
x_vals <- seq(-2, 2, length.out = 200)
y_exact <- exp(x_vals)

# Función de Taylor
taylor_approx <- function(x, coefs) {
  result <- 0
  for (i in 0:(length(coefs) - 1)) {
    result <- result + coefs[i + 1] * x^i
  }
  return(result)
}

y_taylor <- sapply(x_vals, function(x) taylor_approx(x, coefs_exp))

df_taylor <- data.frame(
  x = rep(x_vals, 2),
  y = c(y_exact, y_taylor),
  tipo = rep(c("Función exacta", "Aproximación Taylor"), each = length(x_vals))
)

ggplot(df_taylor, aes(x = x, y = y, color = tipo)) +
  geom_line(size = 1.2) +
  labs(title = "Serie de Taylor de e^x (5 términos)",
       x = "x", y = "y", color = "Tipo") +
  theme_minimal() +
  coord_cartesian(ylim = c(-2, 8))

cat("\nÁLGEBRA SIMBÓLICA\n")
cat("==================\n\n")

# Ejemplo 1: Simplificación de expresiones
cat("1. Simplificación de expresiones:\n")
expr1 <- yacas("Simplify((x+1)^2 - (x-1)^2)")
cat("   Simplificar: (x+1)² - (x-1)² =", as.character(expr1), "\n\n")

# Ejemplo 2: Factorización
cat("2. Factorización:\n")
expr2 <- yacas("Factor(x^2 - 5*x + 6)")
cat("   Factorizar: x² - 5x + 6 =", as.character(expr2), "\n\n")

# Ejemplo 3: Expansión de productos
cat("3. Expansión de productos:\n")
expr3 <- yacas("Expand((x+2)*(x-3)*(x+1))")
cat("   Expandir: (x+2)(x-3)(x+1) =", as.character(expr3), "\n\n")

# Ejemplo 4: Resolución de ecuaciones
cat("4. Resolución de ecuaciones:\n")
solucion <- yacas("Solve(x^2 - 4*x + 3 == 0, x)")
cat("   Resolver: x² - 4x + 3 = 0\n")
cat("   Soluciones: x =", as.character(solucion), "\n")

cat("\nAPLICACIÓN: ANÁLISIS DE MOVIMIENTO PARABÓLICO\n")
cat("==============================================\n\n")

# Definir la posición en función del tiempo
# x(t) = v0*cos(θ)*t
# y(t) = v0*sin(θ)*t - (1/2)*g*t²

# Para simplificar, usamos valores específicos
v0 <- 20  # velocidad inicial (m/s)
theta <- pi/4  # ángulo de lanzamiento (45°)
g <- 9.81  # gravedad (m/s²)

# Funciones de posición
x_pos <- function(t) v0 * cos(theta) * t
y_pos <- function(t) v0 * sin(theta) * t - 0.5 * g * t^2

# Calcular velocidades (derivadas de posición)
vx <- Deriv(x_pos, "t")
vy <- Deriv(y_pos, "t")

cat("Ecuaciones de movimiento:\n")
cat("Posición x(t) = v₀·cos(θ)·t\n")
cat("Posición y(t) = v₀·sin(θ)·t - ½gt²\n\n")

cat("Velocidades (derivadas):\n")
cat("vₓ(t) = dx/dt =", capture.output(print(vx))[1], "= constante\n")
cat("vᵧ(t) = dy/dt =", capture.output(print(vy))[1], "\n\n")

# Calcular aceleraciones (segundas derivadas)
ax <- Deriv(vx, "t")
ay <- Deriv(vy, "t")

cat("Aceleraciones (segundas derivadas):\n")
cat("aₓ(t) = d²x/dt² =", capture.output(print(ax))[1], "\n")
cat("aᵧ(t) = d²y/dt² =", capture.output(print(ay))[1], "\n\n")

# Encontrar el tiempo de vuelo (cuando y = 0)
# Resolviendo: v0*sin(θ)*t - 0.5*g*t² = 0
# t(v0*sin(θ) - 0.5*g*t) = 0
# t = 0 o t = 2*v0*sin(θ)/g

t_vuelo <- 2 * v0 * sin(theta) / g
alcance <- x_pos(t_vuelo)

cat(sprintf("Tiempo de vuelo: %.2f segundos\n", t_vuelo))
cat(sprintf("Alcance máximo: %.2f metros\n", alcance))

# Visualización
t_vals <- seq(0, t_vuelo, length.out = 100)
x_vals <- sapply(t_vals, x_pos)
y_vals <- sapply(t_vals, y_pos)

df_trayectoria <- data.frame(x = x_vals, y = y_vals)

ggplot(df_trayectoria, aes(x = x, y = y)) +
  geom_path(size = 1.2, color = "red") +
  geom_point(data = data.frame(x = c(0, alcance), y = c(0, 0)),
             aes(x = x, y = y), color = "blue", size = 3) +
  labs(title = "Trayectoria del proyectil (análisis simbólico)",
       x = "Distancia horizontal (m)",
       y = "Altura (m)") +
  theme_minimal() +
  coord_fixed() +
  annotate("text", x = alcance/2, y = max(y_vals) + 1, 
           label = "Trayectoria parabólica", color = "red")

cat("\nCÁLCULO DE LÍMITES\n")
cat("==================\n\n")

# Ejemplo 1: Límite simple
cat("1. Límite cuando x → 2:\n")
lim1 <- yacas("Limit(x, 2) (x^2 - 4)/(x - 2)")
cat("   lim[x→2] (x² - 4)/(x - 2) =", as.character(lim1), "\n\n")

# Ejemplo 2: Límite al infinito
cat("2. Límite cuando x → ∞:\n")
lim2 <- yacas("Limit(x, Infinity) (3*x^2 + 2*x)/(x^2 - 1)")
cat("   lim[x→∞] (3x² + 2x)/(x² - 1) =", as.character(lim2), "\n\n")

# Ejemplo 3: Límite trigonométrico fundamental
cat("3. Límite trigonométrico:\n")
lim3 <- yacas("Limit(x, 0) Sin(x)/x")
cat("   lim[x→0] sin(x)/x =", as.character(lim3), "\n\n")

# Visualización del comportamiento cerca del límite
f_limit <- function(x) {
  ifelse(x != 0, sin(x)/x, NA)
}

x_vals <- seq(-1, 1, length.out = 200)
x_vals <- x_vals[x_vals != 0]  # Excluir x = 0
y_vals <- f_limit(x_vals)

df_limit <- data.frame(x = x_vals, y = y_vals)

ggplot(df_limit, aes(x = x, y = y)) +
  geom_line(size = 1.2, color = "purple") +
  geom_hline(yintercept = 1, color = "red", linetype = "dashed") +
  labs(title = "Comportamiento de sin(x)/x cerca de x = 0",
       x = "x", y = "sin(x)/x") +
  theme_minimal() +
  annotate("text", x = 0.5, y = 1.02, label = "Límite = 1", color = "red")

cat("\nANÁLISIS COMPLETO DE UNA FUNCIÓN\n")
cat("=================================\n\n")

# Función a analizar
f_analisis <- function(x) x^3 - 3*x^2 - 9*x + 27

cat("Función: f(x) = x³ - 3x² - 9x + 27\n\n")

# 1. Dominio
cat("1. DOMINIO: ℝ (todos los reales)\n\n")

# 2. Intersecciones con los ejes
cat("2. INTERSECCIONES:\n")
# Con eje Y (x = 0)
y_int <- f_analisis(0)
cat(sprintf("   Con eje Y: (0, %.0f)\n", y_int))

# Con eje X (resolver f(x) = 0)
raices <- yacas("Solve(x^3 - 3*x^2 - 9*x + 27 == 0, x)")
cat("   Con eje X:", as.character(raices), "\n\n")

# 3. Primera derivada (puntos críticos)
f_prime_analisis <- Deriv(f_analisis, "x")
cat("3. PRIMERA DERIVADA:\n")
cat("   f'(x) =", capture.output(print(f_prime_analisis))[1], "\n")

# Resolver f'(x) = 0
# 3x² - 6x - 9 = 0  →  x² - 2x - 3 = 0  →  (x-3)(x+1) = 0
puntos_criticos <- c(-1, 3)
cat("   Puntos críticos: x =", paste(puntos_criticos, collapse = ", "), "\n\n")

# 4. Segunda derivada (concavidad)
f_double_prime_analisis <- Deriv(f_prime_analisis, "x")
cat("4. SEGUNDA DERIVADA:\n")
cat("   f''(x) =", capture.output(print(f_double_prime_analisis))[1], "\n")

# Punto de inflexión (f''(x) = 0)
# 6x - 6 = 0  →  x = 1
punto_inflexion <- 1
cat("   Punto de inflexión: x =", punto_inflexion, "\n\n")

# 5. Análisis de monotonía y concavidad
cat("5. ANÁLISIS DE INTERVALOS:\n")
test_points <- c(-2, 0, 2, 4)
for (i in 1:length(test_points)) {
  x_test <- test_points[i]
  f_prime_val <- f_prime_analisis(x_test)
  f_double_val <- f_double_prime_analisis(x_test)
  
  crecimiento <- ifelse(f_prime_val > 0, "Creciente", "Decreciente")
  concavidad <- ifelse(f_double_val > 0, "Cóncava hacia arriba", "Cóncava hacia abajo")
  
  if (i == 1) cat("   x < -1: ")
  else if (i == 2) cat("   -1 < x < 1: ")
  else if (i == 3) cat("   1 < x < 3: ")
  else cat("   x > 3: ")
  
  cat(crecimiento, ",", concavidad, "\n")
}

# 6. Visualización completa
x_vals <- seq(-4, 6, length.out = 300)
y_vals <- sapply(x_vals, f_analisis)

# Puntos importantes
puntos_importantes <- data.frame(
  x = c(puntos_criticos, punto_inflexion, -3, 3),
  y = sapply(c(puntos_criticos, punto_inflexion, -3, 3), f_analisis),
  tipo = c("Máximo local", "Mínimo local", "Punto de inflexión", 
           "Raíz", "Raíz")
)

ggplot(data.frame(x = x_vals, y = y_vals), aes(x = x, y = y)) +
  geom_line(size = 1.2, color = "darkblue") +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed") +
  geom_vline(xintercept = 0, color = "gray", linetype = "dashed") +
  geom_point(data = puntos_importantes, 
             aes(x = x, y = y, color = tipo), size = 4) +
  labs(title = "Análisis completo de f(x) = x³ - 3x² - 9x + 27",
       x = "x", y = "f(x)", color = "Puntos importantes") +
  theme_minimal() +
  scale_color_manual(values = c("Máximo local" = "red", 
                               "Mínimo local" = "blue",
                               "Punto de inflexión" = "green",
                               "Raíz" = "purple"))