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. Este método se caracteriza en que se utiliza en casos en el que el crecimiento es muy rápido hasta llegar a su máximo y a partir de entonces, empieza a disminuir. Este modelo se usa para casos como por ejemplo las ventas de dispositivos electrónicos como móviles o el crecimiento de tumores o virus.

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 :

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

• 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 de una 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.

• 3.Muestra muy bien el significado del modelo de Gompertz,ya que existe un crecimiento incial muy potente y un posterior descenso,causadas por factores adversos a la rentabilidad de nuestra explotación minera una vez transcurrido los años.

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. En la tercera subventana se comparan ambas curvas mencionadas previamente donde se muestra que para este tipo de calculo de crecimientos en el que es rápido al principio y luego disminuye lentamente, el modelo de Gompertz se aproxima mas y queda mejor definido .

1.4 Utilización del método de Euler,RK4 y Heun. Desarrollo teórico.

• 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.

Runge Kutta de orden 4 es un método de resolución para la aproximación de soluciones de problemas de ecuaciones diferenciales del problema del valor inicial. Para poder aplicarlo tendremos que hallar una serie de constantes (K1, K2, K3, K4) que están definidas de manera explícita, y así obtener nuestras aproximaciones. El método de aproximación Runge Kutta que vamos a aplicar es de orden 4 y se designa como "RK4".

Heun es otro método de aproximación que podemos aplicar, y que también hemos de hallar unas constantes (K1,K2), que como en RK4, las podremos obtener de manera explícita.

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

• Análisis de nuestra curva:

LEYENDA:MÉTODO DE EULER

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

• Análisis de nuestras gráficas :

LEYENDA:RK4 Y HEUN

• 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.
• Si aumentamos la imágen se puede apreciar mejor cada método

LEYENDA:RK4 Y HEUN AUMENTADA

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

• Análisis de nuestra curva:

LEYENDA:TENDENCIA DE NUESTRA CANTIDAD Q

• 1.La función de Gompertz nos dice que lim→∞P(t)=K

• 2.En nuestro caso la población P(t) es nuestra cantidad Q(t), por ello si realizamos una curva de nuestra cantidad en un periodo de tiempo largo podremos ver que tendecia llevará al aproximarse al infinito.

• 3.Aquí se ve con claridad la tendencia que sigue con nuestro tiempo elegido ya que se observa una línea horizontal bastante precisa, después de un crecimiento en los primeros t(años).

1.8 Función de P(t) durante la vida útil de la población

• 1. Es fácil de ver que Q depende del tiempo y por lo consiguiente, P también depende del tiempo.
• Para hallar su curva representativa realizaremos un programa en el que obtendremos las soluciones de Q(t) con las que podremos obtener los valores de P(t).

• 2. Al ejecutar el programa nos da el punto de máxima producción y el momento en el que se dá esta.

% Apartado 8. P(t)
clear all, clf
                    %DATOS DEL PROBLEMA:
t0=0;
Q0=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/12
                    %Definimos la vble independiente pero como no conocemos el tiempo en el que se paraliza la excavacíon,por lo que t=t0
                    
t=t0;               %vector t solo definimos el primer elemento
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;            %el primer elemento de Q es Q0
i=1;                %nuestro buble empieza con el elemento 1
while 1             %bucle que rastrea indefinidamente hasta encontrar el elemnto que cumple nuestra condicion if
   P(i)=r*Q(i)*log(K/Q(i)); %definimos nuestra funcion
    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))%establecemos la condicion
                    %si se cumple dicha condicion con break se rompe el bucle y con end se acaba
break
    end
    Q(i+1)=Q(i)+h*(r*Q(i)*log(K/Q(i))); %euler
     t(i+1)=t(i)+h;%tiempo
   i=i+1;           %si no se ha cumplido con i=1 tenemos que continuar haciendo el bucle con el elemnto siguiente y asi sucesivamente
end                 %se cierra el bucle
    [maximo,tiempo]=max(P)% obtenemos el maximo de nuestra funcion y la posicion de este que sera el tiempo
    tiempo=tiempo/12%el tiempo lo manejamos en nuestro ejercicio en meses

