Análisis Numérico. Tarea 5. Splines cúbicos.

Sandoval Lazcano Kenjhy Axcel

Marzo del 2022

Ejercicio 1.

En cada inciso considera la función de valores reales, usa splines cúbicos para encontrar una aproximación en el intervalo dado. Calcula el valor de la función, de la derivada y en cada caso calcula el error real.

  1. \(f(x)=e^{2x}\). Puntos: \(x_0=0, x_1=0.25, x_2=0.5, x_3=0.75\). Aproximar \(f(0.43)\) y \(f'(0.43)\).
#Función
f_1a <- function(x){exp(2*x)}
x1a <- seq(0, 0.75, 0.001)
y1a <- f_1a(x1a)

#Nodos
x_1an <- c(0, 0.25, 0.5, 0.75)
y_1an <- f_1a(x_1an)

#Spline cúbico
sp_1a <- cubicspline(x_1an, y_1an, x1a)

graf_1a <- ggplot()+
  geom_line(aes(x1a, y1a), color="black", size=1)+
  geom_point(aes(x_1an, y_1an), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x1a, sp_1a), color="red", size=1)+
  theme_bw()

ggplotly(graf_1a)

Valor de la función en \(x=0.43\), \(f(0.43)=2.3631607\)

cubicspline(x_1an, y_1an, 0.43)
## [1] 2.347353

Error absoluto:

#Abs sirve para obtener el valor absoluto
abs(f_1a(0.43)-cubicspline(x_1an, y_1an, 0.43))
## [1] 0.01580817

Derivada de la función \(f´(x)=2 e^(2x)\), evaluada en \(x=0.43\)

