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

clear all

%       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

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)
k1=5;
k2=1;

M=[-k1,0;k1,-k2];

%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

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)
k1=5;
k2=1;

M=[-k1,0;k1,-k2];

%Bucle:
%Para el Trapecio:
%y(n+1)=y(n)+h/2*(f(tn,yn)+f(t(n+1),y(n+1))
%En este caso:
%y(n+1)=y(n)+h/2*(M*yn+M*y(n+1))
%Despejando manualmente:
%[1-(h/2)*M]*y(n+1)=y(n)+h/2*(M*yn)
%Llamando Z a:
%Z=[1-(h/2)*M]
%Z*(INV(z))*y(n+1)=(INV(z))*[y(n)+h/2*(M*yn)]
%Y queda finalmente:

%y(n+1)=(INV(z))*[y(n)+h/2*(M*yn)]

for i=1:N
Z=eye(2)-(h/2)*M;
y(:,i+1)=inv(Z)*(y(:,i)+(h/2)*M*y(:,i));%Trapecio
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