plot(t,P)           %dibujamos la curva con tiempo en el eje x y la produccion en funcion de toneladas año
xlabel('tiempo')    %titulo al eje x
ylabel('produccion (tn/año)')% titulo al eje y

• Análisis de nuestra curva:

• Y como podemos comprobar al ejecutarlo, la máxima producción es de 239.9998 toneladas (como ya podíamos adelantar) y sé da a los 41.2500 años (NOTA: hay que tener en cuenta que este valor no es el real ya que hemos supuesto la cantidad de mineral extraible es de 0.1 para poder hacer el problema).

• Con la gráfica adjunta podemos apreciar que los valores son los que nos ha dado el programa al ejecutarlo.

LEYENDA:Producción en función del tiempo

1.9 Cantidad de mineral sin extraer una vez transcurrida la vida útil

• Para obtener la cantidad de mineral sin extraer una vez transcurrida la vida útil utilizaremos el programa que usamos en el apartado 4, en el cual aproximábamos la solución Q(t).
• Basta con obtener el valor de Q cuando transcurre la vida útil de la explotación.
• Como sabemos que la cantidad mineral que se espera extraer es K, simplemente bastará con restar la Q que se ha extraído a la K y obtenemos la cantidad de mineral que queda por extraer. Vamos a ejecutar el siguiente problema que nos sacará a pantalla los datos pedidos.

% Apartado 9. Valor de Q cuando transcurre su vida útil.
clear all, clf
                     %DATOS DEL PROBLEMA
t0=0;
Q0=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/12
                     % Definimos la vble independiente perocomo no conocemos el tiempo en el que se paraliza la excavacíon,por lo quet=t0
                     %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 del vector Q, por lo que solo definimos el primer
                     %elemento del vector
Q(1)=Q0;
i=1;                 %el bucle para empezar necesita un primer elemento i=1en nuesdtra condicion if
while 1              %este bucle realiza un rastreeo elementoa elemento hasta encontrar el buscado 
     Q(i+1)=Q(i)+h*(r*Q(i)*log(K/Q(i))); %euler
     t(i+1)=t(i)+h;  %tiempo
    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))%condicion
      
break                %si se encuentra el emento de la condicion rompemos el bucle y se acaba (end)
    end
   i=i+1;            %si en nuestro 1er elemento no esta el buscado el bucle tiene que continuar con el siguiente y asi hasta encontrarlo
end                  %se cierra el bucle
N=length(Q);         %tamaño del vector N
MineralExtraido=Q(N);%cantidad de mineral extraida
MineralSinExtraer=K-MineralExtraido %cantidad de mineral que queda por extraer
plot(t,Q,'r')        %dibujamos en rojo nuestra curva tiempo eje x Q eje y
legend('Euler','Location','best') %añadimos leyenda

• La cantidad de mineral que queda sin extraer es de 424.4179 toneladas.
• Adjuntamos un zoom en la parte final en la que se puede apreciar la cantidad de mineral final cuando ha transcurrido la vida útil y que la diferencia de esta con la total es de 420 aproximadamente, es decir, el valor que nos ha dado al ejecutar el programa.

LEYENDA:Cantidad de mineral extraido

1.10 Apartado 1.8 y 1.9 con el método logístistico de Verhulst

• Para este apartado vamos a proceder de la misma manera que en los dos apartados anteriores pero con la diferencia de que ahora utilizaremos el método logístico de Verhulst. Resolveremos las dos soluciones con un mismo programa.

% Apartado 10. P(t) con Heun.
clear all, clf
                %DATOS DEL PROBLEMA;
t0=0;
Q0=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/12
                 % Definimos la vble 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*4)/10875; %el valor de r ha sido hallado en el apartado 4
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 para empezar necesita un primer elemento i=1en nuesdtra condicion if 
while 1          %este bucle realiza un rastreeo elementoa elemento hasta encontrar el buscado 
     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)%nuestra condicion
      
break            %si encontramos el elemento que cumple la condicion rompemos el bucle        
    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
