#Ejercicio 1
x=10
y=3
u=x+y;u
## [1] 13
v=x*y;v
## [1] 30
w=x/y;w
## [1] 3.333333
z=sin(x);z
## [1] -0.5440211
r=8*sin(y);r
## [1] 1.12896
s=5*sin(2*y);s
## [1] -1.397077
#Ejercicio 2
x=2; y=5
(y*x^3)/(x-y)
## [1] -13.33333
(3*x)/(2*y)
## [1] 0.6
3/2*x*y
## [1] 15
(x^5)/(x^5-1)
## [1] 1.032258
#Ejercicio 3
x=3
y=4
(1-1/x^5)^(-1)
## [1] 1.004132
3*pi*x^2
## [1] 84.823
3*y/(4*x-8)
## [1] 3
(4*(y-5))/(3*x-6)
## [1] -1.333333
#Ejercicio 4
x=2
y=6*x^3+4/x;y
## [1] 50
x=8
y=x/4*3;y
## [1] 6
x=10
y=(4*x)^2/(25);y
## [1] 64
x=2
y=2*(sin(x)/5);y
## [1] 0.363719
x=20
y=7*(x^(1/3))+4*x^(0.58);y
## [1] 41.73396
#Ejercicio 5
a=c(1.12);b=c(2.34);c=c(0.72);d=c(0.81);f=c(19.83)
x=1+a/b+c/(f^2);x
## [1] 1.480463
s=(b-a)/(d-c);s
## [1] 13.55556
r=1/(1/a+1/b+1/c+1/d);r
## [1] 0.2535713
y=a*b*(1/c)*(f^2/2);y
## [1] 715.6766
#Ejercicio 6
3/4*6*(7^2)+(4^5)/(7^3-145)
## [1] 225.6717
(48.2*55-9^3)/(53+14^2)
## [1] 7.718876
27^2/4+319^(4/5)/5+60*(14^(-3))
## [1] 202.412
#Ejercicio 7
r=5
V=(4*pi*(r^3))/3;V
## [1] 523.5988
V2=V*1.3;V2
## [1] 680.6784
R2=((V2*3)/(4*pi))^(1/3);R2
## [1] 5.456964
#Ejercicio # 8
x=-7-5i
y=4+3i
x+y
## [1] -3-2i
x*y
## [1] -13-41i
x/y
## [1] -1.72+0.04i
R=0.08206
n=1
T=273.2
V=22.4
#Aquí se comprueba el cálculo de la presión de 1mol de un gas ideal bajo las condiciones mencionadas.
P=(n*R*T)/V;P
## [1] 1.000839
R=0.08206
n=1
T=273.2
V=22.4
a=6.49
b=0.0562
P=(n*R*T)/(V-n*b)-(a*(n^2))/(V^2);P
## [1] 0.9904218
#Efecto del volumen molecular
P=(n*R*T)/(V-n*b);P
## [1] 1.003356
#Efecto de la atracción molecular
P=(n*R*T)/(V)-(a*(n^2))/(V^2);P
## [1] 0.9879045
#La principal causa de la desviación del Cl2 como gas ideal se ejerce por el efecto de la atracción molecular, ya que disminuye la presión en 0,012 atm, mientras que solo el volumen molecular aumenta la presión en 0.003 atm
M1=7.3
E1=(10^(4.4))*(10^(1.5*M1));E1
## [1] 2.238721e+15
M2=5.5
E2=(10^(4.4))*(10^(1.5*M2));E2
## [1] 4.466836e+12
D=E1-E2;D
## [1] 2.234254e+15
#Un terremoto de escala 7.3 libera 2.23 x 10^15 joules más que un terremoto de escala 5.5
t= seq(1,3,0.1)
T=6*log(t)-7*(exp(1))^(0.2*t)
plot(t,T,main="Temperatura vs tiempo",sub="Métodos Multivariados 2023-1S",xlab="t (Min)", ylab="T (C)",col="blue",type="o")
#install.packages("ggplot2")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
#install.packages("reshape2")
library(reshape2)
x=seq(0,2,0.05)
u=2*log10(60*x+1)
v=3*cos(6*x)
df=data.frame(x, u, v)
df2=melt(data = df, id.vars = "x")
ggplot(data = df2, aes(x = x, y = value, colour = variable)) +
geom_line()+
xlab("Distancia (mi)")+
ylab("Velocidad (mi/h)")+
ggtitle("Gráfico Velocidad / distancia para 2 funciones")
Opción 2
x=seq(0,2,0.05)
u=2*(log10(60*x+1))
v=3*cos(6*x)
plot(x,v,main="Speed vs Distance", xlab="Distance (mi)", ylab="speed (mi/h)", ylim=c(-4,5), col="green", type="o")
lines(x,u, col="blue")
legend("topright", legend=c("v","u"), fill=c("green", "blue"))
\[A=W\times L+\frac{D^{2}}2\]
\[D=\sqrt{\frac{W^{2}}2{}}\]
\[A=W\times L+\frac{W^{2}}4{}\]
\[L=\frac{A}{W}-\frac{W}4{}\]
W=6
A=80
L=(A/W)-(W/4);L
## [1] 11.83333
#Longitud total de la cerca (LF)
D=sqrt((W^2)/2)
LF=2*L+2*D+W;LF
## [1] 38.15195
#Para generar el requerimiento de que el usuario ingrese el valor de cada variable se puede usar la funcion readline() aunque al hacerlo no corre RPubs.
#W=readline();
#A=readline()
#Y aquí se colocan las funciones deseadas, el usuario debe ingresar el dato en la Consola
#14.a
#Método 1 seq() function
d=(28-5)/99
seq(from=5,to=28,by=d)
## [1] 5.000000 5.232323 5.464646 5.696970 5.929293 6.161616 6.393939
## [8] 6.626263 6.858586 7.090909 7.323232 7.555556 7.787879 8.020202
## [15] 8.252525 8.484848 8.717172 8.949495 9.181818 9.414141 9.646465
## [22] 9.878788 10.111111 10.343434 10.575758 10.808081 11.040404 11.272727
## [29] 11.505051 11.737374 11.969697 12.202020 12.434343 12.666667 12.898990
## [36] 13.131313 13.363636 13.595960 13.828283 14.060606 14.292929 14.525253
## [43] 14.757576 14.989899 15.222222 15.454545 15.686869 15.919192 16.151515
## [50] 16.383838 16.616162 16.848485 17.080808 17.313131 17.545455 17.777778
## [57] 18.010101 18.242424 18.474747 18.707071 18.939394 19.171717 19.404040
## [64] 19.636364 19.868687 20.101010 20.333333 20.565657 20.797980 21.030303
## [71] 21.262626 21.494949 21.727273 21.959596 22.191919 22.424242 22.656566
## [78] 22.888889 23.121212 23.353535 23.585859 23.818182 24.050505 24.282828
## [85] 24.515152 24.747475 24.979798 25.212121 25.444444 25.676768 25.909091
## [92] 26.141414 26.373737 26.606061 26.838384 27.070707 27.303030 27.535354
## [99] 27.767677 28.000000
#Método 2 usando seq() con length.out
seq(5,28,length.out=100)
## [1] 5.000000 5.232323 5.464646 5.696970 5.929293 6.161616 6.393939
## [8] 6.626263 6.858586 7.090909 7.323232 7.555556 7.787879 8.020202
## [15] 8.252525 8.484848 8.717172 8.949495 9.181818 9.414141 9.646465
## [22] 9.878788 10.111111 10.343434 10.575758 10.808081 11.040404 11.272727
## [29] 11.505051 11.737374 11.969697 12.202020 12.434343 12.666667 12.898990
## [36] 13.131313 13.363636 13.595960 13.828283 14.060606 14.292929 14.525253
## [43] 14.757576 14.989899 15.222222 15.454545 15.686869 15.919192 16.151515
## [50] 16.383838 16.616162 16.848485 17.080808 17.313131 17.545455 17.777778
## [57] 18.010101 18.242424 18.474747 18.707071 18.939394 19.171717 19.404040
## [64] 19.636364 19.868687 20.101010 20.333333 20.565657 20.797980 21.030303
## [71] 21.262626 21.494949 21.727273 21.959596 22.191919 22.424242 22.656566
## [78] 22.888889 23.121212 23.353535 23.585859 23.818182 24.050505 24.282828
## [85] 24.515152 24.747475 24.979798 25.212121 25.444444 25.676768 25.909091
## [92] 26.141414 26.373737 26.606061 26.838384 27.070707 27.303030 27.535354
## [99] 27.767677 28.000000
#14.b
#Metodo 1 seq() function
seq(from=2, to=14, by=0.2)
## [1] 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8
## [16] 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8
## [31] 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 10.2 10.4 10.6 10.8
## [46] 11.0 11.2 11.4 11.6 11.8 12.0 12.2 12.4 12.6 12.8 13.0 13.2 13.4 13.6 13.8
## [61] 14.0
#Método 2 seq() function con lenght out
d=(14-2)/0.2+1
seq(2,14,length.out=d)
## [1] 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8
## [16] 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8
## [31] 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0 10.2 10.4 10.6 10.8
## [46] 11.0 11.2 11.4 11.6 11.8 12.0 12.2 12.4 12.6 12.8 13.0 13.2 13.4 13.6 13.8
## [61] 14.0
#14.c
#Método 1 seq () con length
seq(-2,5,length.out=50)
## [1] -2.0000000 -1.8571429 -1.7142857 -1.5714286 -1.4285714 -1.2857143
## [7] -1.1428571 -1.0000000 -0.8571429 -0.7142857 -0.5714286 -0.4285714
## [13] -0.2857143 -0.1428571 0.0000000 0.1428571 0.2857143 0.4285714
## [19] 0.5714286 0.7142857 0.8571429 1.0000000 1.1428571 1.2857143
## [25] 1.4285714 1.5714286 1.7142857 1.8571429 2.0000000 2.1428571
## [31] 2.2857143 2.4285714 2.5714286 2.7142857 2.8571429 3.0000000
## [37] 3.1428571 3.2857143 3.4285714 3.5714286 3.7142857 3.8571429
## [43] 4.0000000 4.1428571 4.2857143 4.4285714 4.5714286 4.7142857
## [49] 4.8571429 5.0000000
#Método 2 seq con los intervalos previamente calculados
d=(5+2)/49
seq(from=-2,to=5,by=d)
## [1] -2.0000000 -1.8571429 -1.7142857 -1.5714286 -1.4285714 -1.2857143
## [7] -1.1428571 -1.0000000 -0.8571429 -0.7142857 -0.5714286 -0.4285714
## [13] -0.2857143 -0.1428571 0.0000000 0.1428571 0.2857143 0.4285714
## [19] 0.5714286 0.7142857 0.8571429 1.0000000 1.1428571 1.2857143
## [25] 1.4285714 1.5714286 1.7142857 1.8571429 2.0000000 2.1428571
## [31] 2.2857143 2.4285714 2.5714286 2.7142857 2.8571429 3.0000000
## [37] 3.1428571 3.2857143 3.4285714 3.5714286 3.7142857 3.8571429
## [43] 4.0000000 4.1428571 4.2857143 4.4285714 4.5714286 4.7142857
## [49] 4.8571429 5.0000000
#15.a
#install.packages("emdbook")
library(emdbook)
lseq(from = 10,to = 1000,length.out = 50)
## [1] 10.00000 10.98541 12.06793 13.25711 14.56348 15.99859
## [7] 17.57511 19.30698 21.20951 23.29952 25.59548 28.11769
## [13] 30.88844 33.93222 37.27594 40.94915 44.98433 49.41713
## [19] 54.28675 59.63623 65.51286 71.96857 79.06043 86.85114
## [25] 95.40955 104.81131 115.13954 126.48552 138.94955 152.64180
## [31] 167.68329 184.20700 202.35896 222.29965 244.20531 268.26958
## [37] 294.70517 323.74575 355.64803 390.69399 429.19343 471.48664
## [43] 517.94747 568.98660 625.05519 686.64885 754.31201 828.64277
## [49] 910.29818 1000.00000
#15.b
lseq(from = 10,to = 1000,length.out = 20)
## [1] 10.00000 12.74275 16.23777 20.69138 26.36651 33.59818
## [7] 42.81332 54.55595 69.51928 88.58668 112.88379 143.84499
## [13] 183.29807 233.57215 297.63514 379.26902 483.29302 615.84821
## [19] 784.75997 1000.00000
vector_a = c(3,-5,6,15)
vector_b = c(7,9,13,5)
vector_c = c(-4,10,8,4)
vector_d = c(12,2,11,1)
A <- cbind(vector_a,vector_b,vector_c,vector_d)
A
## vector_a vector_b vector_c vector_d
## [1,] 3 7 -4 12
## [2,] -5 9 10 2
## [3,] 6 13 8 11
## [4,] 15 5 4 1
class(A)
## [1] "matrix" "array"
#16.a
v=A[,2];v
## [1] 7 9 13 5
#16.b
w=A[2,];w
## vector_a vector_b vector_c vector_d
## -5 9 10 2
#17.a
B=A[,c(2:4)];B
## vector_b vector_c vector_d
## [1,] 7 -4 12
## [2,] 9 10 2
## [3,] 13 8 11
## [4,] 5 4 1
#17.b
C=A[c(2:4),];C
## vector_a vector_b vector_c vector_d
## [1,] -5 9 10 2
## [2,] 6 13 8 11
## [3,] 15 5 4 1
#17.c
D=A[c(1:2),c(2:4)];D
## vector_b vector_c vector_d
## [1,] 7 -4 12
## [2,] 9 10 2
#18.a
x=c(2,4,7)
length(x)
## [1] 3
abs(x)
## [1] 2 4 7
#18.b
y=c(2,-4,7)
length(y)
## [1] 3
abs(y)
## [1] 2 4 7
#19.a
MAXC1=max(A[,1]);MAXC1
## [1] 15
MAXC2=max(A[,2]);MAXC2
## [1] 13
MAXC3=max(A[,3]);MAXC3
## [1] 10
MAXC4=max(A[,4]);MAXC4
## [1] 12
MINC1=min(A[,1]);MINC1
## [1] -5
MINC2=min(A[,2]);MINC2
## [1] 5
MINC3=min(A[,3]);MINC3
## [1] -4
MINC4=min(A[,4]);MINC4
## [1] 1
#19.b
MAXR1=max(A[1,]);MAXR1
## [1] 12
MAXR2=max(A[2,]);MAXR2
## [1] 10
MAXR3=max(A[3,]);MAXR3
## [1] 13
MAXR4=max(A[4,]);MAXR4
## [1] 15
MINR1=min(A[1,]);MINR1
## [1] -4
MINR2=min(A[2,]);MINR2
## [1] -5
MINR3=min(A[3,]);MINR3
## [1] 6
MINR4=min(A[4,]);MINR4
## [1] 1
#20.a
C1=sort(A[,1])
C2=sort(A[,2])
C3=sort(A[,3])
C4=sort(A[,4])
B=cbind(C1,C2,C3,C4);B
## C1 C2 C3 C4
## [1,] -5 5 -4 1
## [2,] 3 7 4 2
## [3,] 6 9 8 11
## [4,] 15 13 10 12
#20.b
R1=sort(A[1,])
R2=sort(A[2,])
R3=sort(A[3,])
R4=sort(A[4,])
C=rbind(R1,R2,R3,R4);C
## vector_c vector_a vector_b vector_d
## R1 -4 3 7 12
## R2 -5 2 9 10
## R3 6 8 11 13
## R4 1 4 5 15
#20.c
D=as.array(colSums(A));D
## vector_a vector_b vector_c vector_d
## 19 34 18 26
E=as.array(rowSums(A));E
## [1] 18 16 38 25
#Ejercio #21
A=matrix(data = c(1,2,7,3,4,4,9,pi,2,100,7,42),nrow = 4,ncol = 3)
B=log(A);B
## [,1] [,2] [,3]
## [1,] 0.0000000 1.386294 0.6931472
## [2,] 0.6931472 1.386294 4.6051702
## [3,] 1.9459101 2.197225 1.9459101
## [4,] 1.0986123 1.144730 3.7376696
class(B)
## [1] "matrix" "array"
#21.a
a=B[2,];a
## [1] 0.6931472 1.3862944 4.6051702
#21.b
sum(a)
## [1] 6.684612
#Duda de si esto es correcto*
#21.c
c=B[,2]*A[,1];c
## [1] 1.386294 2.772589 15.380572 3.434190
#21.d
max(c)
## [1] 15.38057
#21.e
d=c(B[1,3],B[2,3],B[3,3]);d
## [1] 0.6931472 4.6051702 1.9459101
e=A[1,]/d;e
## [1] 1.442695 0.868589 1.027797
sum(e)
## [1] 3.339081
#Ejercicio 22
#22.a
A=matrix(data = c(3,6,7,-2,8,9,1,-5,10),nrow = 3,ncol = 3)
B=matrix(data = c(6,7,-8,9,5,2,-4,3,1),nrow = 3,ncol = 3)
C=matrix(data = c(-7,10,3,-5,6,-9,2,1,8),nrow = 3,ncol = 3)
D=array(c(A,B,C),dim = c(3,3,3));D
## , , 1
##
## [,1] [,2] [,3]
## [1,] 3 -2 1
## [2,] 6 8 -5
## [3,] 7 9 10
##
## , , 2
##
## [,1] [,2] [,3]
## [1,] 6 9 -4
## [2,] 7 5 3
## [3,] -8 2 1
##
## , , 3
##
## [,1] [,2] [,3]
## [1,] -7 -5 2
## [2,] 10 6 1
## [3,] 3 -9 8
dim(D)
## [1] 3 3 3
#Consultar
#22.b
max(D[,,1])
## [1] 10
max(D[,,2])
## [1] 9
max(D[,,3])
## [1] 10
max(D)
## [1] 10
#Ejercicio 23
A=matrix(data = c(-7,4,16,9),nrow = 2,ncol = 2)
B=matrix(data = c(6,12,-5,-2),nrow = 2,ncol = 2)
C=matrix(data = c(-3,6,-9,8),nrow = 2,ncol = 2)
#23.a
A+B+C
## [,1] [,2]
## [1,] -4 2
## [2,] 22 15
#23.b
A-B+C
## [,1] [,2]
## [1,] -16 12
## [2,] -2 19
#23.c
(A+B)+C
## [,1] [,2]
## [1,] -4 2
## [2,] 22 15
A+(B+C)
## [,1] [,2]
## [1,] -4 2
## [2,] 22 15
#Ejercicio 24
F=matrix(data = c(400,550,700,500,600),nrow = 1,ncol = 5);F
## [,1] [,2] [,3] [,4] [,5]
## [1,] 400 550 700 500 600
D=matrix(data = c(2,0.5,0.75,1.5,3),nrow = 1,ncol = 5);D
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2 0.5 0.75 1.5 3
#24.a
W=F*D;W
## [,1] [,2] [,3] [,4] [,5]
## [1,] 800 275 525 750 1800
#24.b
sum(W)
## [1] 4150
#Ahora con la función rowSums
rowSums(W)
## [1] 4150
W=matrix(data = c(5,40,1000,5.50,43,1100,6.50,37,1000,6,50,1200,6.25,45,1100),nrow = 3,ncol = 5);W
## [,1] [,2] [,3] [,4] [,5]
## [1,] 5 5.5 6.5 6 6.25
## [2,] 40 43.0 37.0 50 45.00
## [3,] 1000 1100.0 1000.0 1200 1100.00
#25.a
PT=W[1,]*W[2,];PT
## [1] 200.00 236.50 240.50 300.00 281.25
#25.b
sum(PT)
## [1] 1258.25
#25.c
sum(W[3,])
## [1] 5400
D1=c(-60,-25,30)
D2=c(-30,-55,20)
#26.a
euclidiana <- sqrt(sum((D1)^2));euclidiana
## [1] 71.58911
#26.b
D=rbind(D1,D2);D
## [,1] [,2] [,3]
## D1 -60 -25 30
## D2 -30 -55 20
D1_D2=D[2,]-D[1,];D1_D2
## [1] 30 -30 -10
#El buzo 1 debe moverse 30ft hacia el oeste, 30 ft hacia el sur y 10 ft hacia arriba para alcanzar al buzo 2
#26.c
euclidiana2 <- sqrt(sum((D1_D2)^2));euclidiana2
## [1] 43.58899
SP=matrix(data = c(11,1000,7,800,8,900,10,1200,9,700),nrow = 2,ncol = 5);SP
## [,1] [,2] [,3] [,4] [,5]
## [1,] 11 7 8 10 9
## [2,] 1000 800 900 1200 700
#27.a
x=SP[1,]/SP[2,];x
## [1] 0.011000000 0.008750000 0.008888889 0.008333333 0.012857143
#27.b
EP=(SP[2,]*(x^2))/2;EP
## [1] 0.06050000 0.03062500 0.03555556 0.04166667 0.05785714
#28.a
P=matrix(data = c(300,550,400,250,500),nrow = 5,ncol = 1);P
## [,1]
## [1,] 300
## [2,] 550
## [3,] 400
## [4,] 250
## [5,] 500
M=matrix(data = c(5,3,6,3,2,4,2,5,5,4,6,4,3,4,3),nrow = 5,ncol = 3);M
## [,1] [,2] [,3]
## [1,] 5 4 6
## [2,] 3 2 4
## [3,] 6 5 3
## [4,] 3 5 4
## [5,] 2 4 3
#28.b
GT=c(P*M[,1],P*M[,2],P*M[,3]);GT
## [1] 1500 1650 2400 750 1000 1200 1100 2000 1250 2000 1800 2200 1200 1000 1500
MGT=matrix(data = GT,nrow=5,ncol = 3);MGT
## [,1] [,2] [,3]
## [1,] 1500 1200 1800
## [2,] 1650 1100 2200
## [3,] 2400 2000 1200
## [4,] 750 1250 1000
## [5,] 1000 2000 1500
GMAY=sum(MGT[,1]);GMAY
## [1] 7300
GJUN=sum(MGT[,2]);GJUN
## [1] 7550
GJUL=sum(MGT[,3]);GJUL
## [1] 7700
#28.c
G_M=rowSums(MGT);G_M
## [1] 4500 4950 5600 3000 4500
#28.d
GTMESES=sum(MGT);GTMESES
## [1] 22550
#Ahora se verifica sumando cada mes
GTMESES2=GMAY+GJUN+GJUL;GTMESES2
## [1] 22550
A=matrix(data=c(11,-9,5,-4),nrow = 2,ncol = 2);A
## [,1] [,2]
## [1,] 11 5
## [2,] -9 -4
B=matrix(data=c(-7,6,-8,2),nrow = 2,ncol = 2);B
## [,1] [,2]
## [1,] -7 -8
## [2,] 6 2
#29.a
A%*%B
## [,1] [,2]
## [1,] -47 -78
## [2,] 39 64
#29.b
B%*%A
## [,1] [,2]
## [1,] -5 -3
## [2,] 48 22
A=matrix(data=c(3,6,7,-2,8,9,1,-5,10),nrow = 3,ncol = 3);A
## [,1] [,2] [,3]
## [1,] 3 -2 1
## [2,] 6 8 -5
## [3,] 7 9 10
B=matrix(data=c(6,7,-8,9,5,2,-4,3,1),nrow = 3,ncol = 3);B
## [,1] [,2] [,3]
## [1,] 6 9 -4
## [2,] 7 5 3
## [3,] -8 2 1
C=matrix(data=c(-7,10,3,-5,6,-9,2,1,8),nrow = 3,ncol = 3);C
## [,1] [,2] [,3]
## [1,] -7 -5 2
## [2,] 10 6 1
## [3,] 3 -9 8
#30.a
M1=A%*%(B+C);M1
## [,1] [,2] [,3]
## [1,] -42 -17 -5
## [2,] 155 147 -25
## [3,] 96 57 112
M2=A%*%B+A%*%C;M2
## [,1] [,2] [,3]
## [1,] -42 -17 -5
## [2,] 155 147 -25
## [3,] 96 57 112
ifelse(M1==M2,"Son iguales","No son iguales")
## [,1] [,2] [,3]
## [1,] "Son iguales" "Son iguales" "Son iguales"
## [2,] "Son iguales" "Son iguales" "Son iguales"
## [3,] "Son iguales" "Son iguales" "Son iguales"
#30.b
M3=(A%*%B)%*%C;M3
## [,1] [,2] [,3]
## [1,] 167 287 -125
## [2,] -99 -111 308
## [3,] 1132 562 250
M4=A%*%(B%*%C);M4
## [,1] [,2] [,3]
## [1,] 167 287 -125
## [2,] -99 -111 308
## [3,] 1132 562 250
ifelse(M3==M4,"Son iguales","No son iguales")
## [,1] [,2] [,3]
## [1,] "Son iguales" "Son iguales" "Son iguales"
## [2,] "Son iguales" "Son iguales" "Son iguales"
## [3,] "Son iguales" "Son iguales" "Son iguales"
#31.a
A=matrix(data=c(2,3,1,-9),nrow = 2,ncol = 2);A
## [,1] [,2]
## [1,] 2 1
## [2,] 3 -9
B=matrix(data=c(5,2),nrow=2,ncol = 1);B
## [,1]
## [1,] 5
## [2,] 2
#Usando la función solve
SOL=solve(A,B);SOL
## [,1]
## [1,] 2.2380952
## [2,] 0.5238095
#Sacando la inversa usando la libreria matlib
library(matlib)
INVA=inv(A);INVA
## [,1] [,2]
## [1,] 0.4285714 0.04761905
## [2,] 0.1428571 -0.09523810
INVA%*%B
## [,1]
## [1,] 2.2380952
## [2,] 0.5238095
#31.b
A=matrix(data=c(-8,-5,-2,7),nrow = 2,ncol = 2,byrow = TRUE);A
## [,1] [,2]
## [1,] -8 -5
## [2,] -2 7
B=matrix(data=c(4,10),nrow=2,ncol = 1);B
## [,1]
## [1,] 4
## [2,] 10
#Usando la función solve
SOL=solve(A,B);SOL
## [,1]
## [1,] -1.181818
## [2,] 1.090909
#Sacando la inversa usando la libreria matlib
library(matlib)
INVA=inv(A);INVA
## [,1] [,2]
## [1,] -0.10606061 -0.07575758
## [2,] -0.03030303 0.12121212
INVA%*%B
## [,1]
## [1,] -1.181818
## [2,] 1.090909
#31.c
A=matrix(data=c(12,-5,0,-3,4,7,6,2,3),nrow = 3,ncol = 3,byrow = TRUE);A
## [,1] [,2] [,3]
## [1,] 12 -5 0
## [2,] -3 4 7
## [3,] 6 2 3
B=matrix(data=c(11,-3,22),nrow=3,ncol = 1);B
## [,1]
## [1,] 11
## [2,] -3
## [3,] 22
#Usando la función solve
SOL=solve(A,B);SOL
## [,1]
## [1,] 3
## [2,] 5
## [3,] -2
#Sacando la inversa usando la libreria matlib
library(matlib)
INVA=inv(A);INVA
## [,1] [,2] [,3]
## [1,] 0.00716846 -0.05376344 0.1254480
## [2,] -0.18279570 -0.12903226 0.3010753
## [3,] 0.10752688 0.19354839 -0.1182796
INVA%*%B
## [,1]
## [1,] 3
## [2,] 5
## [3,] -2
#31.d
A=matrix(data=c(6,-3,4,12,5,-7,-5,2,6),nrow = 3,ncol = 3,byrow = TRUE);A
## [,1] [,2] [,3]
## [1,] 6 -3 4
## [2,] 12 5 -7
## [3,] -5 2 6
B=matrix(data=c(41,-26,14),nrow=3,ncol = 1);B
## [,1]
## [1,] 41
## [2,] -26
## [3,] 14
#Usando la función solve
SOL=solve(A,B);SOL
## [,1]
## [1,] 2
## [2,] -3
## [3,] 5
#Sacando la inversa usando la libreria matlib
library(matlib)
INVA=inv(A);INVA
## [,1] [,2] [,3]
## [1,] 0.07705779 0.04553415 0.00175131
## [2,] -0.06479860 0.09807356 0.15761821
## [3,] 0.08581436 0.00525394 0.11558669
INVA%*%B
## [,1]
## [1,] 2
## [2,] -3
## [3,] 5
32.a \[A(BC+A)=B\] Multiplicando por la izquierda la matriz inversa de A en ambos lados de la ecuación Se debe tener mucho cuidado con el orden en que se hace la multiplicación porque en matrices no hay propiedad conmutativa \[A^{-1}A(BC+A)=A^{-1}B\] Siendo I la matriz identidad \[I(BC+A)=A^{-1}B\] \[BC+A=A^{-1}B\] \[BC=A^{-1}B-A\] Hacemos lo mismo con la matriz inversa de B, multiplicando por la izquierda a ambos lados \[B^{-1}BC=B^{-1}(A^{-1}B-A)\] \[IC=B^{-1}(A^{-1}B-A)\] \[C=B^{-1}(A^{-1}B-A)\]
#32.b
A=matrix(data = c(3,-2,9,4),nrow = 2,ncol = 2);A
## [,1] [,2]
## [1,] 3 9
## [2,] -2 4
B=matrix(data = c(2,7,-3,6),nrow = 2,ncol = 2);B
## [,1] [,2]
## [1,] 2 -3
## [2,] 7 6
#Se verifica que el determinante de A y B no sea cero
det(A)
## [1] 30
det(B)
## [1] 33
INVA=solve(A);INVA
## [,1] [,2]
## [1,] 0.13333333 -0.3
## [2,] 0.06666667 0.1
A%*%INVA
## [,1] [,2]
## [1,] 1 1.110223e-16
## [2,] 0 1.000000e+00
INVB=solve(B);INVB
## [,1] [,2]
## [1,] 0.1818182 0.09090909
## [2,] -0.2121212 0.06060606
B%*%INVB
## [,1] [,2]
## [1,] 1 5.551115e-17
## [2,] 0 1.000000e+00
C=INVB%*%((INVA%*%B)-A);C
## [,1] [,2]
## [1,] -0.6212121 -2.363636
## [2,] 1.1969697 2.157576
#Aquí se verifica el resultado
SOL=A%*%(B%*%C+A);SOL
## [,1] [,2]
## [1,] 2 -3
## [2,] 7 6
B
## [,1] [,2]
## [1,] 2 -3
## [2,] 7 6
#Sol es igual a B por tanto se comprueba el resultado
#33.a
library(matlib)
#Aquí se observan sistemas de ecuaciones sin solución
A=matrix(data=c(-2,1,-2,1),nrow = 2,ncol = 2,byrow = TRUE);A
## [,1] [,2]
## [1,] -2 1
## [2,] -2 1
B=matrix(data=c(-5,3),nrow=2,ncol = 1);B
## [,1]
## [1,] -5
## [2,] 3
#Usando la función solve de la siguiente manera se obtiene el siguiente reporte de R
#SOL=solve(A,B);SOL
#Error in solve.default(A, B) :
#Lapack routine dgesv: system is exactly singular: U[2,2] = 0
#Se coloca a continuación la función de matlib eliminación gausiana para mostrar el paso a paso
gaussianElimination(A = A,B = B,verbose = TRUE)
##
## Initial matrix:
## [,1] [,2] [,3]
## [1,] -2 1 -5
## [2,] -2 1 3
##
## row: 1
##
## multiply row 1 by -0.5
## [,1] [,2] [,3]
## [1,] 1 -0.5 2.5
## [2,] -2 1.0 3.0
##
## multiply row 1 by 2 and add to row 2
## [,1] [,2] [,3]
## [1,] 1 -0.5 2.5
## [2,] 0 0.0 8.0
##
## row: 2
#Se observa que se llega a una solución inconsistente con 0=8
#33.b
#Aquí se encuentra un sistema con infinitas soluciones ya que la ecuación 2 es la misma 1 pero multiplicada por 4
A=matrix(data=c(-2,1,-8,4),nrow = 2,ncol = 2,byrow = TRUE);A
## [,1] [,2]
## [1,] -2 1
## [2,] -8 4
B=matrix(data=c(3,12),nrow=2,ncol = 1);B
## [,1]
## [1,] 3
## [2,] 12
#SOL=solve(A,B);SOL
#Error in solve.default(A, B) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0
#Se coloca a continuación la función de matlib eliminación gausiana para mostrar el paso a paso
gaussianElimination(A = A,B = B,verbose = TRUE)
##
## Initial matrix:
## [,1] [,2] [,3]
## [1,] -2 1 3
## [2,] -8 4 12
##
## row: 1
##
## exchange rows 1 and 2
## [,1] [,2] [,3]
## [1,] -8 4 12
## [2,] -2 1 3
##
## multiply row 1 by -0.125
## [,1] [,2] [,3]
## [1,] 1 -0.5 -1.5
## [2,] -2 1.0 3.0
##
## multiply row 1 by 2 and add to row 2
## [,1] [,2] [,3]
## [1,] 1 -0.5 -1.5
## [2,] 0 0.0 0.0
##
## row: 2
#Con este método se observan soluciones infinitas, al quedar una fila con valor de 0
#33.c
#Aquí vemos un sistema sin solución
A=matrix(data=c(-2,1,-2,1),nrow = 2,ncol = 2,byrow = TRUE);A
## [,1] [,2]
## [1,] -2 1
## [2,] -2 1
B=matrix(data=c(-5,-5.00001),nrow=2,ncol = 1);B
## [,1]
## [1,] -5.00000
## [2,] -5.00001
#SOL=solve(A,B);SOL
#Error in solve.default(A, B) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0
#Se coloca a continuación la función de matlib eliminación gausiana para mostrar el paso a paso
gaussianElimination(A = A,B = B,verbose = TRUE)
##
## Initial matrix:
## [,1] [,2] [,3]
## [1,] -2 1 -5.00000
## [2,] -2 1 -5.00001
##
## row: 1
##
## multiply row 1 by -0.5
## [,1] [,2] [,3]
## [1,] 1 -0.5 2.50000
## [2,] -2 1.0 -5.00001
##
## multiply row 1 by 2 and add to row 2
## [,1] [,2] [,3]
## [1,] 1 -0.5 2.50000
## [2,] 0 0.0 -0.00001
##
## row: 2
#Esto llega a una solución inconsistente al tener una fila con resultado 0 = -0.00001
ft_i=1
m_i=65
N_i=18
lbf_i=5
sl_i=6
kg_i=15
U_i=c(ft_i,m_i,N_i,lbf_i,sl_i,kg_i)
m=ft_i*0.3048
ft=m_i/0.3048
lbf=N_i/4.4482216152605
N=lbf_i*4.4482216152605
kg=sl_i*14.594
sl=kg_i/14.594
CONV=c(m,ft,lbf,N,kg,sl)
dataframe=data.frame(U_i,U=c("ft","m","N_i","lbf","slug","kg"),CONV,U_f=c("m","ft","lbf","N","kg","slug"));dataframe
## U_i U CONV U_f
## 1 1 ft 0.304800 m
## 2 65 m 213.254593 ft
## 3 18 N_i 4.046561 lbf
## 4 5 lbf 22.241108 N
## 5 6 slug 87.564000 kg
## 6 15 kg 1.027820 slug
class(dataframe)
## [1] "data.frame"
#MC=matrix(ft_i,m_i,N_i,lbf_i,sl_i,kg_i,ft,m,N,lbf,sl,kg,byrow = TRUE,nrow=2,ncol=6);MC
X=as.matrix(dataframe)
class(X)
## [1] "matrix" "array"
#35.a
x=5
y=8
RAD=atan2(y,x)
GRAD=RAD*57.2958;GRAD
## [1] 57.99464
#35.b
x=-5
y=8
RAD=atan2(y,x)
GRAD=RAD*57.2958;GRAD
## [1] 122.0054
#35.c
x=5
y=-8
RAD=atan2(y,x)
GRAD=RAD*57.2958;GRAD
## [1] -57.99464
#35.d
x=-5
y=-8
RAD=atan2(y,x)
GRAD=RAD*57.2958;GRAD
## [1] -122.0054
Fa=80
C=5/9*(Fa-32);C
## [1] 26.66667
Grados_C=function(Fa)
{5/9*(Fa-32); return(C)}
Grados_C(80)
## [1] 26.66667
t1=function(vo,h,g){t=ifelse((-vo+sqrt(vo^2+2*g*h))/g>=0,(-vo+sqrt(vo^2+2*g*h))/g,(-vo-sqrt(vo^2+2*g*h))/g);return(t)}
t1(vo = 50,h = 100,g = 9.81)
## [1] 1.712355
#t2=function(vo,h,g){t=(-vo-sqrt(vo^2+2*g*h))/g;return(t)}
#t2(vo = 50,h = 100,g = 9.81)
F3=function(x1,y1,x2,y2,x3,y3,x4,y4)
{{A=matrix(data = c(x1^3,x1^2,x1,1,x2^3,x2^2,x2,1,x3^3,x3^2,x3,1,x4^3,x4^2,x4,1),nrow = 4,ncol = 4,byrow = TRUE)}
{B=matrix(data = c(y1,y2,y3,y4),nrow = 4,ncol = 1)}
{S=solve(a = A,b = B)};
return(S)}
F3(x1 = -2,y1 = -20,x2 = 0,y2 = 4,x3 = 2,y3 = 68,x4 = 4,y4 = 508)
## [,1]
## [1,] 7
## [2,] 5
## [3,] -6
## [4,] 4
i=seq(0,2,0.1);i
## [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## [20] 1.9 2.0
Función=function(x)
{10*exp(-2*x)}
y=Función(i)
plot(i,y,main='Función 10*e^(-2x)',xlab='x',ylab='y',lwd=8,col='green')
x=seq(3,7,0.1);x
## [1] 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8
## [20] 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7
## [39] 6.8 6.9 7.0
FA=function(x)
{20*x^2-200*x+3}
y=FA(x)
plot(x,y,main='20*x^2-200*x+3',xlab='x',ylab='y',lwd=1,col='brown',ylim=c(-500,-450))
#axis(1,at=seq(0,10,0.1,lwd=1))
#axis(2,at=seq(-550,-400,10,lwd=1))
#DE acuerdo a este gráfico la ubicación aproximada del mínimo es en x=5, y=-500
x=data.frame(x,y)
min=min(x);min
## [1] -497
x=6
z=(x<10);z
## [1] TRUE
z=(x==10);z
## [1] FALSE
z=(x>=4);z
## [1] TRUE
z=(x!=7);z
## [1] TRUE
#En R no es igual es !=
x=matrix(data = c(-3,0,0,2,6,8),nrow = 1,ncol = 6,byrow = TRUE)
y=matrix(data = c(-5,-2,0,3,4,10),nrow = 1,ncol = 6,byrow = TRUE)
z=x>y;z
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] TRUE TRUE FALSE FALSE TRUE FALSE
Price=matrix(data = c(19,18,22,21,25,19,17,21,27,29),nrow = 1,ncol = 10,byrow = TRUE)
y=matrix(data = c(1,2,3,4,5,6,7,8,9,10),nrow = 1,ncol = 10,byrow = TRUE)
z=Price>20;z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE
sum(z)
## [1] 6
PriceA=matrix(data = c(19,18,22,21,25,19,17,21,27,29),nrow = 1,ncol = 10,byrow = TRUE)
PriceB=matrix(data = c(22,17,20,19,24,18,16,25,28,27),nrow = 1,ncol = 10,byrow = TRUE)
y=matrix(data = c(1,2,3,4,5,6,7,8,9,10),nrow = 1,ncol = 10,byrow = TRUE)
z=PriceA>PriceB;z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE
sum(z)
## [1] 7
PriceA=matrix(data = c(19,18,22,21,25,19,17,21,27,29),nrow = 1,ncol = 10,byrow = TRUE)
PriceB=matrix(data = c(22,17,20,19,24,18,16,25,28,27),nrow = 1,ncol = 10,byrow = TRUE)
PriceC=matrix(data = c(17,13,22,23,19,17,20,21,24,28),nrow = 1,ncol = 10,byrow = TRUE)
y=matrix(data = c(1,2,3,4,5,6,7,8,9,10),nrow = 1,ncol = 10,byrow = TRUE)
#45.a
z=PriceA>PriceB&PriceA>PriceC;z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE
sum(z)
## [1] 4
#45.b
z=PriceA>PriceB|PriceA>PriceC;z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
sum(z)
## [1] 9
#45.c
z=(PriceA>PriceB&PriceA<=PriceC)|(PriceA<=PriceB&PriceA>PriceC);z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE
sum(z)
## [1] 5
#45.a
a=6
b=6
c=2
z1=(a==b)&((b==c)|(a==c));z1
## [1] FALSE
z2=(a==b)|((b==c)&(a==c));z2
## [1] TRUE
#Por tanto no son equivalentes
#45.b
a=2
b=3
c=4
d=0
z1=(a<b)&((a>c)|(a>d));z1
## [1] TRUE
z2=((a<b)&(a>c))|((a<b)&(a>d));z2
## [1] TRUE
#Se evalua con 2,3,1,0 y 2,3,4,0, se encuentra que las dis expresiones son equivalentes
FX=function(x){y=ifelse(x<(-1),(exp(x+1)),(ifelse(x<5,(2+cos(pi*x)),no = 10*(x-5))));return(y)}
FX(x = -5)
## [1] 0.01831564
FX(x = 3)
## [1] 1
FX(x = 15)
## [1] 100
x=2
y=4
z=5
w=ifelse((x<y)&(z<10),x*y*z,'error');w
## [1] 40
k=seq(1,10,1)
for(i in k) {y=(5*k^3)}
sum(y)
## [1] 15125
#50.a
x = function(t) 5*t - 10
y = function(t) 25*t^2 - 120*t + 144
t = seq(0, 4, length.out = 1000)
#Inicializar las variables
dist_min = Inf
t_min = 0
for (i in 1:length(t)) {dist = sqrt(x(t[i])^2 + y(t[i])^2)
if (dist < dist_min) {dist_min <- dist
t_min = t[i]
}
}
cat("El objeto está más cerca del origen en el tiempo:", t_min, "\n")
## El objeto está más cerca del origen en el tiempo: 2.234234
cat("La distancia mínima al origen es:", dist_min, "\n")
## La distancia mínima al origen es: 1.357775
#50.b
x = function(t) 5*t - 10
y = function(t) 25*t^2 - 120*t + 144
t = seq(0, 4, length.out = 1000)
dist = sqrt(x(t)^2 + y(t)^2)
indice_min = which.min(dist)
tiempo_min = t[indice_min]
cat("El objeto está más cerca del origen en el tiempo:", tiempo_min, "\n")
## El objeto está más cerca del origen en el tiempo: 2.234234
cat("La distancia mínima al origen es:", dist[indice_min], "\n")
## La distancia mínima al origen es: 1.357775
A = matrix(data = c(3,-8,-17,5,-1,6,-4,33,-9), nrow = 3, ncol = 3)
#Se debe inicializar la matriz B
B = matrix(0, nrow = nrow(A), ncol = ncol(A))
# Calcular la matriz B
for (i in 1:nrow(A)) {
for (j in 1:ncol(A)) {
if (A[i,j] >= 1) {
B[i,j] = log(A[i,j])
} else {
B[i,j] = A[i,j] + 20
}
}
}
B
## [,1] [,2] [,3]
## [1,] 1.098612 1.609438 16.000000
## [2,] 12.000000 19.000000 3.496508
## [3,] 3.000000 1.791759 11.000000
# Inicializar las variables
k = 1
suma = 0
# Bucle while para encontrar cuántos términos se requieren para que la suma exceda 2000
while (suma <= 2000) {
suma = suma + 2^k
k = k + 1
}
# Restar uno al valor de k para obtener el número de términos requeridos
num_term = k - 1
# Calcular la suma de los términos requeridos
suma_term = sum(2^(1:num_term))
# Imprimir los resultados
cat("Se requieren", num_term, "términos para que la suma exceda 2000.\n")
## Se requieren 10 términos para que la suma exceda 2000.
cat("La suma de los", num_term, "términos requeridos es", suma_term, ".\n")
## La suma de los 10 términos requeridos es 2046 .
#Inicializar las variables
tB1=0
tB2=0
i_B1=5.5
i_B2=4.5
deposito_anual=1000
valor_anualB1=1000
valor_anualB2=1000
#Banco # 1
while (valor_anualB1<=50000) {
valor_anualB1=deposito_anual+valor_anualB1*(1+i_B1/100)
tB1=tB1+1
}
#Tiempo para llegar a 50000 en el banco 1
Tiempo_B1=tB1;Tiempo_B1
## [1] 24
#Banco # 2
while (valor_anualB2<=50000) {
valor_anualB2=deposito_anual+valor_anualB2*(1+i_B2/100)
tB2=tB2+1
}
#Tiempo para llegar a 50000 en el banco 2
Tiempo_B2=tB2;Tiempo_B2
## [1] 26
DifTiempo=Tiempo_B2-Tiempo_B1
# Se muestra el resultado
cat("Se requiere dejar el dinero", Tiempo_B1, "años en el banco 1 para alcanzar 50000.\n")
## Se requiere dejar el dinero 24 años en el banco 1 para alcanzar 50000.
cat("Se requiere dejar el dinero", Tiempo_B2, "años en el banco 2 para alcanzar 50000.\n")
## Se requiere dejar el dinero 26 años en el banco 2 para alcanzar 50000.
cat("Se requiere", DifTiempo, "años más en el banco 2 para alcanzar 50000 con respecto al banco 1.\n")
## Se requiere 2 años más en el banco 2 para alcanzar 50000 con respecto al banco 1.
#Se define la función V(t)
t=seq(1,100,by=1)
V= 10^9+10^8*(1-exp(-t/100))-(10^7)*t
plot(t,V,main="Volumen vs Tiempo",sub="Métodos Multivariados 2023-1S",xlab="t (días)", ylab="V(L)",col="blue",type="o")
abline(h = ((10^9)/2),col='red',lwd=3)
abline(v = 54.5,col='green',lwd=3)
#Aproximadamente 54.5
V_t=function(t){V=10^9+10^8*(1-exp(-t/100))-(10^7)*t;return(V)}
V_t(54.5)
## [1] 497015822
#Fue una buena aproximación
V_t(54.18)
## [1] 500029975
S = function(k) {sum((-1)^seq(0,k)/(2*seq(0,k) + 1))}
k = seq(0, 200, by = 1)
Dif = pi/4 - sapply(k,S);Dif
## [1] -0.214601837 0.118731497 -0.081268503 0.061588640 -0.049522472
## [6] 0.041386619 -0.035536458 0.031130209 -0.027693320 0.024938259
## [11] -0.022680789 0.020797472 -0.019202528 0.017834509 -0.016648250
## [16] 0.015609815 -0.014693215 0.013878213 -0.013148814 0.012492212
## [21] -0.011898032 0.011357782 -0.010864440 0.010412155 -0.009996008
## [26] 0.009611835 -0.009256089 0.008925729 -0.008618131 0.008331022
## [31] -0.008062421 0.007810595 -0.007574020 0.007351353 -0.007141401
## [36] 0.006943106 -0.006755524 0.006577809 -0.006409204 0.006249024
## [41] -0.006096655 0.005951538 -0.005813168 0.005681085 -0.005554870
## [46] 0.005434141 -0.005318547 0.005207768 -0.005101510 0.004999500
## [51] -0.004901490 0.004807248 -0.004716562 0.004629233 -0.004545079
## [56] 0.004463930 -0.004385628 0.004310025 -0.004236984 0.004166377
## [61] -0.004098085 0.004031996 -0.003968004 0.003906012 -0.003845926
## [66] 0.003787661 -0.003731136 0.003676272 -0.003622998 0.003571246
## [71] -0.003520952 0.003472055 -0.003424497 0.003378224 -0.003333185
## [76] 0.003289331 -0.003246616 0.003204997 -0.003164430 0.003124878
## [81] -0.003086302 0.003048667 -0.003011939 0.002976085 -0.002941075
## [86] 0.002906878 -0.002873468 0.002840817 -0.002808900 0.002777692
## [91] -0.002747170 0.002717311 -0.002688094 0.002659499 -0.002631506
## [96] 0.002604096 -0.002577251 0.002550954 -0.002525188 0.002499938
## [101] -0.002475187 0.002450922 -0.002427127 0.002403791 -0.002380898
## [106] 0.002358438 -0.002336398 0.002314765 -0.002293530 0.002272680
## [111] -0.002252207 0.002232098 -0.002212346 0.002192940 -0.002173872
## [116] 0.002155132 -0.002136713 0.002118606 -0.002100803 0.002083297
## [121] -0.002066080 0.002049146 -0.002032487 0.002016096 -0.001999968
## [126] 0.001984096 -0.001968473 0.001953095 -0.001937955 0.001923048
## [131] -0.001908369 0.001893912 -0.001879673 0.001865646 -0.001851826
## [136] 0.001838210 -0.001824793 0.001811570 -0.001798538 0.001785692
## [141] -0.001773027 0.001760542 -0.001748230 0.001736090 -0.001724117
## [146] 0.001712309 -0.001700661 0.001689170 -0.001677833 0.001666648
## [151] -0.001655611 0.001644719 -0.001633969 0.001623360 -0.001612886
## [156] 0.001602548 -0.001592341 0.001582263 -0.001572311 0.001562485
## [161] -0.001552780 0.001543195 -0.001533728 0.001524376 -0.001515138
## [166] 0.001506010 -0.001496993 0.001488082 -0.001479277 0.001470576
## [171] -0.001461976 0.001453476 -0.001445075 0.001436770 -0.001428560
## [176] 0.001420443 -0.001412418 0.001404483 -0.001396637 0.001388878
## [181] -0.001381205 0.001373616 -0.001366110 0.001358686 -0.001351341
## [186] 0.001344076 -0.001336889 0.001329778 -0.001322742 0.001315780
## [191] -0.001308892 0.001302075 -0.001295328 0.001288651 -0.001282043
## [196] 0.001275502 -0.001269027 0.001262618 -0.001256273 0.001249992
## [201] -0.001243773
plot(k, Dif, type = "l", xlab = "k", ylab = "pi/4 - S(k)", main = "Diferencia entre pi/4 y la suma S(k) para k entre 0 y 200",sub="Métodos Multivariados 2023-1S")
#Aquí se utiliza la función sumatoria y la función sapply para aplicar la función S a cada valor del vector k
A=matrix(data = c(0,5,10,15,20,-8,-4,-1,1,2,6,3,1,0,-1),nrow = 5,ncol = 3);A
## [,1] [,2] [,3]
## [1,] 0 -8 6
## [2,] 5 -4 3
## [3,] 10 -1 1
## [4,] 15 1 0
## [5,] 20 2 -1
Tiempo=A[,1]
Fuerza1=A[,2]
Fuerza2=A[,3]
plot(Tiempo,Fuerza1, type = "o", xlab = "t(s)", ylab = "F(N)", ylim=c(-10,7), xlim = c(-1,21),col='red', main = "Fuerza vs tiempo",sub="Métodos Multivariados 2023-1S")
lines(Tiempo,Fuerza2,col='blue',type = "o")
legend("topright", legend=c("Fuerza1","Fuerza2"), fill=c("red", "blue"))
#abline(v = 0)
x=seq(0,1,length=100)
Sin_x=sin(x)
X1=x
#Aquí se grafuca Sin x vs x
plot(x,Sin_x, type = "l", xlab = "x", ylab = "f(x)", ylim=c(0,1),main = "Análisis apróximación del pequeño ángulo",sub="Métodos Multivariados 2023-1S",col='red')
lines(x,X1,type = "l",col='blue')
legend("topleft", legend=c("Sin x","x"), fill=c("red", "blue"))
#Aquí se grafica ek error de aproximación sin x - x vs x
Err_Ap=sin(x)-x
plot(x,Err_Ap, type = "l", xlab = "x", ylab = "Sin(x)-x", ylim=c(-0.17,0),main = "Error de la apróximación del ángulo pequeño",sub="Métodos Multivariados 2023-1S")
#Ahora se grafica el error relativo
Err_R=(sin(x)-x)/sin(x)
plot(x,Err_R, type = "l", xlab = "x", ylab = "(Sin (x) - x)/sin(x)", ylim=c(-0.2,0),main = "Error relativo de la aproximación del ángulo pequeño",sub="Métodos Multivariados 2023-1S")
#Ahora se calcula el valor para el que el error relativo es 5%
#Aquí buscamos ahora el valor para el que el error relativo es 5% y usamos la función curve para graficar funciones
#SE crea la función error relativo con valor absoluto
Error_Rel=function(x){abs((sin(x)-x)/sin(x))}
#Se grafica con la función curve
curve(Error_Rel,from = 0,to = 1,xlab = "x",ylab = "(Sin (x) - x)/sin(x)",main = "Error relativo de la aproximación del ángulo pequeño")
#Se inicializa la variable x en 0.0001 (Para evitar un error al dividir entre 0)
x=0.0001
#Con la herramienta "while" se empieza a iterar hasta que el valor iguale o supere el 5%
while (Error_Rel(x)<=0.05) {
x=x+0.00001
}
cat("El valor de x debe ser menor o igual a", x, "para que el error relativo sea menor o igual a 5%.")
## El valor de x debe ser menor o igual a 0.53842 para que el error relativo sea menor o igual a 5%.
x=seq(0,2*pi,length=1000)
F1=tan(2*x)
F2=2*tan(x)/(1-(tan(x))^2)
plot(x = x,y = F1,type="l",xlab="x",ylab = "f(x)",main="f(x)=tan(2*x)",sub="Métodos Multivariados 2023-1S",col='red')
plot(x = x,y = F2,type="l",xlab="x",ylab = "f(x)",main="2*tan(x)/(1-(tan(x))^2)",sub="Métodos Multivariados 2023-1S",col='blue')
#Se observa que las gráficas son equivalentes, se colcoan en gráficos diferentes porque se superponen
t=seq(0,1,length=100)
F_tb1=1-exp(-1*t)
F_tb5=1-exp(-5*t)
F_tb10=1-exp(-10*t)
plot(x = t,y = F_tb1,ylim=c(0,1),col='red',type='l',main="y(t)=1-exp(-b*t) para diferentes valores de b",sub="Métodos Multivariados 2023-1S",lwd=2,ylab = "y(t)",xlab="t")
lines(t,F_tb5,type='l',col='blue',lwd=2)
lines(t,F_tb10,type='l',col='green',lwd=2)
legend("bottomright", legend=c("b=1","b=5","b=10"), fill=c("red", "blue","green"))
#Para encontrar el valor de t en el que la función alcanza el 98% de su valor estacionario se crean como funciones.
t=0
F_tb1=function(t){1-exp(-1*t)}
while (F_tb1(t)<=0.98) {
t=t+0.0001
}
cat("el valor de t para el que la f(t) alcanza el 98% de su valor estático es",t,"cuando b=1 1.\n")
## el valor de t para el que la f(t) alcanza el 98% de su valor estático es 3.9121 cuando b=1 1.
t=0
F_tb5=function(t){1-exp(-5*t)}
while (F_tb5(t)<=0.98) {
t=t+0.0001
}
cat("el valor de t para el que la f(t) alcanza el 98% de su valor estático es",t,"cuando b=5 1.\n")
## el valor de t para el que la f(t) alcanza el 98% de su valor estático es 0.7825 cuando b=5 1.
t=0
F_tb10=function(t){1-exp(-10*t)}
while (F_tb10(t)<=0.98) {
t=t+0.0001
}
cat("el valor de t para el que la f(t) alcanza el 98% de su valor estático es",t,"cuando b=10 1.\n")
## el valor de t para el que la f(t) alcanza el 98% de su valor estático es 0.3913 cuando b=10 1.
#La ecuación de gases ideales despejando p
V=seq(20,100,length=100)
R=286.7
T=293
p1kg=(1*R*T)/V
p3kg=(3*R*T)/V
p7kg=(7*R*T)/V
plot(x = V,y = p1kg,ylim=c(0,30000),type = "l",main="Presión de aire vs Volumen",sub="Métodos Multivariados 2023-1S",xlab="V(m^3)",ylab = "P(N/(m^2))",col='red',lwd=2)
lines(V,p3kg,type='l',col='blue',lwd=2)
lines(V,p7kg,type='l',col='green',lwd=2)
legend("topright",legend = c('m=1kg','m=3kg','m=7kg'),fill=c('red','blue','green'))
#Primero creamos un data frame para los datos
tiempo=c(1,2,3,4,5,6,7,8,10)
velocidad=c(1210,1866,2301,2564,2724,2881,2879,2915,3010)
df1=data.frame(tiempo,velocidad);df1
## tiempo velocidad
## 1 1 1210
## 2 2 1866
## 3 3 2301
## 4 4 2564
## 5 5 2724
## 6 6 2881
## 7 7 2879
## 8 8 2915
## 9 10 3010
colnames(df1)=c("tiempo (s)","velocidad (rpm)")
#Para determinar si la función describe los datos de velocidad en función del tiempo, primero realizamos una gráfica para ver como se comportan los datos
#Se hace una exploración de los datos para t=1sec colocando arbitrariamente c=-0.5 y se obtiene un b= de 3075.2 con estos datos se hace una gráfica para ver el comportamiento de la función y evaluar si se comportan de forma similar.
f1=function(tiempo){3075*(1-exp((-0.5)*tiempo))}
curve(f1,from = 1,to = 10,xlab='tiempo(s)',ylab='velocidad(rpm)',main='Velocidad rotacional vs tiempo en un motor',col='red')
lines(tiempo,velocidad,col='blue')
legend("bottomright", legend=c("s(t)=3075*(1-exp((-0.5)*t ","Datos experimentales"), fill=c("red", "blue"))
#Se ve aquí que las dos curvas se comportan de forma similar y la función describe los datos, ahora buscaremos los datos de los coeficientes b y c
#Se observa un patrón no lineal, ahora para aplicar una regresión no lineal usamos la función nls de RStudio
modelo=nls(velocidad ~ b*(1-exp(c*tiempo)), data = df1, start = list(b=3075.2, c=-0.5))
## Warning in min(x): ningún argumento finito para min; retornando Inf
## Warning in max(x): ningun argumento finito para max; retornando -Inf
summary(modelo)
##
## Formula: velocidad ~ b * (1 - exp(c * tiempo))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b 2996.28857 19.14635 156.49 1.15e-13 ***
## c -0.49291 0.01134 -43.46 8.91e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.9 on 7 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 3.773e-06
#El modelo requirió de 3 iteraciones Ahora extraemos los coeficientes del modelo.
b=coef(modelo)["b"]
c=coef(modelo)["c"]
cat("b =", b, "\n")
## b = 2996.289
cat("c =", c, "\n")
## c = -0.4929142
#Ahora usamos estos coeficientes obtenidos del modelo para comparara las gráficas
f1=function(tiempo){2996.289 *(1-exp((-0.4929142 )*tiempo))}
curve(f1,from = 1,to = 10,xlab='tiempo(s)',ylab='velocidad(rpm)',main='Velocidad rotacional vs tiempo en un motor',col='red')
lines(tiempo,velocidad,col='blue')
legend("bottomright", legend=c("s(t)=3075*(1-exp((-0.5)*t ","Datos experimentales"), fill=c("red", "blue"))
#Por lo cual se observa que estos valores en general muestran un buen comportamiento de la función, pueden variar de acuerdo a los valores iniciales de partida,
library(ggplot2)
#install.packages("plotly")
library(plotly)
## Warning: package 'plotly' was built under R version 4.2.3
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
df<- data.frame(t=c(1990:1994),
T=c(18,19,21,17,20))
#Stairs
st<- ggplot() +
geom_step(data=df, mapping=aes(x=t, y=T)) +
geom_step(data=df, mapping=aes(x=t, y=T), direction="vh", linetype=3) +
geom_point(data=df, mapping=aes(x=t, y=T), color="red")+
xlab("Año")+
ylab("Temperatura (°C)")+
ggtitle("Temperatura promedio en una ciudad Stairs Plot")
ggplotly(st)
#stem
s<- ggplot(data=df, aes(x=t, y=T)) +
geom_hline(aes(yintercept=0)) +
geom_point(aes(t,T),size=3) +
geom_segment(aes(t,T,xend=t,yend=T-T))+
xlab("Año")+
ylab("Temperatura (°C)")+
ggtitle("Temperatura promedio en una ciudad Stem Plot")+
geom_text(aes(label=T),nudge_y=+1)
ggplotly(s)
#bar
b<-ggplot(data=df, aes(x=t, y=T))+geom_bar(stat="identity")+
xlab("Año")+
ylab("Temperatura (°C)")+
ggtitle("Temperatura promedio en una ciudad Bar Plot")+
geom_text(aes(label=T))
ggplotly(b)
#63.a
#Para hacer subplots usamos la librería plotly
library(plotly)
#Incluimos aquí las funciones de volumen y área
V=function(r) {((4/3)*pi*r^3)}
A=function(r) {(4*pi*r^2)}
r=seq(0.1, 100, by = 0.1)
# Graficar la función del volumen de una esfera
plot_ly(x = r, y = V(r), type = "scatter", mode = "lines", name = "V") %>%
layout(title = "Volumen de una Esfera",
xaxis = list(title = "Radio r(m)"), yaxis = list(title = "Volumen (m^3)"))
# Graficar la función del área de una esfera
plot_ly(x = r, y = A(r), type = "scatter", mode = "lines", name = "A") %>%
layout(title = "Área de una Esfera",
xaxis = list(title = "Radio r(m)"), yaxis = list(title = "Área (m^2)"))
#Ahora para poder permitir que estas funciones resulten en líneas rectas debemos transformar logarítmicamente los datos, para esto
log_V=log10(V(r))
log_A=log10(A(r))
# Ahora graficamos el logaritmo del volumen de la esfera y le damos un nombre
grV=plot_ly(x = log10(r), y = log_V, type = "scatter", mode = "lines", name = "V") %>%
layout(title = "Volumen de una Esfera escala logarítmica",
xaxis = list(title = "log10 (Radio r(m))"), yaxis = list(title = "log10(Volumen (m^3))"))
# Ahora hacemos lo mismo con el área de la esfera
grA=plot_ly(x = log10(r), y = log_A, type = "scatter", mode = "lines", name = "A") %>%
layout(title = "Área de una Esfera escala logarítmica",
xaxis = list(title = "log10 (Radio r(m))"), yaxis = list(title = "log10(Área (m^2))"))
#Ahora hacemos los subplots
fig=subplot(grV,grA,nrows = 2) %>%
layout(title="Subplots V y A")
fig
#63.b
#Definimos ahora las funciones de V y r vs A
r=function(A) {(A/(4*pi))^(0.5)}
V=function(A) {4/3*pi*((A/(4*pi))^(0.5))^3}
A=seq(1, 10^4, by = 1)
# Graficar la función del volumen de una esfera en función de su área
plot_ly(x = A, y = V(A), type = "scatter", mode = "lines", name = "V") %>%
layout(title = "Volumen de una Esfera en función del área",
xaxis = list(title = "área(m^2)"), yaxis = list(title = "Volumen (m^3)"))
# Graficamos la función del radio de una esfera en función de su área
plot_ly(x = A, y = r(A), type = "scatter", mode = "lines", name = "r") %>%
layout(title = "Radio de una Esfera en función del área",
xaxis = list(title = "área(m^2)"), yaxis = list(title = "Radio (m)"))
#Ahora para poder permitir que estas funciones resulten en líneas rectas debemos transformar logarítmicamente los datos, para esto
log_r=log10(r(A))
log_V=log10(V(A))
# Ahora graficamos el logaritmo del volumen de la esfera
grV=plot_ly(x = log10(A), y = log_V, type = "scatter", mode = "lines", name = "V") %>%
layout(title = "log Volumen de una Esfera en función de logaritmo Área",
xaxis = list(title = "log10 (Área)"), yaxis = list(title = "log10(Volumen)"))
# Ahora graficamos el logaritmo del radio de la esfera
grr=plot_ly(x = log10(A), y = log_r, type = "scatter", mode = "lines", name = "r") %>%
layout(title = "log Radio de una Esfera en función de logaritmo Área",
xaxis = list(title = "log10 (Área)"), yaxis = list(title = "log10(Radio)"))
#Ahora hacemos los subplots
fig=subplot(grV,grr,nrows = 2) %>%
layout(title="Subplots V y r")
fig
#Ejercicio 64
#EJERCICIO 64
#64A Sofía
library(ggplot2)
x<-c(25,30,35,40,45)
y<-c(5,260,480,745,1100)
df <- data.frame(x,y)
lineal.model <-lm(y~x,df)
log.model <-lm(log(y) ~ x, df)
exp.model <-lm(y ~ exp(x), df)
pot.model <-lm (log (y) ~ log (x))
log.model.df <- data.frame(x = df$x,
y = exp(fitted(log.model)))
pot.model.df <- data.frame(x = df$x,
y = exp(fitted(pot.model)))
ggplot(df, aes(x=x, y=y)) +
geom_point() +
geom_smooth(method="lm", aes(color="lineal"), formula= (y ~ x), se=FALSE, linetype = 1) + geom_smooth(method="lm", aes(color="Exp Model"), formula= (y ~ exp(x)), se=FALSE, linetype = 1)+
geom_line(data = log.model.df, aes(x, y, color = "Log Model"), size = 1, linetype = 2) +
geom_line(data = pot.model.df, aes(x, y, color = "Pot Model"), size = 1, linetype = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
guides(color = guide_legend("Model Type"))
## $colour
## $title
## [1] "Model Type"
##
## $title.position
## NULL
##
## $title.theme
## NULL
##
## $title.hjust
## NULL
##
## $title.vjust
## NULL
##
## $label
## [1] TRUE
##
## $label.position
## NULL
##
## $label.theme
## NULL
##
## $label.hjust
## NULL
##
## $label.vjust
## NULL
##
## $keywidth
## NULL
##
## $keyheight
## NULL
##
## $direction
## NULL
##
## $override.aes
## named list()
##
## $nrow
## NULL
##
## $ncol
## NULL
##
## $byrow
## [1] FALSE
##
## $reverse
## [1] FALSE
##
## $order
## [1] 0
##
## $available_aes
## [1] "any"
##
## $name
## [1] "legend"
##
## attr(,"class")
## [1] "guide" "legend"
##
## attr(,"class")
## [1] "guides"
#El modelo que mejor se ajusta en este caso es el lineal.
#64B
library(ggplot2)
xb<-c(2.5,3,3.5,4,4.5,5,5.5,6,7,8,9,10)
yb<-c(1500,1220,1050,915,810,745,690,620,520,480,410,390)
dfb <- data.frame(xb,yb)
lineal.modelb <-lm(yb~xb,dfb)
log.modelb <-lm(log(yb) ~ xb, dfb)
exp.modelb <-lm(yb ~ exp(xb), dfb)
pot.modelb <-lm (log (yb) ~ log (xb))
log.model.dfb <- data.frame(x = dfb$xb,
y = exp(fitted(log.modelb)))
pot.model.dfb <- data.frame(x = dfb$xb,
y = exp(fitted(pot.modelb)))
ggplot(dfb, aes(x=xb, y=yb)) +
geom_point() +
geom_smooth(method="lm", aes(color="lineal"), formula= (y~x),
se=FALSE, linetype = 1) +
geom_smooth(method="lm", aes(color="Exp Model"), formula= (y ~ exp(x)),
se=FALSE, linetype = 1)+
geom_line(data = log.model.dfb, aes(x, y, color = "Log Model"),
size = 1, linetype = 2) +
geom_line(data = pot.model.dfb, aes(x, y, color = "Pot Model"),
size = 1, linetype = 2)
guides(color = guide_legend("Model Type"))
## $colour
## $title
## [1] "Model Type"
##
## $title.position
## NULL
##
## $title.theme
## NULL
##
## $title.hjust
## NULL
##
## $title.vjust
## NULL
##
## $label
## [1] TRUE
##
## $label.position
## NULL
##
## $label.theme
## NULL
##
## $label.hjust
## NULL
##
## $label.vjust
## NULL
##
## $keywidth
## NULL
##
## $keyheight
## NULL
##
## $direction
## NULL
##
## $override.aes
## named list()
##
## $nrow
## NULL
##
## $ncol
## NULL
##
## $byrow
## [1] FALSE
##
## $reverse
## [1] FALSE
##
## $order
## [1] 0
##
## $available_aes
## [1] "any"
##
## $name
## [1] "legend"
##
## attr(,"class")
## [1] "guide" "legend"
##
## attr(,"class")
## [1] "guides"
#El modelo potencial en este caso es el que mejor se ajusta en este caso
#64C
library(ggplot2)
xc<-c(550,600,650,700,750)
yc<-c(41.2,18.62,8.62,3.92,1.86)
dfc <- data.frame(xc,yc)
lineal.modelc <-lm(yc~xc,dfc)
log.modelc <-lm(log(yc) ~ xc, dfc)
exp.modelc <- nls(yc ~ a*exp(b*xc), start = list(a = 5, b = 0.01), dfc)
coef_expo=coef(exp.modelc);coef_expo
## a b
## 2.354840e+05 -1.573054e-02
pot.modelc <-lm (log (yc) ~ log (xc))
log.model.dfc <- data.frame(x = dfc$xc,
y = exp(fitted(log.modelc)))
pot.model.dfc <- data.frame(x = dfc$xc,
y = exp(fitted(pot.modelc)))
#Aquí se usa la función predict para generar valores de y en el modelo exponencial a lo largo del rango del eje x, aunque explorando los datos en excel se observa que el modelo exponencial no se ajusta, el modelo logarítmico parece tener un mejor comportamiento de los datos
exp.model.dfc <- data.frame(x = dfc$xc,
y = predict(exp.modelc, newdata = data.frame(xc = xc)))
ggplot(dfc, aes(x=xc, y=yc)) +
geom_point() +
geom_smooth(method="lm", aes(color="lineal"), formula= (y~x),
se=FALSE, linetype = 1) +
geom_line(data = exp.model.dfc, aes(x,y, color = "Exp Model"),
size = 1, linetype = 2)+
geom_line(data = log.model.dfc, aes(x, y, color = "Log Model"),
size = 1, linetype = 2) +
geom_line(data = pot.model.dfc, aes(x, y, color = "Pot Model"),
size = 1, linetype = 2)+
scale_y_continuous(limits = c(0, max(yc) * 1.1))
## Warning: Removed 9 rows containing missing values (`geom_smooth()`).
guides(color = guide_legend("Model Type"))
## $colour
## $title
## [1] "Model Type"
##
## $title.position
## NULL
##
## $title.theme
## NULL
##
## $title.hjust
## NULL
##
## $title.vjust
## NULL
##
## $label
## [1] TRUE
##
## $label.position
## NULL
##
## $label.theme
## NULL
##
## $label.hjust
## NULL
##
## $label.vjust
## NULL
##
## $keywidth
## NULL
##
## $keyheight
## NULL
##
## $direction
## NULL
##
## $override.aes
## named list()
##
## $nrow
## NULL
##
## $ncol
## NULL
##
## $byrow
## [1] FALSE
##
## $reverse
## [1] FALSE
##
## $order
## [1] 0
##
## $available_aes
## [1] "any"
##
## $name
## [1] "legend"
##
## attr(,"class")
## [1] "guide" "legend"
##
## attr(,"class")
## [1] "guides"
#El modelo logarítmico en este caso es el que mejor se ajusta. El modelo exponencial genera datos muy por fuera de la gráfica.
Year <- c(2002,2003,2004,2005,2006,2007)
Pop <-c(10,10.8,11.7,12.7,13.8,14.9)
df65 <-data.frame(Year,Pop)
plot(Year,Pop)
cor.test(Year, Pop)
##
## Pearson's product-moment correlation
##
## data: Year and Pop
## t = 32.407, df = 4, p-value = 5.406e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9818926 0.9998023
## sample estimates:
## cor
## 0.9981011
model1 <- lm(Year~Pop)
summary(model1)
##
## Call:
## lm(formula = Year ~ Pop)
##
## Residuals:
## 1 2 3 4 5 6
## -0.1586759 0.0328093 0.1232302 0.1125867 0.0008788 -0.1108290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.992e+03 3.877e-01 5138.25 8.61e-15 ***
## Pop 1.011e+00 3.119e-02 32.41 5.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1288 on 4 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9953
## F-statistic: 1050 on 1 and 4 DF, p-value: 5.406e-06
model1$coefficients
## (Intercept) Pop
## 1992.052241 1.010643
y <- function(x) {
(x-1992.05)/(1.016)
}
curve(y, from = min(Year), to = max(Year),
type="l",
col="red",
lwd=3,
add=TRUE,
main="Crecimiento poblacional",
xlab="Population (Millions)",
ylab="Year"
)
Población=1992.05+1.0106*20
Población
## [1] 2012.262
#Según este modelo la población se duplicará en algún momento del año 2012
#66.a
#Para encontrar el coeficiente b, es necesario despejar b de la ecuación cuando la función es 0.5 y el tiempo es 5500 años
\[x=e^{-bt}\] \[ln(x)=-bt\] \[b=-\frac{ln(x))}{t}\]
#Por tanto, si definimos t=5500 y x como 0.5
t=5500
x=0.5
b=-log(x)/t;b
## [1] 0.0001260268
cat('El valor del coeficiente b es',b,'\n')
## El valor del coeficiente b es 0.0001260268
#Ahora dibujamos la función
FC14=function(t){exp(-0.0001260268 *t)}
curve(FC14,from = 0,to = 11000,type="l",col="blue",lwd=3,main="Fracción Carbono 14",xlab="Tiempo (Años)", ylab="Fracción de carbono 14")
#66.b
#Si el 90 porciento del carbono 14 permanece entonces, usamos la función uniroot
C90 = uniroot(function(t) FC14(t) - 0.9, interval = c(0, 10000))
C90
## $root
## [1] 836.0167
##
## $f.root
## [1] -6.850076e-14
##
## $iter
## [1] 6
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-05
#De acuerdo con esto han transcurrido 836 años desde que el organismo falleció
#66.c
#Si el estimado de b varía en +/-1% entonces los valores de be pueden estar entre
b=0.0001260268
bmin=b-b*0.01;bmin
## [1] 0.0001247665
bmax=b+b*0.01;bmax
## [1] 0.0001272871
#Para estos valores vamos a graficar las funciones para determinar visualmente
FC14min=function(t){exp(-0.0001247665 *t)}
FC14max=function(t){exp(-0.0001272871 *t)}
curve(FC14min,from = 0,to = 11000,type="l",col="blue",lwd=1,main="Fracción Carbono 14",xlab="Tiempo (Años)", ylab="Fracción de carbono 14")
curve(FC14max,from = 0,to = 11000,type="l",col="red",lwd=1,add = TRUE)
#SE observa graficamente que las mayores diferencias ocurren en los tiempos mayores, así que estimamos la diferencia de las funciones en el tiempo 11000
FC14min(t = 11000);FC14min
## [1] 0.2534898
## function(t){exp(-0.0001247665 *t)}
## <bytecode: 0x000001fb9e1c7628>
FC14max(t = 11000);FC14max
## [1] 0.246558
## function(t){exp(-0.0001272871 *t)}
## <bytecode: 0x000001fb9ae1d930>
#Aquí vemos que para una fracción muy baja (25% el tiempo son 11000 años), calcularemos ahora las diferencias de edad para una fracción de 25% para estimar el error, la ecuación sería la siguiente
\[t=-\frac{ln(x))}{b}\]
#Por tanto,
tmax=-log(0.25)/0.0001247665;tmax
## [1] 11111.11
tmin=-log(0.25)/0.0001272871;tmin
## [1] 10891.08
Dift=tmax-tmin
cat('El error en estimación para una fracción de C14 de 0.25, 11000 años de edad aprox desde que fallece un organismo es de',Dift,'años\n')
## El error en estimación para una fracción de C14 de 0.25, 11000 años de edad aprox desde que fallece un organismo es de 220.0275 años
#Vamos inicialmente a graficar los datos para hacer una explotración
Time=seq(0,6,by=1);Time
## [1] 0 1 2 3 4 5 6
Temperature=c(300,150,75,35,12,5,2)
plot(Time,Temperature,xlab = 'Tiempo(s)',ylab='Temperature(°C)',main = 'Quenching temperature of a cupper sphere')
#Se observa que tiene una correlación no lineal por lo que se va a buscar la mejor correlación posible, la cual podría ser una tendencia exponencial para esto usamos la función lm para construir el modelo
logT=log(Temperature)
Modelo=lm(logT~Time)
summary(Modelo)
##
## Call:
## lm(formula = logT ~ Time)
##
## Residuals:
## 1 2 3 4 5 6 7
## -0.171204 -0.019106 0.132993 0.216099 -0.009097 -0.039320 -0.110365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.87499 0.10008 58.70 2.71e-08 ***
## Time -0.84525 0.02776 -30.45 7.17e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1469 on 5 degrees of freedom
## Multiple R-squared: 0.9946, Adjusted R-squared: 0.9936
## F-statistic: 927.2 on 1 and 5 DF, p-value: 7.167e-07
#Ahora vamos a graficar la función en la misma gráfica, aquí usamos la función predict para encontrar las predicciones del modelo ajustado para los valores de x, y se eleva a un exponente ya que la función lm obtiene los datos de una regresión lineal y se deben pasar a una función exponencial.
lines(Time,exp(predict(Modelo)),col='red')
#68.a
A=seq(0,9,by=1)
t=c(130,115,110,90,89,89,95,100,110,125)
#Para encontrar los polinomios que se ajustan a estos datos usamos la función lm y con la función poly creamos una variable que representa x (en este caso la cantidad de aditivo) elevada a diferentes potencias
#Polinomio de primer grado (línea recta)
Modelo1=lm(t~A)
#Polinomio de segundo grado, creamos la variable A2, con la aclaración de raw=True
A2=poly(x = A,degree = 2,raw = TRUE)
Modelo2=lm(t~A2)
#Polinomio de tercer grado, creamos la variable A3
A3=poly(x = A,degree = 3,raw = TRUE)
Modelo3=lm(t~A3)
#Polinomio de CUARTO grado, creamos la variable A4
A4=poly(x = A,degree = 4,raw = TRUE)
Modelo4=lm(t~A4)
#Ahora graficamos los datos y los polinomios con la función plot, y utilizamos la función "lines" para graficar los otros polinomios
plot(A,t,xlab = "A(oz)",ylab = 't(Min)',main='Tiempos de secado de una pintura en función de la cantidad de aditivo A, ajustes con polinomios de de grados 1,2,3 y 4',cex.main=0.7,sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='blue',pch=16)
lines(A,predict(Modelo1),col='red')
lines(A,predict(Modelo2),col='orange')
lines(A,predict(Modelo3),col='green')
lines(A,predict(Modelo4),col='purple')
legend("top", legend=c("Datos","P1","P2","P3","P4"), fill=c("blue","red","orange","green","purple"))
#Para obtener los parámetros J, S y r^2 utilizamosla función "summary" sobre cada modelo
summary(Modelo1)
##
## Call:
## lm(formula = t ~ A)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.65 -14.28 -0.30 7.23 22.84
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.436 9.284 11.680 2.63e-06 ***
## A -0.697 1.739 -0.401 0.699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.8 on 8 degrees of freedom
## Multiple R-squared: 0.01968, Adjusted R-squared: -0.1029
## F-statistic: 0.1606 on 1 and 8 DF, p-value: 0.6991
summary(Modelo2)
##
## Call:
## lm(formula = t ~ A2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9136 -1.1523 -0.4212 0.1830 6.7682
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 131.3000 2.6462 49.62 3.54e-10 ***
## A21 -17.8447 1.3693 -13.03 3.65e-06 ***
## A22 1.9053 0.1465 13.01 3.69e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.366 on 7 degrees of freedom
## Multiple R-squared: 0.9611, Adjusted R-squared: 0.9499
## F-statistic: 86.38 on 2 and 7 DF, p-value: 1.165e-05
summary(Modelo3)
##
## Call:
## lm(formula = t ~ A3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1124 -0.9311 -0.4212 0.2839 6.5438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 131.03077 3.29207 39.802 1.68e-08 ***
## A31 -17.35218 3.35107 -5.178 0.00206 **
## A32 1.76107 0.89508 1.967 0.09669 .
## A33 0.01068 0.06526 0.164 0.87534
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.627 on 6 degrees of freedom
## Multiple R-squared: 0.9612, Adjusted R-squared: 0.9418
## F-statistic: 49.59 on 3 and 6 DF, p-value: 0.0001256
summary(Modelo4)
##
## Call:
## lm(formula = t ~ A4)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 0.04545 -1.76573 5.52739 -4.93298 -0.40676 0.44406 2.55536 -0.53904
## 9 10
## -1.70746 0.77972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 129.95455 3.58854 36.214 3.02e-07 ***
## A41 -12.86791 6.22566 -2.067 0.0936 .
## A42 -0.75510 3.05724 -0.247 0.8147
## A43 0.45911 0.52415 0.876 0.4212
## A44 -0.02491 0.02888 -0.863 0.4278
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.707 on 5 degrees of freedom
## Multiple R-squared: 0.9663, Adjusted R-squared: 0.9393
## F-statistic: 35.79 on 4 and 5 DF, p-value: 0.0007146
#J es el error cuadrático medio
#S es la desviación estándar residual
#r^2 es el coeficiente de determinación
#Extraemos los parámetros de cada modelo
#Modelo 1
R2_1=summary(Modelo1)$r.squared;R2_1
## [1] 0.01968261
J_1=summary(Modelo1)$sigma^2;J_1
## [1] 249.503
S_1=sqrt(J_1);S_1
## [1] 15.79566
#Modelo 2
R2_2=summary(Modelo2)$r.squared;R2_2
## [1] 0.9610582
J_2=summary(Modelo2)$sigma^2;J_2
## [1] 11.32706
S_2=sqrt(J_2);S_2
## [1] 3.365569
#Modelo 3
R2_3=summary(Modelo3)$r.squared;R2_3
## [1] 0.9612314
J_3=summary(Modelo3)$sigma^2;J_3
## [1] 13.15614
S_3=sqrt(J_3);S_3
## [1] 3.627139
#Modelo 4
R2_4=summary(Modelo4)$r.squared;R2_4
## [1] 0.9662528
J_4=summary(Modelo4)$sigma^2;J_4
## [1] 13.74254
S_4=sqrt(J_4);S_4
## [1] 3.707093
#El modelo 4 presenta un mejor ajuste de los datos ya que tiene el menor coeficiente de determinación.
#68.b
#La ecuación del modelo 4, que es la que mejor se austó a los datos es (extrayendo los datos del modelo)
tc=129.95455-12.86791*A-0.75510*A^2+0.45911*A^3-0.02491*A^4
A=seq(0,9,by=1)
to=c(130,115,110,90,89,89,95,100,110,125)
tc
## [1] 129.95455 116.76574 104.47265 94.93318 89.40739 88.55750 92.44789
## [8] 100.54510 111.71783 124.23694
Datos=data.frame(A,tc,to);Datos
## A tc to
## 1 0 129.95455 130
## 2 1 116.76574 115
## 3 2 104.47265 110
## 4 3 94.93318 90
## 5 4 89.40739 89
## 6 5 88.55750 89
## 7 6 92.44789 95
## 8 7 100.54510 100
## 9 8 111.71783 110
## 10 9 124.23694 125
#Aquí se muestran los tc (tiempos calculados) contra los to (tiempos observados) en un dataframe
#Para poder encontrar el mínimo de t podemos derivar la función y encontrar el valor donde la derivada es igual a cero (tangente a la curva)
#Al derivar la función obtenemos
dt=-12.86791-1.5102*A+1.37733*A^2-0.09964*A^3
#Para encontrar el punto enel cual la derivada es cero podemos utilizar la función uniroot aplicada a una función y buscamos la solución en el intervalo desde 0 a 9 (rango de A)
fdt=function(A){-12.86791-1.5102*A+1.37733*A^2-0.09964*A^3}
AMínimo=uniroot(fdt,c(0,9))
AMínimo$root
## [1] 4.676351
ft=function(A){129.95455-12.86791*A-0.75510*A^2+0.45911*A^3-0.02491*A^4}
tMínimo=ft(AMínimo$root)
tMínimo
## [1] 88.30475
cat('Según esto, la cantidad de aditivo óptima es ',round(AMínimo$root,1),'oz, la cual obtendrá un tiempo de secado de ',round(tMínimo,1),'min \n')
## Según esto, la cantidad de aditivo óptima es 4.7 oz, la cual obtendrá un tiempo de secado de 88.3 min
#Los datos no aparecen en la tabla, por lo que se busca en internet algunos datos experimentales de presión de vapor de agua como función de la temperatura
#Temperatura en grados Celsius
T=seq(0,100,by=10);T
## [1] 0 10 20 30 40 50 60 70 80 90 100
#Presión de vapor en kPa
P=c(0.611,1.228,2.338,4.246,7.376,12.336,19.92,31.14,47.39,70.39,101.3)
#Primero graficamos los datos para analizar el comportamiento
plot(x = T,y = P,xlab = "Temperatura °C",ylab = 'Presion de vapor (kPa)',main='Presión de vapor del agua a diferentes temperaturas',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
#Construimos un data frame
datos = data.frame(
temperatura = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100) + 273.15, # Temperatura en Kelvin
presion = c(0.611, 1.228, 2.338, 4.246, 7.376, 12.336, 19.92, 31.14, 47.39, 70.39, 101.3) # Presión en kPa
)
#Transformamos los datos
datos_transformados = data.frame(
inverso_temperatura = 1/datos$temperatura,
log_presion = log(datos$presion)
);datos_transformados
## inverso_temperatura log_presion
## 1 0.003660992 -0.4926583
## 2 0.003531697 0.2053868
## 3 0.003411223 0.8492959
## 4 0.003298697 1.4459774
## 5 0.003193358 1.9982315
## 6 0.003094538 2.5125218
## 7 0.003001651 2.9917243
## 8 0.002914177 3.4384932
## 9 0.002831658 3.8584112
## 10 0.002753683 4.2540512
## 11 0.002679887 4.6180864
# Hacemos el ajuste a un modelo lineal con la función lm
modelo = lm(log_presion ~ inverso_temperatura, data = datos_transformados)
# Graficación de los datos y la recta ajustada
plot(1/datos$temperatura, log(datos$presion), xlab = "1/T", ylab = "ln(P)",main='Datos transformados de P y T', sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
abline(modelo, col = "red", lwd = 2)
summary(modelo)
##
## Call:
## lm(formula = log_presion ~ inverso_temperatura, data = datos_transformados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.034176 -0.011029 0.007022 0.017078 0.021353
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.861e+01 6.526e-02 285.1 <2e-16 ***
## inverso_temperatura -5.208e+03 2.078e+01 -250.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02132 on 9 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9998
## F-statistic: 6.279e+04 on 1 and 9 DF, p-value: < 2.2e-16
#De este modelo se obtienen los coeficientes de la ecuación.1/T=m*ln(P)+b
#La ecuación es del tipo ln(P) = m*1/T+b
#m= -5.208e+03
#Intercepto= 1.861e+01
#P=e^(m*1/T+b)
#P=e^(-5.208e+03*1/T+1.861e+01)
P1=exp(-5.208e+03*1/285+1.861e+01);P1
## [1] 1.399781
P2=exp(-5.208e+03*1/300+1.861e+01);P2
## [1] 3.490343
#Este modelo entonces obtiene estos valores estimados para las temperaturas de 285K y 300K
#Hacemos una gráfica para explorar los datos
#Construimos un dataframe
datos = data.frame(
Temperatura = c(10, 20, 30, 40, 50, 60, 70, 80, 90), # Temperatura en Celsius
Solubilidad_NaCl = c(35,35.6,36.25,36.9,37.5,38.1,38.8,39.4,40) # Solubilidad en g NaCl / 100g H2O
)
plot(datos,main='Solubilidad NaCl en H2O',xlab='T(°C)',ylab='Sol (gNaCl/100gH2O)',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
#Inicialmente se observa una distribución lineal, se observará el coeficiente de correlación de varios modelos para encontrar el mejor.
Modelo_Lineal= lm(Solubilidad_NaCl ~ Temperatura, data = datos)
summary(Modelo_Lineal)
##
## Call:
## lm(formula = Solubilidad_NaCl ~ Temperatura, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033889 -0.018889 0.001111 0.009444 0.037778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.436e+01 1.749e-02 1965.2 < 2e-16 ***
## Temperatura 6.283e-02 3.107e-04 202.2 1.91e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02407 on 7 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 4.089e+04 on 1 and 7 DF, p-value: 1.91e-14
#El MOdelo lineal tiene una correlación lo suficientemente buena como para definir que es un modelo correcto
Solubilidad_NaCl=function(Temperatura){Temperatura*0.06283+34.36}
S25=Solubilidad_NaCl(25);S25
## [1] 35.93075
cat('La solubilidad de NaCl a 25°C es',S25,'gNaCl/100gH2O\n')
## La solubilidad de NaCl a 25°C es 35.93075 gNaCl/100gH2O
#Construimos un dataframe
datos = data.frame(
Temperatura = seq(5,45,by=5), # Temperatura en Celsius
Solubilidad_O2 = c(1.95,1.7,1.55,1.40,1.30,1.15,1.05,1.00,0.95) # Solubilidad en milimoles O2/LH2O)
)
plot(datos,main='Solubilidad O2 en H2O',xlab='T(°C)',ylab='Sol (mmO2/LH2O)',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
#Ahora miramos cual es el modelo que mejor se ajusta a la curva
#MOdelo lineal
Modelo_Lineal= lm(Solubilidad_O2 ~ Temperatura, data = datos)
summary(Modelo_Lineal)
##
## Call:
## lm(formula = Solubilidad_O2 ~ Temperatura, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06639 -0.04389 -0.03389 0.02861 0.12111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.951389 0.053887 36.21 3.18e-09 ***
## Temperatura -0.024500 0.001915 -12.79 4.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07418 on 7 degrees of freedom
## Multiple R-squared: 0.959, Adjusted R-squared: 0.9531
## F-statistic: 163.6 on 1 and 7 DF, p-value: 4.135e-06
#El MOdelo lineal tiene un R2 de 0.9531
#Ahora miramos un modelo exponencial, dada la forma de los datos
datos_exp = data.frame(
Temperatura = datos$Temperatura,
log_S = log(datos$Solubilidad_O2)
);datos_exp
## Temperatura log_S
## 1 5 0.66782937
## 2 10 0.53062825
## 3 15 0.43825493
## 4 20 0.33647224
## 5 25 0.26236426
## 6 30 0.13976194
## 7 35 0.04879016
## 8 40 0.00000000
## 9 45 -0.05129329
Modelo_Exponencial= lm(log_S ~ Temperatura, data = datos_exp)
summary(Modelo_Exponencial)
##
## Call:
## lm(formula = log_S ~ Temperatura, data = datos_exp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033388 -0.017907 -0.005218 0.008555 0.047996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7173133 0.0224565 31.94 7.62e-09 ***
## Temperatura -0.0181467 0.0007981 -22.74 8.06e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03091 on 7 degrees of freedom
## Multiple R-squared: 0.9866, Adjusted R-squared: 0.9847
## F-statistic: 517 on 1 and 7 DF, p-value: 8.063e-08
#El MOdelo Exponencial tiene un R2 de 0.9847
#Con esta correlación la ecuación para obtener la solubilidad es
Solubilidad_O2=function(Temperatura){exp(Temperatura*-0.0181467+0.7173133)}
S8=Solubilidad_O2(8);S8
## [1] 1.772055
S50=Solubilidad_O2(50);S50
## [1] 0.8269412
cat('La solubilidad de O2 a 8°C es',S8,'mmO2/LH2O\n')
## La solubilidad de O2 a 8°C es 1.772055 mmO2/LH2O
cat('La solubilidad de O2 a 50°C es',S50,'mmO2/LH2O\n')
## La solubilidad de O2 a 50°C es 0.8269412 mmO2/LH2O
#Construimos un dataframe
datos = data.frame(
x = seq(1,10,by=1), # x
y = c(10,14,16,18,19,20,21,22,23,23) # y
)
#Transformamos los datos ahora
datos_transf = data.frame(
lnx = log(datos$x), # lnx
y = datos$y # y
)
plot(datos,main='y vs lnx',xlab='lnx',ylab='y',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
Modelo=lm(y~lnx,data = datos_transf)
summary(Modelo)
##
## Call:
## lm(formula = y ~ lnx, data = datos_transf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23126 -0.16611 -0.00851 0.11077 0.44979
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.9123 0.1721 57.59 9.17e-12 ***
## lnx 5.7518 0.1035 55.57 1.22e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2276 on 8 degrees of freedom
## Multiple R-squared: 0.9974, Adjusted R-squared: 0.9971
## F-statistic: 3088 on 1 and 8 DF, p-value: 1.22e-11
#Con esto obtenemos la función
#Los parámetros a1 y a2 son 9.9123 y 5.7518
Fy=function(x){9.9123+5.7518*log(x)}
F1.5=Fy(1.5);F1.5
## [1] 12.24445
F11=Fy(11);F11
## [1] 23.70451
cat('El valor de y para x=1.5 es',F1.5,'\n')
## El valor de y para x=1.5 es 12.24445
cat('El valor de y para x=11 es',F11,'\n')
## El valor de y para x=11 es 23.70451
#Ejercicio73
x=seq(-10,10,length=11)
y=seq(-10,10,length=11)
func=function(x,y) x^2-(2*x*y)+(4*y^2)
z=outer(x,y,func)
library(plotly)
fig <- plot_ly(
type='surface',
contours=list(
x=list(
show=TRUE,
start=0,
end=1,
size=0.1,
color='white'),
y=list(
show=TRUE,
start=0,
end=1,
size=0.1)),
x = ~x,
y = ~y,
z = ~z)
fig
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 5),
zaxis = list(nticks =5),
camera = list(eye = list(x = 0, y = -1, z = 0.5)),
aspectratio = list(x = 1, y = 1, z = 1)))
fig
#Ejercicio74
x=seq(-10,10,length=100)
y=seq(-10,10,length=100)
func=function(x,y) -(x^2)+(2*x*y)+(3*y^2)
z=outer(x,y,func)
library(plotly)
fig <- plot_ly(
type='surface',
contours=list(
x=list(
show=TRUE,
start=0,
end=1,
size=0.1,
color='white'),
y=list(
show=TRUE,
start=0,
end=1,
size=0.1)),
x = ~x,
y = ~y,
z = ~z)
fig
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 5),
zaxis = list(nticks =5),
camera = list(eye = list(x = 0, y = -1, z = 0.5)),
aspectratio = list(x = .9, y = .8, z = .9)))
fig
x=seq(-10,10,length=100)
y=seq(-10,10,length=100)
func=function(x,y) (x-y^2)*(x-(3*y^2))
z=outer(x,y,func)
library(plotly)
fig <- plot_ly(
type='surface',
contours=list(
x=list(
show=TRUE,
start=0,
end=1,
size=0.1,
color='white'),
y=list(
show=TRUE,
start=0,
end=1,
size=0.1)),
x = ~x,
y = ~y,
z = ~z)
fig
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 5),
zaxis = list(nticks =5),
camera = list(eye = list(x = 0, y = -1, z = 0.5)),
aspectratio = list(x = 1, y = 1, z = 1)))
fig
#Ejercicio76
x=seq(0,1,length=11)
y=seq(0,1,length=11)
func=function(x,y) 80*exp(-(x-1)^2)*exp(-3*(y-1)^2)
T=outer(x,y,func)
library(plotly)
fig <- plot_ly(
type='surface',
contours=list(
x=list(
show=TRUE,
start=0,
end=1,
size=0.1,
color='white'),
y=list(
show=TRUE,
start=0,
end=1,
size=0.1)),
x = ~x,
y = ~y,
z = ~T)
fig
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 11),
zaxis = list(nticks =11),
camera = list(eye = list(x = 0, y = -1, z = 0.5)),
aspectratio = list(x = 1, y = 1, z = 1)))
fig
#Ejercicio 77
library(ggplot2)
gm<-c(23,25,26,25,27,25,24,22,23,25,26,26,24,24,22,25,26,24,24,24,27,23)
dt <- data.frame(gm)
ggplot(dt, aes (x=gm))+xlab('Gas miliage in miles per gallon')+ylab('Number of cars')+ggtitle('Absolute frequency histogram')+geom_histogram(bins=8)
ggplot(dt, aes (x=gm))+xlab('Gas miliage in miles per gallon')+ylab('Percent of cars')+ggtitle('Relative frequency histogram')+geom_histogram(aes(y=after_stat(count/sum(count))),bins=8)+scale_y_continuous(labels=scales::percent)
#Ejercicio 77
x<-c(23,25,26,25,27,25,24,22,23,25,26,26,24,24,22,25,26,24,24,24,27,23)
hist(x,freq = TRUE,main="histograma de frecuencia", xlab = "miles / galon")
#Ejercicio 78
x<-c(243,236,389,628,143,417,205,404,464,605,137,123,372,439,497,500,535,577,441,231,675,132,196,217,660,569,865,725,457,347)
intervals1<-seq(min(x),max(x)+50,by=50)
intervals2<-seq(min(x),max(x)+100,by=100)
intervals3<-seq(min(x),max(x)+200,by=200)
hist(x, breaks=5, freq = TRUE, main = "Histograma de frecuencia absoluta", xlab="Fuerza en Libras")
hist(x, breaks=intervals1, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 50 lb", xlab="Fuerza en Libras", ylab="Frecuencia")
hist(x, breaks=intervals2, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 100 lb", xlab="Fuerza en Libras", ylab="Frecuencia")
hist(x, breaks=intervals3, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 200 lb", xlab="Fuerza en Libras", ylab="Frecuencia")
#DE LOS 3 VALORES DE ANCHO DE BLOQUE SOLICITADOS POR EL EJERCICIO, EL MEJOR ES EL DE 100 LB, YA QUE EN ANCHO DE 50 LB, ALGUNOS BLOQUES NO TIENEN VALORES, Y EN ANCHO DE 200 LB SE AGRUPAN MUCHOS VALORES EN CADA BLOQUE, LO QUE LIMITA EL ANALISIS.
intervals4<-seq(min(x),max(x)+95,by=95)
hist(x, breaks=intervals4, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 95 LB", xlab="Fuerza en Newtons", ylab="Frecuencia")
#CON ANCHO DE 95 LB PARA LOS BLOQUES SE TIENE MAYOR DISTRIBUCION DE LOS VALORES, MIENTRAS NINGUN BLOQUE QUEDA SIN DATOS.
#Ejercicio 79
x<-c(311,138,340,199,270,255,332,279,231,296,198,269,257,236,313,281,288,225,216,250,259,323,280,205,279,159,276,354,278,221,192,281,204,361,321,282,254,273,334,172,240,327,261,282,208,213,299,318,356,269,355,232,275,234,267,240,331,222,370,226)
hist(x,freq=TRUE,main="Histograma de frecuencia absoluta", xlab="Fuerza en Newtons", ylab="Frecuencia")
intervals1<-seq(min(x),max(x)+10,by=10)
intervals2<-seq(min(x),max(x)+30,by=30)
intervals3<-seq(min(x),max(x)+50,by=50)
hist(x, breaks=intervals1, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 10 N", xlab="Fuerza en Newtons", ylab="Frecuencia")
hist(x, breaks=intervals2, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 30 N", xlab="Fuerza en Newtons", ylab="Frecuencia")
hist(x, breaks=intervals3, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 50 N", xlab="Fuerza en Newtons", ylab="Frecuencia")
#DE LOS 3 VALORES DE ANCHO DE BLOQUE SOLICITADOS POR EL EJERCICIO, EL MEJOR ES EL DE 30 N, YA QUE EN ANCHO DE 10 N, ALGUNOS BLOQUES NO TIENEN VALORES, Y EN ANCHO DE 50 SE AGRUPAN MUCHOS VALORES EN CADA BLOQUE, LO QUE LIMITA EL ANALISIS.
intervals4<-seq(min(x),max(x)+15,by=15)
hist(x, breaks=intervals4, freq = TRUE, main = "Histograma de frecuencia absoluta ancho de 15 N", xlab="Fuerza en Newtons", ylab="Frecuencia")
#CON ANCHO DE 15 N PARA LOS BLOQUES SE TIENE MAYOR DISTRIBUCION DE LOS VALORES, MIENTRAS NINGUN BLOQUE QUEDA SIN DATOS.
#Ejercicio 80
x<-c(23,25,26,25,27,25,24,22,23,25,26,26,24,24,22,25,26,24,24,24,27,23)
mean_x <- mean(x)
mean_x
## [1] 24.54545
sd_x <- sd(x)
sd_x
## [1] 1.438494
hist(x,freq = FALSE,probability = TRUE, main="histograma de frecuencia escalada", xlab = "miles / galon")
abline(v=mean(x), col="violet")
#DEFINIR EL LIMITE SUPERIOR E INFERIOR PARA EL 68% DE LOS DATOS, DE ACUERDO A LA REGLA EMPIRICA CORRESPONDE A UNA DESVIACION ESTANDAR ARRIBA Y ABAJO DE LA MEDIA.
lower_limit <- mean_x - sd_x
lower_limit
## [1] 23.10696
upper_limit <- mean_x + sd_x
upper_limit
## [1] 25.98395
sum(x >= lower_limit & x <= upper_limit)
## [1] 11
#Ejercicio 81
x<-c(243,236,389,628,143,417,205,404,464,605,137,123,372,439,497,500,535,577,441,231,675,132,196,217,660,569,865,725,457,347)
mean_x <- mean(x)
mean_x
## [1] 414.3
sd_x <- sd(x)
sd_x
## [1] 198.4462
hist(x,freq = FALSE,probability = TRUE, main="histograma de frecuencia escalada", xlab = "FUERZA EN LIBRAS")
abline(v=mean(x), col="violet")
#DEFINIR EL LIMITE SUPERIOR E INFERIOR PARA EL 68% DE LOS DATOS, DE ACUERDO A LA REGLA EMPIRICA CORRESPONDE A UNA DESVIACION ESTANDAR ARRIBA Y ABAJO DE LA MEDIA.
lower_limit68 <- mean_x - sd_x
lower_limit68
## [1] 215.8538
upper_limit68 <- mean_x + sd_x
upper_limit68
## [1] 612.7462
#DEFINIR EL LIMITE SUPERIOR E INFERIOR PARA EL 96% DE LOS DATOS, DE ACUERDO A LA REGLA EMPIRICA CORRESPONDE A DOS DESVIACIONES ESTANDAR ARRIBA Y ABAJO DE LA MEDIA.
lower_limit96 <- mean_x - (sd_x*2)
lower_limit96
## [1] 17.40754
upper_limit96 <- mean_x + (sd_x*2)
upper_limit96
## [1] 811.1925
sum(x >= lower_limit68 & x <= upper_limit68)
## [1] 19
sum(x >= lower_limit96 & x <= upper_limit96)
## [1] 29
#Ejercicio 82
x<-c(311,138,340,199,270,255,332,279,231,296,198,269,257,236,313,281,288,225,216,250,259,323,280,205,279,159,276,354,278,221,192,281,204,361,321,282,254,273,334,172,240,327,261,282,208,213,299,318,356,269,355,232,275,234,267,240,331,222,370,226)
mean_x <- mean(x)
mean_x
## [1] 266.95
sd_x <- sd(x)
sd_x
## [1] 52.82123
hist(x,freq = FALSE,probability = TRUE, main="histograma de frecuencia escalada", xlab = "FUERZA EN NEWTONS")
abline(v=mean(x), col="violet")
#DEFINIR EL LIMITE SUPERIOR E INFERIOR PARA EL 68% DE LOS DATOS, DE ACUERDO A LA REGLA EMPIRICA CORRESPONDE A UNA DESVIACION ESTANDAR ARRIBA Y ABAJO DE LA MEDIA.
lower_limit68 <- mean_x - sd_x
lower_limit68
## [1] 214.1288
upper_limit68 <- mean_x + sd_x
upper_limit68
## [1] 319.7712
#DEFINIR EL LIMITE SUPERIOR E INFERIOR PARA EL 96% DE LOS DATOS, DE ACUERDO A LA REGLA EMPIRICA CORRESPONDE A DOS DESVIACIONES ESTANDAR ARRIBA Y ABAJO DE LA MEDIA.
lower_limit96 <- mean_x - (sd_x*2)
lower_limit96
## [1] 161.3075
upper_limit96 <- mean_x + (sd_x*2)
upper_limit96
## [1] 372.5925
sum(x >= lower_limit68 & x <= upper_limit68)
## [1] 38
sum(x >= lower_limit96 & x <= upper_limit96)
## [1] 58
#Ejercicio 83
#83a
1-pnorm(194, mean = 200, sd = sqrt(9))
## [1] 0.9772499
#83b
pnorm(203, mean = 200, sd = sqrt(9)) - pnorm(197, mean = 200, sd = sqrt(9))
## [1] 0.6826895
#Ejercicio 84
z = (60 - 50) / 5
p = pnorm(60, mean = 50, sd = 5)
1 - p
## [1] 0.02275013
#Probabilidad de 2.28%
#Ejercicio 85
# Definiendo el promedio y la desviación de la distribución normal
mu <- 5.007
sigma <- 0.005
p_up <- pnorm(5.01, mean = mu, sd = sigma)
p_low <- pnorm(4.99, mean = mu, sd = sigma)
# En la tolerancia
p_within_tolerance <-(p_up-p_low) * 100
# Resultado
cat(sprintf("Porcentaje de conexiones dentro de la tolerancia: %.2f%%", p_within_tolerance))
## Porcentaje de conexiones dentro de la tolerancia: 72.54%
#Ejercicio 86
#86a
mean_c <- 3 - 2.96
var_c <- 0.0064 + 0.0036
cat("Promedio: ", mean_c, "\n")
## Promedio: 0.04
cat("Varianza: ", var_c, "\n")
## Varianza: 0.01
#86b
#P(c < 0) = P(Z < -10*0.04) = P(Z < -0.4)
pnorm(-0.4)
## [1] 0.3445783
#Ejercicio 87
#87a
parts_per_box <- 300
mean_part_weight <- 1
sd_part_weight <- 0.2
box_weight <- parts_per_box * mean_part_weight
box_weight
## [1] 300
variance_per_part <- 0.2^2
variance_per_part
## [1] 0.04
boxes_per_pallet <- 10
variance_per_box <- variance_per_part * 300
variance_per_box
## [1] 12
sd_per_box <- sqrt(variance_per_box)
sd_per_box
## [1] 3.464102
mean_weight_per_pallet <- 10 * box_weight
mean_weight_per_pallet
## [1] 3000
sd_per_pallet <- sqrt(10) * sd_per_box
sd_per_pallet
## [1] 10.95445
#87b
pnorm_cum<-pnorm(3015,mean_weight_per_pallet,sd_per_pallet)
1-pnorm_cum
## [1] 0.08545176
#Ejercicio 88
#88a
mean <- 1 + 2 + 1.5
mean
## [1] 4.5
variance <- 0.00014 + 0.0002 + 0.0003
variance
## [1] 0.00064
#88b
z1 <- (4.48 - mean) / sqrt(variance)
z2 <- (4.52 - mean) / sqrt(variance)
p <- pnorm(z2) - pnorm(z1)
p
## [1] 0.5708047
#Ejercicio 89
set.seed(123) # Set the seed for reproducibility
n <- 1000
min <- 2
max <- 18
x <- runif(n, min, max)
mean(x)
## [1] 9.956445
hist(x, breaks = seq(2, 18, by = 0.5), col = "lightblue", xlab = "Random Numbers", main = "Histogram of Random Numbers")
# La distribución obtenida se observa relativamente uniforme con números generados entre 2 y 18 con media cercana a 10, visualmente parece que los números se distribuyen uniformemente con la media deseada, ya que no hay agrupación o asimetría evidentes en el histograma.
#Ejercicio 90
set.seed(123)
x <- rnorm(1000, mean = 20, sd = 2)
mean(x)
## [1] 20.03226
hist(x, breaks = 30)
qqnorm(x)
qqline(x)
#Se puede observar que los datos siguen una distribución normal en ambos gráficos puesto que en el histograma se puede observar la forma de campana y en el gráfico normal se observa el la linealidad.
#Ejercicio 91
set.seed(123)
mean_x <- 10
sd_x <- 2
mean_y <- 15
sd_y <- 3
trials <- c(100, 1000, 5000)
simulate_sum <- function(n) {
x <- rnorm(n, mean_x, sd_x)
y <- rnorm(n, mean_y, sd_y)
z <- x + y
mean_z <- mean(z)
var_z <- var(z)
list(mean = mean_z, variance = var_z)
}
results <- lapply(trials, simulate_sum)
for (i in 1:length(trials)) {
cat("Number of trials:", trials[i], "\n")
cat("Simulated mean of z:", results[[i]]$mean, "\n")
cat("Theoretical mean of z:", mean_x + mean_y, "\n")
cat("Simulated variance of z:", results[[i]]$variance, "\n")
cat("Theoretical variance of z:", sd_x^2 + sd_y^2, "\n\n")
}
## Number of trials: 100
## Simulated mean of z: 24.85817
## Theoretical mean of z: 25
## Simulated variance of z: 11.22385
## Theoretical variance of z: 13
##
## Number of trials: 1000
## Simulated mean of z: 25.07761
## Theoretical mean of z: 25
## Simulated variance of z: 13.75186
## Theoretical variance of z: 13
##
## Number of trials: 5000
## Simulated mean of z: 24.99773
## Theoretical mean of z: 25
## Simulated variance of z: 12.99906
## Theoretical variance of z: 13
x <- function(n){
rnorm(n,10,sqrt(2))
}
y <- function(n){
rnorm(n,15,sqrt(3))
}
z <- function(n){
x(n)*y(n)
}
resultados <- data.frame(
"Trials" = c(100,1000,5000),
"Media x por y)" = c(mean(x(100))*mean(y(100)) , mean(x(1000))*mean(y(1000)) , mean(x(5000))*mean(y(5000))),
"Media(z)" = c(mean(z(100)),mean(z(1000)),mean(z(5000))),
"Varianza x por y)" = c(var(x(100))*var(y(100)) , var(x(1000))*var(y(1000)) , var(x(5000))*var(y(5000))),
"Varianza(z)" = c(var(z(100)),var(z(1000)),var(z(5000)))
)
resultados
## Trials Media.x.por.y. Media.z. Varianza.x.por.y. Varianza.z.
## 1 100 146.9607 148.6575 9.857421 840.5323
## 2 1000 150.2524 148.5724 6.007688 865.8661
## 3 5000 149.3809 150.0606 5.965283 767.8845
temp_df <- data.frame(
tiempo = c(1,2,3,4,6,7,8,10,11,12),
temperatura = c(10, 12, 18, 24, 21, 20, 18, 15, 13, 8)
)
interp_temp <- approx(temp_df$tiempo, temp_df$temperatura, xout = c(5,9))
print(interp_temp)
## $x
## [1] 5 9
##
## $y
## [1] 22.5 16.5
temp_df <- data.frame(
hora = c(1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5),
dia = c("Lunes","Lunes","Lunes","Lunes","Lunes","Martes","Martes","Martes","Martes","Martes","Miercoles","Miercoles","Miercoles","Miercoles","Miercoles","Jueves","Jueves","Jueves","Jueves","Jueves","Viernes","Viernes","Viernes","Viernes","Viernes"),
temperatura = c(17,13,14,17,23,15,NA,14,15,18,12,8,9,14,17,16,11,NA,15,20,16,12,15,19,24)
)
interp_temp <- approx(temp_df$hora, temp_df$temperatura, xout = c(2,4))
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
x <- c(0, 0.25, 0.75, 1.25, 1.5, 1.75, 1.875, 2, 2.125, 2.25)
y <- c(1.2, 1.18, 1.1, 1, 0.92, 0.8, 0.7, 0.55, 0.35, 0)
spline_fit <- spline(x, y, n = 100)
plot(x, y, pch = 16, col = "blue", xlim = c(0, 2.5), ylim = c(0, 1.5), xlab = "x(ft)", ylab = "y(ft)", main = "Cubic Splines Fit")
lines(spline_fit$x, spline_fit$y, col = "red", lwd = 2)
legend("topright", legend = c("Coordenadas", "Cubic Splines"), col = c("blue", "red"), pch = c(16, NA), lwd = c(NA, 2))
tiempo <- 0:10
temperatura <- c(72.5, 78.1, 86.4, 92.3, 110.6, 111.5, 109.3, 110.2, 110.5, 109.9, 110.2)
plot(tiempo, temperatura, type = "o", col = "VIOLET", xlab = "Tiempo (minutos)", ylab = "Temperatura (°F)")
plot(tiempo, temperatura, type = "o", col = "blue", xlab = "Tiempo (minutos)", ylab = "Temperatura (°F)")
lines(spline(tiempo, temperatura, n = 100), col = "red")
# Interpolación lineal
t_lin <- c(0.6, 2.5, 4.7, 8.9)
T_lin <- approx(tiempo, temperatura, t_lin)$y
T_lin
## [1] 75.86 89.35 111.23 109.96
# Interpolación inversa lineal
T_objetivo <- c(75, 85, 90, 105)
t_objetivo_lin <- approx(temperatura, tiempo, T_objetivo)$y
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
t_objetivo_lin
## [1] 0.4464286 1.8313253 2.6101695 5.2411765
v <- function(t) {
5 + 7*t^2
}
x <- function(t) {
5 + integrate(v, 2, t)$value
}
x(10)
## [1] 2359.667
tasa_bombeo <- c(0, 80, 130, 150, 150, 160, 165, 170, 160, 140, 120)
volumen_total <- cumsum(tasa_bombeo)
area_tanque <- 100 # ft^2
altura_agua <- volumen_total[11] / area_tanque
altura_agua
## [1] 14.25
Para calcular la integral indefinida seguimos el siguiente cálculo
\[\int p(x)dx=\int \left ( 5x^{2}-6x+8 \right )dx\] \[\int p(x)dx=5\int x^{2}\, dx-6\int x\, dx+8\int 1\, dx\] \[\int p(x)dx=5\left ( \frac{x^{3}}{3} \right )+6\left ( \frac{x^{2}}{2} \right )+8\, x+C\]
#La integral indefinida de la función es
INtp=function(x,C){5*x^3/3-3*x^2+8*x+C}
#Para integrar en R Primero definimos la función
p=function(x){5*x^2-6*x+8}
#Con la función integrate podemos calcular la integral definida de esta función entre 2 valores dados, por ejemplo usaremos los valores 0 y 10 como intervalo para encontrar la integral.
integral=integrate(p,lower = 0,upper = 10);integral
## 1446.667 with absolute error < 1.6e-11
#Esta da un valor de 1446,667
#Para estimar la derivada dy/dx de los siguientes datos primero debemos buscar la función que más se aproxime a la curva.Primero hacemos un gráfico exploratorio.
datos=data.frame(x=seq(0,10,by=1),y=c(0,2,5,7,9,12,15,18,22,20,17))
plot(datos,main='Datos y vs x Ejercicio 102',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
#Intenté buscar una función en las librerías pracma y numDeriv para poder obtener las derivadas hacia adelante, hacia atrás y centrales pero no lo encontré.
#La derivada hacia adelante f'(x)≈[f(x+h)-f(x)]/h' donde h=(xi+1)-(xi)
#Aquí colocamos los vectores
x=seq(0,10,by=1)
y=c(0,2,5,7,9,12,15,18,22,20,17)
#h es la diferencia del vector x
h=diff(x)
#Para calcular la derivada hacia adelante, donde requerimos un vector de longitud y-1. la función numeric crea un vector de longitud y-1, ya que no se puede calcular la derivada hacia adelante del último valor
dy_dx_adelante= c(length = (length(y) - 1))
#Ahora usamos la función for, donde para cada valor i en un rango desde 1 hasta y-1 calcula cada valor del vector dy_dx_adelante
for (i in 1:(length(y) - 1)) {
dy_dx_adelante[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
}
dy_dx_adelante
## length
## 2 3 2 2 3 3 3 4 -2 -3
#La derivada hacia atrás es f'(x)≈[f(x)-f(x-h)]/h' donde h=(xi+1)-(xi)
#Hacemos los mismos pasos que la anterior, en este caso la longitud del vector es y-1 porque no se pued calcular el primer dato
dy_dx_atras = c(length = (length(y) - 1))
for (i in 2:length(y)) {
dy_dx_atras[i] = (y[i] - y[i-1]) / (x[i] - x[i-1])
}
dy_dx_atras
## length
## 10 2 3 2 2 3 3 3 4 -2 -3
#Este valor es igual que las diferencias hacia adelante solo que arrancan desde el punto i=2
#La derivada central es f'(x)≈[f(x+h)-f(x-h)]/2h' donde h=((xi+1)-(xi).
#Aquí la longitud es menor ya que no se puede calcular el dato para el primer ni el último valor de y
dy_dx_central = c(length = length(y) - 2)
for (i in 2:(length(y)-1)) {
dy_dx_central[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1])
}
dy_dx_central
## length
## 9.0 2.5 2.5 2.0 2.5 3.0 3.0 3.5 1.0 -2.5
#Aquí se obtienen las derivadas centrales, las cuales son el promedio de las derivadas hacia adelante y hacia atrás, cabe anotar que este vector tiene solamente 9 datos, mientras que los 2 anteriores tenían 10.
#Ahora dibujamos las derivadas
plot(x,y,xlab = "x",ylab = 'y',main='Datos y vs x Ejercicio 102',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='blue',pch=16,ylim=c(-4,23),xlim=c(-0.5,11))
lines(dy_dx_adelante,col='red',lwd=2)
lines(dy_dx_atras,col='orange',lwd=2)
lines(dy_dx_central,col='green',lwd=2)
legend("topleft", legend=c("Datos","deriv. adelante","deriv. atrás","deriv. central"), fill=c("blue","red","orange","green"))
#Se observa como la pendiente es casi una línea recta en los primeros datos y luego pasa a un valor casi similar pero negativo, la derivada central puede mostrar unos datos suavizados más reales.
#Primero graficamos los datos para hacer una exploración inicial
x=seq(0,10,by=1)
y=c(0,2,5,7,9,10,8,7,6,8,10)
plot(x,y,xlab = "x",ylab = 'y',main='Datos y vs x Ejercicio 103',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='blue',pch=16)
#Calculamos ahora las derivadas hacia adelante y centrales y las graficamos para encontrar los puntos donde la derivada es cero
dy_dx_atras = c(length = (length(y) - 1))
for (i in 2:length(y)) {
dy_dx_atras[i] = (y[i] - y[i-1]) / (x[i] - x[i-1])
}
dy_dx_atras
## length
## 10 2 3 2 2 1 -2 -1 -1 2 2
dy_dx_adelante= c(length = (length(y) - 1))
for (i in 1:(length(y) - 1)) {
dy_dx_adelante[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
}
dy_dx_adelante
## length
## 2 3 2 2 1 -2 -1 -1 2 2
dy_dx_central = c(length = length(y) - 2)
for (i in 2:(length(y)-1)) {
dy_dx_central[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1])
}
dy_dx_central
## length
## 9.0 2.5 2.5 2.0 1.5 -0.5 -1.5 -1.0 0.5 2.0
plot(x = x,y = y,xlab = "x",ylab = 'dy/dx',main='dy/dx vs x Ejercicio 103',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='green',pch=16,ylim=c(-3,4))
lines(dy_dx_adelante,col='blue',lwd=2)
lines(dy_dx_central,col='red',lwd=2)
lines(dy_dx_atras,col='orange',lwd=2)
#Al trazar una línea horizontal podemos explorar el punto en el que los datos tienen une pendiente de cero
legend("bottomleft", legend=c("Datos","deriv. adelante","deriv. central","derv.atras"), fill=c("green","blue","red","orange"))
abline(h = 0,col='green',lwd=2)
#Ahora buscamos el intercepto buscando la abline correspondiente vertical
abline(v = 5.35,col='green',lwd=2)
#Por derivada hacia adelante el punto está cerca a 5.35 en x
abline(v = 5.75,col='green',lwd=2)
#Por derivadas centrales el punto está cerca a 5.75 en x
abline(v = 6.34,col='green',lwd=2)
#Por derivadas hacia atrás el punto está cerca a 6.34 en x
#Para confirmar este valor volvemos a la gráfica, pero esta vez trazamos en líneas
plot(x,y,xlab = "x",ylab = 'y',main='Datos y vs x Ejercicio 103',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,type="o",col='blue',pch=16)
abline(v = 5.35,col='green',lwd=2)
abline(v = 5.75,col='green',lwd=2)
abline(v = 6.34,col='green',lwd=2)
#El valor más cercano usando derivadas hacia adelante es el de 5.35, haciendo una exploración de los datos, para este valor de x, el y es 9.28
abline(h = 9.28,col='green',lwd=2)
#Por tanto se concluye que para los datos dados el valor máximo de ye es 10, el cual se obtiene para valores de x y y.
#Para tener mayor precisión se podría buscar una ecuación de una curva con un ajuste bastante cercano, que podría ser un polinomio de mayor grado, (4 o más), y buscar el valor máximo posible en su derivada = 0, que por la inspección de los datos podría estar por debajo de 5 un poco.
#Colocamos los datos xi = x inicial (sin el error aditivo)
xi=seq(0,4,length=101)
#Creamos ahora el error aditivo aleatorio usando una distribución normal, con 101 datos y desviación estándar de 0.01
error=rnorm(length(xi), mean = 0, sd = 0.01)
x=xi + error
f_x=function(x){exp(-x)*sin(3*x)}
y=f_x(x)
#Ahora calculamos las derivadas hacia atrás, hacia adelante y central usando el algoritmo que ya habíamos creado.
dy_dx_atras = c(length = (length(y) - 1))
for (i in 2:length(y)) {
dy_dx_atras[i] = (y[i] - y[i-1]) / (x[i] - x[i-1])
}
dy_dx_atras
## length
## 100.000000000 2.844277523 2.589862338 2.292916729 2.008831727
##
## 1.714782203 1.350276019 1.026281996 0.721366946 0.476834265
##
## 0.266310483 0.022397467 -0.276155443 -0.502539560 -0.652437104
##
## -0.803502320 -0.966377133 -1.090511198 -1.175162994 -1.246439331
##
## -1.287315625 -1.301274597 -1.300421642 -1.265644410 -1.214815005
##
## -1.161149065 -1.100324771 -1.006698916 -0.896277830 -0.826281254
##
## -0.746610795 -0.634812139 -0.504753162 -0.396024210 -0.300269980
##
## -0.177152770 -0.096140614 -0.050700473 0.052642690 0.184062949
##
## 0.246842591 0.271646541 0.326480445 0.377445377 0.402267948
##
## 0.430638197 0.450420339 0.456435611 0.457177904 0.448541458
##
## 0.437475412 0.421072704 0.393441883 0.367315662 0.329394333
##
## 0.295119118 0.265658020 0.232676542 0.190168367 0.149502288
##
## 0.113141755 0.073533981 0.036227697 0.005949042 -0.021572769
##
## -0.051499908 -0.076049252 -0.096672686 -0.112907073 -0.128636709
##
## -0.140346402 -0.149750935 -0.156752856 -0.159923041 -0.160264198
##
## -0.158298350 -0.153788845 -0.148986276 -0.141031728 -0.129332539
##
## -0.119784251 -0.107041958 -0.095270057 -0.080118856 -0.063731115
##
## -0.048148468 -0.037639712 -0.027910879 -0.014233183 -0.004759362
##
## 0.002489296 0.013486083 0.024834182 0.032113937 0.039539388
##
## 0.044825769 0.048456213 0.052127625 0.054678457 0.055972171
##
## 0.056292906
dy_dx_adelante= c(length = (length(y) - 1))
for (i in 1:(length(y) - 1)) {
dy_dx_adelante[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
}
dy_dx_adelante
## length
## 2.844277523 2.589862338 2.292916729 2.008831727 1.714782203 1.350276019
##
## 1.026281996 0.721366946 0.476834265 0.266310483 0.022397467 -0.276155443
##
## -0.502539560 -0.652437104 -0.803502320 -0.966377133 -1.090511198 -1.175162994
##
## -1.246439331 -1.287315625 -1.301274597 -1.300421642 -1.265644410 -1.214815005
##
## -1.161149065 -1.100324771 -1.006698916 -0.896277830 -0.826281254 -0.746610795
##
## -0.634812139 -0.504753162 -0.396024210 -0.300269980 -0.177152770 -0.096140614
##
## -0.050700473 0.052642690 0.184062949 0.246842591 0.271646541 0.326480445
##
## 0.377445377 0.402267948 0.430638197 0.450420339 0.456435611 0.457177904
##
## 0.448541458 0.437475412 0.421072704 0.393441883 0.367315662 0.329394333
##
## 0.295119118 0.265658020 0.232676542 0.190168367 0.149502288 0.113141755
##
## 0.073533981 0.036227697 0.005949042 -0.021572769 -0.051499908 -0.076049252
##
## -0.096672686 -0.112907073 -0.128636709 -0.140346402 -0.149750935 -0.156752856
##
## -0.159923041 -0.160264198 -0.158298350 -0.153788845 -0.148986276 -0.141031728
##
## -0.129332539 -0.119784251 -0.107041958 -0.095270057 -0.080118856 -0.063731115
##
## -0.048148468 -0.037639712 -0.027910879 -0.014233183 -0.004759362 0.002489296
##
## 0.013486083 0.024834182 0.032113937 0.039539388 0.044825769 0.048456213
##
## 0.052127625 0.054678457 0.055972171 0.056292906
dy_dx_central = c(length = length(y) - 2)
for (i in 2:(length(y)-1)) {
dy_dx_central[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1])
}
dy_dx_central
## length
## 99.000000000 2.739893065 2.405215115 2.201192458 1.809062372 1.551439989
##
## 1.187957314 0.879467392 0.618428258 0.364838170 0.128254000 -0.149905196
##
## -0.347226370 -0.589323865 -0.725951473 -0.904591706 -1.007500462 -1.143923137
##
## -1.207400727 -1.268166069 -1.289502507 -1.300513596 -1.285020173 -1.246144023
##
## -1.187398541 -1.132693635 -1.039657335 -0.963675596 -0.872346689 -0.771511934
##
## -0.691036801 -0.561151378 -0.469382574 -0.334899629 -0.236061466 -0.157468955
##
## -0.069504444 0.027819380 0.114803535 0.193003706 0.263675602 0.311678899
##
## 0.344236409 0.388227780 0.422738872 0.437107720 0.453507109 0.456825228
##
## 0.451760914 0.447376194 0.422814374 0.411388500 0.379036597 0.345672394
##
## 0.318597303 0.277060370 0.251240113 0.204846468 0.177338678 0.125741499
##
## 0.096585500 0.051920926 0.026363658 -0.012459805 -0.034615943 -0.064482483
##
## -0.085238223 -0.104976129 -0.121702070 -0.133207549 -0.146374657 -0.152645918
##
## -0.158485361 -0.160097110 -0.159528602 -0.155338056 -0.152965723 -0.142122806
##
## -0.137396037 -0.123990531 -0.112850294 -0.102618031 -0.084733697 -0.074447079
##
## -0.053796475 -0.047422448 -0.028613093 -0.022934976 -0.009445066 -0.002058510
##
## 0.011011963 0.017155835 0.028954159 0.036011198 0.041574911 0.046800279
##
## 0.050610078 0.053188829 0.055399256 0.056137651
plot(x = x,y = y,type="l",xlab = "x",ylab = 'dy/dx',main='dy/dx vs x Ejercicio 104',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='black',ylim=c(-1.5,3))
lines(x[-1],dy_dx_atras[-1],type="l",col='orange')
#Aquí utilizamos el término [-1] indicando que la gráfica comienza en el segundo punto de la secuencia, debido a que la derivada hacia atrás solo puede calcularse desde el segundo punto hasta el último
lines(x[-length(x)],dy_dx_adelante,type="l",col='blue')
#Aquí utilizamos el término [-length(x)] indicando que la gráfica comienza en el primer punto de la secuencia pero debe elminar el último, debido a que la derivada hacia adelante solo puede calcularse hasta el penúltimo punto de la secuencia.
lines(x[-length(x)], dy_dx_central, type="l", col='red')
#Al trazar una línea horizontal podemos explorar el punto en el que los datos tienen une pendiente de cero
legend("topright", legend=c("Datos","deriv.atras","deriv. adelante","deriv. central"), fill=c("black","orange","blue","red"))
#No se logró retirar el primer punto de la gráfica, se podría crear un dataframe para eliminar los datos, pero de esta manera podemos ver que la derivada central tiene un comportamiento "medio" entre las otras y dos, el error aditivo aleatorio genera este comportamiento un poco oscilatorio que es más notorio en la derivada central.
\[\frac{d\left (\frac{p_{2}}{p_{1}}\right)}{dx}=\frac{\frac{d}{dx}(5x^2-6x+8)(3x^2+7)-(5x^2-6x+8)\frac{d}{dx}(3x^2+7)}{\left(3x^2+7\right)^{2}}\] \[\frac{d\left (\frac{p_{2}}{p_{1}}\right)}{dx}=\frac{(10x-6)(3x^2+7)-(5x^2-6x+8)(6x)}{\left(3x^2+7\right)^{2}}\]
\[\frac{d\left (\frac{p_{2}}{p_{1}}\right)}{dx}=\frac{(30x^3+70x-18x^2-42)-(30x^3-36x^2+48x)}{\left(3x^2+7\right)^{2}}\] \[\frac{d\left (\frac{p_{2}}{p_{1}}\right)}{dx}=\frac{(18x^2+12x-42)}{\left(9x^4+42x^2+49\right)}\]f=function(x,y) {-x^2 + 2*x*y + 3*y^2}
#Creamos 2 vectores de x y y
x=seq(-100,100, length.out = 100)
y=seq(-100,100, length.out = 100)
#Creamos una matriz de valores para x y y con la función 'expand grid' y obtenemos todas las combinaciones posibles de x y y en el rango establecido
xy=expand.grid(x, y)
#Ahora con la función apply calculamos los valores de la función, el valor 1 indica que se aplicara a cada fila de la matriz xy, la función que se aplica en los valores x[1] y x[2] de cada fila
z=apply(xy, 1, function(x) f(x[1], x[2]))
#Creamos un vector de colores para diferenciar un poco más los niveles de la gráfica de contorno
colors=c("red", "red", "red", "red", "red", "green", "green", "green", "green", "blue", "blue", "blue", "blue", "blue", "blue")
#En la función contorno, creamos una matriz con el argumento "matrix", que es el que permite organizar los datos de manera bidimensional.
contour(x, y, matrix(z, ncol = length(x)), xlab = "x", ylab = "y",col=colors,main='Gráfica de contorno y vectores gradiente para f(x.y) -x^2 + 2*x*y + 3*y^2 Ejercicio 106',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5)
#El gradiente es un vector que muestra el sentido de máxima variación de una función, por lo que en este caso se representa por el valor de las derivadas parciales en x y y, para esto utilizamos la librería pracma
library(pracma)
## Warning: package 'pracma' was built under R version 4.2.3
##
## Attaching package: 'pracma'
## The following objects are masked from 'package:matlib':
##
## angle, inv
#Y aquí con la función grad, obtenemos las derivadas parciales en x y y de la función.
grad=cbind(-2*x + 2*y, 2*x + 6*y)
#La función quiver permite colocar los vectores con base en la función grad dada anteriormente, y permite mostrar la dirección de máxima variación de la función, la cual coincide con los valores más cercanos de z.
quiver(x, y, grad[,1], grad[,2], length = 0.1, col = "blue")
Para resolver este ejercicio primero debemos obtener una la temperatura del objeto com función del tiempo, por tanto procedemos a encontrar esta función.
\[10\frac{dT}{dt}+T=T_{b}\] Aquí procedemos a ubicar los términos correspondientes a T y t en lados contrarios \[\frac{10}{T_{b}-T}dT=dt\] Ahora integramos para encontrar la función \[\int \frac{10}{T_{b}-T}dT=\int{dt}\] Para poder hacer esta integral utilizamos la técnica de cambio de variable para integrales logarítmicas \[u=T_{b}-T\] \[\frac{du}{dT}=-1\] \[du=-dT\] Ahora después de hacer el cambio de variable tenemos \[\int \frac{10}{u}dT=\int{dt}\] \[\int -\frac{10}{u}du=\int{dt}\] \[-10\int \frac{1}{u}du=\int{dt}\] \[-10\,ln\left|u\right|=t+C\] \[-10\,ln\left|T_{b}-T\right|=t+C\] Donde C es una constante de integración que se debe determinar \[ln\left|T_{b}-T\right|=-\frac{t+C}{10}\] \[\left|T_{b}-T\right|=e^{-\left(\frac{t+C}{10}\right )}\] Como Tb siempre será >= a T se puede omitir el valor absoluto, y C/10 se cambia por C ya que sigue siendo una constante. \[T_{b}-T=e^{-\left(\frac{t}{10}+C\right )}\]
Y aquí entonces obtenemos la función \[T=T_{b}-e^{-\left(\frac{t}{10}+C\right )}\] Para hallar C usamos los datos suministrados en el problema con Tb=170, T=70 en t=0 \[70=170-e^{-\left(\frac{0}{10}+C\right )}\] \[e^{-C}=100\]
\[-C=ln(100)=4.605\] \[C=-4.605\] Por tanto la ecuación sería en este caso \[T=170-e^{-\left(\frac{t}{10}-4.605\right)}\]
t en funciòn de T sería \[ln(170-T)=-\left(\frac{t}{10}-4.605\right)\] \[10\,\left(-ln(170-T)+4.605\right)=t\]
#b.
t=function(T){10*(-log(170-T)+4.605)}
print(t(168))
## [1] 39.11853
#De acuerdo a esto, aunque el problema no menciona las unidades de t para la función, se alcanzaría 168°F en un tiempo de 39.1 (u.t) (unidades de tiempo)
#c.
#Para la gráfica elegimos una escala de tiempo entre 0 y 100 u.t. (unidades de tiempo)
t=seq(0,100,by=0.1)
Temp=function(t){170-exp(-(t/10-4.605))}
#Con la función sapply, se da un valor de Temp para cada tiempo t
y=sapply(t,Temp)
plot(x = t,y = y,type='l',xlab = "t(u.t.)",ylab = 'T(°F)',main='Temperatura vs tiempo de un objeto sumergido en baño a 170°F',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex.sub=0.5,col='black')
abline(h = 170,lty=2,col = "gray80")
legend("bottomright", legend = c("T de 170°F"), lty = c(2), col = c("gray80"), bty = "o")
Se entiende que A es el área transversal del orificio de salida. a. Para poder apriximar cuanto tiempo tarda el tanque en desocuparse, necesitamos encontrar la ecuación de altura (h) en función del tiempo (t), por lo que debemos integrar la ecuación. \[\pi (2rh-h^{2})\frac{dh}{dt}=-C_{d}A\sqrt{2gh}\] \[\frac{\pi (2rh-h^{2})}{\sqrt{2gh}}dh=-C_{d}A\, dt\] \[\int \frac{\pi (2rh-h^{2})}{\sqrt{2gh}}dh=\int -C_{d}A\, dt\] \[\frac{\pi}{\sqrt{2g}}\int \frac{ (2rh-h^{2})}{\sqrt{h}}dh=\int -C_{d}A\, dt\] \[\frac{\pi}{\sqrt{2g}}\int \left ( \frac{ (2rh)}{\sqrt{h}}-\frac{ (h^{2})}{\sqrt{h}}\right )dh=\int -C_{d}A\, dt \] \[\frac{\pi}{\sqrt{2g}}\int \left ( 2rh^{\left ( 0.5 \right )}-h^{1.5}\right )dh=\int -C_{d}A\, dt\] \[\frac{\pi}{\sqrt{2g}}\left (2r\frac{h^{\left ( 1.5 \right )}}{1.5}-\frac{h^{2.5}}{2.5}\right )=-C_{d}At+C \] \[4r\frac{h^{\left ( 1.5 \right )}}{3}-\frac{2h^{2.5}}{5}=\frac{\sqrt{2g}(-C_{d}At+C)}{\pi} \] \[4r\frac{h^{\left ( 1.5 \right )}}{3}-\frac{2h^{2.5}}{5}=\frac{-\sqrt{2g}C_{d}At}{\pi}+\frac{\sqrt{2g}C}{\pi} \] \[4r\frac{h^{\left ( 1.5 \right )}}{3}-\frac{2h^{2.5}}{5}=\frac{-\sqrt{2g}C_{d}At}{\pi}+C_{2}\]
Para encontrar el valor de C2 reemplazamos los valores por la altura en el tiempo 0 \[4(3m)\frac{(5m)^{\left ( 1.5 \right )}}{3}-\frac{2(5m)^{2.5}}{5}=\frac{-\sqrt{2g}C_{d}A(0s)}{\pi}+C_{2} \] \[12m\frac{(5m)^{\left ( 1.5 \right )}}{3}-\frac{2(5m)^{2.5}}{5}=C_{2} \] \[C_{2}=22.36 m^{2.5}\] \[C_{2}=\frac{\sqrt{2g}C}{\pi}\] \[C=\frac{C_{2}\pi}{\sqrt{2g}}\] \[C=\frac{22.36m^{2.5}\pi}{\sqrt{2\times 9.81m/s^2}}\] \[C=15.86\,m^{2}s\] Ahora reemplazamos la constante C2 en la acuación
\[4r\frac{h^{\left ( 1.5 \right )}}{3}-\frac{2h^{2.5}}{5}=\frac{-\sqrt{2g}C_{d}At}{\pi}+22.36m^{2.5}\] Para encontrar el tiempo en el que el tanque se desocupa reemplazamos h por 0: \[0=\frac{-\sqrt{2\times 9.81 m/s^2}\times0.5\times (\pi\times(0.02m)^2)\,t}{\pi}+22.36m^{2.5}\] \[22.36m^{2.5}=\frac{\sqrt{2\times 9.81 m/s^2}\times0.5\times (\pi\times(0.02m)^2)\,t}{\pi}\] \[t=\frac{22.36m^{2.5}\ \times\pi}{\sqrt{2\times 9.81 m/s^2}\,\times 0.5\times (\pi\times(0.02m)^2)}\] \[t=25240s\] Equivale a 421 minutos aproximadamente (7 horas) \[4r\frac{h^{\left ( 1.5 \right )}}{3}-\frac{2h^{2.5}}{5}=\frac{-\sqrt{2g}C_{d}At}{\pi}+22.36m^{2.5}\]
#Ahora en R
#a.
r=3
g=9.81
Cd=0.5
rd=0.02
t_llenado=function(h){(((4*r*(h^(1.5))/3)-(2/5*(h^(2.5)))-22.36)*pi)/(-(2*g)^(0.5)*Cd*pi*rd^2)}
t_llenado(0)
## [1] 25240.17
#El tanque se desocupará en 25240 segundos.
#b
#Ahora graficaremos la altura en función del tiempo
#curve(t_llenado,from = 0,to = 5)
h_vals = seq(0, 5, length.out = 100)
t_vals = sapply(h_vals,t_llenado)
df = data.frame(h = h_vals, t = t_vals)
df
## h t
## 1 0.00000000 25240.1715311
## 2 0.05050505 25189.1816799
## 3 0.10101010 25096.6825396
## 4 0.15151515 24977.9103505
## 5 0.20202020 24838.4646830
## 6 0.25252525 24681.6629708
## 7 0.30303030 24509.7969962
## 8 0.35353535 24324.5880559
## 9 0.40404040 24127.3997599
## 10 0.45454545 23919.3532027
## 11 0.50505051 23701.3956598
## 12 0.55555556 23474.3445322
## 13 0.60606061 23238.9169836
## 14 0.65656566 22995.7507801
## 15 0.70707071 22745.4194418
## 16 0.75757576 22488.4435703
## 17 0.80808081 22225.2995167
## 18 0.85858586 21956.4261462
## 19 0.90909091 21682.2302101
## 20 0.95959596 21403.0906748
## 21 1.01010101 21119.3622568
## 22 1.06060606 20831.3783435
## 23 1.11111111 20539.4534309
## 24 1.16161616 20243.8851782
## 25 1.21212121 19944.9561530
## 26 1.26262626 19642.9353264
## 27 1.31313131 19338.0793614
## 28 1.36363636 19030.6337316
## 29 1.41414141 18720.8336973
## 30 1.46464646 18408.9051620
## 31 1.51515152 18095.0654288
## 32 1.56565657 17779.5238693
## 33 1.61616162 17462.4825199
## 34 1.66666667 17144.1366146
## 35 1.71717172 16824.6750622
## 36 1.76767677 16504.2808765
## 37 1.81818182 16183.1315646
## 38 1.86868687 15861.3994787
## 39 1.91919192 15539.2521356
## 40 1.96969697 15216.8525078
## 41 2.02020202 14894.3592901
## 42 2.07070707 14571.9271423
## 43 2.12121212 14249.7069135
## 44 2.17171717 13927.8458478
## 45 2.22222222 13606.4877733
## 46 2.27272727 13285.7732783
## 47 2.32323232 12965.8398727
## 48 2.37373737 12646.8221386
## 49 2.42424242 12328.8518701
## 50 2.47474747 12012.0582033
## 51 2.52525253 11696.5677372
## 52 2.57575758 11382.5046471
## 53 2.62626263 11069.9907899
## 54 2.67676768 10759.1458035
## 55 2.72727273 10450.0871992
## 56 2.77777778 10142.9304488
## 57 2.82828283 9837.7890666
## 58 2.87878788 9534.7746862
## 59 2.92929293 9233.9971332
## 60 2.97979798 8935.5644934
## 61 3.03030303 8639.5831774
## 62 3.08080808 8346.1579817
## 63 3.13131313 8055.3921466
## 64 3.18181818 7767.3874103
## 65 3.23232323 7482.2440614
## 66 3.28282828 7200.0609876
## 67 3.33333333 6920.9357226
## 68 3.38383838 6644.9644904
## 69 3.43434343 6372.2422475
## 70 3.48484848 6102.8627231
## 71 3.53535354 5836.9184572
## 72 3.58585859 5574.5008374
## 73 3.63636364 5315.7001336
## 74 3.68686869 5060.6055308
## 75 3.73737374 4809.3051618
## 76 3.78787879 4561.8861365
## 77 3.83838384 4318.4345717
## 78 3.88888889 4079.0356185
## 79 3.93939394 3843.7734893
## 80 3.98989899 3612.7314829
## 81 4.04040404 3385.9920092
## 82 4.09090909 3163.6366128
## 83 4.14141414 2945.7459953
## 84 4.19191919 2732.4000373
## 85 4.24242424 2523.6778188
## 86 4.29292929 2319.6576397
## 87 4.34343434 2120.4170385
## 88 4.39393939 1926.0328115
## 89 4.44444444 1736.5810302
## 90 4.49494949 1552.1370587
## 91 4.54545455 1372.7755700
## 92 4.59595960 1198.5705626
## 93 4.64646465 1029.5953753
## 94 4.69696970 865.9227023
## 95 4.74747475 707.6246080
## 96 4.79797980 554.7725403
## 97 4.84848485 407.4373442
## 98 4.89898990 265.6892749
## 99 4.94949495 129.5980105
## 100 5.00000000 -0.7673362
plot(df$t, df$h, xlab = "Tiempo (s)", ylab = "h(m)",main='Altura de agua de un tanque (h) en función del tiempo',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex=0.5,cex.sub=0.5,col='blue',pch=16)
\[\frac{dy}{dt}+\frac{2}{10+2t}y=4\] Esta ecuación es de la forma
\[\frac{dy}{dt}+P(x)y=f(x)\] Que se resuelve mediante factor integrante y es de la forma
\[y=y_{c}+y_{p}\] yc es la solución del sistema homogéneo que es de la forma.
\[y_{c}=Ce^{-\int P(x)dx}\]
\[y_{c}=Ce^{-\int \frac{2}{10+2t}dt}\] Aquí aplicamos la técnica del cambio de variable
\[u=10+2t\] \[\frac{du}{dt}=2\] \[dt=\frac{du}{2}\] \[\int \frac{2}{10+2t}dt\] \[\int \frac{2}{u}\frac{du}{2}\]
\[\int \frac{1}{u}du=ln|u|\]
Reemplazamos la variable u
\[ln|10+2t|\] Ahora volvemos a la solución del sistema homogéneo
\[y_{c}=Ce^{-ln(10+2t)}\] \[y_{c}=-C(10+2t)\] Como C es una constante podemos factorizar sacando el 2 como factor común y tener una nueva C \[y_{c}=-C(5+t)\] Ahora calculamos la solución del sistema no homogéneo que es de la forma
\[y_{p}=\frac{1}{e^{\int P(x)dx}}\int e^{\int P(x)dx}f(x)dx\] \[y_{p}=\frac{1}{(10+2t)}\int (10+2t)4dt\] \[y_{p}=\frac{1}{(10+2t)}\int (40+8t)dt\] \[y_{p}=\frac{1}{(10+2t)}(40t+\frac{8t^{2}}{2})\] \[y_{p}=\frac{1}{(10+2t)}(40t+4t^{2})\] \[y_{p}=\frac{4(10t+t^{2})}{2(5+t)}\] \[y_{p}=\frac{2(10t+t^{2})}{(5+t)}\] Ahora obtenemos la solución de la ecuación diferencial \[y=-C(5+t)+\frac{2(10t+t^{2})}{(5+t)}\] Calculamos el factor C en el tiempo 0 donde y = 0 \[0=-C(5)+\frac{0}{(5)}\] \[0=-C(5)\] Por tanto \[C=0\] La ecuación es \[y=\frac{2(10t+t^{2})}{(5+t)}\]
#Ahora en R
y=function(t){(2*(10*t+t^2))/(5+t)}
t = seq(0, 10, length.out = 100)
y_app = sapply(t,y)
df = data.frame(y = y_app, t = t)
df
## y t
## 1 0.0000000 0.0000000
## 2 0.4000400 0.1010101
## 3 0.7923899 0.2020202
## 4 1.1774892 0.3030303
## 5 1.5557444 0.4040404
## 6 1.9275322 0.5050505
## 7 2.2932023 0.6060606
## 8 2.6530795 0.7070707
## 9 3.0074660 0.8080808
## 10 3.3566434 0.9090909
## 11 3.7008743 1.0101010
## 12 4.0404040 1.1111111
## 13 4.3754619 1.2121212
## 14 4.7062626 1.3131313
## 15 5.0330072 1.4141414
## 16 5.3558844 1.5151515
## 17 5.6750713 1.6161616
## 18 5.9907344 1.7171717
## 19 6.3030303 1.8181818
## 20 6.6121065 1.9191919
## 21 6.9181019 2.0202020
## 22 7.2211476 2.1212121
## 23 7.5213675 2.2222222
## 24 7.8188784 2.3232323
## 25 8.1137910 2.4242424
## 26 8.4062097 2.5252525
## 27 8.6962339 2.6262626
## 28 8.9839572 2.7272727
## 29 9.2694689 2.8282828
## 30 9.5528534 2.9292929
## 31 9.8341910 3.0303030
## 32 10.1135579 3.1313131
## 33 10.3910268 3.2323232
## 34 10.6666667 3.3333333
## 35 10.9405432 3.4343434
## 36 11.2127189 3.5353535
## 37 11.4832536 3.6363636
## 38 11.7522041 3.7373737
## 39 12.0196248 3.8383838
## 40 12.2855675 3.9393939
## 41 12.5500818 4.0404040
## 42 12.8132150 4.1414141
## 43 13.0750124 4.2424242
## 44 13.3355173 4.3434343
## 45 13.5947712 4.4444444
## 46 13.8528139 4.5454545
## 47 14.1096832 4.6464646
## 48 14.3654158 4.7474747
## 49 14.6200466 4.8484848
## 50 14.8736092 4.9494949
## 51 15.1261357 5.0505051
## 52 15.3776572 5.1515152
## 53 15.6282032 5.2525253
## 54 15.8778024 5.3535354
## 55 16.1264822 5.4545455
## 56 16.3742690 5.5555556
## 57 16.6211882 5.6565657
## 58 16.8672642 5.7575758
## 59 17.1125206 5.8585859
## 60 17.3569799 5.9595960
## 61 17.6006642 6.0606061
## 62 17.8435943 6.1616162
## 63 18.0857906 6.2626263
## 64 18.3272727 6.3636364
## 65 18.5680594 6.4646465
## 66 18.8081690 6.5656566
## 67 19.0476190 6.6666667
## 68 19.2864265 6.7676768
## 69 19.5246078 6.8686869
## 70 19.7621787 6.9696970
## 71 19.9991547 7.0707071
## 72 20.2355505 7.1717172
## 73 20.4713805 7.2727273
## 74 20.7066584 7.3737374
## 75 20.9413978 7.4747475
## 76 21.1756115 7.5757576
## 77 21.4093122 7.6767677
## 78 21.6425121 7.7777778
## 79 21.8752228 7.8787879
## 80 22.1074559 7.9797980
## 81 22.3392223 8.0808081
## 82 22.5705329 8.1818182
## 83 22.8013980 8.2828283
## 84 23.0318277 8.3838384
## 85 23.2618318 8.4848485
## 86 23.4914198 8.5858586
## 87 23.7206008 8.6868687
## 88 23.9493839 8.7878788
## 89 24.1777778 8.8888889
## 90 24.4057908 8.9898990
## 91 24.6334311 9.0909091
## 92 24.8607067 9.1919192
## 93 25.0876254 9.2929293
## 94 25.3141946 9.3939394
## 95 25.5404216 9.4949495
## 96 25.7663137 9.5959596
## 97 25.9918775 9.6969697
## 98 26.2171200 9.7979798
## 99 26.4420476 9.8989899
## 100 26.6666667 10.0000000
plot(df$t, df$y, xlab = "Tiempo", ylab = "y(t)",main='Concentración de salmuera en función del tiempo',sub='Métodos Multivariados 2023-1S Sofía Arévalo, David Acosta, Juan Ramón Vilaseca',cex=0.5,cex.sub=0.5,col='blue',pch=16)