Diferencia entre revisiones de «Desintegracion Radioactiva (Grupo 19C)»
De MateWiki
(→Euler) |
(→Euler) |
||
| Línea 92: | Línea 92: | ||
}} | }} | ||
| − | [[Image:Compmat.jpg|800px|thumb|center|Evolucion de la desintegración para distinto valores iniciales de materia ]] | + | [[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= | ||
Revisión del 13:08 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
3 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.
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.
4 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)
5 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.
6 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 k_{1} y k_{2} . 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 offA continuación proponemos una visión comparativa de ambas soluciones:
