Desintegración radiactiva (Grupo 13A)

Trabajo realizado por estudiantes
Título	Desintegración radiactiva. Grupo 13A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2015-16
Autores	Manuel Bécares Martín, 1077 Óscar Lázaro González, 993 Jorge Barrena Blázquez, 1818 Damiann Cruz Rodriguez, 1861
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El estudio del presente trabajo es el análisis matemático de la desintegración radiactiva de ciertos materiales y conceptos relacionados con el mismo y su aplicación a la datación de restos arqueológicos.

1 Interpretación del fenómeno físico, la función M(t) y la constante k

El fenómeno físico de la desintegración radiactiva se basa en la observación de que materiales como el plutonio, el radio o el isótopo C¹⁴ se desintegran naturalmente para formar otro elemento o isótopo del mismo elemento con una rapidez proporcional a la cantidad de material radiactivo presente.

El proceso puede modelizarse con la siguiente ecuación diferencial:

       [math]M'(t)=\frac{dM(t)}{dt}={-k\cdot M(t)}\,[/math] , con [math]k\gt0[/math]

donde [math]M'(t)[/math] representa la variación en la cantidad en masa del elemento, [math]M(t)[/math] es la cantidad del elemento en cada instante de tiempo y la constante [math]k[/math] recibe el nombre de contante de desintegración radiactiva o constante de semidesintegración y representa la probabilidad por unidad de tiempo de que los núcleos pertenecientes al material se desintegren. También puede interpretarse como la rapidez con que se produce la desintegración (véase gráfico).

Código MATLAB empleado para obtener la gráfica: clear all; clf; %Datos conocidos t0 = 0; M0 = 10; h = 0.1; k1 = 1.24e-4; k2 = 3*1.24e-4; lim = 0.08*M0; %Metodo de Euler para k1 %%Iteracion hasta obtener 0.08*M0 t1(1) = t0; M1(1) = M0; i = 1; while M1(i) > lim M1(i+1) = M1(i)+h*(-k1*M1(i)); t1(i+1) = t1(i)+h; i = i+1; end %Metodo de Euler para k2 %%Iteracion hasta obtenr 0.08*M0 t2(1) = t0; M2(1) = M0; i = 1; while M2(i) > lim M2(i+1) = M2(i)+h*(-k2*M2(i)); t2(i+1) = t2(i)+h; i = i+1; end %Representacion grafica hold on plot(t1,M1,'r') plot(t2,M2,'k') grid on title('Efecto de la constante k en el fenómeno de desintegración') xlabel('Tiempo t (años)') ylabel('Cantidad de Carbono-14 M(t) (kg)') legend('Desintegración con k=1.24e-4','Desintegración con k=3*1.24e-4') hold off

La solución de la ecuación diferencial es del tipo [math]M(t) = C \cdot e^{-kt}[/math], siendo [math]C[/math] una constante.

Si la cantidad de materia radiactiva en el instante incial ([math]t = 0[/math]) es [math]M_0[/math], se puede determinar el valor de [math]C[/math] tomando el valor de la solución en [math]t = 0[/math]:

       [math]M(0) = M_0 = C \cdot e^{-k \cdot 0}\, = C[/math]

Por tanto, la ecuación que representa el valor de la cantidad de materia para cualquier instante de tiempo es:

       [math]M(t) = M_0 \cdot e^{-kt}\,[/math]

2 Planteamiento y resolución del problema de valor inicial que modeliza el fenómeno

Habiendo determinado anteriormente la ecuación diferencial que rige el fenómeno y su solución, para obtener el problema de valor inicial sólo se necesita determinar la condición inicial del mismo, la cual será [math]M(0) = M_0[/math].

Por tanto, el problema se valor inicial será:

       [math]\begin{cases}M'(t) = {-k\cdot M(t)}\, \\ M(0)=M_0 \end{cases}[/math]

donde [math]M_0\;[/math] se determinará en función de la masa de material radiactivo presente en instante considerado como inicial.

