Desintegración Radiactiva 11-A

Trabajo realizado por estudiantes
Título	DESINTEGRACIÓN RADIACTIVA. GRUPO 10-A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Alejandro Carrillo del Aguila (1400) Antonio Carrillo del Aguila (80) Humberto del Castillo Montes de Oca (1281) Laura de la Morena Mendez (1421) Miguel Coello Guijarro (1408)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

El trabajo propuesto nos plantea el cálculo de la desintegración de un material radiactivo a lo largo del tiempo, sabiendo que estos materiales se van desintegrando proporcionalmente a la cantidad restante. Por lo que analíticamente, un material radiactivo se desintegra en función de la siguiente ecuación diferencial:

[math]M'(t) = −k·M(t)[/math]

2 Interpretación de M(t) y K

M(t) es la cantidad de material radiactivo restante respecto del tiempo y K es una constante de desintegración que variará dependiendo del material. Por lo que cuanto mas alta sea la constante de desintegración, mas alto será el valor absoluto de la velocidad, es decir, el material se desintegrará en menos tiempo.

3 Planteamiento del Problema de Valor Inicial (P.V.I)

Tomando M₀ como la cantidad inicial y se plantea el siguiente problema de valor inicial:

[math] P.V.I. \left\{\begin{matrix}\ M(t_{0}) = −kM(t)\\ M(t_{0})=M_{0}\end{matrix}\right. \Longrightarrow M'(t)=\frac{\operatorname dM(t)}{\operatorname dt}= -K·M(t) \Longrightarrow \frac{1}{M(t)}·\operatorname dM(t)= -k·\operatorname dt [/math]

Integrando e imponiendo la condición inicial se obtiene la siguiente solución:

[math]M(t)=M_{0}·e^{-kt}[/math]

Como se puede observar, cuando el tiempo es cero aun quedará todo el material inicial y según transcurra el tiempo la cantidad de materia ira decreciendo exponencialmente.

4 Interpretación y Resolución Numérica

El trabajo propuesto nos plantea concretamente el caso del [math]C^{14}[/math], informándonos que este tiene una constante de desistegración [math]K=1,24·10^{-4}[/math] y unos los huesos encontrados por un arqueólogo contienen un 8% de [math]C^{14}[/math].

4.1 Resolución por el método de Euler

La cantidad inicial del material es indiferente para cualquier tipo de cálculos porcentuales o temporales como se demuestra a continuación:

[math] \left\{\begin{matrix}\ M(t_{0}) = −kM(t)\\ M(t^{0.08})=0.08·M_{0}\end{matrix}\right. \Longrightarrow 0.08M_{0}=M_{0}·e^{-kt^{0.08}} \Longrightarrow 0.08=e^{-kt^{0.08}} \Longrightarrow t^{0.08}=\frac{-Ln(0.08)}{K} [/math]

Por lo que podemos afirmar que el problema es muy estable, ya que al cambiar la condición inicial no varia nada.

Pero para su resolución numérica con Octave (o Matlab), debemos elegir una cantidad inicial. Nosotros lo hemos planteado con un valor de 100 para que represente el porcentaje inicial de material en vez de la cantidad material, pero repito, seria indiferente este valor para los cálculos porcentuales y temporales.

%Variables del Problema
h=0.1; % o h=0.01 dependiendo de la gráfica
t0=0;
k=1.24*10^(-4);
m0=100;
tn=40000;
t=0:h:tn;
m=t;
solucion=1;
m(1)=m0;
N=(tn-t0)/h; %si quisiera el numero de subintervalos
  
%con linespace(t0,tn,n+1)
for i=1:N;
  m(i+1)=m(i)+h*(-k*m(i));
  
  if m(1+i)<=0.08*m(1) & length(solucion)==1
    solucion(2)=t(i+1)/h;
  end
  
end

hold on
  plot(t,m)
  anyosDeDesintegracion=solucion(2)
  plot(solucion(2)*h,m(solucion(2)),'+')
hold off

4.1.1 Gráfica Euler h=01

En la gráfica se observa lo dicho anteriormente. La cantidad de materia ira disminuyendo exponencialmente. Debido a que inicialmente hay la máxima cantidad de material radiactivo, disminuirá mas rápidamente. A medida que pasa el tiempo, el material disminuye, y por ende la velocidad de desintegración también. Por lo que como se observa en la gráfica, cada vez la pendiente se va haciendo mas suave.

4.1.2 Gráfica Euler h=001

4.2 Trapecio

Para la resolución del problema mediante el método del trapecio ha sido necesario despejar la cantidad de materia [math]M_{(i+1)}[/math] de la ecuación implícita, para poder meterla posteriormente en un bucle en el cual impondremos que una vez alcanzada cierta cantidad [math]M=0.8·M_{0}[/math] pare, tal y como nos plantea el problema.

%Variables del problema
t(1)=0;
h=0.1;
tn=40000;
k=1.24*(10^(-4));
m(1)=100;
t=t(1):h:tn;
n=(tn-t(1))/h;

for i=1:n;
  m(i+1)=(m(i)*(1-k*h/2))/(1+h*k/2);
  if m(i+1)<8;
    tmaximo=(1+i)
    break
  end
end

t=0:h:(length(m)-1)*h;
plot(t,m,'r')

Como se observa en la gráfica, la cantidad de material radiactivo llegará al 8% cuando casi hayan pasado 20369 años

4.3 Runge-Kutta

%DATOS DEL PROBLEMA 
t0=0; 
tN=40000; 
m0=100;
k=1.24*(10^(-4));
h=0.1;
N=(tN-t0)/h; 
% Definimos la variable independiente 
t=t0:h:tN; 
m(1)=m0;
solucion=1;
% Solución exacta 
  for i=1:N ;
  k1=-k*m(i);
  k2=-k*(m(i)+1/2*k1*h);
  k3=-k*(m(i)+1/2*k2*h);
  k4=-k*(m(i)+k3*h);
  m(i+1)=m(i)+h/6*(k1+2*k2+2*k3+k4);
  if m(1+i)<=50 & length(solucion)==1
  solucion(2)=t(i+1)/h;
  VidaMedia=solucion(2)
  end
end
plot(t,m)
hold on
plot(VidaMedia*h,m(VidaMedia),'+')

Como se observa en la gráfica, el elemento radiactivo llegará a su vida media (tiempo que tarda en reducirse a la mitad) casi llegados los 5590 años.

5 Desintegración de un compuesto A en C a través de otro B

5.1 Planteamiento

Se considera la descomposición de un elemento radiactivo A para dar C, a través de un isótopo intermedio B. Donde [math]k_1[/math] es la constante de desintegración de A y [math]k_2[/math] es la constante de desintegración de B.

A' será la velocidad de desintegración de dicho compuesto. Como hemos visto en apartados anteriores, un material radiactivo se descompone proporcionalmente a la cantidad de materia restante, por lo que [math]M_A(t)'=-K_1·M_A(t)[/math] (con [math]k_1\gt0[/math]). Debido a que A decrece, A' es negativa. Por lo que C aumentará hasta alcanzar la cantidad inicial de A mas la cantidad inicial de C y B. Como A se transforma en C a través de B, la velocidad de C será [math]M_C(t)'= K_2·M_B(t)[/math] (con [math]k_2\gt0[/math]). [math]M_C(t)'[/math] es positiva dado que crece.
Finalmente, podremos componer la velocidad del compuesto B, que por un lado aumentará lo que decrezca A y por el otro disminuirá lo que aumente C, por lo que [math]M_B(t)'= K_1·M_A(t) - K_2·M_B(t)[/math].

De este modo, el sistema obtenido será:

[math]\\\left\{\begin{matrix} M_A(t)'=-K_1·M_A(t) \\ M_B(t)'= K_1·M_A(t) - K_2·M_B(t) \\ M_C(t)'= K_2·M_B(t) \\ \end{matrix}\right. [/math]

5.2 Resolución numérica

Suponiendo que [math]K_1=5[/math] y [math]K_2=1[/math], la velocidad de desintegración de A será cinco veces mayor que la velocidad de desintegración de B. Además suponiendo que la cantidad inicial de B y de C es nula, explica que la masa de A decrecezca rápidamente hasta anularse mientras que la de C es creciente hasta llegar a la masa inicial de A, que suponemos que es 1.

La masa de B crece hasta el instante en que lo que recibe de A es igual a lo que se transforma en C, es decir hasta el instante en que hay 5 veces mas masa de B que de A, a partir de este momento la masa decrece hasta anularse.(Es decir, [math] M_B(t)[/math] alcanzara su máximo cuando su derivada se haga cero. Eso ocurrira cuando [math]K_1·M_A=K_2·M_B)[/math]
Se observan soluciones ilógicas ya que no tiene sentido que la cantidad final del elemento C sea mayor que la de A si la cantidad inicial de C es nula. Pero es facilmente demostrable que es debido al método utilizado, pues la h elegida por el enunciado no sería la mas correcta. Cuanto mas pequeña fuera esta, mas se aproxima la gráfica a 1.

5.2.1 Euler

%Variables
h=0.1;
t0=0;
tn=10;
t=0:h:tn;
y0=[0;0];
A=[-1,0;1,0];

%si quisiera el numero de subintervalos
N=round((tn-t0)/h);
y=zeros(2,N+1);
y(:,1)=y0;

%con linespace(t0,tn,n+1)
for i=1:N;
  y(:,i+1)=y(:,i)+h*(A*y(:,i)+[5*exp(-5*t(i));0]);
end
hold on
  plot(t,y(1,:));
  plot(t,y(2,:),'r');
  z=exp((-5)*t);
  plot(t,z,'g');
hold off

Como ya se ha comentado anteriormente. El material radiactivo A disminuye exponencialmente, mientras que C aumenta al mismo tiempo exponencialmente. Por otra parte, B tendrá un periodo inicial de crecimiento, y un periodo final de decrecimiento.

5.2.2 Trapecio

%Variables
h=0.1
t0=0
tn=10
t=0:h:tn;
y0=[0;0];
A=[-1,0;1,0];
%si quisiera el numero de subintervalos
N=round((tn-t0)/h);
y=zeros(2,N+1);
y(:,1)=y0;

%con linespace(t0,tn,n+1)
for i=1:N;
  y(:,i+1)=(eye(2)-h/2*A)\(y(:,i)+h/2*(A*y(:,i)+[5*exp(-5*t(i));0]+[5*exp(-5*t(i+1));0]));
end

hold on
  plot(t,y(1,:))
  plot(t,y(2,:),'r')
  z=exp((-5)*t);
  plot(t,z,'g')
hold off

5.2.3 Cambio de constantes

Al intercambiar las constantes de desintegración, la masa de A decrece a menor velocidad. El punto de derivada nula de la función que describe la masa de B se encuentra en este caso en el instante en que tenemos 5 veces más masa de A que de B. La masa de B comienza a decrecer hasta el instante en que hay 5 veces mas masa de A que de B.

5.2.3.1 Euler

5.2.3.2 Trapecio