Ver el código de «Desintegración radiactiva (Grupo 13-C)»

 {{ TrabajoED | DESINTEGRACIÓN RADIACTIVA. GRUPO 13-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Beatriz Oliva Manzanero

Rubén Peláez Moreno

José María Pérez Doval

Álvaro Luis Pérez Martín

Ignacio Nieto Peña

María Pablos Romero }}

La radiactividad es un fenómeno físico por el cual los núcleos de algunos elementos químicos que son inestables son capaces de transformarse, o decaer, espontáneamente, en núcleos atómicos de otros elementos mas estables. Lo cual se produce mediante emisiones de partículas con una determinada energía cinética.

Esta propiedad, característica de los isótopos inestables, es útil para estimar la edad de variedad de muestras naturales, como rocas y materia orgánica. Esto es posible, siempre y cuando, se conozca el ritmo de desintegración de determinado isótopo, en relación a los que ya han decaído.

Un ejemplo de ello es el carbono 14, el cual, es comúnmente utilizado para, en función del grado de desintegración en que se encuentre, calcular     la edad de restos arqueológicos.

Este proceso puede simularse resolviendo la siguiente ecuación diferencial para la concentración de elemento radiactivo:
'''M'(t) = - kM(t)'''


==.-Interpretación==
Resolviendo la ecuación anterior obtenemos la siguiente expresión:

[[Archivo:Sistema1.png|250px|thumb|centro|]]

A partir esta podemos deducir que M(t) representa la cantidad de masa del elemento radiactivo en cada instante del tiempo. La cual depende de la masa inicial, M0, y de la constante k, constante de desintegración, que es característica de cada elemento y determina la velocidad de desintegración de este, como podemos observar en la ecuación diferencial.
Según avanza el tiempo la masa sigue una ley decreciente de forma exponencial con respecto al tiempo, siendo k el grado de proporcionalidad. 


==.-Problema de valor inicial(Método de Euler)==
Para calcular la edad de los restos arqueológicos comenzamos con el cálculo de un problema de valor inicial, eligiendo una condición inicial y aplicando el método de Euler para los intervalos h=0.1 y h=0.01. 
                                       
                              [[Archivo:Sistema2.png|300px|thumb|centro|    Resolución numérica del problema de valor inicial]]

{{matlab|codigo=

t0=0; y0=1; h1=0.1; h2=0.01;
y1(1)=y0; y2(1)=y0;
t1(1)=t0; t2(1)=t0;
i=1; n=1; k=1.24e-04;
cond=y0*0.08;
while y1(i)>cond
    y1(i+1)=y1(i)+h1*(-k*y1(i));
    t1(i+1)=t1(i)+h1;
    i=i+1;
end
while y2(n)>cond
    t2(n+1)=t2(n)+h2;
    y2(n+1)=y2(n)+h2*(-k*y2(n));
    n=n+1;
end
disp('Años de antigüedad de los huesos:')
disp(t2(end))
disp(i)
hold on
  plot(t1,y1)
  plot(t2,y2,'--g')
  legend('euler h=0.1','euler h=0.01','location','best')
hold off 
}}
                       
                        [[Archivo:Edad.png|300px|centro|thumb|]]
                        [[Archivo:Grafica32.png|600px|centro|thumb|]]


Observando la gráfica anterior, vemos que tanto para un paso de 0.1 como para uno de 0.01, el método de Euler mantiene un error similar a lo largo de todo el intervalo, por lo tanto podemos afirmar que es estable.
En el problema planteado la masa está expresada como porcentaje de la masa inicial, por lo que esto, sumado a que el error se mantiene estable, hará que el tiempo no dependa de la masa inicial.

==.-Problema de valor inicial (Método del Trapecio)==
Resolución del problema anterior por el Método del Trapecio para h=0.1 de forma que se detiene al alcanzar el 8% de la cantidad inicial.
{{matlab|codigo=


t0=0;
h1=0.1;
y0=1; 
z(1)=y0;       
t1(1)=t0;
i=1; n=1; k=1.24*(10^(-4));
cond=y0*0.08
while z(i)>cond
    z(i+1)=((1-h1*k/2)/(1+h1*k/2))*z(i);
    t1(i+1)=t1(i)+h1;
    i=i+1;
end
disp('tiempo final:')
disp(t1(end))
disp(i)
hold on
  plot(t1,z,'m')
  legend('trapecio','location','best')
hold off
}}

