Revisión actual del 22:40 6 mar 2015

1 Introducción

1.1 Relación entre Producción-Material extraido

Dada la ecuación diferencial proporcionada en el problema, donde se tiene la derivada de Q con respecto al tiempo, la función de producción es la propia derivada, ya que Q es la función de distribución y P es la densidad, con lo cual la relación entre ambas es la derivada, esto es:

dQ/dt=rQ log(K/Q)=P

1.2 Valor de la tasa intrínseca de crecimiento (r)

Como ya se ha comentado anteriormente la función P(producción) es la derivada de la cantidad extraída de material,

dx/dt=rQ log(K/Q)=P

sabemos que la producción máxima es de 240, por tanto, derivando la función P e igualando a 0 obtenemos donde se produce el máximo

dP/dt=0    t= c/r

que introducido en la ecuación P e igualando a 240 se obtiene el valor de la tase de crecimiento r pedido

r= 240e/k =0.06

1.3 Relación entre la producción y el volumen de toneladas extraidas según el modelo de Gompertz

Mediante el modelo de Gompertz, podemos expresar la producción de minerales (P en función de Q), como hemos visto en el apartado anterior. En la gráfica se observa cómo dicha producción aumenta hasta llegar a su máximo y a continuación disminuye. Para representar la curva necesitamos el siguiente código matlab:

b=240
k=10875;
r=b*exp(1)/k;
q=0:1:10875;
N=length(q);
P=zeros(1,N);
for i=1:N
    P(i)=r*q(i)*log(k/q(i));
end
plot(q,P)

1.4 Relación entre la producción y el volumen de toneladas extraídas según el modelo de Verhulst

Volvemos a dibujar la curva del modelo de Gompertz descrita en el apartado anterior, y en una segunda ventana la curva de Verhulst. Para ello, necesitamos definir una nueva función, así como una nueva tasa intrínseca de crecimiento.

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

Añadimos una tercera ventana en la que se dibujan las gráficas de ambos modelos, para poder compararlas. Necesitamos el siguiente código matlab para obtener los resultados:

b=240;
k=10875;
r1=b*exp(1)/k;%Tasa intrínseca de crecimiento según el modelo de Gompertz 
r2=(240*4)/10875; %Tasa intrinseca de creciemiento según el modelo de Verhulst
q=0:1:10875;
N=length(q);
P=zeros(1,N);
for i=1:N
P(i)=r1*q(i)*log(k/q(i)); %Gompertz
PV(i)=r2*q(i)*(1-q(i)/k);  %Verhulst
end
subplot(1,3,1)
plot(q,P,'k') %Curva para la función de Gompertz
xlabel('cantidad')            
ylabel('produccion')         
subplot(1,3,2)
plot(q,PV,'r') %Curva para la función de Verhulst
xlabel('cantidad (tn)')      
ylabel('produccion (tn/año)') 
subplot(1,3,3) %Comparación entre las curvas de ambos modelos
plot(q,P,'k')                   
xlabel('cantidad')             
ylabel('produccion')         
hold on %Para superponer gráficas                   
plot(q,PV,'r')              
xlabel('cantidad (tn)')        
ylabel('produccion (tn/año)')  
legend('Modelo Gompertz','Modelo Verhulst','Location','best') %Leyenda
hold off

1.5 Problema de valor inicial y aproximación de la cantidad de material extraído por el método de Euler

Sea el problema de valor inicial [math]Q'=rQ(1-\frac{Q}{k})[/math] tal que [math]Q(0)=0.1[/math] , la gráfica representa la variación de cantidad de mineral extraible a lo largo del tiempo. Como se puede observar, en un primer momento podemos extraer mucha cantidad de mineral hasta que, al llegar a un tiempo t y debido a que los recursos minerales limitados, la tasa de rendimiento de extracción del mineral es inferior a 25 toneladas/año (estamos alcanzando la máxima cantidad extraible). Por lo tanto, ya no es rentable seguir con la extracción.

