Desintegración Radioactiva (G.2A)

Trabajo realizado por estudiantes
Título	Desintegración Radiactiva. (Grupo 5-C).
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Juan Raúl Ruiz Méndez (531) Jaime Enrech Martínez (532) Jose Manuel Alonso de Caso Gilsanz (618) Guillermo Díaz Rivera (649) Iago Rodríguez Romero (824).
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se conoce la propensión de algunos compuestos e isótopos a transformarse en otros mas estables con el paso del tiempo, como es el caso de Carbono 14, objeto de este trabajo. Partiendo de una concentración inicial de dicho isotopo con cierto grado de inestabilidad, sabemos entonces que se producirá una perdida de cantidad del material con el paso del tiempo, llegando un punto en el que el material se haya convertido en su totalidad en compuestos mas estables.

El Carbono14 es un isótopo presente en estado de desintegración en muestras arqueológicas, y es comúnmente utilizado para, en función del grado de disociación en el que se encuentre, fechar la muestra que lo porta.

De esta manera, sabiendo que M(t) representa la cantidad en función del tiempo, expresaremos la velocidad de desintegración como la derivada con respecto a [math]t[/math] de dicha función; [math]\frac{\operatorname dM(t)}{\operatorname dt} = M'(t) [/math] y la expresaremos como [math]M'(t) = -k M(t)[/math] donde [math]k=1,24e^-4[/math] es la constante de desintegración característica del Carbono 14.

1 Datación de Muestras Arqueológicas

Habiéndose realizado la prueba del Carbono 14 en una muestra, se obtiene que ésta posee una concentración M de este elemento del 8% del total que es natural en un ser vivo, habiendose desintegrado por lo tanto un 92% de la cantidad de Carbono 14 original. Conocida la concentración del isotopo presente en la muestra, se procede a calcular su edad (el tiempo que ha tardado en desintegrarse el 92% del contenido de C14). Para ello se utiliza la ecuación diferencial ya expresada y se resuelve mediante el método numérico de Euler.

derecha

clear all
%Datos iniciales del problema
y0=0.08;
h=input('Introducir tamaño del paso: ');
h=-h;
t0=0;
i=1;
y(1)=y0;
t(1)=t0;
k=1.24*10^-4;
%Bucle generación vector 't' y método Euler
while y(i)<1
   tN=t0+i*h;
   t=t0:h:tN;
   y(i+1)=y(i)+h*(-k*y(i));
   i=i+1;
end
   a=length(y);
  
   y(a)
   z=t(a-1);
   z=-z
   plot(t,y,'r')
   legend('Cantidad de C14 (Método de Euler)','Location','Best')

Tras ejecutar nuestro programa observamos que el carbono 14 ha llegado al 8% tras un periodo de 20.369 años. Esto queda reflejado en la gráfica adjunta.

2 Previsión de desintegración de Carbono 14

Con la ecuación diferencial ya comentada, no solo se puede calcular la antiguedad de una muestra, sino que también permite hacer una estimación de cuando se llegará a un nivel de Carbono 14 concreto. El próximo cálculo realizado determinará cuanto tiempo tarda la muestra de Carbono 14 restante (8%) en reducirse hasta el 8% (es decir, el 0,64% del total). Para ello se utiliza el método del trapecio (Método implícito).

derecha

clear all
%Datos iniciales del problema
y0=0.08;
h=input('Introducir tamaño del paso: ');
t0=0;
i=1;
y(1)=y0;
t(1)=t0;
k=1.24*10^-4;
%Bucle generación vector 't' y método trapecio
while y(i)>0.0064
   tN=t0+i*h;
   t=t0:h:tN;
   y(i+1)=(y(i)-(h/2)*k*y(i))/(1+(h/2)*k);
%Expresión obtenida de despejar y sub(n+1)
% de la fórmula del método del trapecio
   i=i+1;
end
   a=length(y);
   y(a)
   z=t(a-1)
   plot(t,y,'r')
   legend('Cantidad de C14','Location','Best')

El resultado de este código nos muestra que el tiempo de descomposición hasta llegar al 0.64% del total es prácticamente el mismo que en el caso anterior (20368)

Ésto ocurre por conceptos de desintegracion atómica, en los que se explica que el tiempo que tarda en descomponerse una cantidad de un elemento, se produce en función de, exclusivamente el porcentaje de materia desintegrada, sin tener en cuenta cual fuese la cantidad de la que se parte para este proceso. Esto significa que si originalmente (a día de hoy) tenemos 8% de la cantidad del animal vivo hace 20.000, ésta ha tardado una cantidad T de años. Y por esa misma razón si de una cantidad 8% desintegramos un nuevo 8% estamos actuando como si la cantidad a dia de hoy fuese un nuevo total a desintegrar. Demostrandose entonces que independientemente de la masa de la que se parta, un porcentaje dado a disociar tarda siempre la misma cantidad de tiempo T.

3 Periodo de Semidesintegración

En física nuclear se conoce como periodo de semidesintegración al espacio de tiempo que se requiere para que se desintegren la mitad de los núcleos de una muestra inicial de radioisotopos. Utilizando el método de Runge-Kutta, y teniendo en cuenta que no se conoce el valor final del intervalo de tiempo (puesto que es la incognita que buscamos), se ejecuta un bucle while que se interrumpe cuando la cantidad de C14 presente en la muestra se ha reducido a la mitad

