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
clear all
clf

%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:

%PASO 1
%DATOS ENUNCIADO:
t0=0;
y0=1;
tN=10000;

%NUMERO DE SUBINTERVALOS:N
h=0.1;
%h=(tN-t0)/N
N=(tN-t0)/h;

%VECTOR DE TIEMPOS:
t=t0:h:tN;
%t=linspace(t0,tN,N+1);

%VECTOR DE SOLUCIONES:
y=zeros(1,N+1);
%relleno el primer valor de "y" con y0.
y(1)=y0;

%Constante de desintegración:
k=1.24*10^(-4);

%Bucle:
i=1;
while y(i)>(0.08*y0)
 y(i+1)=(y(i)-(k*h*y(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 y=0.08:')
disp(t(end))

hold on
plot(t,y);
legend('Trapecio','Location','best');

clear all

clear all, clf
%DATOS DEL PROBLEMA
t0=0;
tN=10000;
y0=1;

%h ES EL ESPACIO QUE OCUPA CADA SUBINTERVALOS (El paso).
%N ES EL NUMERO DE SUBINTERVALOS
%Tendremos N+1 NODOS
h=input('Introduce tamaño del paso h: ')
%h=0.1;
N=(tN-t0)/h;

%La variable independiente t:
t=t0:h:tN;

%Se crea un vector "y" para irlo rellenando posteriormente con las soluciones.
%La longitud de y tiene que ser de la misma longitud que la del vector tiempo.
y=zeros(1,N+1); %Será la solución del Runge Kutta;

%relleno el primer valor de "y" con y0.
y(1)=y0;

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

n=1;
while y(n)>(0.5*y0)
    %K1=f(  tn  ,  yn  )
    K1=-k*y(n);
    %Definicion de variable K2
    %K2=f(  tn+h/2,  yn+K1*h/2 );
    t2=t(n)+(h/2);
    y2=y(n)+(h/2)*K1;
    K2=-k*y2;
    %Definicion de variable K3;
     %K3=f(  tn+h/2,  yn+K2*h/2 );
    t3=t2;
    y3=y(n)+(h/2)*K2;
    K3=-k*y3;
    %Definicion de variable K4;
     %K4=f(  tn+h,  yn+K3*h );
    t4=t(n)+h;
    y4=y(n)+h*K3;
    K4=-k*y4;
    %Funcion de RungeKutta;
    %y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
    y(n+1)=y(n)+(h/6)*(K1+2*K2+2*K3+K4);
    t(n+1)=t(n)+h;
    n=n+1;
end

disp('Tiempo medio: ')
disp(t(end))


hold on
plot(t,y)
legend('Runge Kutta','Location','best')
hold off