[maximo,tiempo]=max(P)                %obtenemos em maximo de nuestra funcion P y la posicion de este que sera nuestro tiempo
tiempo=tiempo/12                      %manejamos el tiempo en meses
N=length(Q);                          %tamaño del vector Q
MineralExtraido=Q(N);                 %cantidad de mineral extraida
MineralSinExtraer=K-MineralExtraido   %cantidad de mineral que queda por extraer

plot(t,P)                             %dibujamos la curva
legend('Heun','Location','best')      %añadimos leyenda
xlabel('tiempo')                      %añadimos titulo al eje x
ylabel('produccion (tn/año)')

• Los resultados obtenidos son:

• 1- La producción máxima es:239.9993 toneladas y se da a los 131.4167 años.

• 2- La cantidad de mineral sin extraer una vez transcurrida la vida útil es de 290.9657 toneladas.
• Estos resultados aparecen en pantalla una vez se ejecutan el programa.

• Como se puede comprobar y como ya hemos hablado en el apartado tres, el método de Gompertz y el método de Verhulst difieren bastante. La solución del primer método es una doble exponencial mientras que en el segundo método, el logístico tradicional, la solución es simplemente una exponencial. Para este tipo de estudios de crecimiento en el que al principio es muy rápido hasta llegar a un punto en el que disminuye lentamente.

Producción con método Verhulst

1.11 Revisión del comportamiento del modelo a los 12 años

• A los doce años de comernzar la explotación se realiza una revisión del modelo. Se comprueba que la cantidad extraída de mineral hasta la fecha es de 2695 toneladas. Además nos indica que la cantidad de mineral que queda por extraer es de 9075 toneladas (K-Q(12)=9075). Por lo que la cantidad total extraible esperada K será K=9075+2695=11770 toneladas.

• Ahora que ya tenemos las nuevas constantes r y K aproximamos la producción mediante el método de Heun. al ejecutar el programa nos dará la máxima producción y cuando se da esta. También saldrá a pantalla el valor de la cantidad mineral que quedará al final de la vida útil.

• Para hallar la tasa intrinseca de creciemiento (r), vamos a calcular su valor, para ello tiene que cumplir que la cantidad de mineral extraído en 12 años sea 2695.

% Apartado 11. P(t) revisada
clear all, clf
%DATOS DEL PROBLEMA ANTIGUO
t0=0;
Q0=100;  %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/12
         %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 del vector Q, por lo que solo definimos el primer
                    %elemento del vector
Q(1)=Q0;
i=1;
while 1
   P(i)=r*Q(i)*log(K/Q(i));
    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
    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
    [maximo,tiempo]=max(P);
    tiempo=tiempo/12;
subplot(1,3,1)
plot(t,P)
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo viejo','Location','best') 
%DATOS DEL PROBLEMA NUEVO
t0=0;
Q0=100;  %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/12
         %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
t1=t0;
r2=0.08; %el valor de r revisado
K=11770;            %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
Q2(1)=Q0;
c=1;
while 1
   P1(c)=r2*Q2(c)*log(K/Q2(c));
    if c>1&&abs((r2*Q2(c)*log(K/Q2(c)))-25)<0.1&r2*Q2(c-1)*log(K/Q2(c-1))>r2*Q2(c)*log(K/Q2(c))
      
break
    end
    K1=r2*Q2(c)*(1-Q2(c)/K);
    K2=r2*(Q2(c)+K1*h)*(1-(Q2(c)+K1*h)/K);
    Q2(c+1)=Q2(c)+(h/2)*(K1+K2);       %Heun
    t1(c+1)=t1(c)+h;
   c=c+1;
end
    [maximo,tiempo]=max(P1)
    tiempo=tiempo/12

subplot(1,3,2)
plot(t1,P1,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo nuevo','Location','best') 
subplot(1,3,3)
plot(t,P)
hold on
plot(t1,P1,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo antiguo','Modelo nuevo','Location','best') 
hold off

Comparación de los dos métodos

Como podemos comprobar en la gráfica, con los nuevos datos el crecimiento va a ser mucho más rápido, además de que el máximo sea mayor que para los datos iniciales. C