df_1a <- function(x){2*exp(2*x)}
2*exp(2*0.43)
## [1] 4.726321
  1. \(f(x)=x\, log(x)\), \(x\in [2,12]\), \(h=2\). Aproximar \(f(8.4)\) y \(f'(8.4)\).
f_1b <- function(x){x*log(x)}
x1b <- seq(2, 12, 0.001)
y1b <- f_1b(x1b)

x_1bn <- c(2, 4, 6, 8, 10, 12)
y_1bn <- f_1b(x_1bn)

sp_1b <- cubicspline(x_1bn, y_1bn, x1b)

graf_1b <- ggplot()+
  geom_line(aes(x1b, y1b), color="black", size=1)+
  geom_point(aes(x_1bn, y_1bn), color="dodgerblue3", size=3)+
 
  geom_line(aes(x1b, sp_1b), color="red", size=1)+
  theme_bw()

ggplotly(graf_1b)

Valor de la función en \(x=8.4\), \(f(8.4)=17.8771463\)

cubicspline(x_1bn, y_1bn, 8.4)
## [1] 17.87406

Error absoluto:

abs(f_1b(8.4)-cubicspline(x_1bn, y_1bn, 8.4))
## [1] 0.003082143

Derivada de la función \(f´(x)=log(x)+1\), evaluada en \(x=8.4\)

df_1b <- function(x){log(x)+1}
log(8.4)+1
## [1] 3.128232
  1. \(f(x)=sen(e^x-2)\), \(x\in [0,5]\), \(h=1\). Aproximar \(f(0.9)\) y \(f'(0.9)\).
f_1c <- function(x){sin(exp(x)-2)}
x1c <- seq(0, 5, 0.001)
y1c <- f_1c(x1c)

x_1cn <- c(0, 1, 2, 3, 4, 5)
y_1cn <- f_1c(x_1cn)

sp_1c <- cubicspline(x_1cn, y_1cn, x1c)

graf_1c <- ggplot()+
  geom_line(aes(x1c, y1c), color="black", size=1)+
  geom_point(aes(x_1cn, y_1cn), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1c, sp_1c), color="red", size=1)+
  theme_bw()

ggplotly(graf_1c)

Valor de la función en \(x=0.9\), \(f(0.9)=0.4435924\)

cubicspline(x_1cn, y_1cn, 0.9)
## [1] 0.656121

Error absoluto:

abs(f_1c(0.9)-cubicspline(x_1cn, y_1cn, 0.9))
## [1] 0.2125286

Derivada de la función \(f´(x)=cos(e^x-2)\), evaluada en \(x=0.9\)

df_1c <- function(x){cos(exp(x)-2)}
cos(exp(0.9)-2)
## [1] 0.8962286
  1. \(f(x)=x\, cos\,x-2x^2+3x-1\). \(x\in [0,2]\), \(h=0.5\). Aproximar \(f(0.25)\) y \(f'(0.25)\).
f_1d <- function(x){x*cos(x)-2*x^2+3*x-1}
x1d <- seq(0, 2, 0.001)
y1d <- f_1d(x1d)

x_1dn <- c(0, 0.5, 1, 1.5, 2)
y_1dn <- f_1d(x_1dn)

sp_1d <- cubicspline(x_1dn, y_1dn, x1d)

graf_1d <- ggplot()+
  geom_line(aes(x1d, y1d), color="black", size=1)+
  geom_point(aes(x_1dn, y_1dn), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1d, sp_1d), color="red", size=1)+
  theme_bw()

ggplotly(graf_1d)

Valor de la función en \(x=0.25\), \(f(0.25)=-0.1327719\)

cubicspline(x_1dn, y_1dn, 0.25)
## [1] -0.1773416

Error absoluto:

abs(f_1d(0.25)-cubicspline(x_1dn, y_1dn, 0.25))
## [1] 0.04456967

Derivada de la función \(f´(x)=-xsen(x)-4x+cos(x)+3\), evaluada en \(x=0.25\)

df_1d <- function(x){-x*sin(x)-4*x+cos(x)+3}
-0.25*sin(0.25)-4*0.25+cos(0.25)+3
## [1] 2.907061
  1. \(f(x)=x\,cos\,x-3x\). Puntos: \(x_0=0.1, x_1=0.2, x_2=0.3, x_3=0.4\). Aproximar \(f(0.18)\) y \(f'(0.18)\).
f_1e <- function(x){x*cos(x)-3*x}
x1e <- seq(0, 0.4, 0.001)
y1e <- f_1e(x1e)

x_1en <- c(0, 0.1, 0.2, 0.3, 0.4)
y_1en <- f_1e(x_1en)

sp_1e <- cubicspline(x_1en, y_1en, x1e)

graf_1e <- ggplot()+
  geom_line(aes(x1e, y1e), color="black", size=1)+
  geom_point(aes(x_1en, y_1en), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1e, sp_1e), color="red", size=1)+
  theme_bw()

ggplotly(graf_1e)

Valor de la función en \(x=0.18\), \(f(0.18)=-0.3629081\)

cubicspline(x_1en, y_1en, 0.18)
## [1] -0.362941

Error absoluto:

abs(f_1e(0.18)-cubicspline(x_1en, y_1en, 0.18))
## [1] 3.28817e-05

Derivada de la función \(f´(x)=cos(x)-xsen(x)-3\), evaluada en \(x=0.18\)

df_1e <- function(x){cos(x)-x*sin(x)-3}
cos(0.18)-0.18*sin(0.18)-3
## [1] -2.048382

Ejercicio 2

Encuentra los splines cúbicos condicionados para las funciones del ejercicio anterior.

1.a

#Función
f_1a <- function(x){exp(2*x)}
x1a <- seq(0, 0.75, 0.001)
y1a <- f_1a(x1a)

#Nodos
x_1an <- c(0, 0.25, 0.5, 0.75)
y_1an <- f_1a(x_1an)

#Spline cúbico condicionado
sp_2a <- cubicspline(x_1an, y_1an, x1a, endp2nd = TRUE, der=c(df_1a(0), df_1a(0.75)))

graf_2a <- ggplot()+
  geom_line(aes(x1a, y1a), color="black", size=1)+
  geom_point(aes(x_1an, y_1an), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x1a, sp_2a), color="red", size=1)+
  theme_bw()

ggplotly(graf_2a)

1.b

f_1b <- function(x){x*log(x)}
x1b <- seq(2, 12, 0.001)
y1b <- f_1b(x1b)

x_1bn <- c(2, 4, 6, 8, 10, 12)
y_1bn <- f_1b(x_1bn)

sp_2b <- cubicspline(x_1bn, y_1bn, x1b, endp2nd = TRUE, der=c(df_1b(2), df_1b(12)))

graf_2b <- ggplot()+
  geom_line(aes(x1b, y1b), color="black", size=1)+
  geom_point(aes(x_1bn, y_1bn), color="dodgerblue3", size=3)+
 
  geom_line(aes(x1b, sp_2b), color="red", size=1)+
  theme_bw()

ggplotly(graf_2b)

1.c

f_1c <- function(x){sin(exp(x)-2)}
x1c <- seq(0, 5, 0.001)
y1c <- f_1c(x1c)

x_1cn <- c(0, 1, 2, 3, 4, 5)
y_1cn <- f_1c(x_1cn)

sp_2c <- cubicspline(x_1cn, y_1cn, x1c, endp2nd = TRUE, der=c(df_1c(0), df_1c(5)))