La resolución del problema de valor inicial se aplicará al caso particular de la datación de restos arqueológicos por la cantidad de carbono-14 presente en los restos de un ser vivo. Se supone que la cantidad de carbono-14 en el momento de analizar los restos es del 8% de la cantidad que se encuentra en un ser vivo. Esto implica que [math]M(T_F) = {0.08 \cdot M_0}\,[/math], siendo [math]T_F[/math] la edad de los restos. Se supone también que la cantidad de carbono-14 no ha variado, pudiendo tomar la diferencia en el contenido de carbono-14 como la debida únicamente a su desintegración. Se considera que la constante de desintegración del carbono-14 es [math]k = 1.24\ \times\ 10^{-4}[/math] por año. El objetivo es obtener la antigüedad de los restos analizados.

A continuación se procede a su resolución por diferentes métodos:

2.1 Método de Euler

Para la resolución por el método de Euler, se emplea la siguiente fórmula, que aproxima los valores de la función [math]M(t)[/math] en los diferentes puntos de [math]t[/math]:

       [math]M_{i+1} = M_i + h \cdot (-k \cdot M_i)[/math]

Además, se han considerado los siguientes datos: masa inicial [math]M_0=10\;[/math] kg y pasos [math]h_1 = 0.1\;[/math] y [math]h_2 = 0.01\;[/math].

El método de Euler con paso [math]h_1 = 0.1\;[/math] genera como solución la siguiente gráfica:

La gráfica se ha obtenido por medio de un programa de MATLAB, al cual se le ha especificado que finalice la representación de valores cuando se alcance la cantidad de carbono-14 que existe en el momento de datación (es decir, [math]{0.08 \cdot M_0}\,[/math]). El código empleado es el siguiente:

clear all; clf;
%Datos conocidos
t0 = 0;
M0 = 10;
h1 = 0.1;
k = 1.24e-4;
lim = 0.08*M0;

%Metodo de Euler para h = 0.1
    %%Iteracion hasta obtener 0.08*M0
    t1(1) = t0; M1(1) = M0;
    i = 1;
    while M1(i) > lim
        M1(i+1) = M1(i)+h1*(-k*M1(i));
        t1(i+1) = t1(i)+h1;
        i = i+1;
    end

    %%Representacion grafica
    plot(t1,M1,'r')
    grid on
    title('Evolución de la cantidad de Carbono-14 (Método de Euler con paso h = 0.1)')
    xlabel('Tiempo   t    (años)')
    ylabel('Cantidad de Carbono-14   M(t)    (kg)')
    resp = sprintf('La edad de los restos arquelógicos para un paso h = 0.1 es: %d años',t1(end));
    disp(resp)

A través de este código, que muestra el resultado por pantalla, y como puede apreciarse en la gráfica, la datación de los restos otorga a los mismos una edad de [math]2.036870\ \times\ 10^{4}[/math] años.

A continuación, se repite la operación con paso [math]h_2 = 0.01\;[/math] para verificar el resultado. Se obtiene como solución la siguiente gráfica:

Al igual que en el caso anterior, la gráfica se obtiene por medio de un programa de Matlab que detiene la representación de valores al alcanzarse el valor [math]{0.08 \cdot M_0}\,[/math] de masa. El código utilizado es el siguiente:

clear all; clf;
%Datos conocidos
t0 = 0;
M0 = 10;
h2 = 0.01;
k = 1.24e-4;
lim = 0.08*M0;

%Metodo de Euler para h = 0.01
    %%Iteracion hasta obtenr 0.08*M0
    t2(1) = t0; M2(1) = M0;
    i = 1;
    while M2(i) > lim
        M2(i+1) = M2(i)+h2*(-k*M2(i));
        t2(i+1) = t2(i)+h2;
        i = i+1;
    end

    %%Representacion grafica
    plot(t2,M2,'r')
    grid on
    title('Evolución de la cantidad de Carbono-14 (Método de Euler con paso h = 0.01)')
    xlabel('Tiempo   t    (años)')
    ylabel('Cantidad de Carbono-14   M(t)    (kg)')
    resp = sprintf('La edad de los restos arquelógicos para un paso h = 0.01 es: %d años',t2(end));
    disp(resp)

Con este código, se obtiene como resultado de la datación un valor de [math]2.036877\ \times\ 10^{4}[/math] años.

