Explotación minera (Grupo21-C)

Trabajo realizado por estudiantes
Título	Explotación minera. Grupo 21-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Jesús Infestas Robles, Pablo Medina Higueras , Alejandro Perales Juidías, Jaime Delage Ramírez, Mairena Pérez López
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

• La idea principal se basa en el modelo de Gompertz ,el cual muestra las tasas de crecimiento en un periodo de tiempo y su disminución de forma exponencial con el paso del mismo.Aplicado a nuestro caso se trata de una explotación de un yacimiento mineral, en un principio se decide explotar por su posible rentabilidad,el cual en un periodo de tiempo dado -25 años- será durante el cual se extrae la cantidad mayor de mineral y a partir de dicho tiempo,empecerá a descender esta producción, las causas son ajenas a nuestro estudio, pero seguramente ya no exista el mismo rendimiento.
  
• Con lo cual nuestro problema muestra una ecuación diferencial que sigue esta forma:
  

[math] \frac{dQ}{dt}=rQlog(\frac{K}{Q}) [/math]

1.1 Relación de cantidad y produccion

[math] P(Q) = \frac{dQ}{dt}=rQlog(\frac{K}{Q}) [/math]

• 1.Necesitamos saber para que valor Q alcanza el máximo, para ello derivamos nuestra función respecto a Q :

[math] P'(Q) =r(log(\frac{K}{Q})+Q\frac{\frac{-K}{Q^2}}{\frac{K}{Q}})=r(log(\frac{K}{Q})-1)⇐⇒P′(Q) = 0⇐⇒log(\frac{K}{Q})=1⇐⇒(\frac{K}{Q})=e⇐⇒Q=(\frac{K}{e}) [/math]

• 2. Tenemos en cuenta que nuestra producción máxima son 240 toneladas/año, entonces :

• 3. Y así obtendremos nuestra tasa intrínseca de crecimiento(r) para nuestra función Gompertz :

[math] 240=r\frac{K}{e}⇐⇒r=\frac{240e}{K}=\frac{240e}{10875} [/math]

1.2 Relación de la producción y el volumen de toneladas extraídas modelo de Gompertz

• Realizamos el programa teniendo en cuenta nuestra tasa intrínseca de crecimiento anteriormente hallada:

r=240*exp(1)/10785;           %tasa intrínseca de crecimiento
k=10875;                      %la cantidad total (toneladas) extraible
Q=0:1:10875;                  %vector con la cantidad de toneladas desde 0 hasta el maximo que se extraen
N=length(Q);                  %Tamaño vector Q
P=zeros(1,N);                 %vector de ceros deuna fila y N columnas
for i=1:N                     %realizo el bucle 
    P(i)=r*Q(i)*log(k/Q(i));  %Defino la funcion P(Q) 
end 
plot(Q,P)                     %Dibujo mi grafica como una curva con Q(abcisas) y P(ordenadas)
xlabel('cantidad (tn)')       %añado un pequeño titulo a mi Q y P
ylabel('produccion (tn/año)')

• Análisis de nuestra curva:

LEYENDA:GRÁFICO FUNCIÓN P(Q)

1.Como se puede apreciar,nuestra función realiza una curva que transcurre desde 0 hasta la cantidad máxima de toneladas que se pueden extraer.

2.Además;la curva que describre nuestra función alcanza su máximo relativo en 240 toneladas por año, que es la producción máxima de nuestra función.

1.3 Relación de la producción y el volumen de toneladas extraídas modelo de Verhulst

• Realizamos nuestro programa con nuestra nueva función y tasa intrínseca de crecimiento para el modelo de Verhulst :

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

r1=240*exp(1)/10875;           %Tasa intrinseca de creciemiento modelo de Gompertz
r2=(240*4)/10875;              %Tasa intrinseca de creciemiento modelo de Verhulst
K=10875;                       %Cantidad maxima de toneladas extraibles
Q=0:1:10875;                   %Vector con la cantidad en toneladas desde 0 hasta K
N=length(Q);                   %Tamaño de nuestro vector Q
for i=1:N                      %Realizamos un bucle desde uno hasta el numero de elementos(N) de nuestro vector Q
    PG(i)=r1*Q(i)*log(K/Q(i)); %Gompertz
    PV(i)=r2*Q(i)*(1-Q(i)/K);  %Verhulst
end                            %cerramos el bucle
subplot(1,3,1)
plot(Q,PG)                     %Dibujamos una curva en nuestro eje x(Q)y nuestro eje y(funcion de Gompertz)
xlabel('cantidad')             %añadimos titulo a las abcisas
ylabel('produccion')           %añadimos titulo a las ordenadas
subplot(1,3,2)
plot(Q,PV,'g')                 %dibujamos nuestra curva con eje x(Q) y eje y(funcion del modelo de Verhulst)
xlabel('cantidad (tn)')        %añadimos titulo a las abcisas d 
ylabel('produccion (tn/año)')  %añadimos titulo a las ordenadas
subplot(1,3,3)
plot(Q,PG)                     %Dibujamos una curva en nuestro eje x(Q)y nuestro eje y(funcion de Gompertz)
xlabel('cantidad')             %añadimos titulo a las abcisas
ylabel('produccion')           %añadimos titulo a las ordenadas
hold on                        %sobre nuestra curva anterior le superponemos la siguiente
plot(Q,PV,'g')                 %dibujamos nuestra curva con eje x(Q) y eje y(funcion del modelo de Verhulst)
xlabel('cantidad (tn)')        %añadimos titulo a las abcisas d 
ylabel('produccion (tn/año)')  %añadimos titulo a las ordenadas
legend('Modelo Gompertz','Modelo Verhulst','Location','best') 
hold off