t0=0;
h=1/12;
k=10875;
b=240;
r=b*exp(1)/k;
t(1)=t0;
q(1)=0.1;
c=log(log(k/q(1)));
Q=k/(exp(exp(-r*t+c)));
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))));
      
        break
    end
    i=i+1;
end
[t',q']
plot(t,q)

1.6 Aproximación de la cantidad de material extraído con los métodos de Runge kutta y Heun

En este apartado se ha aproximado la cantidad de mineral que se extrae a través de dos métodos mss precisos que el de Euler del apartado anterior, Runge Kutta y Heun. Se aprecia ligeramente la diferencia en cuanto a precisión se refiere entre la gráfica del presente apartado, y la obtenida por el método de Euler.

clear all
%Condiciones iniciales
t0=0;
q0=0.1;
h=1/12;
k=10875;
b=240;
r=b*exp(1)/k;
%Variable dependiente
t=t0;
q(1)=q0;
z(1)=q0;
i=1;
while 1
K1=r*q(i)*(log(k)-log(q(i)));
K2=r*(q(i)+1/2*K1*h)*(log(k)-log(q(i)+1/2*K1*h));
K3=r*(q(i)+1/2*K2*h)*(log(k)-log(q(i)+1/2*K2*h));
K4=r*(q(i)+K3*h)*(log(k)-log(q(i)+K3*h));
q(i+1)=q(i)+h/6*(K1+2*K2+2*K3+K4);
t(i+1)=t(i)+h;
    if i>1&&abs(r*q(i)*(log(k)-log(q(i)))-25)<0.1&&r*q(i-1)*(log(k)-log(q(i-1)))>r*q(i)*(log(k)-log(q(i)))
        break
    end
   i=i+1;
end
i=1;
while 1
    K1=r*z(i)*(log(k)-log(z(i)));
    K2=r*(z(i)+K1*h)*(log(k)-log(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)-log(z(i)))-25)<0.1&&r*z(i-1)*(log(k)-log(z(i-1)))>r*z(i)*(log(k)-log(z(i)))
        break
    end
i=i+1;
end
hold on
plot(t,q,'b*')
plot(t,z,'r+')
legend('RK4','HEUN','location','best')
hold off

1.7 Cantidad de mineral cuando el tiempo tiende a infinito

Al buscar la cantidad de mineral cuando t (tiempo) tiende a infinito, se ha creado un programa matlab con la opción de introducir un tiempo grande que se acerque al infinito, de tal forma que la funcion de Q (cantidad de mineral) dependerá de dicho tiempo. La gráfica que se obtiene, presentará, en un tiempo elevado, un valor prácticamente constante hacia el final de la mencionada gráfica, dejando ver que la cantidad de mineral crece durante unos años que coinciden con el principio de la explotación para luego mantenerse constante.

tN=input('Introduce un valor del tiempo muy grande:');
t=0:1:tN;                             
N=length(t);                             
Q=zeros(1,N);                            
r=240*exp(1)/10875;                      
Q0=0.1;                                  
K=10875;                                 
Q(1)=Q0;                                 
Q=K*exp(exp(-r*t)*(log(Q0/K)));          
for i=1:N                                
    Q(i)=K*exp(exp(-r*t(i))*(log(Q0/K)));
end                                      
Q(tN)                                                
plot(t,Q)

1.8 Representación de la función P(t)

Como se puede ver, en un principio la producción crece rápidamente. Esto es debido a que en un principio la cantidad de material que es posible extraer de la explotación minera es mucho mayor que la cantidad extraída. Pasado un tiempo (42 años ) se alcanza el máximo de producción de 239.9906 toneladas y, a partir de ahí, la producción desciende ya que según pasa el tiempo cada vez queda menos material por extraer y poco a poco la explotación minera se agota. Se adjunta el código matlab que da la función de producción en función del tiempo.

clear all
t=0:1:200;
n=length(t);
Q0=0.1;
K=10875;
r=240*exp(1)/K;
C=log(log(K/Q0));
for i=1:n
    P(i)=(K*r*exp(-r*t(i)+C))/(exp(exp(-r*t(i)+C)));