Habiendo tomado dos pasos diferentes, se puede analizar si el método es estable. La estabilidad consiste en observar si el error entre las aproximaciones no se hace muy grande cuando la distancia entre los puntos de la malla se hace más pequeña.

Código MATLAB empleado para la obtención de las gráficas: clear all; clf; %Datos conocidos t0 = 0; M0 = 10; h1 = 0.1; h2 = 0.01; k = 1.24e-4; lim = 0.08*M0; %Metodo de Euler para h = 0.1 %%Iteracion hasta obtener 0.08*M0 t1(1) = t0; M1(1) = M0; i = 1; while M1(i) > lim M1(i+1) = M1(i)+h1*(-k*M1(i)); t1(i+1) = t1(i)+h1; i = i+1; end %Metodo de Euler para h = 0.01 %%Iteracion hasta obtenr 0.08*M0 t2(1) = t0; M2(1) = M0; i = 1; while M2(i) > 0.08*M0 M2(i+1) = M2(i)+h2*(-k*M2(i)); t2(i+1) = t2(i)+h2; i = i+1; end %Comparativa del método de Euler con pasos h=0.1 y h=0.01 hold on plot(t1,M1,'r') plot(t2,M2,'k--') grid on title('Comparativa del método de Euler con los pasos h = 0.1 y h = 0.01') xlabel('Tiempo t (años)') ylabel('Cantidad de Carbono-14 M(t) (kg)') legend('Método de Euler con paso h = 0.1','Método de Euler con paso h = 0.01') hold off

Como puede apreciarse en las gráficas, el error entre ambas aproximaciones es muy pequeño, produciéndose una variación entre las gráficas de [math]10^{-5}[/math] en el eje de la masa. Las soluciones obtenidas para la datación tienen una diferencia del orden de [math]10^{-2}[/math], en relación a la unidad tomada (años). Por tanto, el método es estable en este caso.

Al ser estable, se analiza si el resultado de la datación se ve influenciado por la condición inicial tomada.

Código MATLAB empleado para obtener la gráfica: clear all; clf; %Comprobacion de que el resultado no depende de la condicion inicial %Datos conocidos t0 = 0; h = 0.1; k = 1.24e-4; %Metodo de Euler para h = 0.1 M01 = 10; lim1 = 0.08*M01; t1(1) = t0; M1(1) = M01; i = 1; while M1(i) > lim1 M1(i+1) = M1(i)+h*(-k*M1(i)); t1(i+1) = t1(i)+h; i = i+1; end M02 = 80; lim2 = 0.08*M02; t2(1) = t0; M2(1) = M02; i = 1; while M2(i) > lim2 M2(i+1) = M2(i)+h*(-k*M2(i)); t2(i+1) = t2(i)+h; i = i+1; end M03 = 35; lim3 = 0.08*M03; t3(1) = t0; M3(1) = M03; i = 1; while M3(i) > lim3 M3(i+1) = M3(i)+h*(-k*M3(i)); t3(i+1) = t3(i)+h; i = i+1; end %Representacion grafica hold on plot(t1,M1,'g') plot(t2,M2,'r') plot(t3,M3,'k') grid on title('Comparativa del resultado con diferentes masas iniciales') xlabel('Tiempo t (años)') ylabel('Cantidad de Carbono-14 M(t) (kg)') legend('Cantidad inicial M0 = 10','Cantidad inicial M0=80','Cantidad inicial M0=35')

Habiendo considerado como masas iniciales 10, 35 y 80 kg, se aprecia en la gráfica que en todos los casos el tiempo obtenido para la datación es el mismo.

2.2 Método del trapecio

En el caso del método del trapecio, para resolver el problema de valor inicial se emplea la siguiente fórmula:

       [math]M_{i+1} = M_i + \frac{h}{2} \cdot (-k \cdot M_i -k \cdot M_{i+1})[/math]

De la cual, despejando obtenemos el valor de la aproximación:

       [math]M_{i+1} = \frac{M_i \cdot (1 - \frac{k \cdot h}{2})}{1 + \frac{k \cdot h}{2}}[/math]

Para la resolución por este método, se han tomado los siguientes datos: masa inicial [math]M_0=10\;[/math] y paso [math]h = 0.1\;[/math].

