Desintegración Radiactiva.Grupo 5

Trabajo realizado por estudiantes
Título	Desintegración Radiactiva. Grupo 5.
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Álvaro Ramón López, Diego García Vaquero, Noemí Palomino Bustos, Mercedes Ruiz Barrajón, Teresa Quintana Romero, Araceli Martín Candilejo.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Apartado 1

clear all

%% TRABAJO 3 APARTADO 1
%       M'(t)=-1,24*10^(-4)*M(t)
%       M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
  M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
  t(i+1)=t(i)+h;
  i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');

2 Resolución por el método del Trapecio

clear all
%Trapecio
 
%y´=f(t,y);
%y(t0)=y0;
 
%y(n+1)=yn+(h/2)*(f(tn,yn)+f(t(n+1),y(n+1)))
%Hay que despejar manualmente para tener SOLO y(n+1) a un lado de la ecuación.
 
% Si representamos la concentración del elemento radiactivo M(t) como y(t):
%M'(t)=-kM(t)---> y'(t)=-ky(t)
 
%En este caso, cuando se despeja manualmente:
%y(n+1)=y(n)+(h/2)(-ky(n)-ky(n+1))
%(1+kh/2)y(n+1)=y(n)-khy(n)/2
%y(n+1)=(y(n)-khy(n)/2)/(1+kh/2)
 
%SOLUCIÓN:
 
t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;
 
%Constante de desintegración:
k=1.24*10^(-4);
 
%Bucle:
i=1;
while M(i)>(0.08*M0)
 M(i+1)=(M(i)-(k*h*M(i))/2)/(1+(k*h)/2); %Trapecio
 t(i+1)=t(i)+h;
 t(i+1)=t(i)+h;
 i=i+1;
end
 
disp('Tiempo para que quede un 8% de la cantidad inicial: ')
disp(t(end))
 
hold on
plot(t,M);
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Trapecio','Location','best');

Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 .

Resultado del valor del tiempo para el cual el contenido de C14 es el 8% del contenido inicial, según el método del Trapecio con un paso h=0.1 .

3 Resolución por el método de Runge-Kutta

clear all

clear all
 
t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;
 
%Se inicializa "M" con M0 para irlo rellenando posteriormente con las soluciones obtenidas
%para cada instante mediante el bucle while. Simultáneamente se va ampliando el vector de tiempos
%hasta que se cumple la condición deseada.

 
%RUNGE-KUTTA:
%K1=f(  tn  ,  yn  )
%K2=f(  tn+h/2,  yn+K1*h/2 );
%K3=f(  tn+h/2,  yn+K2*h/2 );
%K4=f(  tn+h,  yn+K3*h );
%y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
 
%En nuestro caso: y'=f(t,y)=-ky

k=1.24*10^(-4);
 
i=1;
while M(i)>(0.5*M0)
    %K1=f(  tn  ,  yn  )
    K1=-k*M(i);
    %Definicion de variable K2
    %K2=f(  tn+h/2,  yn+K1*h/2 );
    t2=t(i)+(h/2);
    M2=M(i)+(h/2)*K1;
    K2=-k*M2;
    %Definicion de variable K3;
     %K3=f(  tn+h/2,  yn+K2*h/2 );
    t3=t2;
    M3=M(i)+(h/2)*K2;
    K3=-k*M3;
    %Definicion de variable K4;
     %K4=f(  tn+h,  yn+K3*h );
    t4=t(i)+h;
    M4=M(i)+h*K3;
    K4=-k*M4;
    %Funcion de RungeKutta;
    %y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
    M(i+1)=M(i)+(h/6)*(K1+2*K2+2*K3+K4);
    t(i+1)=t(i)+h;
    i=i+1;
end
 
disp('Tiempo medio: ')
disp(t(end))
 
 
hold on
plot(t,y)
legend('Runge Kutta','Location','best')
hold off

Evolución del contenido de C14 (hasta que éste es el 50% del contenido inicial) según el método de Runge Kutta con un paso h=0.1 .

La vida media que da el programa es de 5592,51 años.

4 Estudio de la evolución de concentraciones de compuestos interrelacionados

4.1 Sistemas de ecuaciones lineales: Método de Euler.

clc
clf
clear all

%Datos iniciales:
%Tiempo:
t0=0;
tN=10;
%Paso:
h=0.1;
%Número de subintervalos
N=(tN-t0)/h; 

%Definimos la variable independiente: El vector tiempo
t=t0:h:tN;
y=zeros(2,N+1);

%Valores de las concentraciones iniciales:
A0=1;
B0=0;
C0=0;

%Las soluciones se recogerán en un vector "y" que irá autoformándose.
%Inicialmente, sólo se le asignará los valores de las concentraciones iniciales desde las que parte.

y0=[A0,B0]';
y(:,1)=y0;

%Un sistema Lineal:
%y'=My+T
%En este problema en particular:
% A=(-k1    0)(A)+(0)
% B=( k1  -k2)(B)+(0)

%M=(-k1    0)
%  ( k1  -k2)

M=[-5,0;5,-1];

%Bucle:

for i=1:N
y(:,i+1)=y(:,i)+h*(M*y(:,i));%Euler
end

%Se reasigna cada parte del vector "y" que recoge las soluciones con las concentraciones que representan:
A=y(1,:);
B=y(2,:);
C=C0+A0-A-B;

%Dibujo:
hold on
plot(t,A)
plot(t,B,'r')
plot(t,C,'g')
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off