derecha

clear all
%Datos iniciales del problema
y0=1;
h=input('Introducir tamaño del paso: ');
t0=0;
i=1;
y(1)=y0;
t(1)=t0;
k=1.24*10^-4;
%Bucle generación vector 't' y método RK-4
while y(i)>0.5
   tN=t0+i*h;
   t=t0:h:tN;
   K1=-k*y(i);
   K2=-k*(y(i)+(1/2)*K1*h);
   K3=-k*(y(i)+(1/2)*K2*h);
   K4=-k*(y(i)+K3*h);
   y(i+1)=y(i)+(h/6)*(K1+2*K2+2*K3+K4);
   i=i+1;
end
   a=length(y);
   y(a)
   z=t(a-1);
   z
   plot(t,y,'r')
   legend('Cantidad de C14','Location','Best')

La vida media del elemento radiactivo es 5589 años. Con este resultado apreciamos que el tiempo de desintegración no es lineal, ya que en los primeros años se reduce mucho más rápido.

4 Descomposición Compuesta

Existen casos en los que la reacción entre un elemento A y otro C se produce a través de un elemento intermedio B. Esta reacción compuesta se puede asimilar a un sistema de ecuaciones diferenciales planteado a partir de la ecuación de desintegración ya conocida. El sistema expresado de forma matricial será:

[math]M'(x) = \begin{Bmatrix} M'(x)_A\\ M'(x)_B \end{Bmatrix} = \begin{bmatrix} -k_1 & 0 \\ -k_1 & -k_2 \end{bmatrix} \begin{Bmatrix} M_A\\ M_B \end{Bmatrix} ;\qquad \mathbf{X}(0) = \begin{Bmatrix} 1\\ 0 \end{Bmatrix}[/math]

Siendo [math]k_1[/math] y [math]k_2[/math] las constantes de desintegración de las correspondientes reacciones.

4.1 euler

derecha

%Datos iniciales
t0=0;
tN=10;
h=0.1;
N=round((tN-t0)/h);
y0=[1;0];
%Generación del vector de tiempo 't' e 'y'
t=t0:h:tN;
y=zeros(2,N+1);
y(:,1)=y0;
c(1)=0;
A=[-5 0;5 -1]; %Matriz de coeficientes del sistema
%Bucle Euler
for i=1:N
    y(:,i+1)=y(:,i)+h*(A*y(:,i));
    c(i+1)=1-y(1,i)-y(2,i);
end
plot(t,y(1,:))
hold on
plot(t,y(2,:),'r')
plot(t,c,'g')
legend('[A]','[B]','[C]','location','best')

4.2 trapecio

derecha

%Datos iniciales
t0=0;
tN=10;
h=0.1;
N=round((tN-t0)/h);
y0=[1;0];
%Generación del vector de tiempo 't' e 'y'
t=t0:h:tN;
y=zeros(2,N+1);
y(:,1)=y0;
c(1)=0;
A=[-5 0;5 -1]; %Matriz de coeficientes del sistema
%Bucle trapecio
for i=1:N
    Z=eye(2)-(h/2)*A;
    y(:,i+1)=inv(Z)*(y(:,i)+(h/2)*A*y(:,i));
    c(i+1)=1-y(1,i)-y(2,i);
end
plot(t,y(1,:))
hold on
plot(t,y(2,:),'r')
plot(t,c,'g')
legend('[A]','[B]','[C]','location','best')

Expuestas nuestras gráficas observamos que el proceso comienza con A desintegrándose, creando B y también apareciendo poco a poco C. En t=1 vemos que A se ha reducido a cero y que B continúa desintegrándose hacia el compuesto C. Estudiando las dos gráficas apreciamos una pequeña variación de resultados, por ejemplo, el máximo del compuesto B es mayor con el método de Euler.

5 Cambio de las k1 y k2

5.1 Euler

MAS COSAS

derecha

%Datos iniciales
t0=0;
tN=10;
h=0.1;
N=round((tN-t0)/h);
y0=[1;0];
%Generación del vector de tiempo 't' e 'y'
t=t0:h:tN;
y=zeros(2,N+1);
y(:,1)=y0;
c(1)=0;
A=[-1 0;1 -5]; %Matriz de coeficientes del sistema
%Bucle Euler
for i=1:N
    y(:,i+1)=y(:,i)+h*(A*y(:,i));
    c(i+1)=1-y(1,i)-y(2,i);
end
plot(t,y(1,:))
hold on
plot(t,y(2,:),'r')
plot(t,c,'g')
legend('[A]','[B]','[C]','location','best')

5.2 Trapecio

MAS COSAS

derecha

%Datos iniciales
t0=0;
tN=10;
h=0.1;
N=round((tN-t0)/h);
y0=[1;0];
%Generación del vector de tiempo 't' e 'y'
t=t0:h:tN;
y=zeros(2,N+1);
y(:,1)=y0;
c(1)=0;
A=[-1 0;1 -5]; %Matriz de coeficientes del sistema
%Bucle trapecio
for i=1:N
    Z=eye(2)-(h/2)*A;
    y(:,i+1)=inv(Z)*(y(:,i)+(h/2)*A*y(:,i));
    c(i+1)=1-y(1,i)-y(2,i);
end
plot(t,y(1,:))
hold on
plot(t,y(2,:),'r')
plot(t,c,'g')
legend('[A]','[B]','[C]','location','best')