MateWiki - Contribuciones del usuario [es]

Explotación Minera (G12-A)

2015-03-12T22:47:53Z

Grupo12A: /* Estudio del modelo a los 12 años */

{{ TrabajoED | Explotación minera. Grupo 12-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio }}

=''' 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.

='''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 :'''

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP12.jpg|400px|thumb|centre|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.

='''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:

{{matlab|codigo=
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}}

[[Archivo:ProduccionP220.jpg|800px|thumb|centre|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.

='''Método de Euler'''=

{{matlab|codigo=
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; %Utilizamos 0.1 en vez de 0, sino no el programa no funcionaría correctamente al detectar una indeterminación
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}}

[[Archivo:ProduccionP3.jpg|800px|thumb|centre|Método de Euler]]

='''Modelo de Runge Kutta (de cuarto orden) y Heun'''=

{{matlab|codigo=
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}}

[[Archivo:ProduccionP4.jpg|800px|thumb|centre|Método de Runge Kutta (de 4 orden) y Heun]]
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|centre|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.

='''Cantidad extraída Q cuando lim t→∞'''=

Se hace que Q tienda a infinito, y analizamos los resultados

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP5.jpg|400px|thumb|centre|Curva cuando t tiende a infinito]]

Podemos observar que transcurrido un tiempo, la función tiende a ser constante en el valor de la cantidad máxima extraíble, 10875 toneladas de mineral.

='''Función de prodicción P(t)'''=

Procedemos a analizar la función P(t), con el correspondiente código MATLAB y su gráfica:

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP6.jpg|800px|thumb|centre|Gráfica de la Producción en función del tiempo]]

Como podíamos predecir, el máximo se produce en 239.9998 toneladas, ya que tomamos 0.1 de valor inicial, tenemos cierto error. El máximo se produce para t=41.2500 años.

='''Cantidad de mineral sin extraer'''=

Utilizando el método de euler, obtenido previamente, procedemos a obtener la cantidad de mineral que quedará sin extraer al final de la vida útil de la explotación. Con el valor anteriormente calculado de Q, y restándole a éste el valor de K, obtenemos la cantidad sin extraer.
El código MATLAB es el siguiente:

{{matlab|codigo=
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')}}

[[Archivo:ProduccionP7.jpg|800px|thumb|centre|Cantidad de mineral sin extraer]]

La cantidad de mineral que queda sin extraer será 424.4179 toneladas.

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

A continuación vamos a realizar los apartados 7 y 8 utilizando el modelo logístico , aproximados por el método de Heun.

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP8.jpg|800px|thumb|centre|Método logístico]]

La producción máxima obtenida es de 239.9993 toneladas, en un tiempo t=131.4167 años. La cantidad de mineral sin extraer es de 288.8542 toneladas.
Viendo los resultados, con sus errores, podemos concluir que el modelo de Gompertz se ajusta más a la realidad.

='''Estudio del modelo a los 12 años'''=

Los datos de campo obtenidos nos indican que a los 12 años del comienzo de la explotación la cantidad total extraída es de 2695 toneladas. La cantidad de mineral que queda por extraer debido a las nuevas técnicas utilizadas es de 9075 toneladas, por lo que la cantidad total extraíble K será 9075+2695=11770 toneladas.

