Desintegracion Radioactiva (Grupo 19C)

Trabajo realizado por estudiantes
Título	Desintegración Radioactiva Grupo 19C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Ciro Rodriguez Matamoros, David Fariña Berjon, Agustín Laja Santiago, Gonzalo Pizarro Cuervo-Arango, Jesus Perez Fernandez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

La desintegración radiactiva es un fenómeno que sse produce de manera natural mediante el cuál un elemento (como el plutonio o el radio) se descompone formando un isótopo u otro elemento. Dicho proceso se produce con una rapidez proporcional a la cantidad de material radioactivo. Es especialmente llamativo el caso del isótopo C14 que es usado para datar objetos de alto valor histórico.

En este trabajo analizaremos la desintegracion radiactiva desde el punto de vista del cálculo de ecuaciones diferenciales por medio de métodos numéricos de aproximación de la solución.

2 Metodos numéricos de resolución

EULER

Evolucion de la desintegracion del C14 hata alcanzar el 8% de la cantidad inicial 'M0'

clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
z(1)=M0;
i=1;
j=1;
k=1.24*10^-4;
while M(i)>0.08*M0;
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);
    z(i+1)=z(i)+h*((-1.24*10^(-4))*z(i));
    t(i+1)=t(i)+h;
    i=i+1;
end
[t',M',z']
hold on
plot(t,M)
plot(t,z,'--r')
hold off
disp(t(end))%Muestra en pantalla el valor del vector "t" de tiempo

Tiempo en el que se desintegra el 92% del C14

• Este Programa data la antiguedad de los huesos en 20369 años.Posteriormente ejecutamos el programa con valores de paso h=0.1 y h=0.01 dando una gráfica similar siendo necesario ampliarla para poder ver la comparativa. El resultado final no experimenta cambios al cambiar los valores de paso.

clear all, clf
%M'(t)=-1,24*10^-4*M(t)
% A continuación exponemos una comparación de las dos gráficas que resultan
% de variar el paso de 0.1 a 0.01
h1=0.1;
h2=0.01;
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t1(1)=t0;
t2(1)=t0
M1(1)=M0;
M2(1)=M0;
i=1;
n=1;
%El primer Bucle utiliza un paso de 0.1
while M1(i)>(0.08*M0);
    M1(i+1)=M1(i)+h1*((-1.24*10^(-4))*M1(i));
    t1(i+1)=t1(i)+h1;
    i=i+1;
end
%El segundo bucle utiliza un paso de 0.01
while M2(n)>(0.08*M0);
    M2(n+1)=M2(n)+h2*((-1.24*10^(-4))*M2(n));
    t2(n+1)=t2(n)+h2;
    n=n+1;
end
%expone en pantalla la antigüedad de los restos arqueológicos
hold on
plot(t1,M1,'r')
plot(t2,M2,'--b')
hold off

Gráfica comparativa para los valores de paso h=0.1 y h=0.01.

• En el siguiente programa observaremos como la cantidad inicial de materia no afecta a la hora de calcular el tiempo de desintegración.

clear all, clf
%M'(t)=-1,24*10^-4*M(t)
%Método de Euler
h=0.01;
M0=1;
N0=100;
Z0=50;
t0=0;
t1(1)=t0;%Utilizamos el vector t1 para calcular la edad de un resto arqueológico con una supuesta masa de 1
t2(1)=t0;%Utilizamos el vector t1 para calcular la edad de un resto arqueológico con una supuesta masa de 100
t3(1)=t0;
M(1)=M0;
N(1)=N0;
Z(1)=Z0;
i=1;
%utilizamos un bucle "while" para evitar que la función siga adoptando
%valores una vez superado llegue a ese 8% que se encuentra en el hueso
while M(i)>(0.08*M0);
    M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
    t1(i+1)=t1(i)+h;
    i=i+1;
end
j=1;
while N(j)>(0.08*N0);
    N(j+1)=N(j)+h*((-1.24*10^(-4))*N(j));
    t2(j+1)=t2(j)+h;
    j=j+1;
end
n=1
while Z(n)>(0.08*Z0);
    Z(n+1)=Z(n)+h*((-1.24*10^(-4))*Z(n));
    t3(n+1)=t3(n)+h;
    n=n+1;
end
hold on
plot(t2,N,'g')
plot(t1,M)
plot(t3,Z,'r')
hold off
disp(t1(end))
disp(t2(end))
disp(t3(end))
%Sacamos en pantalla el valor de t1 y t2 para comparar la edad de ambos.
Error=t1-t2;
disp(Error(end))%Expone en pantalla la diferencia entre t1 y t2 y queda 
%demostrado que no depende de la cantidad inicial de M (materia).

Evolucion de la desintegración para distinto valores iniciales de materia. En la imagen de queda demostrado que para valores iniciales distintos como lo son 1, 50 y 100, el tiempo del cual datarían los restos arqueológicos no cambiaría. En otras palabras, la cantidad inicial de materia no afecta a la hora de calcular el tiempo de desintegración.

Trapecio

• A continuación procederemos a la resolucion del mismo problema por el metodo del trapecio como muestra el siguiente codigo:

%El método del trapecio es 'indirecto' lo que quiere decir que previamente
%debemos despejar analíticamente el valor de f(t(n+1),y(n+1)
clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
i=1;
k=1.24*10^-4;
%Igual que en el caso anterior, Utilizaremos un bucle "while" para que el
%programa se detenga automáticamente una vez reducida la cantidad de C14
while M(i)>0.08*M0;
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);
    t(i+1)=t(i)+h;
    i=i+1;
end
disp(t(end))
plot(t,M)

Gráfica que muestra la evolución de la desintegracion resuelta por el metodo del trapecio.

• En el siguiente programa se comparan los dos métodos, Euler y Trapecio. Ambos métodos numericos son muy exactos y su error es pequeño, de esta manera y debido al orden de resolución de ambos sistemas, debemos hacer un zoom profundo en la gráfica para observar la diferencia entre ambas soluciones.

Gráfica comparativa ampliada de los metodos de Euler y Trapecio

clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
z(1)=M0;
i=1;
j=1;
k=1.24*10^-4;
%La "M" representa los valores obtenidos por el método de Euler, y la "z" el
%vector de volores obtenidos mediante Trapecio
while M(i)>0.08*M0;
    %Euler
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);
    %Trapecio
    z(i+1)=z(i)+h*((-1.24*10^(-4))*z(i));
    t(i+1)=t(i)+h;
    i=i+1;
end
%A continuación se muestra una tabla comparativa para los distintos valores
%de "t" comparando los obtenidos mediante "Euler" y mediante "Trapecio"
[t',M',z'];
hold on
plot(t,M,'linewidth',1.5)
plot(t,z,'--r','linewidth',1.5)
hold off

3 Rounge Kutta

%Rounge Kutta
clear all, clf
%M'(t)=-1,24*10^-4*M(t)
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
i=1;
k=1.24*10^-4;
%Los métodos de Runge-Kutta (RK) son un conjuntos de métodos iterativos
%para la aproximación de ecuaciones diferenciales ordinarias.
while M(i)>0.5*M0;
    %K1=f(t(n),M(n))
    K1=-k*M(i);
    %K2=f(t(n)+1/2*h, M(n)+1/2K1*h)
    K2=-k*M(i)-k/2*K1*h;
    %K3=f(t(n)+1/2*h,M(n)+1/2*K2*h)
    K3=-k*M(i)-k/2*K2*h;
    %K4=f(t(n)+h,M(n)+K3*h)
    K4=-k*(M(i)+K3*h);
    %M(n+1)=M(n)+h/6*(K1+2K2+2K3+K4)
    M(i+1)=M(i)+h/6*(K1+2*K2+2*K3+K4);
    t(i+1)=t(i)+h;
    i=i+1;
end
plot(t,M,'g')
disp(t(end))

La vida media es de 5589,9 años.

4 Desintegración Isotópica

Consideramos el estudio de una descomposión de un elemento A hasta uno C , con una fase intermedia de descomposición B, teniendo dos parámetros de descomposición denominados k1 y k2 . De esta manera podemos deducir un sistema de ecuaciones que equipara la descomposición de estos elementos. \begin{cases}MA'=-k_{1}MA\\MB'=k_{1}MA-k_{2}MB\\MC'=-k_{2}MB\\\end{cases} Analizamos el estudios de este sistema por métodos numéricos, en este caso el método de Euler y del Trapecio.

%Condiciones iniciales.
t0=0;
tf=10;
k1=5;
k2=1;
%Tamaño de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%Matriz para la resolución
y=zeros(3,N+1);
%Concentraciones iniciales
A=1;B=0;C=0;
y0=[A B C]';
y(:,1)=y0;
%Resolvemos por euler partiendo de la matriz de coeficientes sabiendo que
%y'=M*y+t;
M=[-k1 0 0;k1 -k2 0;0 k2 0];
for i=1:N
    y(:,i+1)=y(:,i)+h*(M*y(:,i));
end
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
plot(t,y(3,:))
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off

%Condiciones iniciales.
t0=0;
tf=10;
k1=5;
k2=1;
%Tamaño de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%Matriz para la resolución
y=zeros(3,N+1);
I=eye(3);
%Concentraciones iniciales
A=1;B=0;C=0;
y0=[A B C]';
y(:,1)=y0;
%Resolvemos por euler partiendo de la matriz de coeficientes sabiendo que
%y'=M*y+t;
M=[-k1 0 0;k1 -k2 0;0 k2 0];
for i=1:N
    y(:,i+1)=(I-(h/2)*M)\(y(:,i)+(h/2)*(M*y(:,i)));%Método del trapecio
end
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
plot(t,y(3,:))
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off

Soluciones de ambos métodos

En esta gráfica podemos comparar ambas soluciones numéricas, de manera que ambas son correctas pero con un orden distinto, demanera que Trapecio tiene un error menor que Euler, por lo que tomaremos como correcta la solución de Trapecio. Procedemos a interpretar un cambio de las valores k1 y k2, sustituyendo las nuevas condiciones en el sistema. De manera que optenemos:

Nuevos graficos variando k1 y k2

De esta manera variando los valores de conversion de un isótopo a otro, podemos deducir que la fase B de descomposición del sistema es una fase totalmente efímera, es un estado mas transicional ya que toda la materia B se transforma de manera mucho mas veloz a materia C.