graf_2c <- ggplot()+
  geom_line(aes(x1c, y1c), color="black", size=1)+
  geom_point(aes(x_1cn, y_1cn), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1c, sp_2c), color="red", size=1)+
  theme_bw()

ggplotly(graf_2c)

1.d

f_1d <- function(x){x*cos(x)-2*x^2+3*x-1}
x1d <- seq(0, 2, 0.001)
y1d <- f_1d(x1d)

x_1dn <- c(0, 0.5, 1, 1.5, 2)
y_1dn <- f_1d(x_1dn)

sp_2d <- cubicspline(x_1dn, y_1dn, x1d, endp2nd = TRUE, der=c(df_1d(0), df_1d(2)))

graf_2d <- ggplot()+
  geom_line(aes(x1d, y1d), color="black", size=1)+
  geom_point(aes(x_1dn, y_1dn), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1d, sp_2d), color="red", size=1)+
  theme_bw()

ggplotly(graf_2d)

1.e

f_1e <- function(x){x*cos(x)-3*x}
x1e <- seq(0, 0.4, 0.001)
y1e <- f_1e(x1e)

x_1en <- c(0, 0.1, 0.2, 0.3, 0.4)
y_1en <- f_1e(x_1en)

sp_2e <- cubicspline(x_1en, y_1en, x1e, endp2nd = TRUE, der=c(df_1e(0), df_1e(0.4))) 

graf_2e <- ggplot()+
  geom_line(aes(x1e, y1e), color="black", size=1)+
  geom_point(aes(x_1en, y_1en), color="dodgerblue3", size=3)+
  
  geom_line(aes(x1e, sp_2e), color="red", size=1)+
  theme_bw()

ggplotly(graf_2e)

Ejercicio 3

Se sospecha que las elevadas concentraciones de tanina en las hojas de los robles maduros inhiben el crecimiento de las larvas de la polilla invernal (Operophtera bromata L. Geometridae) que tanto dañan a los árboles en algunos años. La tabla anexa contiene el peso promedio de dos muestras de larva, tomadas en los primeros 28 días después de nacimiento. La primera muestra se crió en hojas de robles jóvenes, mientras que la segunda lo hizo en hojas maduras del mismo árbol.

  1. Usa splines cúbicos para aproximar la curva del peso promedio de las muestras.

  2. Para calcular un peso promedio máximo aproximado de cada muestra, determina el máximo del polinomio interpolante.

\[\begin{equation} \begin{array}{l|c|c|c|c|c|c|r} \text{Día} & 0 & 6 & 10 & 13 & 17 & 20 & 28 \\ \hline \text{Peso promedio muestra 1 (mg)} & 6.67 & 17.33 & 42.67 & 37.33 & 30.10 & 29.31 & 28.74 \\ \text{Peso promedio muestra 2 (mg)} & 6.67 & 16.11 & 18.89 & 15.00 & 10.56 & 9.44 & 8.89 \end{array} \end{equation}\]

Debido a que no sé construir un cubic spline con estos datos, es decir, sin contar con la función, lo haré por medio de polinomios interpolantes:

dias <- c(0, 6, 10, 13, 17, 20, 28)
dias_seq <- seq(from=0, to=28, by=0.01)

muestra1 <- c(6.67, 17.33, 42.67, 37.33, 30.10, 29.31, 28.74)
pol_muestra1 <- as.function(poly.calc(dias, muestra1))
muestra_seq <- pol_muestra1(dias_seq)

muestra2 <- c(6.67, 16.11, 18.89, 15.00, 10.56, 9.44, 8.89)
pol_muestra2 <- as.function(poly.calc(dias, muestra2))
muestra_seq2 <- pol_muestra2(dias_seq)

graf_2 <- ggplot()+
  geom_vline(xintercept = 0, linetype="dashed")+
  geom_hline(yintercept = 0, linetype="dashed")+
  #Muestra 1
  geom_line(aes(x=dias_seq, y=muestra_seq), color="green", size=1)+
  geom_point(aes(x=dias, y=muestra1), color="dodgerblue3", size=3)+
  #Muestra 2
  geom_line(aes(x=dias_seq, y=muestra_seq2), color="red", size=1)+
  geom_point(aes(x=dias, y=muestra2), color="orange", size=3)+
  
  labs(x="dias", y="Peso Muestra", title="Interpolación  2, inciso a)")+
  theme_bw()