{{matlab|codigo=
clear all

%DATOS DEL PROBLEMA ANTIGUO
t0=0;
Q0=100; %La cantidad de mineral extraido en el comienzo es 0
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado.
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t=t0;
r=240*exp(1)/10875; %Mismo valor que el hallado en el apartado 2
K=10875; %Cantidad total extraible
Q(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer

%elemento del vector Q.
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); %Método de Heun
t(i+1)=t(i)+h;
i=i+1;
end

[maximo,tiempo]=max(P);
tiempo=tiempo/12;
plot(t,P,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo viejo','Location','best')

%NUEVOS DATOS
t0=0;
Q0=100; %Cantidad de mineral extraido inicialmente
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t1=t0;
r2=0.08; %el valor de r revisado
K=11770; %Cantidad total extraible
Q2(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer
%elemento del vector Q.
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); %Metodo Heun
t1(c+1)=t1(c)+h;
c=c+1;
end

[maximo,tiempo]=max(P1)
tiempo=tiempo/12

hold on
plot(t1,P1)
xlabel('tiempo')
ylabel('produccion tonelada/año')
legend('Modelo nuevo','modelo antiguo')}}

[[Archivo:graficaP3.jpg|800px|thumb|centre|Modelo nuevo]]
[[Archivo:graficaP2.jpg|800px|thumb|centre|Modelo antiguo]]
[[Archivo:graf11.jpg|800px|thumb|centre|Comparación de los modelos]]

La máxima producción la obtenemos en t=52.8333 y es de 346.3951 toneladas.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

Archivo:Graf11.jpg

2015-03-12T22:45:38Z

Grupo12A:

Explotación Minera (G12-A)

2015-03-11T15:40:45Z

Grupo12A: /* Estudio del modelo a los 12 años */

{{ TrabajoED | Explotación minera. Grupo 12-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio }}

=''' 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.

='''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 :'''

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP12.jpg|400px|thumb|centre|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.

='''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:

{{matlab|codigo=
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}}

[[Archivo:ProduccionP220.jpg|800px|thumb|centre|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.

='''Método de Euler'''=

{{matlab|codigo=
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; %Utilizamos 0.1 en vez de 0, sino no el programa no funcionaría correctamente al detectar una indeterminación
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}}

[[Archivo:ProduccionP3.jpg|800px|thumb|centre|Método de Euler]]

='''Modelo de Runge Kutta (de cuarto orden) y Heun'''=

{{matlab|codigo=
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}}

[[Archivo:ProduccionP4.jpg|800px|thumb|centre|Método de Runge Kutta (de 4 orden) y Heun]]
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|centre|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.

='''Cantidad extraída Q cuando lim t→∞'''=

Se hace que Q tienda a infinito, y analizamos los resultados

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP5.jpg|400px|thumb|centre|Curva cuando t tiende a infinito]]

Podemos observar que transcurrido un tiempo, la función tiende a ser constante en el valor de la cantidad máxima extraíble, 10875 toneladas de mineral.

='''Función de prodicción P(t)'''=

Procedemos a analizar la función P(t), con el correspondiente código MATLAB y su gráfica:

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP6.jpg|800px|thumb|centre|Gráfica de la Producción en función del tiempo]]

Como podíamos predecir, el máximo se produce en 239.9998 toneladas, ya que tomamos 0.1 de valor inicial, tenemos cierto error. El máximo se produce para t=41.2500 años.

='''Cantidad de mineral sin extraer'''=

Utilizando el método de euler, obtenido previamente, procedemos a obtener la cantidad de mineral que quedará sin extraer al final de la vida útil de la explotación. Con el valor anteriormente calculado de Q, y restándole a éste el valor de K, obtenemos la cantidad sin extraer.
El código MATLAB es el siguiente:

{{matlab|codigo=
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')}}

[[Archivo:ProduccionP7.jpg|800px|thumb|centre|Cantidad de mineral sin extraer]]

La cantidad de mineral que queda sin extraer será 424.4179 toneladas.

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

A continuación vamos a realizar los apartados 7 y 8 utilizando el modelo logístico , aproximados por el método de Heun.

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP8.jpg|800px|thumb|centre|Método logístico]]

La producción máxima obtenida es de 239.9993 toneladas, en un tiempo t=131.4167 años. La cantidad de mineral sin extraer es de 288.8542 toneladas.
Viendo los resultados, con sus errores, podemos concluir que el modelo de Gompertz se ajusta más a la realidad.

='''Estudio del modelo a los 12 años'''=