[[Archivo:Grafica33.png|600px|thumb|centro|]]

==.-Vida media==
Para cada sustancia radiactiva hay un intervalo t fijo, denominado vida media, durante el cual el número de núcleos que había al comienzo se reduce a la mitad. Poniendo en la ecuación M=M0/2 se obtiene,

[[Archivo:Vida_media.png|250px|thumb|centro|]]


que relaciona la vida media y la constante de desintegración.
En este ejercicio además de utilizar el método de Runge-Kutta, hemos realizado una comparativa con los métodos de los apartados anteriores además de la solución exacta.
{{matlab|codigo=
t0=0;
h1=0.1;
h2=0.01;
y0=1; 
y1(1)=y0;
y2(1)=y0;
trap(1)=y0; 
RK(1)=y0;
ye(1)=y0;
t1(1)=t0;
t2(1)=t0;
t3(1)=t0;
i=1;
n=1;
m=1;
k=1.24e-04;
cond1=y0*0.08;
cond2=y0/2;

while trap(i)>cond1
    
    y1(i+1)=y1(i)+h1*(-k*y1(i));
    trap(i+1)=((1-h1*k/2)/(1+h1*k/2))*trap(i);
    t1(i+1)=t1(i)+h1;
    ye(i+1)=y0*exp(-k*t1(i+1));
    i=i+1;
    
end
while y2(n)>cond1
    
t2(n+1)=t2(n)+h2;
y2(n+1)=y2(n)+h2*(-k*y2(n));
n=n+1;

end
while RK(m)>cond2
    
    K1=-k*RK(m);
    K2=-k*RK(m)-k/2*K1*h1;
    K3=-k*RK(m)-k/2*K2*h1;
    K4=-k*RK(m)-k*K3*h1;
    
    RK(m+1)=RK(m)+h1/6*(K1+2*K2+2*K3+K4);
    t3(m+1)=t3(m)+h1;
    m=m+1;
end

disp('años de antigüedad de los huesos:')
disp(t2(end))
disp('vida media del carbono 14:')
disp(t3(end))
hold on
plot(t1,y1,'r')
plot(t2,y2,'g')
plot(t1,trap)
plot(t3,RK,'*k')
plot(t1,ye,'m')
legend('Euler h=0.1','Euler h=0.01','trapecio','Runge-Kutta','exacta','location','best')
hold off
}}

La solución de este problema sería:
[[Archivo:Vmedia.png|350px|thumb|centro|]]
En la siguiente gráfica hacemos una comparativa de los ejercicios anteriores y la solución exacta. En esta, podemos comprobar que el método de Runge-Kutta es el que mejor se aproxima a la solución exacta, siendo los de Euler los que menos (a mayor intervalo, en este caso 0.1,peor aproximación)
[[Archivo:Ejerfinal4.png|600px|thumb|centro|Hemos hecho zoom en las dos gráficas inferiores para poder verlo con mayor claridad. La de la derecha tiene mayor zoom que la de la izquierda.]]

==.-Serie de desintegración radiactiva==
Cuando una sustancia se desintegra radiactivamente, por lo general no solo se transmuta en en una sustancia estable (con lo que se detiene el proceso), sino que la primera sustancia se desintegra, forma otra sustancia radiactiva y ésta, a su vez, decaer y forma una tercera sustancia, etc. 

[[Archivo:Sistema22.png|300px|thumb|centro|Sistema de ecuaciones diferenciales que permite conocer la cantidad de masa de cada elemento en cada instante t.]]

Siendo K1 Y K2 las constantes de desintegración de MA y MB y MC un elemento estable. Supongamos también que MA(t),MB(t)y MC(t) representan las cantidades de los elementos MA, MB y MC que quedan en cualquier instante.
[[Archivo:Gráfica_desintegración_radiactiva.png|300px|thumb|centro|Serie de desintegración radiactiva]]

