Diferencia entre revisiones de «Desintegracion Radioactiva (Grupo 19C)»
(→Metodos numéricos de resolución) |
(→Metodos numéricos de resolución) |
||
| Línea 143: | Línea 143: | ||
[[Image:Compmat.jpg|800px|thumb|center|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. ]] | [[Image:Compmat.jpg|800px|thumb|center|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 | + | ==Trapecio== |
{{matlab|codigo= | {{matlab|codigo= | ||
Revisión del 13:53 5 mar 2015
| 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 | |
Contenido
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
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
- 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
- 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).
3 Trapecio
%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)
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
4 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.
5 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 offEn 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:
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.