Los datos de campo obtenidos nos indican que a los 12 años del comienzo de la explotación la cantidad total extraída es de 2695 toneladas. La cantidad de mineral que queda por extraer debido a las nuevas técnicas utilizadas es de 9075 toneladas, por lo que la cantidad total extraíble K será 9075+2695=11770 toneladas.

{{matlab|codigo=
clear all

%DATOS DEL PROBLEMA ANTIGUO
t0=0;
Q0=100; %La cantidad de mineral extraido en el comienzo es 0
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado.
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t=t0;
r=240*exp(1)/10875; %Mismo valor que el hallado en el apartado 2
K=10875; %Cantidad total extraible
Q(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer

%elemento del vector Q.
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); %Método de Heun
t(i+1)=t(i)+h;
i=i+1;
end

[maximo,tiempo]=max(P);
tiempo=tiempo/12;
plot(t,P,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo viejo','Location','best')

%NUEVOS DATOS
t0=0;
Q0=100; %Cantidad de mineral extraido inicialmente
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t1=t0;
r2=0.08; %el valor de r revisado
K=11770; %Cantidad total extraible
Q2(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer
%elemento del vector Q.
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); %Metodo Heun
t1(c+1)=t1(c)+h;
c=c+1;
end

[maximo,tiempo]=max(P1)
tiempo=tiempo/12

hold on
plot(t1,P1)
xlabel('tiempo')
ylabel('produccion tonelada/año')
legend('Modelo nuevo','modelo antiguo')}}

[[Archivo:graficaP3.jpg|800px|thumb|centre|Modelo nuevo]]
[[Archivo:graficaP2.jpg|800px|thumb|centre|Modelo viejo]]
[[Archivo:graficaP1.jpg|800px|thumb|centre|Comparación de los modelos]]

La máxima producción la obtenemos en t=52.8333 y es de 346.3951 toneladas.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

Explotación Minera (G12-A)

2015-03-11T15:39:57Z

Grupo12A: /* Estudio del modelo a los 12 años */

{{ TrabajoED | Explotación minera. Grupo 12-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio }}

=''' 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.

='''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 :'''

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP12.jpg|400px|thumb|centre|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.

='''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:

{{matlab|codigo=
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}}

[[Archivo:ProduccionP220.jpg|800px|thumb|centre|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.

='''Método de Euler'''=

{{matlab|codigo=
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; %Utilizamos 0.1 en vez de 0, sino no el programa no funcionaría correctamente al detectar una indeterminación
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}}

[[Archivo:ProduccionP3.jpg|800px|thumb|centre|Método de Euler]]

='''Modelo de Runge Kutta (de cuarto orden) y Heun'''=

{{matlab|codigo=
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}}

[[Archivo:ProduccionP4.jpg|800px|thumb|centre|Método de Runge Kutta (de 4 orden) y Heun]]
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|centre|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.

='''Cantidad extraída Q cuando lim t→∞'''=

Se hace que Q tienda a infinito, y analizamos los resultados

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP5.jpg|400px|thumb|centre|Curva cuando t tiende a infinito]]

Podemos observar que transcurrido un tiempo, la función tiende a ser constante en el valor de la cantidad máxima extraíble, 10875 toneladas de mineral.

='''Función de prodicción P(t)'''=

Procedemos a analizar la función P(t), con el correspondiente código MATLAB y su gráfica:

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP6.jpg|800px|thumb|centre|Gráfica de la Producción en función del tiempo]]

Como podíamos predecir, el máximo se produce en 239.9998 toneladas, ya que tomamos 0.1 de valor inicial, tenemos cierto error. El máximo se produce para t=41.2500 años.

='''Cantidad de mineral sin extraer'''=

Utilizando el método de euler, obtenido previamente, procedemos a obtener la cantidad de mineral que quedará sin extraer al final de la vida útil de la explotación. Con el valor anteriormente calculado de Q, y restándole a éste el valor de K, obtenemos la cantidad sin extraer.
El código MATLAB es el siguiente:

{{matlab|codigo=
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')}}

[[Archivo:ProduccionP7.jpg|800px|thumb|centre|Cantidad de mineral sin extraer]]

La cantidad de mineral que queda sin extraer será 424.4179 toneladas.

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

A continuación vamos a realizar los apartados 7 y 8 utilizando el modelo logístico , aproximados por el método de Heun.

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP8.jpg|800px|thumb|centre|Método logístico]]

La producción máxima obtenida es de 239.9993 toneladas, en un tiempo t=131.4167 años. La cantidad de mineral sin extraer es de 288.8542 toneladas.
Viendo los resultados, con sus errores, podemos concluir que el modelo de Gompertz se ajusta más a la realidad.

='''Estudio del modelo a los 12 años'''=

Los datos de campo obtenidos nos indican que a los 12 años del comienzo de la explotación la cantidad total extraída es de 2695 toneladas. La cantidad de mineral que queda por extraer debido a las nuevas técnicas utilizadas es de 9075 toneladas, por lo que la cantidad total extraíble K será 9075+2695=11770 toneladas.

{{matlab|codigo=
clear all

%DATOS DEL PROBLEMA ANTIGUO
t0=0;
Q0=100; %La cantidad de mineral extraido en el comienzo es 0
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado.
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t=t0;
r=240*exp(1)/10875; %Mismo valor que el hallado en el apartado 2
K=10875; %Cantidad total extraible
Q(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer

%elemento del vector Q.
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); %Método de Heun
t(i+1)=t(i)+h;
i=i+1;
end

[maximo,tiempo]=max(P);
tiempo=tiempo/12;
plot(t,P,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo viejo','Location','best')

%NUEVOS DATOS
t0=0;
Q0=100; %Cantidad de mineral extraido inicialmente
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t1=t0;
r2=0.08; %el valor de r revisado
K=11770; %Cantidad total extraible
Q2(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer
%elemento del vector Q.
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); %Metodo Heun
t1(c+1)=t1(c)+h;
c=c+1;
end

[maximo,tiempo]=max(P1)
tiempo=tiempo/12

hold on
plot(t1,P1)
xlabel('tiempo')
ylabel('produccion tonelada/año')
legend('Modelo nuevo','modelo antiguo')}}

[[Archivo:graficaP3.jpg|800px|thumb|centre|Modelo viejo]]
[[Archivo:graficaP2.jpg|800px|thumb|centre|Modelo nuevo]]
[[Archivo:graficaP1.jpg|800px|thumb|centre|Comparación de los modelos]]

La máxima producción la obtenemos en t=52.8333 y es de 346.3951 toneladas.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

Explotación Minera (G12-A)

2015-03-10T18:11:55Z

Grupo12A:

{{ TrabajoED | Explotación minera. Grupo 12-A | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio }}

=''' 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.

='''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 :'''

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP12.jpg|400px|thumb|centre|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.

='''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:

{{matlab|codigo=
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}}

[[Archivo:ProduccionP220.jpg|800px|thumb|centre|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.

='''Método de Euler'''=

{{matlab|codigo=
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; %Utilizamos 0.1 en vez de 0, sino no el programa no funcionaría correctamente al detectar una indeterminación
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}}

[[Archivo:ProduccionP3.jpg|800px|thumb|centre|Método de Euler]]

='''Modelo de Runge Kutta (de cuarto orden) y Heun'''=

{{matlab|codigo=
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}}

[[Archivo:ProduccionP4.jpg|800px|thumb|centre|Método de Runge Kutta (de 4 orden) y Heun]]
[[Archivo:ProduccionP4ZOOM.jpg|800px|thumb|centre|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.

='''Cantidad extraída Q cuando lim t→∞'''=

Se hace que Q tienda a infinito, y analizamos los resultados

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP5.jpg|400px|thumb|centre|Curva cuando t tiende a infinito]]

Podemos observar que transcurrido un tiempo, la función tiende a ser constante en el valor de la cantidad máxima extraíble, 10875 toneladas de mineral.

='''Función de prodicción P(t)'''=

Procedemos a analizar la función P(t), con el correspondiente código MATLAB y su gráfica:

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP6.jpg|800px|thumb|centre|Gráfica de la Producción en función del tiempo]]

Como podíamos predecir, el máximo se produce en 239.9998 toneladas, ya que tomamos 0.1 de valor inicial, tenemos cierto error. El máximo se produce para t=41.2500 años.

='''Cantidad de mineral sin extraer'''=

Utilizando el método de euler, obtenido previamente, procedemos a obtener la cantidad de mineral que quedará sin extraer al final de la vida útil de la explotación. Con el valor anteriormente calculado de Q, y restándole a éste el valor de K, obtenemos la cantidad sin extraer.
El código MATLAB es el siguiente:

{{matlab|codigo=
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')}}

[[Archivo:ProduccionP7.jpg|800px|thumb|centre|Cantidad de mineral sin extraer]]

La cantidad de mineral que queda sin extraer será 424.4179 toneladas.

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

A continuación vamos a realizar los apartados 7 y 8 utilizando el modelo logístico , aproximados por el método de Heun.

{{matlab|codigo=
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)')}}

[[Archivo:ProduccionP8.jpg|800px|thumb|centre|Método logístico]]

La producción máxima obtenida es de 239.9993 toneladas, en un tiempo t=131.4167 años. La cantidad de mineral sin extraer es de 288.8542 toneladas.
Viendo los resultados, con sus errores, podemos concluir que el modelo de Gompertz se ajusta más a la realidad.

='''Estudio del modelo a los 12 años'''=

Los datos de campo obtenidos nos indican que a los 12 años del comienzo de la explotación la cantidad total extraída es de 2695 toneladas. La cantidad de mineral que queda por extraer debido a las nuevas técnicas utilizadas es de 9075 toneladas, por lo que la cantidad total extraíble K será 9075+2695=11770 toneladas.

{{matlab|codigo=
clear all

%DATOS DEL PROBLEMA ANTIGUO
t0=0;
Q0=100; %La cantidad de mineral extraido en el comienzo es 0
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado.
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t=t0;
r=240*exp(1)/10875; %Mismo valor que el hallado en el apartado 2
K=10875; %Cantidad total extraible
Q(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer

%elemento del vector Q.
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); %Método de Heun
t(i+1)=t(i)+h;
i=i+1;
end

[maximo,tiempo]=max(P);
tiempo=tiempo/12;
plot(t,P,'r')
xlabel('tiempo')
ylabel('produccion (tn/año)')
legend('Modelo viejo','Location','best')

%NUEVOS DATOS
t0=0;
Q0=100; %Cantidad de mineral extraido inicialmente
h=1/12; %Adquirimos una unidad de mes respecto a los doce que tiene un año como dice el enunciado
%Como no sabemos cuando se paraliza la excavación definiremos el
%primer elemento de dicho vector.
t1=t0;
r2=0.08; %el valor de r revisado
K=11770; %Cantidad total extraible
Q2(1)=Q0; %Como no sabemos el valor del tiempo final solo definiermos el primer
%elemento del vector Q.
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); %Metodo Heun
t1(c+1)=t1(c)+h;
c=c+1;
end

[maximo,tiempo]=max(P1)
tiempo=tiempo/12

hold on
plot(t1,P1)
xlabel('tiempo')
ylabel('produccion tonelada/año')
legend('Modelo nuevo','modelo antiguo')}}

[[Archivo:graficaP3.jpg|800px|thumb|centre|Modelo viejo]]
[[Archivo:graficaP2.jpg|800px|thumb|centre|Modelo nuevo]]
[[Archivo:graficaP1.jpg|800px|thumb|centre|Comparación de los modelos]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

Archivo:GraficaP1.jpg

2015-03-10T18:06:25Z

Grupo12A:

Archivo:GraficaP2.jpg

2015-03-10T18:06:14Z

Grupo12A:

2015-03-06T19:06:02Z

Grupo12A:

Explotación Minera (G12-A)