Diferencia entre revisiones de «Explotación minera (Grupo 5C)»
De MateWiki
(→Método de Runge-Kutta) |
|||
| (No se muestran 31 ediciones intermedias del mismo usuario) | |||
| Línea 12: | Línea 12: | ||
: <math> {dQ \over dt}= P </math><br /> | : <math> {dQ \over dt}= P </math><br /> | ||
o lo que es lo mismo:<math> P=Q' </math><br /> | o lo que es lo mismo:<math> P=Q' </math><br /> | ||
| − | + | Para ello, vamos a ver la representación gráfica de cada función. | |
| − | Archivo: | + | |
| − | Archivo: | + | De esta forma, la evolución de Q respecto al tiempo tendría esta forma: |
| − | + | ||
| + | |||
| + | [[Archivo:Graficacomparativaapartado2Q.jpg|400px|Evolución de Q respecto al tiempo]] | ||
| + | |||
| + | Y la evolución de P respecto al tiempo es la siguiente: | ||
| + | |||
| + | |||
| + | [[Archivo:Graficacomparativaapartado2P.jpg|400px|Evolución de P respecto al tiempo]] | ||
| + | |||
| + | Podemos observar a simple vista que el valor de P es la pendiente de la función Q(t), ya que vemos que el punto de inflexión de la función Q coincide con el máximo de la función P. | ||
==Modelo de Gompertz: Valor de la tasa de crecimiento (r) y resolución de la función de producción== | ==Modelo de Gompertz: Valor de la tasa de crecimiento (r) y resolución de la función de producción== | ||
| Línea 38: | Línea 47: | ||
Por lo tanto, nuestra ecuación del modelo de Gompertz quedará de la siguiente forma:<br /> | Por lo tanto, nuestra ecuación del modelo de Gompertz quedará de la siguiente forma:<br /> | ||
: <math> Q(t)=10875e^{-4,47e^{-0,6t}} </math><br /> | : <math> Q(t)=10875e^{-4,47e^{-0,6t}} </math><br /> | ||
| + | |||
| + | Y su representación gráfica: | ||
| + | |||
| + | [[Archivo:Graficaapartado42.jpg|400px|Evolución de Q respecto al tiempo con el modelo de Gompertz]] | ||
===Función de producción P(Q)=== | ===Función de producción P(Q)=== | ||
Para obtener la función P(Q), obtendremos por un lado P(t) y por otro Q(t), para después representarlos en la forma P(Q).<br /> | Para obtener la función P(Q), obtendremos por un lado P(t) y por otro Q(t), para después representarlos en la forma P(Q).<br /> | ||
| − | [[Archivo: | + | [[Archivo:Graficapdeqgompertz.jpg|400px|Evolución de P(Q) con el método de Gompertz]] |
La gráfica hace referencia a la producción que se tiene en función de la cantidad de mineral extraída. Se puede ver que el máximo se alcanza para una P=240 toneladas/año cuando la cantidad extraída es Q=4000,69 toneladas.<br /> | La gráfica hace referencia a la producción que se tiene en función de la cantidad de mineral extraída. Se puede ver que el máximo se alcanza para una P=240 toneladas/año cuando la cantidad extraída es Q=4000,69 toneladas.<br /> | ||
| Línea 81: | Línea 94: | ||
Así, la ecuación del modelo tendrá la siguiente forma: | Así, la ecuación del modelo tendrá la siguiente forma: | ||
: <math> Q'=0,0883Q(1-{Q \over 10875}) </math> | : <math> Q'=0,0883Q(1-{Q \over 10875}) </math> | ||
| + | |||
| + | Y la representación gráfica queda de la siguiente forma: | ||
| + | |||
| + | [[Archivo:Graficaapartado411.jpg|400px|Evolución de Q respecto al tiempo con el modelo de Verhlust]] | ||
===Función de producción P(Q)=== | ===Función de producción P(Q)=== | ||
| Línea 88: | Línea 105: | ||
: <math> Q(t)={1196,25e^{0.088276t} \over 1+11e^{0.088276t}} </math><br /> | : <math> Q(t)={1196,25e^{0.088276t} \over 1+11e^{0.088276t}} </math><br /> | ||
Al igual que en el apartado anterior, obtenemos P(Q) dibujando ambas funciones en la misma gráfica:<br /> | Al igual que en el apartado anterior, obtenemos P(Q) dibujando ambas funciones en la misma gráfica:<br /> | ||
| − | + | ||
| − | Archivo: | + | [[Archivo:Graficapdeqverlhust.jpg|400px|Evolución de P(Q) con el método de Verlhust]] |
| − | + | ||
| − | + | Vemos que la gráfica empieza en 1000 toneladas porque en caso contrario tendríamos una indeterminación a la hora de calcular los valores. | |
| + | |||
CODIGO MATLAB<br /> | CODIGO MATLAB<br /> | ||
{{matlab|codigo= | {{matlab|codigo= | ||
| − | % | + | %Lo que haremos será calcular valores de Q para distintos valores de t, y |
| − | + | %luego calcular valores de P para esos mismo valores de t. Una vez tengamos | |
| − | + | %los valores de P y Q en sendos vectores, los representaremos | |
| − | + | K=10875; | |
| − | t= | + | C=0.11; |
| − | + | r=0.0883; | |
| − | + | t=0:1:100; | |
| − | + | Q=zeros(1,101); | |
| − | + | P=zeros(1,101); | |
| − | + | for i=1:101 | |
| − | + | Q(1,i)=(C*K*exp(r*t(i)))/(1+C*exp(r*t(i))); | |
| − | + | end | |
| − | + | for j=1:101 | |
| − | + | P(j)=r*Q(j)*log(K/Q(j)); | |
| − | plot( | + | end |
| − | + | plot(Q,P,'-b') | |
| − | + | xlabel('Cantidad de mineral extraída en toneladas') | |
| − | xlabel(' | + | ylabel('Producción de mineral en toneladas/año') |
| − | ylabel(' | + | |
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
| − | + | ||
}} | }} | ||
| − | |||
==Comparación de P(Q) según modelos de Gompertz y Verhlust== | ==Comparación de P(Q) según modelos de Gompertz y Verhlust== | ||
| Línea 159: | Línea 170: | ||
[[Archivo:Q(t) hasta fin de vida util_EULER.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | [[Archivo:Q(t) hasta fin de vida util_EULER.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | ||
| + | |||
| + | En la siguiente gráfica se puede apreciar el desarrollo de la función de cantidad de mineral hasta que llega a una producción menor de 25 toneladas/año, momento en el cual deja de ser rentable la explotación.<br /> | ||
CODIGO MATLAB | CODIGO MATLAB | ||
| Línea 168: | Línea 181: | ||
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | ||
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | ||
| − | t=t0:h:tN; %t tiene 1201 elementos. 1201 | + | t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos. |
| − | %1200 elementos entre 100 años que hemos representado = 12 | + | %1200 elementos entre 100 años que hemos representado = 12 elementos. |
%12 elementos por año, es decir, uno por cada mes. | %12 elementos por año, es decir, uno por cada mes. | ||
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | ||
| Línea 216: | Línea 229: | ||
[[Archivo:Q(t)_hasta_fin_de_vida_util_RUNGE_KUTTA.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | [[Archivo:Q(t)_hasta_fin_de_vida_util_RUNGE_KUTTA.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | ||
| + | |||
| + | La siguiente gráfica representa lo mismo que el apartado anterior, pero utilizando el método de aproximación de Runge Kutta, que es más preciso.<br /> | ||
| + | |||
CODIGO MATLAB | CODIGO MATLAB | ||
| Línea 226: | Línea 242: | ||
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | ||
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | ||
| − | t=t0:h:tN; %t tiene 1201 elementos. 1201 | + | t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos. |
| − | %1200 elementos entre 100 años que hemos representado = 12 | + | %1200 elementos entre 100 años que hemos representado = 12 elementos. 12 elementos por año, es decir, uno por cada mes. |
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | ||
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica | %del problema sustituida en la ec. Gompertz programada segun la sol. analitica | ||
| Línea 274: | Línea 290: | ||
[[Archivo:Q(t)_hasta_fin_de_vida_util_HEUN.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | [[Archivo:Q(t)_hasta_fin_de_vida_util_HEUN.jpg|400px|Comparación de modelo de Verhlust y Gompertz]]<br /> | ||
| + | |||
| + | La gráfica vuelve a mostrar lo mismo que en apartados anteriores, esta vez con el método de Heun que también es más preciso que el método de Euler.<br /> | ||
| + | |||
CODIGO MATLAB | CODIGO MATLAB | ||
| Línea 284: | Línea 303: | ||
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente | ||
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año | ||
| − | t=t0:h:tN; %t tiene 1201 elementos. 1201 | + | t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos. |
| − | %12 elementos por año, es decir, uno por cada mes. | + | %1200 elementos entre 100 años que hemos representado = 12 elementos. |
| + | %12 elementos por año, es decir, uno por cada mes. | ||
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años | ||
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica | %del problema sustituida en la ec. Gompertz programada segun la sol. analitica | ||
| Línea 326: | Línea 346: | ||
}} | }} | ||
| − | ===Comparación de los | + | ===Comparación de los métodos de aproximación=== |
| − | [[Archivo: | + | [[Archivo:P(t)_corte_con_P=25_fin_de_vida_util_RUNGE_KUTTA_vs_HEUN_(no_se_aprecia_la_diferencia).jpg|400px|Comparación de los métodos de Heun y Runge Kutta]]<br /> |
| + | [[Archivo:Q(t)_hasta_fin_de_vida_util_RUNGE_KUTTA_vs_HEUN_(ZOOM).jpg|400px|Ampliación de la gráfica]]<br /> | ||
| + | |||
| + | Con estas dos imágenes podemos ver que los métodos de Runge Kutta y Heun tienen una magnitud de error muy similar, por lo que quedan superpuestas. Al ampliar la imagen vemos que hay una ligera diferencia entre una y otra.<br /> | ||
CODIGO MATLAB | CODIGO MATLAB | ||
| Línea 342: | Línea 365: | ||
legend('Cantidad extraida (Gompertz / Runge Kutta)','Cantidad extraida (Gompertz / Heun)','Location','Best') | legend('Cantidad extraida (Gompertz / Runge Kutta)','Cantidad extraida (Gompertz / Heun)','Location','Best') | ||
hold off | hold off | ||
| − | |||
| − | |||
figure | figure | ||
hold on | hold on | ||
| Línea 354: | Línea 375: | ||
}} | }} | ||
| − | === | + | ===Tendencia de Q en el largo plazo=== |
| + | |||
| + | Ahora aproximaremos el valor al que tiende Q cuando tomamos un valor del tiempo muy grande. | ||
| + | |||
| + | Como sabemos por la teoría, en el método de Gompertz (que es el que vamos a utilizar para calcularlo), cuando aproximamos el valor de la función a estudiar a largo plazo, ésta tiende al valor de K, que en este caso es 10875 toneladas. Lo veremos más claro con la siguiente gráfica: | ||
| + | |||
| + | [[Archivo:Apartado7tinfinito1.jpg|400px|Valor al que tiende Q cuando t tiende a infinito]] | ||
| + | |||
| + | CODIGO MATLAB<br /> | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | %Lo que haremos será calcular valores de Q para distintos valores de t, y | ||
| + | %luego calcular valores de P para esos mismo valores de t. Una vez tengamos | ||
| + | %los valores de P y Q en sendos vectores, los representaremos | ||
| + | K=10875; | ||
| + | A=4.47; | ||
| + | r=0.06; | ||
| + | t=0:1:400; | ||
| + | Q=zeros(1,401); | ||
| + | P=zeros(1,401); | ||
| + | for i=1:401 | ||
| + | Q(1,i)=K*exp(-A*exp(-r*t(i))); | ||
| + | end | ||
| + | plot(t,Q,'-g') | ||
| + | xlabel('Cantidad de mineral extraída en toneladas') | ||
| + | ylabel('Producción de mineral en toneladas/año') | ||
| + | }} | ||
==Función de producción P(t): Modelo de Gompertz == | ==Función de producción P(t): Modelo de Gompertz == | ||
* Utilizando el método de Heun aproximamos la función de producción P(t).<br /> | * Utilizando el método de Heun aproximamos la función de producción P(t).<br /> | ||
* A continuación obtenemos el punto de máxima producción, que es P(t)=240,18 para t=24,18 años.<br /> | * A continuación obtenemos el punto de máxima producción, que es P(t)=240,18 para t=24,18 años.<br /> | ||
| − | + | [[Archivo:P(t)_GOMPERTZ.jpg|400px|Función de producción con respecto al tiempo]]<br /> | |
| − | + | ||
| − | + | Esta gráfica representa la función P(t) según el modelo de Gompertz, donde se aprecia el aumento de la producción en su primera fase, y el descenso posterior, tras alcanzar la máxima producción en el instante t que hemos calculado.<br /> | |
| − | + | ||
| − | + | ||
CODIGO MATLAB<br /> | CODIGO MATLAB<br /> | ||
| Línea 400: | Línea 446: | ||
y=g*tN/N; | y=g*tN/N; | ||
}} | }} | ||
| + | |||
| + | Puntos de crecimiento y decrecimiento:<br /> | ||
| + | |||
| + | [[Archivo:Derivada_de_la_produccion_GOMPERTZ.jpg|400px|Función de producción con respecto al tiempo]]<br /> | ||
| + | |||
| + | La grafica representa la derivada de la función de producción P(t), es decir, P'(t). El máximo y mínimo de esta función son los puntos de mayor crecimiento y decrecimiento de la función P(t).<br /> | ||
===Cantidad de mineral que queda por extraer=== | ===Cantidad de mineral que queda por extraer=== | ||
| Línea 407: | Línea 459: | ||
* En este caso volvemos a utilizar el método de Heun para aproximar la función.<br /> | * En este caso volvemos a utilizar el método de Heun para aproximar la función.<br /> | ||
* El punto de máxima producción que obtenemos es P(t)=240,18 para un tiempo t=50,38 años<br /> | * El punto de máxima producción que obtenemos es P(t)=240,18 para un tiempo t=50,38 años<br /> | ||
| − | + | [[Archivo:P(t)_LOGISTICA.jpg|400px|Función de producción con respecto al tiempo]]<br /> | |
| − | + | ||
| − | Archivo: | + | Esta gráfica representa la función P(t) según el modelo de Verhlust, donde se aprecia el aumento de la producción en su primera fase, y el descenso posterior, tras alcanzar la máxima producción en el instante t que hemos calculado.<br /> |
| − | + | ||
| − | </ | + | |
CODIGO MATLAB:<br /> | CODIGO MATLAB:<br /> | ||
{{matlab|codigo= | {{matlab|codigo= | ||
| Línea 446: | Línea 498: | ||
y=g*tN/N; | y=g*tN/N; | ||
}} | }} | ||
| + | |||
| + | Puntos de crecimiento y decrecimiento:<br /> | ||
| + | |||
| + | [[Archivo:Derivada_de_la_produccion_LOGISTICA.jpg|400px|Derivada de la función de producción]]<br /> | ||
| + | |||
| + | La grafica representa la derivada de la función de producción P(t), es decir, P'(t). El máximo y mínimo de esta función son los puntos de mayor crecimiento y decrecimiento de la función P(t).<br /> | ||
===Cantidad de mineral que queda sin extraer=== | ===Cantidad de mineral que queda sin extraer=== | ||
==Comparación de P(t) según Gompertz y Verhlust== | ==Comparación de P(t) según Gompertz y Verhlust== | ||
| + | |||
| + | [[Archivo:P(t)_LOGISTICAvsGOMPERTZ.jpg|Comparación del modelo logístico y Gompertz]]<br /> | ||
| + | |||
==Aproximación del modelo con nuevos datos== | ==Aproximación del modelo con nuevos datos== | ||
| Línea 459: | Línea 520: | ||
* Para obtener K, simplemente sumamos la cantidad de mineral extraída hasta el momento y la cantidad que queda por extraer, y nos queda <math>K=2695+9075=11770 </math>toneladas. | * Para obtener K, simplemente sumamos la cantidad de mineral extraída hasta el momento y la cantidad que queda por extraer, y nos queda <math>K=2695+9075=11770 </math>toneladas. | ||
| + | Para obtener r, tenemos que programar un bucle que nos vaya asignando valores de r hasta que la diferencia entre el Q a los 12 años estimado y el real sea muy pequeña. | ||
| + | |||
| + | CODIGO MATLAB | ||
| + | |||
| + | {{matlab|codigo= | ||
| + | clear all | ||
| + | t0=1; | ||
| + | tN=12; | ||
| + | q0=124.48; | ||
| + | A=4.47; | ||
| + | t=t0:1:tN; | ||
| + | Q=zeros(1,tN); | ||
| + | Q(1,1)=q0; | ||
| + | K=11770; | ||
| + | r=0.01:0.001:1; | ||
| + | m=length(r); | ||
| + | for i=1:m | ||
| + | Q(1,12)=K*exp(-A*exp(-r(i)*t(12))); | ||
| + | P=2695-Q(1,12); | ||
| + | if P<=0.01 | ||
| + | break | ||
| + | end | ||
| + | end | ||
| + | r(i) | ||
| + | }} | ||
| + | |||
| + | De esta forma, tenemos que el valor de r que mejor se adecua a nuestras necesidades es r=0,093. | ||
===Producción=== | ===Producción=== | ||
* Volvemos a aproximar el modelo utilizando el método de Heun.<br /> | * Volvemos a aproximar el modelo utilizando el método de Heun.<br /> | ||
| Línea 464: | Línea 552: | ||
* La cantidad de mineral que queda sin extraer la obtenemos...??? | * La cantidad de mineral que queda sin extraer la obtenemos...??? | ||
| − | + | [[Archivo:P(t)_11.jpg|Comparación del modelo logístico y Gompertz]]<br /> | |
| − | Archivo: | + | |
| − | </ | + | |
| − | + | ||
CODIGO MATLAB | CODIGO MATLAB | ||
Revisión actual del 16:53 6 mar 2015
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Explotación minera |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | • Jaume Martorell Cerdá • Miguel Angel Serrano Leo • Carla Vázquez Gómara • Pablo Alonso Medina • Joaquín Sánchez Molina • Fernando Millán Cobo |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Relación entre cantidad y producción
- 3 Modelo de Gompertz: Valor de la tasa de crecimiento (r) y resolución de la función de producción
- 4 Modelo de Verhlust: Valor de la tasa de crecimiento (r) y resolución de la función de producción
- 5 Comparación de P(Q) según modelos de Gompertz y Verhlust
- 6 Problema de valor inicial utilizando distintos métodos
- 7 Función de producción P(t): Modelo de Gompertz
- 8 Función de producción P(t): Modelo de Verhlust
- 9 Comparación de P(t) según Gompertz y Verhlust
- 10 Aproximación del modelo con nuevos datos
1 Introducción
- El problema nos pide el análisis de la explotación de un yacimiento de mineral. Dicha explotación sigue un modelo logístico de Gompertz, cuya ecuación tiene la siguiente forma:
:[math] {dQ \over dt}=rQlog ({K \over Q}) [/math]
donde Q(t) es la cantidad de mineral extraído, K la cantidad total extraíble y r la tasa intrínseca de crecimiento.
- En nuestro caso, sabemos que la extracción de mineral tendrá un crecimiento muy rápido de producción durante los primeros 25 años, momento a partir del cual descenderá lentamente debido a diversos factores. Además de esto, conocemos la cantidad total extraíble del yacimiento, por lo que nuestra ecuación inicial quedará de la siguiente forma:
:[math] {dQ \over dt}=rQlog ({10875 \over Q}) [/math]
2 Relación entre cantidad y producción
Mediante las siguientes gráficas, se demuestra que la relación entre Q y P es la siguiente:
- [math] {dQ \over dt}= P [/math]
o lo que es lo mismo:[math] P=Q' [/math]
Para ello, vamos a ver la representación gráfica de cada función.
De esta forma, la evolución de Q respecto al tiempo tendría esta forma:
Y la evolución de P respecto al tiempo es la siguiente:
Podemos observar a simple vista que el valor de P es la pendiente de la función Q(t), ya que vemos que el punto de inflexión de la función Q coincide con el máximo de la función P.
3 Modelo de Gompertz: Valor de la tasa de crecimiento (r) y resolución de la función de producción
3.1 Tasa de crecimiento r
- Partimos del siguiente dato inicial:
- K= 10875 toneladas
- La ecuación del modelo de Gompertz es:
:[math]{dQ \over dt}=rQlog({K \over Q})[/math]
- Condición:
- P(t=25 años)= 240 t/año [math] \rightarrow [/math] P(25)=240 [math] \rightarrow [/math] Q'(25)=240
- Resolución:
La solución de la ecuación de Gompertz es la siguiente:
: [math] Q(t)=Ke^{-Ae^{-rt}} [/math]
El enunciado dice que la producción es máxima a los 25 años, por lo que si Q' tiene un máximo en un punto, Q"=0 en dicho punto.
Derivamos la ecuación inicial del modelo y la igualamos con la condición que hemos obtenido anteriormente.
: [math] Q"=r^2Qln({K \over Q})ln({K \over Q}-1)=0 [/math]
Despejando esta ecuación llegamos a que Q(25 años)= 4000,69 t.
Conociendo Q(25) podemos despejar "r" de nuestro problema:
: [math] Q'= rQ(25)ln({10375 \over Q(25)})=240 [/math]
Así obtenemos r=0,05999 [math] \simeq [/math] 0,06
Con estos datos, ya podemos despejar A de [math] Q(t)=Ke^{-Ae^{-rt}} [/math], siendo A=4,4705
Por lo tanto, nuestra ecuación del modelo de Gompertz quedará de la siguiente forma:
: [math] Q(t)=10875e^{-4,47e^{-0,6t}} [/math]
Y su representación gráfica:
3.2 Función de producción P(Q)
Para obtener la función P(Q), obtendremos por un lado P(t) y por otro Q(t), para después representarlos en la forma P(Q).
La gráfica hace referencia a la producción que se tiene en función de la cantidad de mineral extraída. Se puede ver que el máximo se alcanza para una P=240 toneladas/año cuando la cantidad extraída es Q=4000,69 toneladas.
CODIGO MATLAB
%Lo que haremos será calcular valores de Q para distintos valores de t, y
%luego calcular valores de P para esos mismo valores de t. Una vez tengamos
%los valores de P y Q en sendos vectores, los representaremos
K=10875;
A=4.47;
r=0.06;
t=0:1:100;
Q=zeros(1,101);
P=zeros(1,101);
for i=1:101
Q(1,i)=K*exp(-A*exp(-r*t(i)));
end
for j=1:101
P(j)=r*Q(j)*log(K/Q(j));
end
plot(Q,P,'*g')
xlabel('Cantidad de mineral extraída en toneladas')
ylabel('Producción de mineral en toneladas/año')
4 Modelo de Verhlust: Valor de la tasa de crecimiento (r) y resolución de la función de producción
4.1 Tasa de crecimiento r
- Partimos del siguiente dato inicial:
- K= 10875 toneladas
- La ecuación del modelo de Verhlust es:
- [math]{dQ \over dt}=rQ(1-{Q \over K})[/math]
- Condición:
- P(t=25 años)= 240 t/año [math] \rightarrow [/math] P(25)=240 [math] \rightarrow [/math] Q'(25)=240
- Resolución:
Para la resolución, utilizamos el mismo planteamiento que en el apartado anterior del modelo de Gompertz.
: [math] Q"=rQ'(1-{Q \over 5437,5})=0 [/math]
En el caso del modelo de Verhlust obtenemos r=0,0883 y Q(25 años)=5437,5 t. Así, la ecuación del modelo tendrá la siguiente forma:
: [math] Q'=0,0883Q(1-{Q \over 10875}) [/math]
Y la representación gráfica queda de la siguiente forma:
4.2 Función de producción P(Q)
Como sabemos que P=Q', la ecuación de P(t) es la siguiente:
: [math] P(t)=0,0883Q(t)(1-{Q(t) \over 10875})=0 [/math]
La ecuación de Q(t) la hemos obtenido integrando la del modelo obtenida anteriormente, y tiene la forma:
: [math] Q(t)={1196,25e^{0.088276t} \over 1+11e^{0.088276t}} [/math]
Al igual que en el apartado anterior, obtenemos P(Q) dibujando ambas funciones en la misma gráfica:
Vemos que la gráfica empieza en 1000 toneladas porque en caso contrario tendríamos una indeterminación a la hora de calcular los valores.
CODIGO MATLAB
%Lo que haremos será calcular valores de Q para distintos valores de t, y
%luego calcular valores de P para esos mismo valores de t. Una vez tengamos
%los valores de P y Q en sendos vectores, los representaremos
K=10875;
C=0.11;
r=0.0883;
t=0:1:100;
Q=zeros(1,101);
P=zeros(1,101);
for i=1:101
Q(1,i)=(C*K*exp(r*t(i)))/(1+C*exp(r*t(i)));
end
for j=1:101
P(j)=r*Q(j)*log(K/Q(j));
end
plot(Q,P,'-b')
xlabel('Cantidad de mineral extraída en toneladas')
ylabel('Producción de mineral en toneladas/año')
5 Comparación de P(Q) según modelos de Gompertz y Verhlust
_LOGISTICAvsGOMPERTZ.jpg)
La gráfica muestra la comparación entre el modelo logístico de Verhlust y el modelo de Gompertz. Nos quedaremos con el modelo de Gompertz ya que como podemos observar en las gráficas se ajusta más a nuestro modelo, ya que la producción inicial se aproxima más al valor nulo al comenzar el modelo.
CODIGO MATLAB
%Comparativa LOGISTICO vs GOMPERTZ
t0=0;
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente
h=0.01;
t=t0:h:tN;
%MODELO LOGISTICO
Q_log=(1196.25.*exp(0.088276.*t))./(1+0.11.*exp(0.088276.*t));
P_log=0.088276.*Q_log.*(1-Q_log./10875);
%MODELO GOMPERTZ
Q_gom=10875.*exp(-4.47.*exp(-0.06.*t));
P_gom=0.05999.*Q_gom.*log(10875./Q_gom);
%Comparativa de cantidades extraidas
figure
hold on
plot(t,Q_log,'g')
plot(t,Q_gom)
plot(t,t.*0+10875,'r')
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad extraida (modelo logistico)','Cantidad extraida (modelo Gompertz)','Cantidad total extraible (limite)','Location','Best')
hold off
%Comparativa de producciones
figure
hold on
6 Problema de valor inicial utilizando distintos métodos
6.1 Método de Euler
_hasta_fin_de_vida_util_EULER.jpg)
En la siguiente gráfica se puede apreciar el desarrollo de la función de cantidad de mineral hasta que llega a una producción menor de 25 toneladas/año, momento en el cual deja de ser rentable la explotación.
CODIGO MATLAB
r=0.06;
k=10875;
t0=0;
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año
t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos.
%1200 elementos entre 100 años que hemos representado = 12 elementos.
%12 elementos por año, es decir, uno por cada mes.
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica
Q=zeros(1,length(t));
P=zeros(1,length(t));
Q(1,1)=Q0;
P(1,1)=r*Q(1)*log(k/Q(1));
for i=1:length(t)-1
Q(1,i+1)=Q(1,i)+h*r*Q(1,i)*log(k/Q(1,i));
P(1,i+1)=r*Q(1,i+1)*log(k/Q(1,i+1));
end
for n=1:length(t)-1
P(1,n);
if P(1,n)>=25
Plim=P(1,n);
m=n;
end
end
Plim %Plim es el ultimo valor de la produccion que esta por encima de las 25t/año
m %Paso correspondiente al valor de Plim, se puede comprobar que P(m)=Plim=25t/año
%Ahora que ya conocemos el paso de tiempo en el que finaliza la vida util
%de la explotacion (m=787), se puede reescribir la funcion Q(t), que valdra
%0 a partir de dicho paso de tiempo
for j=m:length(t)-1
Q(1,j+1)=0;
end
figure
plot(t,Q(1,:))
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad de material extraido (Gompertz / Euler)','Location','Best')
figure
hold on
plot(t,P(1,:))
plot(t,t.*0+25,'r')
xlabel('Tiempo en años')
ylabel('Produccion en toneladas/año')
legend('Produccion (Gompertz / Euler)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')
hold off
6.2 Método de Runge-Kutta
_hasta_fin_de_vida_util_RUNGE_KUTTA.jpg)
La siguiente gráfica representa lo mismo que el apartado anterior, pero utilizando el método de aproximación de Runge Kutta, que es más preciso.
CODIGO MATLAB
clear all
r=0.06;
k=10875;
t0=0;
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año
t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos.
%1200 elementos entre 100 años que hemos representado = 12 elementos. 12 elementos por año, es decir, uno por cada mes.
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica
%Runge Kutta
Q_rk=zeros(1,length(t));
P_rk=zeros(1,length(t));
Q_rk(1,1)=Q0;
P_rk(1,1)=r*Q_rk(1)*log(k/Q_rk(1));
%Aproximacion por Runge Kutta
for i=1:length(t)-1
K1=r*Q_rk(1,i)*log(k/Q_rk(1,i));
K2=r*(Q_rk(1,i)+0.5*K1*h)*log(k/(Q_rk(1,i)+0.5*K1*h));
K3=r*(Q_rk(1,i)+0.5*K2*h)*log(k/(Q_rk(1,i)+0.5*K2*h));
K4=r*(Q_rk(1,i)+K3*h)*log(k/(Q_rk(1,i)+K3*h));
Q_rk(1,i+1)=Q_rk(1,i)+h/6*(K1+2*K2+2*K3+K4);
P_rk(1,i+1)=r*Q_rk(1,i+1)*log(k/Q_rk(1,i+1));
end
for n=1:length(t)-1
P_rk(1,n);
if P_rk(1,n)>=25
Plim=P_rk(1,n);
m=n;
end
end
for j=m:length(t)-1
Q_rk(1,j+1)=0;
end
%Representacion Runge Kutta
figure
plot(t,Q_rk(1,:),'c')
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad extraida (Gompertz / Runge Kutta)','Location','Best')
figure
hold on
plot(t,P_rk(1,:),'c')
plot(t,t.*0+25,'r')
xlabel('Tiempo en años')
ylabel('Produccion en toneladas/año')
legend('Produccion (Gompertz / Runge Kutta)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')
hold off
6.3 Método de Heun
_hasta_fin_de_vida_util_HEUN.jpg)
La gráfica vuelve a mostrar lo mismo que en apartados anteriores, esta vez con el método de Heun que también es más preciso que el método de Euler.
CODIGO MATLAB
clear all
r=0.06;
k=10875;
t0=0;
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año
t=t0:h:tN; %t tiene 1201 elementos. 1201 elementos - 1 elemento inicial = 1200 elementos.
%1200 elementos entre 100 años que hemos representado = 12 elementos.
%12 elementos por año, es decir, uno por cada mes.
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica
%Heun
Q_heun=zeros(1,length(t));
P_heun=zeros(1,length(t));
Q_heun(1,1)=Q0;
P_heun(1,1)=r*Q_heun(1)*log(k/Q_heun(1));
%Aproximacion por Heun
for i=1:length(t)-1
k1=r*Q_heun(1,i)*log(k/Q_heun(1,i));
k2=r*(Q_heun(1,i)+k1*h)*log(k/(Q_heun(1,i)+k1*h));
Q_heun(1,i+1)=Q_heun(1,i)+h/2*(k1+k2);
P_heun(1,i+1)=r*Q_heun(1,i+1)*log(k/Q_heun(1,i+1));
end
for n=1:length(t)-1
P_heun(1,n);
if P_heun(1,n)>=25
Plim=P_heun(1,n);
m=n;
end
end
for j=m:length(t)-1
Q_heun(1,j+1)=0;
end
%Representacion Heun
figure
plot(t,Q_heun(1,:),'k')
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad extraida (Gompertz / Heun)','Location','Best')
figure
hold on
plot(t,P_heun(1,:),'k')
plot(t,t.*0+25,'r')
xlabel('Tiempo en años')
ylabel('Produccion en toneladas/año')
legend('Produccion (Gompertz / Heun)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')
hold off
6.4 Comparación de los métodos de aproximación
_corte_con_P%3D25_fin_de_vida_util_RUNGE_KUTTA_vs_HEUN_(no_se_aprecia_la_diferencia).jpg)
_hasta_fin_de_vida_util_RUNGE_KUTTA_vs_HEUN_(ZOOM).jpg)
Con estas dos imágenes podemos ver que los métodos de Runge Kutta y Heun tienen una magnitud de error muy similar, por lo que quedan superpuestas. Al ampliar la imagen vemos que hay una ligera diferencia entre una y otra.
CODIGO MATLAB
%Representacion simultanea de Runge Kutta y Heun
figure
hold on
plot(t,Q_rk(1,:),'c')
plot(t,Q_heun(1,:),'k')
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad extraida (Gompertz / Runge Kutta)','Cantidad extraida (Gompertz / Heun)','Location','Best')
hold off
figure
hold on
plot(t,P_rk(1,:),'c')
plot(t,P_heun(1,:),'k')
plot(t,t.*0+25,'r')
xlabel('Tiempo en años')
ylabel('Produccion en toneladas/año')
legend('Produccion (Gompertz / Runge Kutta)','Produccion (Gompertz / Heun)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')
6.5 Tendencia de Q en el largo plazo
Ahora aproximaremos el valor al que tiende Q cuando tomamos un valor del tiempo muy grande.
Como sabemos por la teoría, en el método de Gompertz (que es el que vamos a utilizar para calcularlo), cuando aproximamos el valor de la función a estudiar a largo plazo, ésta tiende al valor de K, que en este caso es 10875 toneladas. Lo veremos más claro con la siguiente gráfica:
CODIGO MATLAB
%Lo que haremos será calcular valores de Q para distintos valores de t, y
%luego calcular valores de P para esos mismo valores de t. Una vez tengamos
%los valores de P y Q en sendos vectores, los representaremos
K=10875;
A=4.47;
r=0.06;
t=0:1:400;
Q=zeros(1,401);
P=zeros(1,401);
for i=1:401
Q(1,i)=K*exp(-A*exp(-r*t(i)));
end
plot(t,Q,'-g')
xlabel('Cantidad de mineral extraída en toneladas')
ylabel('Producción de mineral en toneladas/año')
7 Función de producción P(t): Modelo de Gompertz
- Utilizando el método de Heun aproximamos la función de producción P(t).
- A continuación obtenemos el punto de máxima producción, que es P(t)=240,18 para t=24,18 años.
_GOMPERTZ.jpg)
Esta gráfica representa la función P(t) según el modelo de Gompertz, donde se aprecia el aumento de la producción en su primera fase, y el descenso posterior, tras alcanzar la máxima producción en el instante t que hemos calculado.
CODIGO MATLAB
clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=10875;
r=0.06;
q(1)=q0;
p=zeros(1,N+1);
p(1)=r*q0*log(k/q0);
for j=1:N
q(j+1)=q(j)+h*r*q(j)*log(k/q(j));
end
for i=1:N
k1=r^2*q(i)*log(k/q(i))*(log(k/q(i))-1);
k2=r^2*(q(i)+k1*h)*log(k/(q(i)+k1*h))*(log(k/(q(i)+k1*h))-1);
p(i+1)=p(i)+h/2*(k1+k2);
end
[t',p'];
plot(t,p)
%Obtenemos el año de máxima producción
max(p)
for g=1:N
if max(p)-p(g)<=0.01
break
end
end
y=g*tN/N;
Puntos de crecimiento y decrecimiento:
La grafica representa la derivada de la función de producción P(t), es decir, P'(t). El máximo y mínimo de esta función son los puntos de mayor crecimiento y decrecimiento de la función P(t).
7.1 Cantidad de mineral que queda por extraer
8 Función de producción P(t): Modelo de Verhlust
- En este caso volvemos a utilizar el método de Heun para aproximar la función.
- El punto de máxima producción que obtenemos es P(t)=240,18 para un tiempo t=50,38 años
_LOGISTICA.jpg)
Esta gráfica representa la función P(t) según el modelo de Verhlust, donde se aprecia el aumento de la producción en su primera fase, y el descenso posterior, tras alcanzar la máxima producción en el instante t que hemos calculado.
CODIGO MATLAB:
clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=10875;
r=0.0883;
q(1)=q0;
p=zeros(1,N+1);
p(1)=r*q0*(1-(q0/k));
for j=1:N
q(j+1)=q(j)+h*r*q(j)*(1-(q(j)/k));
end
for i=1:N
k1=r*p(i)*(1-(2*q(i)/k));
k2=r*(p(i)*(1-(2*(q(i)+k1*h)/k)));
p(i+1)=p(i)+h/2*(k1+k2);
end
[t',p'];
plot(t,p)
%Obtenemos el año de máxima producción
max(p)
for g=1:N
if max(p)-p(g)<=0.01
break
end
end
y=g*tN/N;
Puntos de crecimiento y decrecimiento:
La grafica representa la derivada de la función de producción P(t), es decir, P'(t). El máximo y mínimo de esta función son los puntos de mayor crecimiento y decrecimiento de la función P(t).
8.1 Cantidad de mineral que queda sin extraer
9 Comparación de P(t) según Gompertz y Verhlust
_LOGISTICAvsGOMPERTZ.jpg)
10 Aproximación del modelo con nuevos datos
- El enunciado nos dice que transcurridos 12 años del inicio de la explotación, se hace una revisión del comportamiento del modelo y se ajustan los datos de forma más real.
- Ahora, la cantidad extraída hasta ese año 12 es Q(12)=2695 toneladas, y se estima que quedan por extraer en la mina 9075 toneladas más.
- Con estos nuevos datos, calculamos los parámetros r y K, y aproximamos la función P(t).
10.1 Obtención de los parámetros r y K
- Para obtener K, simplemente sumamos la cantidad de mineral extraída hasta el momento y la cantidad que queda por extraer, y nos queda [math]K=2695+9075=11770 [/math]toneladas.
Para obtener r, tenemos que programar un bucle que nos vaya asignando valores de r hasta que la diferencia entre el Q a los 12 años estimado y el real sea muy pequeña.
CODIGO MATLAB
clear all
t0=1;
tN=12;
q0=124.48;
A=4.47;
t=t0:1:tN;
Q=zeros(1,tN);
Q(1,1)=q0;
K=11770;
r=0.01:0.001:1;
m=length(r);
for i=1:m
Q(1,12)=K*exp(-A*exp(-r(i)*t(12)));
P=2695-Q(1,12);
if P<=0.01
break
end
end
r(i)
De esta forma, tenemos que el valor de r que mejor se adecua a nuestras necesidades es r=0,093.
10.2 Producción
- Volvemos a aproximar el modelo utilizando el método de Heun.
- El punto de máxima producción que obtenemos es P(t)=400,87 toneladas para el tiempo t=16,33 años.
- La cantidad de mineral que queda sin extraer la obtenemos...???
_11.jpg)
CODIGO MATLAB
clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=11770;
r=0.0925;
q(1)=q0;
p=zeros(1,N+1);
p(1)=r*q0*log(k/q0);
for j=1:N
q(j+1)=q(j)+h*r*q(j)*log(k/q(j));
end
for i=1:N
k1=r^2*q(i)*log(k/q(i))*(log(k/q(i))-1);
k2=r^2*(q(i)+k1*h)*log(k/(q(i)+k1*h))*(log(k/(q(i)+k1*h))-1);
p(i+1)=p(i)+h/2*(k1+k2);
end
[t',p'];
figure
plot(t,p)
figure
plot(t,q)
%Obtenemos el año de máxima producción
max(p)
for g=1:N
if max(p)-p(g)<=0.01
break
end
end
x=g*tN/N;
10.3 Comparación del nuevo modelo con la previsión inicial
Gráfica comparativa
