#1. Repaso Código Integrales (aproximaciones de Área)

\[\int_{0}^{2e}ln(x)dx = lim_{b \to 0}\int_{b}^{2e}ln(x)dx \]

inte <- function(x) x^3
x1 <- seq(0,1,by = 0.01)
plot(x1, inte(x1),type='l', col = 'blue',lwd=2)
polygon(c(0,x1,1),c(0,inte(x1),0), col ='green')

#2. Suma trapezoidal

n <- 20
a <- 0
b <- 1
h <- (b-a)/n
xi <- a+(0:n)*h;xi
##  [1] 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70
## [16] 0.75 0.80 0.85 0.90 0.95 1.00
plot(x1, inte(x1),type='l', col = 'black',lwd=2, main = 'suma trapezoidal')
for(j in 1:n){
  segments(xi[j],0,xi[j+1],0,col = 'blue',lwd = 2)
  segments(xi[j+1],0,xi[j+1],inte(xi[j+1]),col = 'blue')
  segments(xi[j+1],inte(xi[j+1]),xi[j],inte(xi[j]),col = 'blue')
}

pot1 <- (h/2)*(inte(a)+inte(b))
pot2 <- (h/2)*(sum(2*inte(a+(1:(n-1)*h))))
protra <- pot1+pot2
protra
## [1] 0.250625

#3. Sumas de Riemann por la derecha

plot(x1, inte(x1),type='l', col = 'black',lwd=2, main = 'Riemann por la derecha')
for (l in 1:n) {
  segments(xi[l],0,xi[l+1],0,col = 'red',lwd=2)
  segments(xi[l+1],0,xi[l+1],inte(xi[l+1]),col = 'red',lwd=2)
  segments(xi[l+1],inte(xi[l+1]),xi[l],inte(xi[l+1]),col = 'red',lwd=2)
  segments(xi[l],inte(xi[l+1]),xi[l],0,col = 'red',lwd=2)
}

prode<-sum(inte(xi[-1]))*h
prode
## [1] 0.275625

Por la izquierda

plot(x1, inte(x1),type='l', col = 'black',lwd=2, main = 'Riemann por la izquierda')
for (l in 1:n) {
  segments(xi[l],0,xi[l+1],0,col = 'green',lwd=2)
  segments(xi[l],0,xi[l],inte(xi[l]),col = 'green',lwd=2)
  segments(xi[l+1],0,xi[l+1],inte(xi[l]),col = 'green',lwd=2)
  segments(xi[l+1],inte(xi[l]),xi[l],inte(xi[l]),col = 'green',lwd=2)
}

proiz<-sum(inte(xi[-(n+1)]))*h
proiz
## [1] 0.225625

Función de Montecarlo

plot(x1, inte(x1),type='l', col = 'black',lwd=2, main = 'Montecarlo')

mc.in <- function(n1,a,b){
  x <- runif(n1,a,b)
  y <- runif(n1,0,inte(b))
  co <- ifelse(y<=inte(x),'red','black')
  points(x,y,col = ifelse(y<=inte(x),'red','black'))
}
n1=5000
mc.in(n1,a,b)

x2 <- runif(n1,a,b)
y2 <- runif(n1,0,inte(b))
co <- ifelse(y2<=inte(x2),'red','black')
table(co)/n1
## co
##  black    red 
## 0.7568 0.2432
mc.3 <- function(n3,a,b){
  x <- runif(n3,a,b)
  f <- (b-a)*mean(inte(x))
  return(f)
}
n3=1000
mc.3(n3,a,b)
## [1] 0.2528943

Tarea

Resolver por simulación la siguiente integral

\[\int_{0}^{2e}ln(x)dx = lim_{b \to 0}\int_{b}^{2e}ln(x)dx \]

integralT <- function(x) log(x)
x4 <- seq(0,2*exp(1),by = 0.01)
plot(x4, integralT(x4),type='l', col = 'blue',lwd=2,ylab = 'y',xlab = 'x')
abline(v=1,h=0)

plot(x4, integralT(x4),type='l',ylim = c(0,2), col = 'black',lwd=2,ylab = 'y',xlab = 'x')
a <- 0.0001
b <- 2*exp(1)
puntos <- 10000
SimuInte <- function(puntos,a,b){
  x <- runif(puntos,a,b)
  y <- runif(puntos,0,integralT(b))
  co <- ifelse(y<=integralT(x),'yellow','black')
  points(x,y,col = ifelse(y<=integralT(x),'yellow','black'))
}
SimuInte(puntos,a,b)
rect(1,0,b,integralT(b), border = 'red')

x5 <- runif(puntos,a,b)
y5 <- runif(puntos,0,integralT(b))
co <- ifelse(y5<=integralT(x5),'yellow','black')
table(co)/puntos
## co
##  black yellow 
## 0.4843 0.5157
(table(co)/puntos)*(b-1)*(integralT(b))
## co
##    black   yellow 
## 3.637943 3.873812
plot(x4, integralT(x4),type='l',ylim = c(-10,2),xlim = c(0,1), col = 'black',lwd=2,ylab = 'y',xlab = 'x')
a <- 0.0001
b <- 1
puntos <- 50000
SimuInte <- function(puntos,a,b){
  x <- runif(puntos,a,b)
  y <- runif(puntos,-10,0)
  co <- ifelse(y<=integralT(x),'black','yellow')
  points(x,y,col = ifelse(y<=integralT(x),'black','yellow'))
}
SimuInte(puntos,a,b)
rect(0,-10,1,0,border = 'red')

x5 <- runif(puntos,a,b)
y5 <- runif(puntos,-10,0)
co <- ifelse(y5<=integralT(x5),'black','yellow')
(table(co)/puntos)*-4
## co
##    black   yellow 
## -3.60288 -0.39712