Desintegración Radiactiva.Grupo 5
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Desintegración Radiactiva. (Grupo 5-C). |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Diego García Vaquero
Araceli Martín Candilejo Noemí Palomino Bustos Teresa Quintana Romero Alvaro Ramón López Mercedes Ruiz Barrajón. |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Interpretación de las funciones M(t) y la constante K
Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función [math]m=m(t)[/math], la velocidad de descomposición viene dada por la derivada de [math]M(t)[/math] respecto de [math]t[/math]. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que [math]M' = -k.M[/math] ([math]k\gt0[/math] coeficiente de proporcionalidad).
La constante [math]k[/math] se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.
2 Método Euler
Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del [math]8[/math]% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es [math]1.24·10^{-4}[/math] por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos [math]h=0.1[/math] y [math]h=0.01[/math].
clear all
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1
% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.
t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;
% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end
disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Tiempo (años)');
ylabel('Cantidad de Carbono 14');
legend('Euler explícito','Location','best');
El método de Euler suele funcionar muy bien con pasos pequeños, sin embargo no nos asegura que el método sea convergente, dependerá de los factores consistencia y estabilidad.
La consistencia trata en verificar que el error que damos en un paso (calcular [math]y_{n+1}[/math] a partir de [math]y_n[/math]) sea pequeño suponiendo que [math]y_n[/math] es exacto.
La estabilidad consiste de controlar la posible amplificación del error por el hecho de que en cada paso del método el valor de [math]y_n[/math] no es exacto sino una aproximación de [math]y(t_n)[/math]
Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial[math](C_0)[/math] de la muestra de C14:
El método de Euler es convergente para las ecuaciones que verifiquen que la función es diferenciable. En nuestro caso,tenemos que [math]f(t,M)[/math] es diferenciable,por lo tanto el método es convergente.
3 Resolución por el método del Trapecio
Resolveremos el apartado anterior con el método del trapecio utilizando un paso [math]h=0.1[/math] y con la condición de que el programa se detenga cuando la masa alcance un [math]8[/math]% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.
clear all
%% Trapecio
% y´=f(t,y);
% y(t0)=y0;
% y(n+1)=yn+(h/2)*(f(tn,yn)+f(t(n+1),y(n+1)))
% Hay que despejar manualmente para tener SOLO y(n+1) a un lado de la ecuación.
% Si representamos la concentración del elemento radiactivo M(t) como y(t):
% M'(t)=-kM(t)---> y'(t)=-ky(t)
% En este caso, cuando se despeja manualmente:
% y(n+1)=y(n)+(h/2)(-ky(n)-ky(n+1))
% (1+kh/2)y(n+1)=y(n)-khy(n)/2
% y(n+1)=(y(n)-khy(n)/2)/(1+kh/2)
%% SOLUCIÓN:
t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;
% Constante de desintegración:
k=1.24*10^(-4);
% Bucle:
i=1;
while M(i)>(0.08*M0)
M(i+1)=(M(i)-(k*h*M(i))/2)/(1+(k*h)/2); %Trapecio
t(i+1)=t(i)+h;
t(i+1)=t(i)+h;
i=i+1;
end
disp('Tiempo para que quede un 8% de la cantidad inicial: ')
disp(t(end))
hold on
plot(t,M);
xlabel('Tiempo (años)');
ylabel('Cantidad de Carbono 14');
legend('Trapecio','Location','best');
El tiempo para el cual queda el 8% del contenido inicial de Carbono 14 en ambos métodos es [math]20369[/math] años (redondeando, quedan soluciones similares).
El método del trapecio nos garantiza estabilidad incondicional, este método es más preciso que el de Euler, sin embargo es más complicado trabajar con él ya que debemos despejar [math]y_{n+1}[/math].
4 Vida media. Resolución por el método de Runge-Kutta
El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.
clear all
clear all
t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;
%Se inicializa "M" con M0 para irlo rellenando posteriormente con las soluciones obtenidas
%para cada instante mediante el bucle while. Simultáneamente se va ampliando el vector de tiempos
%hasta que se cumple la condición deseada.
%RUNGE-KUTTA:
%K1=f( tn , yn )
%K2=f( tn+h/2, yn+K1*h/2 );
%K3=f( tn+h/2, yn+K2*h/2 );
%K4=f( tn+h, yn+K3*h );
%y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
%En nuestro caso: y'=f(t,y)=-ky
k=1.24*10^(-4);
i=1;
while M(i)>(0.5*M0)
%K1=f( tn , yn )
K1=-k*M(i);
%Definicion de variable K2
%K2=f( tn+h/2, yn+K1*h/2 );
t2=t(i)+(h/2);
M2=M(i)+(h/2)*K1;
K2=-k*M2;
%Definicion de variable K3;
%K3=f( tn+h/2, yn+K2*h/2 );
t3=t2;
M3=M(i)+(h/2)*K2;
K3=-k*M3;
%Definicion de variable K4;
%K4=f( tn+h, yn+K3*h );
t4=t(i)+h;
M4=M(i)+h*K3;
K4=-k*M4;
%Funcion de RungeKutta;
%y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
M(i+1)=M(i)+(h/6)*(K1+2*K2+2*K3+K4);
t(i+1)=t(i)+h;
i=i+1;
end
disp('Tiempo medio: ')
disp(t(end))
hold on
plot(t,M)
xlabel('Tiempo (años)')
ylabel('Cantidad de Carbono 14')
legend('Runge Kutta','Location','best')
hold off
La vida media que da el programa es de [math]5592,51[/math] años.
5 Estudio de la evolución de concentraciones de compuestos interrelacionados
Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[math]A \to _{k1} B \to _{k2} C[/math]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
\begin{array}{c} A'=-K_1A \\ B'=-K_2B+K_1A \end{array}
Con estas ecuaciones y conociendo los valores iniciales, resolveremos este problema mediante Euler y el método del trapecio
5.1 Sistemas de ecuaciones lineales: Método de Euler y trapecio
Le daremos unos valores a las constantes para resolver nuestro sistema, siendo [math]K_1=5[/math], y [math]K_2=1[/math]
Método de Euler
Resolvemos el sistema anterior con el método de Euler:
clc
clf
clear all
%% Datos iniciales:
% Tiempo:
t0=0;
tN=10;
% Paso:
h=0.1;
% Número de subintervalos
N=(tN-t0)/h;
% Definimos la variable independiente: El vector tiempo
t=t0:h:tN;
y=zeros(2,N+1);
% Valores de las concentraciones iniciales:
A0=1;
B0=0;
C0=0;
% Las soluciones se recogerán en un vector "y" que irá autoformándose.
% Inicialmente, sólo se le asignará los valores de las concentraciones iniciales desde las que parte.
y0=[A0,B0]';
y(:,1)=y0;
% Un sistema Lineal:
% y'=My+T
% En este problema en particular:
% A=(-k1 0)(A)+(0)
% B=( k1 -k2)(B)+(0)
% M=(-k1 0)
% ( k1 -k2)
k1=5;
k2=1;
M=[-k1,0;k1,-k2];
% Bucle:
for i=1:N
y(:,i+1)=y(:,i)+h*(M*y(:,i)); % Euler
end
% Se reasigna cada parte del vector "y" que recoge las soluciones con las concentraciones que representan:
A=y(1,:);
B=y(2,:);
C=C0+A0-A-B;
%% Representación gráfica
hold on
plot(t,A)
plot(t,B,'r')
plot(t,C,'g')
xlabel('Tiempo (años)')
ylabel('Cantidad de los elementos')
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off
Método trapecio
Ahora aplicaremos el método del trapecio para resolver el sistema:
clc
clf
clear all
%% Datos iniciales:
% Tiempo:
t0=0;
tN=10;
% Paso:
h=0.1;
% Número de subintervalos
N=(tN-t0)/h;
% Definimos la variable independiente: El vector tiempo
t=t0:h:tN;
y=zeros(2,N+1);
% Valores de las concentraciones iniciales:
A0=1;
B0=0;
C0=0;
% Las soluciones se recogerán en un vector "y" que irá construyéndose.
% Inicialmente, sólo se le asignará los valores de las concentraciones iniciales desde las que parte.
y0=[A0,B0]';
y(:,1)=y0;
% Un sistema Lineal:
% y'=My+T
% En este problema en particular:
% A=(-k1 0)(A)+(0)
% B=( k1 -k2)(B)+(0)
% M=(-k1 0)
% ( k1 -k2)
k1=5;
k2=1;
M=[-k1,0;k1,-k2];
% Bucle:
% Para el Trapecio:
% y(n+1)=y(n)+h/2*(f(tn,yn)+f(t(n+1),y(n+1))
% En este caso:
% y(n+1)=y(n)+h/2*(M*yn+M*y(n+1))
% Despejando manualmente:
% [1-(h/2)*M]*y(n+1)=y(n)+h/2*(M*yn)
% Llamando Z a:
% Z=[1-(h/2)*M]
% Z*(INV(z))*y(n+1)=(INV(z))*[y(n)+h/2*(M*yn)]
% Y queda finalmente:
% y(n+1)=(INV(z))*[y(n)+h/2*(M*yn)]
for i=1:N
Z=eye(2)-(h/2)*M;
y(:,i+1)=inv(Z)*(y(:,i)+(h/2)*M*y(:,i));%Trapecio
end
% Se reasigna cada parte del vector "y" que recoge las soluciones con las concentraciones que representan:
A=y(1,:);
B=y(2,:);
C=C0+A0-A-B;
%% Representación gráfica
hold on
plot(t,A)
plot(t,B,'r')
plot(t,C,'g')
xlabel('Tiempo (años)')
ylabel('Cantidad de los elementos')
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off
Podemos ver que, el compuesto A decrece, el C crece y el B empieza creciendo para más tarde decrecer y desaparecer. Hemos tomado como tiempo el eje de abscisas, y la cantidad del como eje de ordenadas (hemos tomado en tanto por uno la cantidad). A decrece mas rápidamente que el resto, es decir, desaparece antes, B llega a una cantidad de [math]0.6[/math]. En conclusión el compuesto A desaparece de forma rápida y esto da lugar a que haya más cantidad de B.
5.2 Sistemas de ecuaciones lineales: Método de Euler y trapecio (Constantes de integración intercambiadas)
Ahora, resolveremos el sistema pero intercambiando los valores de las constantes, es decir [math]K_1=1[/math] y [math]K_2=5[/math] y aplicando los mismos programas anteriores, Euler y trapecio.
Método de Euler
Método del trapecio
En todas las gráficas, con las K iniciales y con las K intercambiadas y con ambos métodos, podemos observar que el compuesto A decrece, el C crece y el B empieza creciendo para más tarde decrecer y desaparecer. En esas últimas gráficas también hemos tomado como tiempo el eje de abscisas, y en tanto por uno la cantidad del compuesto como eje de ordenadas. Podemos ver que la cantidad en tanto por uno de B sin cambiar las constantes llega a una cantidad de [math]0.6[/math], mientras que con el cambio de constantes no llega a [math]0.2[/math]. En conclusión, al comparar los dos pares de gráficas observamos que ahora con el cambio de constantes A decrece más lentamente, y debido a eso la cantidad de B que se forma no llega a ser tan grande como la del apartado anterior, en el que el compuesto A desaparecía antes y eso daba lugar a más cantidad de B. La creación de C no difiere mucho entre los pares de gráficas, lo que significa que el tiempo que tarda A en descomponerse en C es independiente de las constantes de desintegración.
La solución de las Ecuaciones diferenciales por medio de métodos numéricos involucra varios tipos de errores:
1.Error del Método (Error de Truncamiento Local y Global): Este se debe a que, cómo la aproximación de una curva mediante una línea recta no es exacta, se comete un error propio del método.
2.Local: Es la diferencia que se produce entre el valor real de la función y el aproximado mediante la recta tangente -en lugar de moverse por la curva- suponiendo que el punto desde el que partimos -donde se cruzan la curva real y la recta que la aproxima- no tiene error alguno.
3.Propagado: Acumulación de errores por las aproximaciones producidas durante los pasos previos acumuladas. Es decir, ya no se supone que el punto del cual partimos -donde se cruzan la curva real y la recta que la aproxima- no tenía error sino que asumimos que dicho error existe y que se propaga de paso en paso. Dicha propagación es, en el peor de los casos, lineal.
El método de Euler tiene orden 1, el del trapecio orden 2 y Runge Kutta tiene orden 4.