ggplotly(graf_2)
  1. La muestra 1 está representada por la línea color verde, la muestra 2, por la línea color roja. Una vez que obtuvimos el Polinomio Interpolante para cada una de estas muestras, podemos observar a simple vista que el valor máximo de la 1 se encuentra en el día 10, pero calcularemos con mayor precisión el valor máximo de ambos polinomios mediante el siguiente procedimiento:
  1. Calcular el peso máximo de la muestra 1:

a)Derivar el polinomio interpolante para obtener la ecuación que describe su pendiente y graficarla

#Solicitamos a R que nos muestre el polinomio que describe el comportamiento de la muestra 1
poly.calc(dias, muestra1)
## 6.67 - 42.64348*x + 16.14272*x^2 - 2.094639*x^3 + 0.1269024*x^4 -  
## 0.003671679*x^5 + 4.094576e-05*x^6
deri_pol1 <- function(x){-42.64348+2*16.14272*x-3*2.094639*x^2+4*0.1269024*x^3-5*0.003671679*x^4+6*4.094576e-05*x^5}
x1 <- seq(0,28,0.01)
y1 <- deri_pol1(dias_seq)

graf_muestra_1 <- ggplot()+
  #Ejes
  geom_vline(xintercept = 0, linetype="dashed")+
  geom_hline(yintercept = 0, linetype="dashed")+
  
    #Pendiente
    geom_line(aes(x1, y1), color="dodgerblue3", size=1)+
  
  labs(x="dias", y="Peso Muestra", title="Pendiente de la muestra 1")+
  theme_bw()
ggplotly(graf_muestra_1)
  1. Una vez hecho esto, tenemos que calcular el punto máximo, el cual podemos identificar rápidamente debido a la combinación de ambas gráficas. Por ello, iteraremos usando el método de la secante y encontraremos la raíz de la derivada entre el día 10 y 11.
secant(deri_pol1, 9, 11)
## $root
## [1] 10.18858
## 
## $f.root
## [1] 2.345633e-09
## 
## $iter
## [1] 5
## 
## $estim.prec
## [1] 1.180142e-05
  1. Como podemos observar, el método de la secante ha encontrado una raíz, ahora obtendremos la imagen de este valor en la función que describe el peso de la muestra 1 para finalmente obtener el peso máximo que esta muestra puede tener:
pol_muestra1(10.18858)
## [1] 42.70842
  1. De esta manera, podemos observar que el peso máximo de la muestra 1 ocurre en el día \(10.18858\), en el cual, la muestra alcanza un peso de \(42.70842\)
  1. Calcular el peso máximo de la muestra 2:

a)Derivar el polinomio interpolante para obtener la ecuación que describe su pendiente y graficarla

#Solicitamos a R que nos muestre el polinomio que describe el comportamiento de la muestra 2
poly.calc(dias, muestra2)
## 6.67 - 5.678207*x + 2.912809*x^2 - 0.4137987*x^3 + 0.02584128*x^4 -  
## 0.0007525462*x^5 + 8.361598e-06*x^6
deri_pol2 <- function(x){-5.678207 + 2*2.912809*x - 3*0.4137987*x^2 + 4*0.02584128*x^3 - 5*0.0007525462*x^4 +6*8.361598e-06*x^5}
x2 <- seq(0,28,0.01)
y2 <- deri_pol2(dias_seq)

graf_muestra_2 <- ggplot()+
  #Ejes
  geom_vline(xintercept = 0, linetype="dashed")+
  geom_hline(yintercept = 0, linetype="dashed")+
  
    #Pendiente
    geom_line(aes(x2, y2), color="orange", size=1)+
  
  labs(x="dias", y="Peso Muestra", title="Pendiente de la muestra 2")+
  theme_bw()
ggplotly(graf_muestra_2)
  1. Una vez hecho esto, tenemos que calcular el punto máximo, el cual podemos identificar rápidamente debido a la combinación de ambas gráficas. Por ello, iteraremos usando el método de la secante y encontraremos la raíz de la derivada entre el día 7 y 9.
secant(deri_pol2, 7, 9)
## $root
## [1] 8.769426
## 
## $f.root
## [1] -1.309113e-10
## 
## $iter
## [1] 4
## 
## $estim.prec
## [1] 4.446447e-06
  1. Como podemos observar, el método de la secante ha encontrado una raíz, ahora obtendremos la imagen de este valor en la función que describe el peso de la muestra 2 para finalmente obtener el peso máximo que esta muestra puede tener:
pol_muestra2(8.769426)
## [1] 19.41575
  1. De esta manera, podemos observar que el peso máximo de la muestra 2 ocurre en el día \(8.769426\), en el cual, la muestra alcanza un peso de \(19.41575\)

