Diferencia entre revisiones de «Desintegración Radiactiva.Grupo 5»

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');

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');

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