Diferencia entre revisiones de «Desintegracion Radioactiva (Grupo 19C)»
(→Desintegración Isotópica) |
(→Euler) |
||
| Línea 44: | Línea 44: | ||
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. | 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. | ||
| + | {{matlab|codigo= | ||
| + | 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 | ||
| + | |||
| + | }} | ||
[[Image:G19pasos.jpg|800px|thumb|center|Gráfica comparativa para los valores de paso h=0.1 y h=0.01. ]] | [[Image:G19pasos.jpg|800px|thumb|center|Gráfica comparativa para los valores de paso h=0.1 y h=0.01. ]] | ||
Revisión del 13:19 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.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
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).
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 off
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:
