Diferencia entre revisiones de «Explotación Minera (G12-A)»
(→Apartados 7 y 8 con el método de logístico) |
(→Estudio del modelo a los 12 años) |
||
| (No se muestran 14 ediciones intermedias del mismo usuario) | |||
| Línea 52: | Línea 52: | ||
| − | [[Archivo:ProduccionP12.jpg|400px|thumb| | + | [[Archivo:ProduccionP12.jpg|400px|thumb|centre|Curva de la función P(Q]] |
| Línea 60: | Línea 60: | ||
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. | 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. | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| Línea 126: | Línea 108: | ||
hold off}} | hold off}} | ||
| − | [[Archivo:ProduccionP220.jpg|800px|thumb| | + | [[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. | 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. | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| Línea 170: | Línea 124: | ||
r=240*exp(1)/k; %Tasa intrínseca de crecimiento | r=240*exp(1)/k; %Tasa intrínseca de crecimiento | ||
t(1)=t0; | t(1)=t0; | ||
| − | Q(1)=0.1; | + | 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; | i=1; | ||
while 1 | while 1 | ||
| Línea 185: | Línea 139: | ||
plot(t,Q) %Gráfica de la producción en función del tiempo}} | plot(t,Q) %Gráfica de la producción en función del tiempo}} | ||
| − | [[Archivo:ProduccionP3.jpg|800px|thumb| | + | [[Archivo:ProduccionP3.jpg|800px|thumb|centre|Método de Euler]] |
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
='''Modelo de Runge Kutta (de cuarto orden) y Heun'''= | ='''Modelo de Runge Kutta (de cuarto orden) y Heun'''= | ||
| Línea 275: | Línea 190: | ||
hold off}} | hold off}} | ||
| − | [[Archivo:ProduccionP4.jpg|800px|thumb| | + | [[Archivo:ProduccionP4.jpg|800px|thumb|centre|Método de Runge Kutta (de 4 orden) y Heun]] |
| − | [[Archivo:ProduccionP4ZOOM.jpg|800px|thumb| | + | [[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→∞'''= | ='''Cantidad extraída Q cuando lim t→∞'''= | ||
| + | Se hace que Q tienda a infinito, y analizamos los resultados | ||
{{matlab|codigo= | {{matlab|codigo= | ||
| Línea 373: | Línea 216: | ||
ylabel('Cantidad (ton)')}} | ylabel('Cantidad (ton)')}} | ||
| − | [[Archivo:ProduccionP5.jpg|400px|thumb| | + | [[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= | {{matlab|codigo= | ||
| Línea 431: | Línea 250: | ||
ylabel('produccion (ton/año)')}} | ylabel('produccion (ton/año)')}} | ||
| − | [[Archivo:ProduccionP6.jpg|800px|thumb| | + | [[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= | {{matlab|codigo= | ||
| Línea 497: | Línea 283: | ||
legend('Euler','Location','best')}} | legend('Euler','Location','best')}} | ||
| − | [[Archivo:ProduccionP7.jpg|800px|thumb| | + | [[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= | {{matlab|codigo= | ||
| Línea 571: | Línea 322: | ||
ylabel('produccion (tn/año)')}} | ylabel('produccion (tn/año)')}} | ||
| − | [[Archivo:ProduccionP8.jpg|800px|thumb| | + | [[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:Ecuaciones Diferenciales]] | ||
[[Categoría:ED14/15]] | [[Categoría:ED14/15]] | ||
[[Categoría:Trabajos 2014-15]] | [[Categoría:Trabajos 2014-15]] | ||
Revisión actual del 00:47 13 mar 2015
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Explotación minera. Grupo 12-A |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Javier Abad, José Abad, Jose María Antón-Pacheco, Eduardo Areitio |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Interpretación del problema
- 2 Modelo logístico de Gompertz
- 3 Modelo logístico de Verhulst
- 4 Método de Euler
- 5 Modelo de Runge Kutta (de cuarto orden) y Heun
- 6 Cantidad extraída Q cuando lim t→∞
- 7 Función de prodicción P(t)
- 8 Cantidad de mineral sin extraer
- 9 Apartados 7 y 8 con el método de logístico
- 10 Estudio del modelo a los 12 años
1 Interpretación del problema
Debido a la alta demanda de un determinado mineral, se ha decidido explotar un yacimiento de una region estudiada. Los estudios han concluido que la cantidad total extraíble (K) de dicho mineral es de 10875 toneladas. Se estima un crecimiento muy rápido de la producción (toneladas/año) durante los 25 años, tras los cuales, a causa de dificultades técnicas y la caída de la demanda, decrecerá la producción lentamente. Para estudiar este problema vamos a tomar diferentes modelos matemáticos (aproximaciones numéricas computacionales a partir de los datos obtenidos en el trabajo de campo).
La relación entre la producción (P) y la cantidad extraída (Q) es una relación diferencial. P será la derivada de Q respecto del tiempo.
2 Modelo logístico de Gompertz
Un posible modelo que relaciona la producción con la cantidad extraída es el modelo logístico de Gompertz, basado en la siguiente ecuación::
[math]P(Q) = \frac{dQ}{dt} = rQ\log\left(\frac{K}{Q}\right) [/math]
Tras estudios previos se obtuvo una produccion máxima de 240 toneladas/año (máximo relativo de la función). Derivando nuestra ecuación respecto de Q e igualandola a 0 obtenemos la siguiente ecuación::
[math]P'=0=rlog(\frac{K}{Q})-r=r(log(\frac{K}{Q})-1)[/math]
Despejamos Q y obtenemos::
[math]Q = \frac{K}{e}[/math]
Introduciendo el valor de Q obtenido en los estudios, y el valor de P=240, despejamos la ecuación y obtenemos el coeficiente r::
[math]240=r\frac{K}{e}log(\frac{K}{\frac{K}{e}})=r\frac{K}{e} → r=\frac{240e}{K}=0.0599[/math]
- Modelo computacional de Gompertz en MATLAB :
k=10875; %Cantidad total extraible en toneladas
Q=0:1:10875; %Vector con la cantidad de toneladas extraídas
n=length(Q); %Tamaño del vector Q
P=zeros(1,n); %Vector de ceros de una fila y N columnas
r=240*exp(1)/10785; %coeficiente r
for i=1:n %Realizamos el bucle
P(i)=r*Q(i)*log(k/Q(i)); %Definimos la funcion P(Q)
end
plot(Q,P,'k') %Gráfica de Q (abcisas) y P (ordenadas) en color negro.
xlabel('cantidad (ton)')
ylabel('produccion (ton/año)')
Analizando la gráfica obtenida por MATLAB, podemos observar que la pendiente (en valor absoluto) inicial de la curva es mayor a la del final, como indicaba el estudio previo. La curva muestra un cambio de pendiente en el valor de Q=240 toneladas (máximo) y un fin de producción en el valor de cantidad total extraída de 10875 toneladas.
3 Modelo logístico de Verhulst
Otro posible modelo logístico es el de Verhulst, definido por la siguiente ecuación::
[math]Q'=rQ(1-\frac{Q}{k})[/math]
Para obtener el nuevo coeficiente r, procedemos de la misma forma que en el modelo de Gompertz::
[math]P'=0=r(1-\frac{2Q}{K})[/math]
Despejamos Q y obtenemos::
[math]\frac{2Q}{K}=1→ Q=\frac{K}{2}=5437.5[/math]
Introduciendo el valor de Q obtenido en la ecuación, y el valor de P=240, obtenemos el coeficiente r::
[math]240=r\frac{K}{2}(1-\frac{\frac{K}{2}}{K})=r\frac{K}{4} → r=\frac{960}{K}=0.088[/math]
- Modelo computacional de Verhulst en MATLAB :
Con los datos obtenidos creamos un programa en MATLAB para obtener una gráfica del modelo de Verhulst y poder compararlo con el de Gompertz:
k=10875; %Cantidad maxima extraible
rG=240*exp(1)/k; %Tasa intrinseca de creciemiento (Gompertz)
rV=960/k; %Tasa intrinseca de creciemiento (Verhulst)
Q=0:1:k;
n=length(Q);
for i=1:n
PGom(i)=rG*Q(i)*log(k/Q(i)); %Gompertz
PVer(i)=rV*Q(i)*(1-Q(i)/K); %Verhulst
end
subplot(1,2,1)
plot(Q,PVer,'g') %Gráfica modelo de Verhulst
xlabel('cantidad (ton)')
ylabel('produccion (ton/año)')
subplot(1,2,2)
plot(Q,PGom,'r') %Gráfica modelo de Gompertz
xlabel('cantidad')
ylabel('produccion')
hold on %Superponemos las dos gráficas
plot(Q,PVer,'g') %Gráfica modelo de Verhulst
xlabel('cantidad (ton)')
ylabel('produccion (ton/año)')
legend('Modelo Gompertz','Modelo Verhulst','Location','best')
hold off
La primera curva muestra el gráfico según el modelo logístico de Verhulst, el cual a primera vista parece aproximarse bastante bien al previo estudio de la producción P. Sin embargo al observar la segunda gráfica, en la cual se superponen las curvas de ambos modelos (modelo de Gompertz y modelo de Verhulst) se puede indicar claramente que en este caso el modelo de Gompertz se aproxima mejor a nuestro estudio de producción, indentificable por el elevado crecimiento inicial de dicha producción.
4 Método de Euler
k=10875; %Cantidad máxima extraible
t0=0;
h=1/12; %El paso es de 1 mes por año
r=240*exp(1)/k; %Tasa intrínseca de crecimiento
t(1)=t0;
Q(1)=0.1; %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
5 Modelo de Runge Kutta (de cuarto orden) y Heun
k=10875; %Cantidad máxima extraible
t0=0; %Tiempo inicial = 0
Q0=0.1; %Cantidad inicial extraida = 0
h=1/12; %El paso es de 1 mes por año
t=t0;
r=240*exp(1)/k; %Tasa intrínseca de crecimiento
y(1)=Q0; %RK4
z(1)=Q0; %Heun
i=1;
while 1 %RK4
K1=r*y(i)*log(k/y(i));
K2=r*(y(i)+1/2*K1*h)*log(k/(y(i)+1/2*K1*h));
K3=r*(y(i)+1/2*K2*h)*log(k/(y(i)+1/2*K2*h));
K4=r*(y(i)+K3*h)*log(k/(y(i)+K3*h));
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
t(i+1)=t(i)+h;
if i>1&&abs((r*y(i)*log(k/y(i)))-25)<0.1&r*y(i-1)*log(k/y(i-1))>r*y(i)*log(k/y(i))
break
end
i=i+1;
end
i=1;
while 1 %Heun
k1=r*z(i)*log(k/z(i));
k2=r*(z(i)+k1*h)*log(k/(z(i)+k1*h));
z(i+1)=z(i)+(h/2)*(k1+k2);
t(i+1)=t(i)+h;
if i>1&&abs((r*z(i)*log(k/z(i)))-25)<0.1&r*z(i-1)*log(k/z(i-1))>r*z(i)*log(k/z(i))
break
end
i=i+1;
end
subplot(1,3,1)
plot(t,y,'r'); %Gráfica RK4
legend('RK4','Location','best')
subplot(1,3,2)
plot(t,z,'y'); %Gráfica Heun
legend('Heun','Location','best')
subplot(1,3,3)
hold on %Superponemos las dos gráficas
plot(t,y,'r'); %Gráfica RK4
plot(t,z,'g'); %Gráfica Heun
legend('RK4','Heun','Location','best')
hold off
Tras superponer ambas gráficas se puede observar la gran similitud de los modelos. Sólo aumentando la gráfica hemos sido capaces de poder distinguir las dos curvas.
6 Cantidad extraída Q cuando lim t→∞
Se hace que Q tienda a infinito, y analizamos los resultados
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)')
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.
7 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:
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)')
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.
8 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:
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')
La cantidad de mineral que queda sin extraer será 424.4179 toneladas.
9 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.
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)')
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.
10 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.
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')
La máxima producción la obtenemos en t=52.8333 y es de 346.3951 toneladas.