Ejercicio 4

Construye los splines cúbicos con \(n\) nodos, donde \(n=3,4\) para las siguientes funciones en el intervalo dado.

  1. \(f(x) = e^{2x}\, cos 3x\), \([0,2]\).
#Función
f_4a <- function(x){exp(2*x)*cos(3*x)}
x4a <- seq(0, 2, 0.001)
y4a <- f_4a(x4a)

#Nodos
x_4an <- c(0, 1, 1.5, 2)
y_4an <- f_4a(x_4an)

#Spline cúbico
sp_4a <- cubicspline(x_4an, y_4an, x4a)

graf_4a <- ggplot()+
  geom_line(aes(x4a, y4a), color="black", size=1)+
  geom_point(aes(x_4an, y_4an), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x4a, sp_4a), color="red", size=1)+
  theme_bw()

ggplotly(graf_4a)
  1. \(f(x) = sen(log\,x)\), \([1,3]\).
#Función
f_4b <- function(x){sin(log(x))}
x4b <- seq(1, 3, 0.001)
y4b <- f_4b(x4b)

#Nodos
x_4bn <- c(1, 2, 3)
y_4bn <- f_4b(x_4bn)

#Spline cúbico
sp_4b <- cubicspline(x_4bn, y_4bn, x4b)

graf_4b <- ggplot()+
  geom_line(aes(x4b, y4b), color="black", size=1)+
  geom_point(aes(x_4bn, y_4bn), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x4b, sp_4b), color="red", size=1)+
  theme_bw()

ggplotly(graf_4b)
  1. \(f(x) = e^{x}+e^{-x}\), \([0,2]\).
#Función
f_4c <- function(x){exp(x)+exp(-x)}
x4c <- seq(0, 2, 0.001)
y4c <- f_4c(x4c)

#Nodos
x_4cn <- c(0, 2/3, 4/3, 2)
y_4cn <- f_4c(x_4cn)

#Spline cúbico
sp_4c <- cubicspline(x_4cn, y_4cn, x4c)

graf_4c <- ggplot()+
  geom_line(aes(x4c, y4c), color="black", size=1)+
  geom_point(aes(x_4cn, y_4cn), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x4c, sp_4c), color="red", size=1)+
  theme_bw()

ggplotly(graf_4c)
  1. \(f(x) = cos \,x+sen\,x\), \([0,2\pi]\).
#Función
f_4d <- function(x){cos(x)+sin(x)}
x4d <- seq(0, 2*pi, 0.001)
y4d <- f_4d(x4d)

#Nodos
x_4dn <- c(0, (2/3)*pi, (4/3)*pi, 2*pi)
y_4dn <- f_4d(x_4dn)

#Spline cúbico
sp_4d <- cubicspline(x_4dn, y_4dn, x4d)

graf_4d <- ggplot()+
  geom_line(aes(x4d, y4d), color="black", size=1)+
  geom_point(aes(x_4dn, y_4dn), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básciamente es la aproximación de esta función
  geom_line(aes(x4d, sp_4d), color="red", size=1)+
  theme_bw()

ggplotly(graf_4d)

Podemos observar que, cuando se trata de funciones trigonométricas, la aproximación de los splines cúbicos es muy mala. Esto puede ser explicado debido a la distinta naturaleza entre ambos tipos de funciones.

Ejercicio 5

Dada la partición \(x_0=0, x_1=0.5, x_2=1\), del intervalo \([0,1]\), encuentra el spline cúbico \(S\) para \(f(x)=e^{2x}\). Aproxima \(\int_0^{1} e^{2x}\,dx\) con \(\int_0^{1} S(x)\,dx\) y compara el resultado con el valor real.

#Función
f_5 <- function(x){exp(2*x)}
x5 <- seq(0, 0.75, 0.001)
y5 <- f_5(x5)

#Nodos
x_5n <- c(0, 0.25, 0.5, 0.75)
y_5n <- f_5(x_5n)

#Spline cúbico
sp_5 <- cubicspline(x_5n, y_5n, x5)

graf_5 <- ggplot()+
  geom_line(aes(x5, y5), color="black", size=1)+
  geom_point(aes(x_5n, y_5n), color="dodgerblue3", size=3)+
  #Este es el spline cúbico, y básicamente es la aproximación de esta función
  geom_line(aes(x5, sp_5), color="red", size=1)+
  geom_area(aes(x5, y5), fill="goldenrod1", alpha=0.5)+
  theme_bw()

ggplotly(graf_5)

La integral de la función es:

pracma::integral(f_5,0,1, )
## [1] 3.194528