El resultado obtenido se recoge en la siguiente gráfica:

La gráfica se ha obtenido, al igual que en el caso anterior, por medio de un programa de MATLAB, el cual finaliza la representación de valores cuando se alcanza la cantidad de carbono-14 que existe en el momento de datación (es decir, [math]{0.08 \cdot M_0}\,[/math]). El código empleado es el siguiente:

clear all; clf;
%Datos conocidos
t0 = 0;
M0 = 10;
h = 0.1;
k = 1.24e-4;
lim = 0.08*M0;

%Metodo del trapecio para h = 0.1
%%Iteracion hasta obtener 0.08*M0
t(1) = t0; M(1) = M0;
i = 1;
while M(i) > lim
    M(i+1) = (M(i)*(1-k*h/2))/(1+k*h/2);
    t(i+1) = t(i)+h;
    i = i+1;
end 

%%Representacion grafica
plot(t,M,'r')
grid on
title('Evolución de la cantidad de Carbono-14 (Método del trapecio con paso h = 0.1)')
xlabel('Tiempo   t    (años)')
ylabel('Cantidad de Carbono-14   M(t)    (kg)')
resp = sprintf('La edad de los restos arquelógicos para un paso h = 0.1 es: %d años',t(end));
disp(resp)

Con este codigo, el resultado que se muestra por pantalla es de [math]2.036880\ \times\ 10^{4}[/math] años.

Recopilando los resultados de los métodos empleados obtenemos:

Método empleado	Resultado obtenido para la datación
Euler con paso h = 0.1	[math]2.036870\ \times\ 10^{4}[/math] años
Euler con paso h = 0.01	[math]2.036877\ \times\ 10^{4}[/math] años
Trapecio con paso h = 0.1	[math]2.036880\ \times\ 10^{4}[/math] años

La tabla muestra que los resultados obtenidos por los diferentes métodos son muy similares, con errores muy pequeños, lo que permite considerar todos ellos como buenos.

3 Concepto de vida media y resolución por el método de Runge-Kutta de cuarto orden

Un concepto útil para definir la desintegración de un elemento radiactivo es su vida media, que se define como el tiempo que tarda en reducirse a la mitad. Considerando que el valor de la masa en el tiempo de vida media es:

       [math]M(t_{1/2}) = \frac{M_0}{2}[/math]

Se obtener la relación para obtener el valor de [math]t_{1/2}[/math], es decir, el tiempo de vida media.

       [math]M(t_{1/2}) = \frac{M_0}{2} = M_0 \cdot e^{-kt_{1/2}}[/math]

       [math]\frac{1}{2} = e^{-kt_{1/2}}[/math]

       [math]-kt_{1/2} = \ln \frac{1}{2} = - \ln{2}[/math]

[math]t_{1/2} = \frac{\ln{2}}{k}[/math]

En el caso estudiado, se obtiene un valor para el tiempo de vida media de:

       [math]t_{1/2} = \frac{\ln{2}}{k} = \frac{\ln{2}}{1.24\ \times\ 10^{-4}} = 5589.8966[/math] años.

A continuación, se comparará este valor obtenido de forma analítica con el obtenido a través de la aproximación de la función [math]M(t)[/math] mediante el método de Runge-Kutta de orden 4, el cual se rige por las siguientes ecuaciones:

       [math]M_{i+1} = M_i + \frac{h}{6} \cdot (K_1 + 2K_2 + 2K_3 + K_4)[/math]

Donde:        [math] \begin{cases} K_1 = -k \cdot M_i \\ K_2 = -k \cdot (M_i + \frac {1}{2} \cdot K_1 \cdot h) \\ K_3 = -k \cdot (M_i + \frac {1}{2} \cdot K_2 \cdot h) \\ K_4 = -k \cdot (M_i + K_3 \cdot h) \\ \end{cases} [/math]

Utilizando este método, se obtiene la siguiente gráfica:

El código de MATLAB empleado para resolver el problema por el método de Runge-Kutta y obtener la gráfica de la evolución de la masa, hasta obtener el valor de masa [math]{0.5 \cdot M_0}[/math] , es el siguiente:

clear all; clf;
%Datos conocidos
t0 = 0;
M0 = 10;
h = 0.1;
k = 1.24e-4;
vidamedia = 0.5*M0;

%Metodo de Runge-Kutta para h = 0.1
%%Iteracion hasta obtener la masa equivalente a la vida media
t(1) = t0; M(1) = M0;
i = 1;
while M(i) > vidamedia
    K1 = -k*M(i);
    K2 = -k*(M(i)+1/2*K1*h);
    K3 = -k*(M(i)+1/2*K2*h);
    K4 = -k*(M(i)+K3*h);
    M(i+1) = M(i)+1/6*h*(K1+2*K2+2*K3+K4);
    t(i+1) = t(i)+h;
    i = i+1;
end 

%%Representacion grafica
plot(t,M,'r')
grid on
title('Evolución de la cantidad de Carbono-14 (Método de Runge-Kutta con paso h = 0.1)')
xlabel('Tiempo   t    (años)')
ylabel('Cantidad de Carbono-14   M(t)    (kg)')
resp = sprintf('La vida media del carbono-14 de los restos arquelógicos para un paso h = 0.1 es: %d años',t(end));
disp(resp)

Con el método empleado, se obtiene una valor del tiempo de vida media del carbono-14 de los restos arqueológicos estudiados de [math]5.589900\ \times\ 10^{3}[/math] años. Comparando este valor con el obtenido mediante la fórmula analítica, la diferencia entre ambos valores es del orden de [math]10^{-2}[/math], con lo cual, se puede considerar la aproximación como buena. Este valor permite observar que la desintegración no es lineal, siendo esta de mayor relevancia al comienzo y disminuyendo su velocidad con el tiempo. En valores numéricos, el tiempo que tarda en reducirse a la mitad el carbono-14 es de 5589 años aproximadamente, siendo necesarios 20369 años para reducirse la cantidad hasta el 8% de la masa inicial.

4 Descomposición de un elemento A en otro C formándose un elemento intermedio B

4.1 Interpretación del fenómeno físico y modelización matemática

Se considera la situación en que se dispone de un elemento radiactivo A, el cual al descomponerse se transformará en un elemento C, apareciendo en dicha transformación un elemento B que se formará a partir de la desintegración de A y se desintegrará formando el componente C. Las reacciones de desintegración estarán reguladas por las siguientes constantes de desintegración:  [math]k_1[/math], que es la constante de desintegración de A, y [math]k_2[/math], que es la constante de desintegración de B.

La variación de la cantidad de masa del elemento A se regirá por la ecuación definida en el epígrafe 1, descomponiéndose proporcionalmente a la cantidad de masa en un instante determinado:  [math]M_A'(t) = -k_1 \cdot M_A(t)[/math] , con [math]k_1\gt0[/math]. La masa del elemento B dependerá tanto de la descomposición de A, la cual aporta masa para formarse B, y de la propia descomposición de B para formarse C. Este fenómeno queda recogido en la expresión:  [math]M_B'(t) = k_1 \cdot M_A(t) - k_2 \cdot M_B(t)[/math] , con [math]k_1, k_2\gt0[/math]. Por último, el elemento C se formará en función de la masa de B y de su constante de desintegración ([math]k_2[/math]), por lo que la expresión que determina la variación de C será:  [math]M_C'(t) = k_2 \cdot M_B(t)[/math] , siendo esta positiva dado que se produce la formación del componente.

Con estas ecuaciones se puede determinar el sistema que describe el fenómeno físico completo:

       [math] \begin{cases} M_A'(t) = -k_1 \cdot M_A(t) \\ M_B'(t) = k_1 \cdot M_A(t) - k_2 \cdot M_B(t) \\ M_C'(t) = k_2 \cdot M_B(t) \\ \end{cases} [/math]

4.2 Resolución numérica de un caso particular con constantes [math]k_1=5[/math] y [math]k_2=1[/math]

Con las constantes [math]k_1=5[/math] y [math]k_2=1[/math], se determina que la velocidad de desintegración de A será 5 veces mayor que la velocidad con que se desintegra B. Se supone que en el instante inicial, toda la masa es del elemento A, que se supone de valor 1, no habiendo masa inicial de los elementos B y C. El intervalo de tiempo a estudiar será el [math](0,10)[/math].

