Diferencia entre revisiones de «Explotación Minera (G12-A)»

De MateWiki
Saltar a: navegación, buscar
(→‎Modelo de Runge Kutta (de cuarto orden) y Heun)
Línea 277: Línea 277:
 
[[Archivo:ProduccionP4.jpg|800px|thumb|left|Método de Runge Kutta (de 4 orden) y Heun]]
 
[[Archivo:ProduccionP4.jpg|800px|thumb|left|Método de Runge Kutta (de 4 orden) y Heun]]
 
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|left|Superposición de las gráficas aumentada]]
 
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|left|Superposición de las gráficas aumentada]]
 +
 +
 +
 +
 +
  
  
  
 
Tras superponer ambas gráficas se puede observar la gran similitud de los modelos. Sólo aumentando la gráfica hemos sido capaces de poder distinguir las dos curvas.
 
Tras superponer ambas gráficas se puede observar la gran similitud de los modelos. Sólo aumentando la gráfica hemos sido capaces de poder distinguir las dos curvas.
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
  
 
='''Cantidad extraída Q cuando lim t→∞'''=
 
='''Cantidad extraída Q cuando lim t→∞'''=
Línea 303: Línea 374:
  
 
[[Archivo:ProduccionP5.jpg|400px|thumb|left|Curva cuando t tiende a infinito]]
 
[[Archivo:ProduccionP5.jpg|400px|thumb|left|Curva cuando t tiende a infinito]]
−
 
  
  

Revisión del 21:56 6 mar 2015

Trabajo realizado por estudiantes
Título Explotación minera. Grupo 12-A
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Interpretación del problema

Debido a la alta demanda de un determinado mineral, se ha decidido explotar un yacimiento de una region estudiada. Los estudios han concluido que la cantidad total extraíble (K) de dicho mineral es de 10875 toneladas. Se estima un crecimiento muy rápido de la producción (toneladas/año) durante los 25 años, tras los cuales, a causa de dificultades técnicas y la caída de la demanda, decrecerá la producción lentamente. Para estudiar este problema vamos a tomar diferentes modelos matemáticos (aproximaciones numéricas computacionales a partir de los datos obtenidos en el trabajo de campo).

La relación entre la producción (P) y la cantidad extraída (Q) es una relación diferencial. P será la derivada de Q respecto del tiempo.


2 Modelo logístico de Gompertz

Un posible modelo que relaciona la producción con la cantidad extraída es el modelo logístico de Gompertz, basado en la siguiente ecuación::

[math]P(Q) = \frac{dQ}{dt} = rQ\log\left(\frac{K}{Q}\right) [/math]

Tras estudios previos se obtuvo una produccion máxima de 240 toneladas/año (máximo relativo de la función). Derivando nuestra ecuación respecto de Q e igualandola a 0 obtenemos la siguiente ecuación::

[math]P'=0=rlog(\frac{K}{Q})-r=r(log(\frac{K}{Q})-1)[/math]


Despejamos Q y obtenemos::


[math]Q = \frac{K}{e}[/math]


Introduciendo el valor de Q obtenido en los estudios, y el valor de P=240, despejamos la ecuación y obtenemos el coeficiente r::

[math]240=r\frac{K}{e}log(\frac{K}{\frac{K}{e}})=r\frac{K}{e} → r=\frac{240e}{K}=0.0599[/math]


  • Modelo computacional de Gompertz en MATLAB :


k=10875;                      %Cantidad total extraible en toneladas
Q=0:1:10875;                  %Vector con la cantidad de toneladas extraídas
n=length(Q);                  %Tamaño del vector Q
P=zeros(1,n);                 %Vector de ceros de una fila y N columnas
r=240*exp(1)/10785;           %coeficiente r

for i=1:n                     %Realizamos el bucle 
    P(i)=r*Q(i)*log(k/Q(i));  %Definimos la funcion P(Q) 
end 

plot(Q,P,'k')                     %Gráfica de Q (abcisas) y P (ordenadas) en color negro.
xlabel('cantidad (ton)')       
ylabel('produccion (ton/año)')