==.-Serie de desintegración radiactiva del PVI.(Método de Euler y Trapecio)==

Resolvemos tomando el apartado anterior, el PVI  conocidas las constantes de desintegración y las concentraciones iniciales en un intervalo de h=0.1 para Euler y otro de t∈[0,1] para el método del Trapecio.
{{matlab|codigo=
t0=0;
tN=10;
h=0.1;
k1=5;
k2=1;
N=(tN-t0)/h;
t=t0:h:tN;
y=zeros(3,N+1);
y0=[1,0,0]';
z=zeros(3,N+1);
z0=[1,0,0]';
I=eye(3);
y(:,1)=y0;
z(:,1)=z0;
A=[-k1,0,0;k1,-k2,0;0,k2,0];
 for i=1:N
     y(:,i+1)=y(:,i)+h*(A*y(:,i));%Método de Euler
     z(:,i+1)=(I-(h/2)*A)\(z(:,i)+(h/2)*(A*z(:,i)));%Método del trapecio
 end
 hold on
 plot(t,y(1,:),'--')
 plot(t,y(2,:),'--g')
 plot(t,y(3,:),'--r')
 plot(t,z(1,:))
 plot(t,z(2,:),'g')
 plot(t,z(3,:),'r')
 xlabel('tiempo')
 ylabel('masa')
 hold off
 legend('desMA-Euler','desMB-Euler','desMC-Euler','desMA-Trapecio','desMB-Trapecio','desMC-Trapecio','Location','Best');
}}
[[Archivo:Ejer6.png|500px|thumb|left|]] 
[[Archivo:Ejer64.png|500px|thumb|centro|Cada gráfica por separado para la comparación de los métodos Euler y Trapecio en la desintegración de cada masa]]

Al ser cinco veces mayor la velocidad de desintegración de la masa A que la de B, podemos observar que ésta se desintegra prácticamente entera en un breve periodo de tiempo, comparando con un mayor periodo de desaparición de la masa B. La masa B comienza aumentando y continúa con el paso del tiempo, hasta el punto en el que su relación con la masa A sea igual a la relación entre sus constantes de desintegración, en nuestro caso, cuando la masa B sea 5 veces mayor que la masa A. A partir de este instante la masa B empieza a decrecer.

==.-Serie de desintegración radiactiva del PVI.(Método de Euler y Trapecio).Intercambiadas las constantes== 

{{matlab|codigo=
t0=0;
tN=10;
h=0.1;
k2=5;
k1=1;
N=(tN-t0)/h;
t=t0:h:tN;
y=zeros(3,N+1);
y0=[1,0,0]';
z=zeros(3,N+1);
z0=[1,0,0]';
I=eye(3);
y(:,1)=y0;
z(:,1)=z0;
A=[-k1,0,0;k1,-k2,0;0,k2,0];
 for i=1:N
     y(:,i+1)=y(:,i)+h*(A*y(:,i));%Método de Euler
     z(:,i+1)=(I-(h/2)*A)\(z(:,i)+(h/2)*(A*z(:,i)));%Método del trapecio
 end
 hold on
 plot(t,y(1,:),'--')
 plot(t,y(2,:),'--g')
 plot(t,y(3,:),'--r')
 plot(t,z(1,:))
 plot(t,z(2,:),'g')
 plot(t,z(3,:),'r')
 xlabel('tiempo')
 ylabel('masa')
 hold off
 legend('desMA-Euler','desMB-Euler','desMC-Euler','desMA-Trapecio','desMB-Trapecio','desMC-Trapecio','Location','Best');
}}
[[Archivo:Ejer7.png|500px|thumb|centro|]] 
[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

En este caso, al haber intercambiado las constantes, el máximo de la masa B tendrá lugar cuando sea cinco veces menor que la masa A, al contrario del apartado anterior que se daba cuando era cinco veces mayor, pero se dará en el mismo instante que antes. Por otro lado, vemos que la masa C aumenta del mismo modo en los dos casos, ya que ésta depende tanto de la masa A y la masa B, pero el intercambio entre las constantes de desintegración de A y B no varía la cantidad de C en cada instante.