Para la resolución de este problema se emplearán los métodos de Euler y del trapecio.

4.2.1 Método de Euler

Para la resolución por este método se empleará un paso de valor [math]h = 0.1[/math]. Se obtiene la siguiente gráfica empleando el código de MATLAB de la derecha.

clear all; clf; %Datos conocidos t0 = 0; tN = 10; y0 = [1;0;0]; %Concentraciones iniciales de A;B;C h = 0.1; k1 = 5; %Cte. de desintegración de A k2 = 1; %Cte. de desintegración de B %Funciones que relacionan las cantidades de A, B y C f = inline('[-k1*y(1);-k2*y(2)+(k1*y(1));k2*y(2)]','t','y','k1','k2'); %Discretización en espacio y tiempo N = round((tN-t0)/h); t = linspace(t0,tN,N+1); %Inicialización de la matriz de soluciones y(:,1) = y0; %Iteración por el método de Euler for i = 1:N y(:,i+1) = y(:,i) + h*f(t(i),y(:,i),k1,k2); end %Representación gráfica hold on plot(t,y(1,:),'r') plot(t,y(2,:),'b') plot(t,y(3,:),'g') grid on xlabel('Tiempo t (años)') ylabel('Cantidad de A,B y C M(t) (kg)') title('Descomposicion de un elemento A en otro C, pasando por B, por el Método de Euler') legend('Elemento A','Elemento B','Elemento C') gtext('k_{1}=5 k_{2}=1') hold off

4.2.2 Método del trapecio

Utilizando el mismo paso que en el método de Euler, se obtiene la siguiente gráfica con el código MATLAB de la derecha.

clear all; clf; %Datos conocidos t0 = 0; tN = 10; y0 = [1;0;0]; %Concentraciones iniciales de A;B;C h = 0.1; k1 = 5; %Cte. de desintegración de A k2 = 1; %Cte. de desintegración de B %Funciones que relacionan las cantidades de A, B y C f = inline('[-k1*y(1);-k2*y(2)+(k1*y(1));k2*y(2)]','t','y','k1','k2'); %Discretización en espacio y tiempo N = round((tN-t0)/h); t = linspace(t0,tN,N+1); %Inicialiación de la matriz de soluciones y(:,1) = y0; %Matrices a utilizar I = eye(3); A = [-k1 0 0;k1 -k2 0;0 k2 0]; %Iteración por el método del trapecio for i = 1:N T = [0;0;0]; y(:,i+1) = (I-h/2*A)\(y(:,i)+h/2*(f(t(i),y(:,i),k1,k2)+T)); end %Representación gráfica hold on plot(t,y(1,:),'r') plot(t,y(2,:),'b') plot(t,y(3,:),'g') grid on xlabel('Tiempo t (años)') ylabel('Cantidad de A,B y C M(t) (kg)') title('Descomposicion de un elemento A en otro C, pasando por B, por el Método del trapecio') legend('Elemento A','Elemento B','Elemento C') gtext('k_{1}=5 k_{2}=1') hold off

4.2.3 Interpretación de resultados

Partiendo de la situación inicial en la que toda la masa se encuentra en A, siendo esta igual a 1, al comenzarse a desintegrar A, la masa de B empieza a crecer. Esta situación se mantendrá hasta que la masa que recibe B procedente de la desintegración de A se iguale con la masa que se desintegra de B para formar el compuesto C. Es decir:  [math]k_1 \cdot M_A = k_2 \cdot M_B[/math]. Como [math]k_1=5[/math] y [math]k_2=1[/math], esta situación se dará cuando la masa de B sea 5 veces superior a la masa de A. Una vez superado ese momento, la masa tanto de A como de B decrece en favor de la formación de C, creciendo este último de forma exponencial hasta alcanzar el valor de 1 cuando los otros dos elementos se han descompuesto completamente.

En cuanto a las aproximaciones por los métodos de Euler y del trapecio, estos se diferencian únicamente en el resultado de masa máxima para el elemento B, pero podemos extrapolar la misma información de ambos, pues se obtiene el mismo comportamiento antes descrito.

4.3 Resolución númerica intercambiando el valor de las contantes de desintegración ([math]k_1=1[/math] y [math]k_2=5)[/math]