Curva de la función P(Q




Analizando la gráfica obtenida por MATLAB, podemos observar que la pendiente (en valor absoluto) inicial de la curva es mayor a la del final, como indicaba el estudio previo. La curva muestra un cambio de pendiente en el valor de Q=240 toneladas (máximo) y un fin de producción en el valor de cantidad total extraída de 10875 toneladas.











3 Modelo logístico de Verhulst

Otro posible modelo logístico es el de Verhulst, definido por la siguiente ecuación::

[math]Q'=rQ(1-\frac{Q}{k})[/math]

Para obtener el nuevo coeficiente r, procedemos de la misma forma que en el modelo de Gompertz::

[math]P'=0=r(1-\frac{2Q}{K})[/math]

Despejamos Q y obtenemos::

[math]\frac{2Q}{K}=1→ Q=\frac{K}{2}=5437.5[/math]

Introduciendo el valor de Q obtenido en la ecuación, y el valor de P=240, obtenemos el coeficiente r::

[math]240=r\frac{K}{2}(1-\frac{\frac{K}{2}}{K})=r\frac{K}{4} → r=\frac{960}{K}=0.088[/math]

  • Modelo computacional de Verhulst en MATLAB :

Con los datos obtenidos creamos un programa en MATLAB para obtener una gráfica del modelo de Verhulst y poder compararlo con el de Gompertz: 

k=10875; %Cantidad maxima extraible
rG=240*exp(1)/k; %Tasa intrinseca de creciemiento (Gompertz)
rV=960/k; %Tasa intrinseca de creciemiento (Verhulst)
Q=0:1:k; 
n=length(Q);
for i=1:n 
    PGom(i)=rG*Q(i)*log(k/Q(i)); %Gompertz
    PVer(i)=rV*Q(i)*(1-Q(i)/K); %Verhulst
end                                     
subplot(1,2,1)
plot(Q,PVer,'g') %Gráfica modelo de Verhulst
xlabel('cantidad (ton)')         
ylabel('produccion (ton/año)')  
subplot(1,2,2)
plot(Q,PGom,'r') %Gráfica modelo de Gompertz
xlabel('cantidad')             
ylabel('produccion')          
hold on %Superponemos las dos gráficas
plot(Q,PVer,'g') %Gráfica modelo de Verhulst
xlabel('cantidad (ton)')        
ylabel('produccion (ton/año)')  
legend('Modelo Gompertz','Modelo Verhulst','Location','best') 
hold off


Curva del modelo de Verhulst en comparación con la gráfica de Gompertz




La primera curva muestra el gráfico según el modelo logístico de Verhulst, el cual a primera vista parece aproximarse bastante bien al previo estudio de la producción P. Sin embargo al observar la segunda gráfica, en la cual se superponen las curvas de ambos modelos (modelo de Gompertz y modelo de Verhulst) se puede indicar claramente que en este caso el modelo de Gompertz se aproxima mejor a nuestro estudio de producción, indentificable por el elevado crecimiento inicial de dicha producción.














4 Método de Euler

k=10875; %Cantidad máxima extraible
t0=0;
h=1/12; %El paso es de 1 mes por año
r=240*exp(1)/k; %Tasa intrínseca de crecimiento
t(1)=t0;
Q(1)=0.1;
i=1;
while 1
    Q(i+1)=Q(i)+h*r*Q(i)*log(k/Q(i));
    t(i+1)=t(i)+h;
    if i>1&&abs((r*Q(i)*log(k/Q(i)))-25)<0.1&&abs((r*Q(i-1)*log(k/Q(i-1))))>abs((r*Q(i)*log(k/Q(i))));
 %La producción deja de ser rentable cuando es menor de 25 toneladas
        break
    end
    i=i+1;
end
t=t';
Q=Q';
plot(t,Q) %Gráfica de la producción en función del tiempo


Método de Euler





















5 Modelo de Runge Kutta (de cuarto orden) y Heun

k=10875; %Cantidad máxima extraible
t0=0; %Tiempo inicial = 0
Q0=0.1; %Cantidad inicial extraida = 0
h=1/12; %El paso es de 1 mes por año
t=t0;
r=240*exp(1)/k; %Tasa intrínseca de crecimiento                       
y(1)=Q0; %RK4
z(1)=Q0; %Heun
i=1;
while 1 %RK4
K1=r*y(i)*log(k/y(i)); 
K2=r*(y(i)+1/2*K1*h)*log(k/(y(i)+1/2*K1*h));
K3=r*(y(i)+1/2*K2*h)*log(k/(y(i)+1/2*K2*h));
K4=r*(y(i)+K3*h)*log(k/(y(i)+K3*h));
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
     t(i+1)=t(i)+h;
    if i>1&&abs((r*y(i)*log(k/y(i)))-25)<0.1&r*y(i-1)*log(k/y(i-1))>r*y(i)*log(k/y(i))
break
    end
   i=i+1;
end
i=1;
while 1  %Heun
k1=r*z(i)*log(k/z(i));
k2=r*(z(i)+k1*h)*log(k/(z(i)+k1*h));
z(i+1)=z(i)+(h/2)*(k1+k2);
     t(i+1)=t(i)+h;
    if i>1&&abs((r*z(i)*log(k/z(i)))-25)<0.1&r*z(i-1)*log(k/z(i-1))>r*z(i)*log(k/z(i))
break
    end
   i=i+1;
end
subplot(1,3,1)
plot(t,y,'r'); %Gráfica RK4
legend('RK4','Location','best')
subplot(1,3,2)
plot(t,z,'y'); %Gráfica Heun
legend('Heun','Location','best') 
subplot(1,3,3)
hold on %Superponemos las dos gráficas
plot(t,y,'r'); %Gráfica RK4                                     
plot(t,z,'g'); %Gráfica Heun
legend('RK4','Heun','Location','best') 
hold off


Método de Runge Kutta (de 4 orden) y Heun
Superposición de las gráficas aumentada





Tras superponer ambas gráficas se puede observar la gran similitud de los modelos. Sólo aumentando la gráfica hemos sido capaces de poder distinguir las dos curvas.


































6 Cantidad extraída Q cuando lim t→∞

k=10875; %Cantidad maxima extraible
t=0:1:250; %Tiempo desde 0 hasta 250 (Elegido por nosotros)
N=length(t);
Q=zeros(1,N);
r=240*exp(1)/k; %Tasa intrínseca de crecimiento
Q0=0.1; 
Q(1)=Q0;
Q=k*exp(exp(-r*t)*(log(Q0/k))); %Ecuación de Gompertz
for i=1:N 
    Q(i)=k*exp(exp(-r*t(i))*(log(Q0/k))); %Gompertz
end 
Q(250)                                               
plot(t,Q) %Gráfica de la producción en funcion del tiempo
xlabel('Tiempo(años)')                   
ylabel('Cantidad (ton)')


Curva cuando t tiende a infinito















7 Función de prodicción P(t)

k=10875; %Cantidad máxima extraible       
t0=0; %Tiempo inicial = 0
Q0=0.1; %Cantidad inicial = 0          
h=1/12; %El paso es de 1 mes por año
t=t0;
r=240*exp(1)/k; %Tasa intrínseca de crecimiento           
Q(1)=Q0;
i=1;
while 1 
   P(i)=r*Q(i)*log(k/Q(i)); %Ecuación de Gompertz
    if i>1&&abs((r*Q(i)*log(k/Q(i)))-25)<0.1&r*Q(i-1)*log(k/Q(i-1))>r*Q(i)*log(k/Q(i))              
break
    end
    Q(i+1)=Q(i)+h*(r*Q(i)*log(k/Q(i))); %Euler
     t(i+1)=t(i)+h; 
   i=i+1;          
end
    [max_val,tiem]=max(P) %Maximo de la función y su posición
    tiem=tiem/12
 
plot(t,P,'g') %Gráfica de la producción en funcion del tiempo
xlabel('tiempo (años)')
ylabel('produccion (ton/año)')


Gráfica de la Producción en función del tiempo





















8 Cantidad de mineral sin extraer

k=10875; %Cantidad total extraible
t0=0; %Tiempo inicial = 0
Q0=0.1; %Cantidad inicial = 0
h=1/12; %El paso es de 1 mes por año
t=t0;
r=240*exp(1)/k; %Tasa intrínseca de crecimiento
Q(1)=Q0;
i=1;                
while 1              
     Q(i+1)=Q(i)+h*(r*Q(i)*log(k/Q(i))); %Euler
     t(i+1)=t(i)+h;
    if i>1&&abs((r*Q(i)*log(k/Q(i)))-25)<0.1&r*Q(i-1)*log(k/Q(i-1))>r*Q(i)*log(k/Q(i))
break               
    end
   i=i+1;            
end                 
n=length(Q);
MinEx=Q(n); %Cantidad de mineral extraida
MinSinEx=k-MinEx %Cantidad de mineral sin extraer
plot(t,Q,'g') %Gráfica de la producción en función del tiempo
legend('Euler','Location','best')


Cantidad de mineral sin extraer





















9 Apartados 7 y 8 con el método de logístico

k=10875; %Cantidad máxima extraible
t0=0; %Tiempo inicial = 0
Q0=0.1; %Cantidad inicial = 0
h=1/12; %El paso es de 1 mes por año
t=t0;
r=(240*4)/k; %Tasa intrínseca de crecimiento
Q(1)=Q0;
i=1;
while 1 
   Q(i+1)=Q(i)+h*(r*Q(i)*log(k/Q(i))); %Euler
   P(i)=r*Q(i)*(1-Q(i)/k); %Verhulst
    if i>1&&abs(r*Q(i)*(1-Q(i)/k)-25)<0.1&r*Q(i-1)*(1-Q(i-1)/k)>r*Q(i)*(1-Q(i)/k)
break                   
    end
    k1=r*Q(i)*(1-Q(i)/k);
    k2=r*(Q(i)+k1*h)*(1-(Q(i)+k1*h)/k);
    Q(i+1)=Q(i)+(h/2)*(k1+k2); %Heun
    t(i+1)=t(i)+h;
   i=i+1;       
end
[max_val,tiem]=max(P) %Maximo de la función y su posición
tiem=tiem/12 
n=length(Q); 
MinEx=Q(n); %Cantidad de mineral extraida
MinSinEx=k-MinEx %Cantidad de mineral sin extraer
plot(t,P,'g') %Grafica de la producción en función del tiempo
legend('Heun','Location','best')
xlabel('tiempo')
ylabel('produccion (tn/año)')


Método logístico
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Explotación_Minera_(G12-A)&oldid=28903»