end
[maximo,tiempo]=max(P) %maximo de la funcion P y y tiempo en el que se alcanza en años
plot(t,P)

El máximo alcanzado por la función P(t) es 239.9906 ton., y se alcanza a los 42 años La función resultante es:

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

En apartados anteriores hemos obtenido la función de euler para la cantidad de mineral extraído (Q). Creamos un vector q que varía a lo largo del tiempo, y sabemos que en el último punto del vector la función tiene el valor del total extraído. Sabiendo que k es la cantidad total que se puede extraer, basta con restarle a este valor la q calculada para obtener lo que queda sin extraer. Para ello hemos utilizado el siguiente código matlab:

t0=0;
q0=0.1;             
h=1/12;
t=t0;
k=10875;    
r=240*exp(1)/k; %Tasa intrínseca de crecimiento modelo de Gompertz
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);
extraido=q(N); %cantidad de mineral extraida
sinextraer=k-extraido %cantidad de mineral que queda por extraer

La cantidad de mineral que queda sin extraer es: 424.4179 ton.

1.10 Modelo Logístico. Aproximación de Heun

t0=0;
t=t0;
q0=0.1;
h=1/12;
k=10875;
r=(240*4)/k; %Tasa intrínseca de crecimiento según el modelo de Verhulst
q(1)=q0;
i=1;
while 1 
     q(i+1)=q(i)+h*r*q(i)*log(k/q(i));
     P(i+1)=r*q(i)*(1-q(i)/k);
     t(i+1)=t(i)+h;
    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
    %Heun
    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);       
    t(i+1)=t(i)+h;       
   i=i+1;
end
[maximo,tiempo]=max(P) %maximo de la funcion P y y tiempo en el que se alcanza en años
N=length(q);
extraido=q(N); %cantidad de mineral extraída
sinextraer=k-extraido %cantidad de mineral que queda por extraer
plot(t,P)

Heun

Repitiendo los apartados anteriores por aproximación de Heun, obtenemos los siguientes resultados:

La producción máxima resulta ser 239.9993 ton. y se alcanza en 157 años. Además quedarán sin extraer 288.8542 ton. De esta gráfica podemos extraer la conclusión de que el modelo de logístico no se ajusta a la realidad ya que no solo la producción máxima es, en realidad, mayor sino que el tiempo necesario para alcanzar dicho valor es muy superior a la realidad. Podemos concluir que en caso de elegir uno de los dos modelos es más recomendable es el modelo de Gompertz por ajustarse más a la realidad.

1.11 Reajuste de Datos

En este apartado, se propone la revisión del modelo pasados 12 años. Los datos reales proporcionados son: la cantidad de mineral extraído hasta ese momento es de 2695 toneladas; y la cantidad de mineral que falta por extraer es de 9075 toneladas [math](K-Q(12)=9075)[/math]. Por lo tanto la cantidad total de mineral que se podía extraer es de [math]K= 9075+2695=11770 [/math] toneladas.

Basándonos en la indicación del problema la tasa intrínseca de crecimiento (r) se calcula mediante un bucle que dando valores a r para que se obtenga el valor de la cantidad de mineral extraído a los 12 años.

En primer lugar hemos calculado con el modelo de Gompertz inicial los valores de la producción (P) y la cantidad de mineral extraído. Para ello hemos utilizado el bucle del apartado uno y la solución del problema de valor inicial planteado. Una vez conocido la cantidad de mineral extraído a lo largo del tiempo hemos tomado el valor para t=12 años que ocupa la posición número 13 del vector (q). A continuación hemos reajustado el modelo de forma que la nueva población límite pasa a ser de 9075 toneladas. Realizando un bucle "while" hemos introducido los valores de la producción inicial obteniendo distintos valores para la (r), siendo el que más se aproxima a nuestro valor r=0.0708.

Se observa que la tasa intrínseca de rendimiento (r) es mayor una vez revisado el modelo, lo cual es evidente ya que se han producido mejoras en las técnicas de extracción del mineral.

