PROGRAMAFICIÓN
R-suelve tus dudas
PRÁCTICA 7:
Integración numérica en R
EJERCICIO 1
Primera parte
SOLUCIÓN:
SOLUCIÓN:
par(mfrow=c(2,1))
x=seq(A,B,0.01)
plot(x,d(x[1:length(x)]),xlim=c(A,B),ylim=c(0,8),xlab='x',ylab='rho')
par(new='TRUE')
plot(T,d(T[1:length(T)]),xlim=c(A,B),ylim=c(0,8),type='h',col='green',xlab='',ylab='')
d=function(x){exp(x)*sin(x)}
Simpson=function(h,n,T,A,B){
s1=0
s2=0
i=1
while(i<=(n-1)){
s1=s1+d(T[i])
i=i+1
}
i=2
while(i<=(n-1)){
s2=s2+d((T[i]+T[i+1])/2)
i=i+1
}
Mab=((h/6)*(d(A)+2*s1+4*s2+d(B)))
return(Mab)
}
A=0
B=3.05
n=35
h=(B-A)/(n-1)
T=0
for(j in 1:n){
T[j]=A+(j-1)*h
}
Masa=Simpson(h,n,T,A,B)
Masa
¡¡CUIDADO!!
ERRORES COMUNES:
-Poner length mal, poniendo lenght
-Olvidar poner los paréntesis en el bucle for 1: (n-1)
-El signo de la multiplicación exp(x)*sin(x)
-No dar nombre a la función para que pueda escribirse
-Poner los datos justo antes de llamar a la función
Segunda parte
SOLUCIÓN:
num_puntos=1000
AA=0
BB=7.4
ss=runif(num_puntos,A,B)
ff=runif(num_puntos,AA,BB)
puntos_dentro=0
puntos_fuera=0
color=0
for(i in 1:num_puntos){
if(ff[i]<=d(ss[i])){
puntos_dentro=puntos_dentro+1
color[i]='blue'
} else {
puntos_fuera=puntos_fuera+1
color[i]='orange'
}
}
¡¡CUIDADO!!
ERRORES COMUNES:
-Para dar valores aleatorios se utiliza el comando runif, no utilizar seq
-Confundir la i con 1
SOLUCIÓN:
plot(ss,ff,xlim=c(A,B),ylim=c(AA,BB),col=color,xlab='x',ylab='Densidad')
par(new='TRUE')
ss=seq(A,B,0.001)
plot(ss,d(ss[1:length(ss)]),type='l',col='dark green',xlim=c(A,B),ylim=c(AA,BB),xlab='',ylab='')
SOLUCIÓN:
Area=puntos_dentro/num_puntos*B*BB
Vexact=11.97907159
ErrorM=(Area-Vexact)/Vexact
ErrorS=(Masa-Vexact)/Vexact
Metodo=c('Valor exacto','Simpson','Montecarlo')
Valores=c(Vexact,Masa,Area)
Error_Relativo=c(0,round(ErrorS,7),round(ErrorM,7))
data.frame(Metodo,Valores,Error_Relativo)
Otra idea:
Haciendo unos cambios al código podemos ver como la naturaleza azarosa del método de montecarlo tiene su impacto en el cálculo de error según se aumenta el número de puntos:
d=function(x){
v=exp(x)*sin(x)
return(v)
}
simpson=function(a,b,d,t){
h=(b-a)/(length(t)-1)
s1=0
s2=0
for (i in 2:(length(t)-1)){
s1=s1+d(t[i])
}
for (i in 1:(length(t)-1)){
s2=s2+d((t[i]+t[i+1])/2)
}
masatotal=(h/6)*(d(a)+2*s1+4*s2+d(b))
return(masatotal)
}
num=10
a=0
b=3.05
t=seq(0,3.05,length = 35)
simpson=simpson(0,3.05,d,t)
par(mfrow=c(4,3))
x=seq(a,b,0.001)
plot(x,d(x),xlim=c(a,b),ylim=c(0,8),xlab='x',ylab='Densidad',type = 'l')
par(new='TRUE')
plot(t,d(t),xlim=c(a,b),ylim=c(0,8),type='h',col='green',xlab='',ylab='')
area=c(0)
for (j in 1:num) {
npuntos=1000*(j)^2
aa=0
bb=7.4
ss=runif(npuntos,a,b)
ff=runif(npuntos,aa,bb)
pdentro=0
pfuera=0
col=0
for (i in 1:npuntos) {
if(ff[i]<=d(ss[i])){
pdentro=pdentro+1
col[i]='blue'
}
else{
pfuera=pfuera+1
col[i]='orange'
}
}
plot(ss,ff,xlim=c(a,b),ylim=c(aa,bb),col=col,ylab='Densidad',xlab='x')
par(new='TRUE')
xx=seq(a,b,0.001)
plot(xx,d(xx),type='l',col='dark green',xlim=c(a,b),ylim=c(aa,bb),xlab='',ylab='')
area[j]=(pdentro/(npuntos))*b*bb
}
vexacto=11.97907159
em=((area-vexacto)/vexacto)
es=((simpson-vexacto)/vexacto)
metodo=c('valor exacto','simpson','montecarlo1','montecarlo2','montecarlo3','montecarlo4','montecarlo5','montecarlo6','montecarlo7','montecarlo8','montecarlo9','montecarlo10')
valores=c(vexacto,simpson,area)
erelativo=c(0,es,em)
data.frame(metodo,valores,erelativo)
La tabla quedaría:
metodo valores erelativo
1 valor exacto 11.97907 0.000000e+00
2 simpson 11.97907 -9.019708e-08
3 montecarlo1 12.54892 4.757033e-02
4 montecarlo2 12.08623 8.945886e-03
5 montecarlo3 12.05990 6.747747e-03
6 montecarlo4 11.96492 -1.181255e-03
7 montecarlo5 11.99641 1.447091e-03
8 montecarlo6 11.95959 -1.626117e-03
9 montecarlo7 11.95657 -1.878187e-03
10 montecarlo8 12.01923 3.352407e-03
11 montecarlo9 11.97910 2.134597e-06
12 montecarlo10 12.02123 3.519623e-03
EJERCICIO 2
SOLUCIÓN:
t=seq(0,3,length = 120)
y=c(0)
w=c(0)
y[1]=2
w[1]=2
r=2
k=1
h=(3/length(t))
f=function(r,k,t,y){
v=r*y*(1-(y/k))
return(v)
}
for (i in 1:(length(t)-1)) {
y[i+1]=y[i]+h*f(r,k,t[i],y[i])
z=w[i]+h*f(r,k,t[i],w[i])
w[i+1]=w[i]+(h/2)*(f(r,k,t[i],w[i])+f(r,k,t[i+1],z))
}
exacta=function(t){
(2*k)/(2+exp(-r*t)*k-2*exp(-r*t))
}
tt=seq(0,3, length=10000)
plot(tt,exacta(tt),type='l', col='blue')
par(new='TRUE')
plot(t,y,col='green',xlab='',ylab='')
par(new='TRUE')
plot(t,w,col='red',xlab='',ylab='')
¡¡CUIDADO!!
ERRORES COMUNES:
-No cerrar las llaves y paréntesis
-Poner plot=() y par=() cuando no se pone el igual
-Confundir . y ,
-Olvidarse de poner las componentes del vector [i]
Explicación en vídeo:
Para más información se puede acceder al repositorio a través del enlace