Al intercambiar las contantes de desintegración se invierten las velocidades de desintegración, de forma que B se desintegra 5 veces más rápido que A. Al igual que en el caso anterior, suponemos que en el instante inicial toda la masa es del elemento A e igual a 1, no habiendo masa de los elementos B y C. El intervalo de tiempo a estudiar será el mismo:  [math](0,10)[/math].

Se resolverá por los métiodos de Euler y del trapecio, utilizando como paso el valor [math]h=0.1[/math].

4.3.1 Método de Euler

clear all; clf; %Datos conocidos t0 = 0; tN = 10; y0 = [1;0;0]; %Concentraciones iniciales de A;B;C h = 0.1; k1 = 1; %Cte. de desintegración de A k2 = 5; %Cte. de desintegración de B %Funciones que relacionan las cantidades de A, B y C f = inline('[-k1*y(1);-k2*y(2)+(k1*y(1));k2*y(2)]','t','y','k1','k2'); %Discretización en espacio y tiempo N = round((tN-t0)/h); t = linspace(t0,tN,N+1); %Inicialización de la matriz de soluciones y(:,1) = y0; %Iteración por el método de Euler for i = 1:N y(:,i+1) = y(:,i)+h*f(t(i),y(:,i),k1,k2); end %Representación gráfica hold on plot(t,y(1,:),'r') plot(t,y(2,:),'b') plot(t,y(3,:),'g') grid on xlabel('Tiempo t (años)') ylabel('Cantidad de A,B y C M(t) (kg)') title('Descomposicion de un elemento A en otro C, pasando por B, por el Método de Euler') legend('Elemento A','Elemento B','Elemento C') gtext('k_{1}=1 k_{2}=5') hold off

4.3.2 Método del trapecio

clear all; clf; %Datos conocidos t0 = 0; tN = 10; y0 = [1;0;0]; %Concentraciones iniciales de A;B;C h = 0.1; k1 = 1; %Cte. de desintegración de A k2 = 5; %Cte. de desintegración de B %Funciones que relacionan las cantidades de A, B y C f = inline('[-k1*y(1);-k2*y(2)+(k1*y(1));k2*y(2)]','t','y','k1','k2'); %Discretización en espacio y tiempo N = round((tN-t0)/h); t = linspace(t0,tN,N+1); %Inicialiación de la matriz de soluciones y(:,1) = y0; %Matrices a utilizar I = eye(3); A = [-1 0 0;1 -5 0;0 5 0]; %Iteración por el método del trapecio for i = 1:N T = [0;0;0]; y(:,i+1) = (I-h/2*A)\(y(:,i)+h/2*(f(t(i),y(:,i),k1,k2)+T)); end %Representación gráfica hold on plot(t,y(1,:),'r')% plot(t,y(2,:),'b')% plot(t,y(3,:),'g') grid on xlabel('Tiempo t (años)') ylabel('Cantidad de A,B y C M(t) (kg)') title('Descomposicion de un elemento A en otro C, pasando por B, por el Método del trapecio') legend('Elemento A','Elemento B','Elemento C') gtext('k_{1}=1 k_{2}=5') hold off

4.3.3 Interpretación de resultados

Se parte de la misma situación inicial que el primer caso estudiado, con toda la masa en A, comenzandose a desintegrar A, pero esta vez a menos velocidad que en el caso anterior, al haber disminuido su constante de desintegración. La masa de B crecerá hasta que, igual que anteriormente, la masa procedente de la desintegración de A y la que se desintegra de B sean iguales, lo que se determina con la expresión:  [math]k_1 \cdot M_A = k_2 \cdot M_B[/math]. En el caso actual, esto implica que la masa de B crecerá hasta que la masa de A sea 5 veces superior a la masa de B. Por ello, no se cruzan las curvas de masa de A y B en este caso. Al igual que anteiormente, la masa de A y B comenzarán a decrecer a partir de este instante en favor de la formación de C, creciendo su valor de forma exponencial hasta alcanzar el valor de 1.

En este caso, las aproximaciones del método de Euler y del método del trapecio son más similares y reflejan este proceso descrito de manera similar.