clear all, clc
    b=240
    k=10875;
    r0=b*exp(1)/k;
    t0=0;
    tN=25;
    h=1/12; 
    N=(tN-t0)/h;
    C=log(log(k/0.1));
    q=zeros(1,N);
    q(1)=0.1
    P=zeros(1,N);
    t=t0:h:tN;
    for i=1:N
       P(i)=r0*q(i)*log(k/q(i));
       q(i+1)=k/(exp(exp(-r0*t(i)+C)));
    end
    q0=q(13);
    K=9075;
    n=1;
    while 1
    r(n)=(C-log(log(K/q(n))))/t(n);
    if n>1&&q(n)>q0
        r(n-1);
        break
    end
    n=n+1;
    end
    r(length(r))

Para la nueva r y la nueva K procedemos a realizar la aproximación mediante el método de Heun y compararla para los valores obtenidos en los apartados anteriores.

clear all
%Aproximacion de Heun en el modelo nuevo
%Condiciones iniciales
t0=0;
q0=2695;
h=1/12;
k=9075;
r=0.0708;
%Variable dependiente
t=t0;
q(1)=q0;
i=1;
while 1
K1=r*q(i)*(log(k)-log(q(i)));
K2=r*(q(i)+K1*h)*(log(k)-log(q(i)+K1*h));
q(i+1)=q(i)+h/2*(K1+K2);
t(i+1)=t(i)+h;
if i>1&&abs(r*q(i)*(log(k)-log(q(i)))-25)<0.1&&r*q(i-1)*(log(k)-log(q(i-1)))>r*q(i)*(log(k)-log(q(i)))
break
end
i=i+1;
end
q
P=zeros(size(q));
for m=1:length(q)
    P(m)=r*q(m)*log(k/q(m));
end
P
%Aproximacion de Heun en el modelo antiguo
T0=0;
z0=2695;
K=10875;
b=240;
R=b*exp(1)/K;
T=T0;
z(1)=z0;
w=1;

while 1
    K1=R*z(w)*(log(K)-log(z(w)));
    K2=R*(z(w)+K1*h)*(log(K)-log(z(w)+K1*h));
    z(w+1)=z(w)+h/2*(K1+K2);
    T(w+1)=T(w)+h;
    if w>1&&abs(R*z(w)*(log(K)-log(z(w)))-25)<0.1&&R*z(w-1)*(log(K)-log(z(w-1)))>R*z(w)*(log(K)-log(z(w)))
        
        break
    end
w=w+1;
end
z
P2=zeros(size(z));
for s=1:length(z)
    P2(s)=R*z(s)*log(K/z(s));
end
P2
%GRAFICAS
subplot(3,1,1)
plot(t,q,'g')
legend('Modelo Gompertz Modificado','Location','best') %Leyenda
xlabel('tiempo (años)')
ylabel('Material extraido (tm)')
subplot(3,1,2)
plot(t,P)
legend('Modelo Gompertz Modificado: Producción','Location','best') %Leyenda
xlabel('tiempo (años)')
ylabel('Produccion (tm/años)')
subplot(3,1,3)
plot(z,P2,'r')
legend('Modelo Gompertz Inicial: Producción','Location','best') %Leyenda
xlabel('Material Extraido (tm)')
ylabel('Produccion (tm/años)')
%Maximo
max(P)

Graficos comparados

Podemos observar que el modelo Logístico es el más impreciso y tiene una desviación con la realidad considerable. La principal razón de esto es que la mayoría de los datos utilizados son supuestos y que los métodos utilizados son aproximaciones. Sin embargo se puede observar que no todos los métodos son igual de imprecisos y que algunos de ellos son buenas aproximaciones a la solución real del problema. Es el caso de los métodos de Heun y Runge-Kutta 4 que son los métodos más exactos. Como conclusión al trabajo podemos decir que el modelo de Gompertz es el que mejor se ajusta a la realidad entre los dos modelos estudiados.