• Análisis de nuestras gráficas :

Comparación de gráficos

• 1.En la primera subventana aparece la curva del modelo de Gompertz explicada en el apartado anterior
• 2.En la segunda subventana aparece la curva del modelo de Verhulst , de acuerdo con este podemos decir que ajusta bastante bien cuando está alejado del momento inicial
• 3.

1.4 Utilización del método de Euler,RK4 y Heun

• Previo

Euler es un método de aproximación de problemas de valor inicial,obtendremos la solución aproximación de soluciones en un intervalo definido.Si nuestra h es pequeña la aproximación será mejor, ya que la acumulación de errores en cada solución será menor.

1.5 Función Q(t), obtención por Euler y dibujo de su gráfica

• Aunque nuestra cantidad inicial es 0 ponemos 0.1 porque sino al hacer Euler saldría una indeterminación, y no se podría dibujar nuestra gráfica, así que hemos elegido ese valor pequeño para que funcione correctamente

t0=0;                                      %El tiempo inicial es cero
Q0=0.1;                                    %Al principio la cantidad de mineral extraida es cero
h=1/12;                                    %El enunciado dice que el paso es de un mes, y como nuestra unidad es el año, el paso será 1/12
                                           %Definimos la variable independiente pero como no conocemos el tiempo en el que se paraliza la excavacíón, 
                                           %por lo que del vector t solo definimos el primer elemento
t=t0;                                      %Primer elemento de t
r=240*exp(1)/10875;                        %El valor de r ha sido hallado en el apartado 2
K=10875;                                   %Cantidad total extraible
                                           %Al no conocer el tiempo final, tampoco podemos saber el valor de N, y por
                                           %lo tanto, el tamaño del vector Q, por lo que solo definimos el primer
                                           %elemento del vector
Q(1)=Q0;
i=1;                                       %el bucle empieza con el elemento 1
while 1                                    %Empezamos el bucle que rastrea todos los elementos hasta el buscado  
     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))  
                                           %Vemos cuando la producción deja de ser rentable, que será cuando la producción baje a 25 toneladas
      
break                                      %en ese momento salimos del bucle
    end                                    %cerramos el bucle
   i=i+1;                                  %para que nuestro bucle prosiga elemento a elemento hasta encontrar el de la condicion if
end


plot(t,Q,'r')                              %dibujamos la curva en rojo, tiempo en abcisas y cantidad en ordenadas
legend('Euler','Location','best')          %ponemos una leyenda nuestra curva

1.6 Utilización de Runge Kutta de cuarto orden y Heun con el apartado anterior

• Realizamos el mismo ejercicio que en el apartado anterior, pero esta vez, utilizando los métodos de Runge Kutta de cuarto orden y Heun

t0=0;                                      %El tiempo inicial es cero
y0=0.1;                                    %Al principio la cantidad de mineral extraida es cero
h=1/12;                                    %El paso nos dice que es de un mes, y como nuestra unidad es el año el paso es 1/10
                                           %Definimos la variable independiente pero como no conocemos el tiempo en el que se paraliza 
                                           %la excavacíon, por lo que del vector t solo definimos el primer elemento
t=t0;
r=240*exp(1)/10875;                        %el valor de r ha sido hallado en el apartado 2
K=10875;                                   %Cantidad total extraible
                                           %Al no conocer el tiempo final, tampoco podemos saber el valor de N, y por lo tanto, el 
                                           %tamaño de los vectores y y z, por lo que solo definimos el primer elemento de los vectores.

y(1)=y0;                                   %RK4
z(1)=y0;                                   %Heun
i=1;
while 1                                    %Empezamos el bucle de Runge Kutta
                                           %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                                   %Empezamos el bucle de Heun
                                          %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


hold on
plot(t,y,'r');                           %Dibujamos nuestra primera gráfica                                      
plot(t,z,'g');                           %Dibujamos nuestra segunda gráfica
legend('RK4','Heun','Location','best') 
hold off

• Se puede observar la gran aproximación que hay en nuestras dos gráficas obtenidas por los métodos de Runge Kutta de cuarto orden y por Heun.

1.7 Tendencia de nuestra cantidad Q cuando el tiempo lim→ ∞

• Realizamos el programa usando un tiempo máximo elegido por nosotros,que sea grande para ver la tendencia de nuestra función (tmáximo 10000)

%Q(t)=K*exp(exp(-r*t)*(log(Q0/K)))
%y por lo tanto el limite valdrá:
t=0:1:10000;                             %vector t desde 0 hasta 10000
N=length(t);                             %tamaño del vector t
Q=zeros(1,N);                            %matriz de una fila y N columnas
r=240*exp(1)/10875;                      %valor de nuestra tasa intriseca de crecimiento en Gompertz
Q0=0.1;                                  %valor inicial
K=10875;                                 %cantidad maxima extraible
Q(1)=Q0;                                 %primer valor de nuestra matriz Q es Q0
Q=K*exp(exp(-r*t)*(log(Q0/K)));          %Definimos nuestra funcion de Gompertz
for i=1:N                                %empezamos el bucle desde 1 hasta el tamaño del vector t ,(N)
    Q(i)=K*exp(exp(-r*t(i))*(log(Q0/K)));%Definimos Gompertz en funcion de cada elemento de t
end                                      %cerramos el bucle
Q(10000)                                 %               
plot(t,Q)                                %dibujamos la curva con tiempo en abcisas y Q en ordenadas
xlabel('tiempo(años)')                   %titulo al eje x
ylabel('cantidad (tn)')                  %titulo al eje y

LEYENDA:TENDENCIA DE NUESTRA CANTIDAD Q