MateWiki - Contribuciones del usuario [es]

Desintegración Radiactiva.Grupo 5

2015-03-01T18:13:50Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero

Araceli Martín Candilejo

Noemí Palomino Bustos

Teresa Quintana Romero

Alvaro Ramón López

Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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('Tiempo (años)');
ylabel('Cantidad de Carbono 14');
legend('Euler explícito','Location','best');
}}

 
[[Image:Eulerpaso0.1b.png|500px|thumb|left|Evolución de la concentración del C14 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1 ]] [[Image:Eulerpaso0.0b1.png|500px|thumb|center| Evolución de la concentración del C14 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.01]]

 

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14:

 
[[Image:Euler_paso_0.1_conc20.png|500px|thumb|left|Evolución de la concentración del C14,desde una concentración inicial=20 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1 ]] [[Image:Euler_paso_0.1_conc60.png|500px|thumb|center|Evolución de la concentración del C14,desde una concentración inicial=60 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1]]

 

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia(2).png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

===Sistemas de ecuaciones lineales: Método de Euler y trapecio===
Le daremos unos valores a las constantes para resolver nuestro sistema, siendo K1=5, y K2=1.

'''Método de Euler'''

Resolvemos el sistema anterior con el método de Euler:

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]

'''Método trapecio'''

Ahora aplicaremos el método del trapecio para resolver el sistema:

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

Podemos ver que, el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer.Hemos tomado como tiempo el eje de abscisas, y la cantidad del como eje de ordenadas. Tomando en tanto por uno la cantidad.
A decrece mas rápidamente, es decir, desaparece antes, B llega a una cantidad de 0.6.
En conclusión el compuesto A desaparece de forma rápida y esto da lugar a que haya mas cantidad de B.

===Sistemas de ecuaciones lineales: Método de Euler y trapecio (Constantes de integración intercambiadas)===

Ahora, resolveremos el sistema pero intercambiando los valores de las constantes, es decir K1=1 y K2=5, y aplicando los mismos programas anteriores, Euler y trapecio.

'''Método de Euler'''

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

'''Método del trapecio'''

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

En todas las gráficas, con las K iniciales y con las K intercambiadas y con ambos métodos, podemos observar que el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer. En esas últimas gráficas también hemos tomado como tiempo el eje de abscisas, y en tanto por uno la cantidad del compuesto como eje de ordenadas.
Podemos ver que la cantidad en tanto por uno de B sin cambiar las constantes llega a una cantidad de 0.6, mientras que con el cambio de constantes no llega a 0.2.
En conclusión, al comparar los dos pares de gráficas observamos que ahora con el cambio de constantes A decrece más lentamente, y debido a eso la cantidad de B que se forma no llega a ser tan grande como la del apartado anterior, en el que el compuesto A desaparecía antes y eso daba lugar a más cantidad de B. La creación de C no difiere mucho entre los pares de gráficas, lo que significa que el tiempo que tarda A en descomponerse en C es independiente de las constantes de desintegración.

Desintegración Radiactiva.Grupo 5

2015-03-01T18:12:51Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero

Araceli Martín Candilejo

Noemí Palomino Bustos

Teresa Quintana Romero

Álvaro Ramón López

Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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('Tiempo (años)');
ylabel('Cantidad de Carbono 14');
legend('Euler explícito','Location','best');
}}

 
[[Image:Eulerpaso0.1b.png|500px|thumb|left|Evolución de la concentración del C14 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1 ]] [[Image:Eulerpaso0.0b1.png|500px|thumb|center| Evolución de la concentración del C14 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.01]]

 

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14:

 
[[Image:Euler_paso_0.1_conc20.png|500px|thumb|left|Evolución de la concentración del C14,desde una concentración inicial=20 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1 ]] [[Image:Euler_paso_0.1_conc60.png|500px|thumb|center|Evolución de la concentración del C14,desde una concentración inicial=60 hasta que se alcanza el 8% de la concentración inicial, según el método de Euler para un paso h=0.1]]

 

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia(2).png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

===Sistemas de ecuaciones lineales: Método de Euler y trapecio===
Le daremos unos valores a las constantes para resolver nuestro sistema, siendo K1=5, y K2=1.

'''Método de Euler'''

Resolvemos el sistema anterior con el método de Euler:

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]

'''Método trapecio'''

Ahora aplicaremos el método del trapecio para resolver el sistema:

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

Podemos ver que, el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer.Hemos tomado como tiempo el eje de abscisas, y la cantidad del como eje de ordenadas. Tomando en tanto por uno la cantidad.
A decrece mas rápidamente, es decir, desaparece antes, B llega a una cantidad de 0.6.
En conclusión el compuesto A desaparece de forma rápida y esto da lugar a que haya mas cantidad de B.

===Sistemas de ecuaciones lineales: Método de Euler y trapecio (Constantes de integración intercambiadas)===

Ahora, resolveremos el sistema pero intercambiando los valores de las constantes, es decir K1=1 y K2=5, y aplicando los mismos programas anteriores, Euler y trapecio.

'''Método de Euler'''

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

'''Método del trapecio'''

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

En todas las gráficas, con las K iniciales y con las K intercambiadas y con ambos métodos, podemos observar que el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer. En esas últimas gráficas también hemos tomado como tiempo el eje de abscisas, y en tanto por uno la cantidad del compuesto como eje de ordenadas.
Podemos ver que la cantidad en tanto por uno de B sin cambiar las constantes llega a una cantidad de 0.6, mientras que con el cambio de constantes no llega a 0.2.
En conclusión, al comparar los dos pares de gráficas observamos que ahora con el cambio de constantes A decrece más lentamente, y debido a eso la cantidad de B que se forma no llega a ser tan grande como la del apartado anterior, en el que el compuesto A desaparecía antes y eso daba lugar a más cantidad de B. La creación de C no difiere mucho entre los pares de gráficas, lo que significa que el tiempo que tarda A en descomponerse en C es independiente de las constantes de desintegración.

Desintegración Radiactiva.Grupo 5

2015-03-01T15:36:00Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia(2).png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

===Sistemas de ecuaciones lineales: Método de Euler y trapecio===
Le daremos unos valores a las constantes para resolver nuestro sistema, siendo K1=5, y K2=1.

'''Método de Euler'''

Resolvemos el sistema anterior con el método de Euler:

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]

'''Método trapecio'''

Ahora aplicaremos el método del trapecio para resolver el sistema:

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

Podemos ver que, el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer.Hemos tomado como tiempo el eje de abscisas, y la cantidad del como eje de ordenadas. Tomando en tanto por uno la cantidad.
A decrece mas rápidamente, es decir, desaparece antes, B llega a una cantidad de 0.6.
En conclusión el compuesto A desaparece de forma rápida y esto da lugar a que haya mas cantidad de B.

===Sistemas de ecuaciones lineales: Método de Euler y trapecio (Constantes de integración intercambiadas)===

Ahora, resolveremos el sistema pero intercambiando los valores de las constantes, es decir K1=1 y K2=5, y aplicando los mismos programas anteriores, Euler y trapecio.

'''Método de Euler'''

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

'''Método del trapecio'''

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

En todas las gráficas, con las K iniciales y con las K intercambiadas y con ambos métodos, podemos observar que el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer. En esas últimas gráficas también hemos tomado como tiempo el eje de abscisas, y en tanto por uno la cantidad del compuesto como eje de ordenadas.
Podemos ver que la cantidad en tanto por uno de B sin cambiar las constantes llega a una cantidad de 0.6, mientras que con el cambio de constantes no llega a 0.2.
En conclusión, al comparar los dos pares de gráficas observamos que ahora con el cambio de constantes A decrece más lentamente, y debido a eso la cantidad de B que se forma no llega a ser tan grande como la del apartado anterior, en el que el compuesto A desaparecía antes y eso daba lugar a más cantidad de B. La creación de C no difiere mucho entre los pares de gráficas, lo que significa que el tiempo que tarda A en descomponerse en C es independiente de las constantes de desintegración.

Archivo:Rungekuttavidamedia(2).png

2015-03-01T15:34:53Z

Teresa Quintana Romero:

Desintegración Radiactiva.Grupo 5

2015-03-01T15:28:06Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia.png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

===Sistemas de ecuaciones lineales: Método de Euler y trapecio===
Le daremos unos valores a las constantes para resolver nuestro sistema, siendo K1=5, y K2=1.

'''Método de Euler'''

Resolvemos el sistema anterior con el método de Euler:

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]

'''Método trapecio'''

Ahora aplicaremos el método del trapecio para resolver el sistema:

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

Podemos ver que, el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer.Hemos tomado como tiempo el eje de abscisas, y la cantidad del como eje de ordenadas. Tomando en tanto por uno la cantidad.
A decrece mas rápidamente, es decir, desaparece antes, B llega a una cantidad de 0.6.
En conclusión el compuesto A desaparece de forma rápida y esto da lugar a que haya mas cantidad de B.

===Sistemas de ecuaciones lineales: Método de Euler y trapecio (Constantes de integración intercambiadas)===

Ahora, resolveremos el sistema pero intercambiando los valores de las constantes, es decir K1=1 y K2=5, y aplicando los mismos programas anteriores, Euler y trapecio.

'''Método de Euler'''

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

'''Método del trapecio'''

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

En todas las gráficas, con las K iniciales y con las K intercambiadas y con ambos métodos, podemos observar que el compuesto A decrece, el C crece y el B empieza creciendo para mas tarde decrecer y desaparecer. En esas últimas gráficas también hemos tomado como tiempo el eje de abscisas, y en tanto por uno la cantidad del compuesto como eje de ordenadas.
Podemos ver que la cantidad en tanto por uno de B sin cambiar las constantes llega a una cantidad de 0.6, mientras que con el cambio de constantes no llega a 0.2.
En conclusión, al comparar los dos pares de gráficas observamos que ahora con el cambio de constantes A decrece más lentamente, y debido a eso la cantidad de B que se forma no llega a ser tan grande como la del apartado anterior, en el que el compuesto A desaparecía antes y eso daba lugar a más cantidad de B. La creación de C no difiere mucho entre los pares de gráficas, lo que significa que el tiempo que tarda A en descomponerse en C es independiente de las constantes de desintegración.

Desintegración Radiactiva.Grupo 5

2015-02-27T19:23:20Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia.png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

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

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]
===Sistemas de ecuaciones lineales: Método del Trapecio.===

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

Ahora aplicaremos los mismos programas, pero intercambiando los valores de K, es decir K1=1 y K2=5

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

Desintegración Radiactiva.Grupo 5

2015-02-27T19:21:34Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia.png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

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

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]
===Sistemas de ecuaciones lineales: Método del Trapecio.===

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

==Método de Euler y trapecio (K cambiadas)==
Ahora aplicaremos los mismos programas, pero intercambiando los valores de K, es decir K1=1 y K2=5

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

Desintegración Radiactiva.Grupo 5

2015-02-27T19:05:43Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia.png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
[[Archivo:SistEcuacionC14.jpg|miniaturadeimagen|centro]]

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

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]
===Sistemas de ecuaciones lineales: Método del Trapecio.===

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

Archivo:SistEcuacionC14.jpg

2015-02-27T19:03:27Z

Teresa Quintana Romero:

Desintegración Radiactiva.Grupo 5

2015-02-27T18:52:32Z

Teresa Quintana Romero:

{{ TrabajoED | Desintegración Radiactiva. Grupo 5. | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Diego García Vaquero, Araceli Martín Candilejo, Noemí Palomino Bustos, Teresa Quintana Romero, Álvaro Ramón López, Mercedes Ruiz Barrajón. }}

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Interpretación de las funciones M(t) y la constante K==

Una reacción química se denomina reacción de primer orden si en ella una molécula se descompone en otras espontáneamente, y el número de moléculas que se descompone en una unidad de tiempo es proporcional al número de moléculas existentes. Si se considera una sustancia cuya masa se descompone en función del tiempo según una función m=m(t), la velocidad de descomposición viene dada por la derivada de m(t) respecto de t. Si se supone que esta velocidad es directamente proporcional a la masa, se tiene que; dm/dt = -k.m (k>0 coeficiente de proporcionalidad).

La constante k se denomina constante de rapidez ya que su valor indica una medida de la velocidad a la que se realiza la reacción.

==Método Euler==

Suponenos que un arqueólogo descubre huesos con un contenido en C14 que resulta ser del 8% del que se encuentra en un ser vivo. Si suponemos además que la cantidad de C14 en la atmósfera no ha variado podemos tomar la diferencia de contenido en C14 del hueso antiguo debida únicamente a su desintegración. Conociendo que la constante de desintegración del C14 es 1.24 × 10−4 por año, calcularemos la edad de los restos arqueológicos. Para ello, planteamos un PVI adecuado eligiendo la condición inicial y lo resolveremos por el método de Euler para diferentes pasos h = 0.1 y h = 0.01.

{{matlab|codigo=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');
}}

Aplicando Euler se comprueba que se puede considerar estable la edad de los restos arqueológicos, es decir, el tiempo de desintegración que obtenemos es independiente de la cantidad inicial de la muestra de C14.

==Resolución por el método del Trapecio==

Resolveremos el apartado anterior con el método del trapecio utilizando un paso h = 0.1 y con la condición de que el programa se detenga cuando la masa alcance un 8% de la masa inicial. Usando un bucle "while" conseguimos realizar un programa que cumpla esta condición.

{{matlab|codigo=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');
}}

[[Image:Trapecio3)cuandoy=0.08.png|1000px|thumb|center|Evolución del contenido de C14 según el método del Trapecio con un paso h=0.1 . ]]

[[Image:Resultado_del_tiempo_del_trapecio.png|500px|thumb|center|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 . ]]

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

El conocimiento de la vida media de los elementos radiactivos que hay en la naturaleza se utiliza para asignar fechas a acontecimientos que ocurrieron hace mucho tiempo. Es usual expresar la descomposición de un elemento radiactivo en función de su vida media, es decir, el tiempo necesario para que una cantidad dada del elemento se reduzca a la mitad. Para ello aplicaremos el método de Runge-Kutta.

{{matlab|codigo=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
}}

[[Image:Rungekuttavidamedia.png|1000px|thumb|center|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.

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

Consideraremos ahora una descomposición de un elemento A en otro C a través de un elemento o isótopo intermedio B
[[Archivo:DescomposicionC14.jpg|miniaturadeimagen|centro]]
donde k1 y k2 son las constantes de desintegración respectivas. Es importante considerar que el elemento C empieza a crearse en cuanto hay elemento B, es decir, podemos tener simultáneamente los 3 elementos. Con este concepto determinaremos el sistema de ecuaciones que nos permitirán conocer las cantidades de cada elemento en cada instante de tiempo.
===Sistemas de ecuaciones lineales: Método de Euler.===

{{matlab|codigo=
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
}}

[[Image:Metodo_de_eulerK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler con un paso de h=0.1. ]]
===Sistemas de ecuaciones lineales: Método del Trapecio.===

{{matlab|codigo=
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
}}

[[Image:Método_del_trapecioK.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método del Trapecio con un paso de h=0.1. ]]

[[Image:Eulerconstantesvariadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Euler,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

[[Image:Trapecioconstantescambiadas.png|1000px|thumb|center|Evolución del contenido de los compuestos A, B y C a lo largo del tiempo, según el método de Trapecio,con las constantes de desintegración cambiadas, con un paso de h=0.1. ]]

Desintegración Radiactiva.Grupo 5

2015-02-27T18:37:39Z

Teresa Quintana Romero:

Archivo:DescomposicionC14.jpg

2015-02-27T18:20:49Z

Teresa Quintana Romero:

Desintegración Radiactiva.Grupo 5

2015-02-27T18:18:26Z

Teresa Quintana Romero:

Desintegración Radiactiva.Grupo 5

2015-02-27T18:08:45Z

Teresa Quintana Romero:

Desintegración Radiactiva.Grupo 5

2015-02-27T17:44:14Z

Teresa Quintana Romero:

Prueba1

2014-12-05T15:50:15Z

Teresa Quintana Romero:

Visualización de campos escalares y vectoriales en elasticidad. Grupo C30

2014-12-05T15:49:53Z

Teresa Quintana Romero:

{{ TrabajoED | Visualización de campos escalares y vectoriales en un sólido (Grupo C30) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] |
María García Fernández

Sergio Ortega Pajares

Noemí Palomino Bustos

Diego Paramio Sastre

Álvaro Ramón López
}}

==Introducción==
{{beta}}
En este trabajo estudiaremos la visualización de campos escalares y vectoriales en elasticidad de la placa plana que ocupa la región comprendida entre estas dos parábolas:

<math>P1: 18y-81x^2-1=0 </math>
:
<math>P2: 2y+x^2-1=0 </math>

Utilizaremos un cambio de coordenadas para hacer más sencilla su representación (sistema de coordenadas adaptado a la geometría dada):

<math>x=uv</math>
:
<math>y=1/2(u^2-v^2)</math>

con u,v definidas en el intervalo [1/3,1]x[-1,1]

Comenzaremos dibujando el mallado del sólido, así como sus líneas coordenadas y los vectores de la base natural en cada punto. Conociendo la fórmula de la temperatura de la placa seremos capaces de representar sus curvas de nivel y su gradiente.
Después aplicaremos un campo de desplazamientos a nuestro sólido y observaremos como ha variado tras el desplazamiento.Calcularemos el tensor de tensiones y lo representaremos, utilizando este tensor llegaremos a obtener la tensión de Von Mises.Finalmente calcularemos numéricamente(método de los rectángulos) la masa de la placa.

==Mallado de los puntos interiores del sólido==
Hemos representado los puntos interiores del sólido mediante un mallado,considerando un paso de muestreo h=1/20 para las coordenadas u,v que hemos definido en la introducción.

.[[Archivo:MalladoC30.jpg|325px|miniaturadeimagen|derecha|Aquí vemos el resultado de dibujar el mallado]]
{{matlab|codigo=
h=1/20; %paso de muestreo
u=1/3:h:1; %intervalo de u
v=-1:h:1; %intervalo de v
[uu,vv]=meshgrid(u,v);%mallado de la gráfica
figure(1)
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2); %cambio de coorde
mesh(xx,yy,0*xx)
axis([-1,1,-1,1])
view(2)
}}

===Líneas coordenadas y vectores de la base natural ===
Los vectores de la base natural son aquellos obtenidos al derivar las compentes x e y respecto de u y v.
Las llamaremos <math>\vec{g_u}</math> y <math>\vec{g_u}</math>
:
<math>\vec{g_u}= v\vec{i}+u\vec{j}</math>
:
<math>\vec{g_v}=u\vec{i}-v\vec{j}</math>

Para hacer más sencillo su interpretación hemos dividido <math>\vec{g_u}</math> y <math>\vec{g_u}</math> en partes.
A continuación vemos el código y la representación gráfica de las líneas coordenadas:
.[[Archivo:LineascoordenadasC30.jpg|300px|miniaturadeimagen|derecha|Líneas de gu (en azul) y líneas de gv (en rojo)]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
guu=vv;
guv=uu;
gvu=uu;
gvv=-vv;
hold on
quiver(xx,yy,guu,guv)%representamos los vectores gu en cada punto del mallado
quiver(xx,yy,gvu,gvv)%representamos los vectores gv en cada punto del mallado
hold off
}}

==Temperatura del sólido==
La temperatura del sólido viene dada por la función <math>T(x,y)=e^{(-y)}</math>.Hemos representado sus curvas de nivel, en el gráfico podemos observar que la temperatura es simétrica en la placa respecto al eje OY, siendo la parte más cálida el extremo inferior y las más fría la superior.
Podríamos obtener el punto más cálido de la placa y el punto más frío calculando sus derivadas parciales e igualandolas a 0, si los resultados quedasen fuera de la placa, buscariamos máximos y mínimos en la frontera.
[[Archivo:TemperaturaC30.jpg|300px|miniaturadeimagen|derecha|Distribución de las temperaturas en el sólido)]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1; v=-1:h:1;
%mallado de la gráfica
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
T=exp(-yy); %Función de la Temperatura
figure(1)
surf(xx,yy,T)
axis([-1,1,-1,1])
view(2)
}}

El gradiente de un campo indica la dirección en la cual el campo varía más rápidamente, se denota con <math>\nabla </math>.Las curvas de nivel de un campo son aquellas que conecta los puntos que tienen un mismo valor constante.
De forma geométrica el gradiente es un vector que se encuentra normal(perpendicular)a a las curvas de nivel,como se puede observar en el siguiente gráfico.

El gradiente y las curvas de nivel de dicho campo (T) son los siguientes:
[[Archivo:GradTemperaturaC30.jpg|400px|miniaturadeimagen|derecha|Distribución de las temperaturas en el sólido)]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1; v=-1:h:1;
%mallado de la gráfica
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
T=exp(-yy); %Función de la Temperatura
Tx=-xx.*0;% derivada parcial en x
Ty=-yy.*exp(-yy);% derivada parcial en y
figure(1)
quiver(xx,yy,Tx,Ty) %Gradiente de la temperatura
axis([-1,1,-1,1])
figure(2)
contour(xx,yy,T,20)%dibujo curvas de nivel
axis([-1,1,-1,1])
view(2)

}}

== Desplazamiento del sólido==

Considerando el vector desplazamiento <math>\vec{u}(u,v)=\vec{a}(\vec{b}·\vec{r_0})^2</math> siendo <math>\vec{a}=\frac{\vec{g_u}}{|\vec{g_u}|}</math> <math>\vec{b}=\frac{-4\vec{g_v}}{|\vec{g_v}|}</math> y <math>\vec{r_0}=\vec{x}+\vec{y}</math>

Cambiaremos <math>\vec{u}</math> a coordenadas curvilíneas

<math>\vec{u}=(uv,\frac{1}{2}(v^2+u^2))\begin{pmatrix} v\over(v^2+u^2) & u\over( v^2+u^2) \\u\over( v^2+u^2) & - v\over(v^2+u^2) \end{pmatrix} </math>
:
<math>\vec{u}=4·v^2·(u^2+v^2)·\vec{g_u}</math>

Si aplicamos <math>\vec{u}</math> a nuestra placa sufrirá una deformación producida por una fuerza determinada, el vector posición después de la deformación viene dada por <math>\vec{r}(u,v)=\vec{r_0}(u,v)+\vec{u}(u,v)</math>. Aplicando el campo al sólido inicial obtenemos el sólido deformado.

A continuación tenemos el sólido antes y después del desplazamiento dado por el campo de vectores <math>\vec{u}</math>
:
[[Archivo:DespC30.jpg|500px|miniaturadeimagen|centro|Representación del campo y de su desplazamiento]]

==Divergencia y rotacional de <math>\vec{u}</math>==

La divergencia de un campo vectorial mide la diferencia entre el flujo saliente y el flujo entrante.
Si la divergencia es positiva, es una fuente, si es negativa es un sumidero.

La divergencia del desplazamiento<math>\vec{u}</math> la hemos manualmente a través de la fórmula
<math>\nabla \cdot \vec{u}=\frac{1}{ \sqrt{g} } \frac{\partial}{\partial u^i} \sqrt{g} u^{i} </math>

<math>\nabla\cdot\vec{u}=\frac{1}{u^2+v^2}( \frac{\partial}{\partial u}(u^2+v^2)(4v^2)+ (\frac{\partial}{\partial v}((u^2+v^2)(0))=\frac{12uv^2}{\sqrt{u^2+v^2}}</math>

[[Archivo:DivergenciaC30.jpg|300px|miniaturadeimagen|derecha|Representación de la divergencia del campo]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
div=12.*uu.*vv.^2./sqrt(uu.^2+vv.^2);
surf(xx,yy,div)
axis([-1,1,-1,1])
view(2)
}}
Vemos que la divergencia alcanza sus valores más altos en aquellos puntos en los que la variación de volumen es mayor, en este caso se corresponden con los extremos de la placa.

Vamos a calcular ahora el rotacional de <math>\vec{u}</math>:
<math>\nabla\times\vec{u} = det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ u_1 & u_2 & u_3 \end{pmatrix} = det\begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ 4v^3 \sqrt{u^2+v^2} & 4uv^2 \sqrt{u^2+v^2} & 0 \end{pmatrix} =-4v^2(u^2+3v^2)/\sqrt{u^2+v^2} \vec{k} </math>
Vemos en su representación que los puntos con mayor rotacional (los de color más amarillo) representarán los puntos del sólido en los que se induce un mayor giro debido a la deformación inducida por el campo <math>\vec{u}</math>
[[Archivo:RotacionalC30.jpg|300px|miniaturadeimagen|derecha|Representación del rotacional del campo]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1;
v=-1:h:1;
figure(1)
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
mesh(xx,yy,0*xx)
rot=4*(vv.^2).*sqrt(uu.^2+vv.^2);
surf(xx,yy,rot)
axis([-1,1,-1,1])
view(2)
}}

==Tensiones==
Vamos a ver la representación del tensor de tensiones resultante de aplicar el campo <math>\vec{u} </math> al sólido. Dicho tensor de tensiones viene dado por la fórmula:
<math>\sigma= \lambda \nabla \cdot \vec{u} 1 + 2 \mu \epsilon </math>
donde <math>\epsilon = \frac{(\nabla \vec{u} + \nabla \vec{u}^t)}{2}</math> . <math> \lambda </math> y <math> \mu </math> son los parámetros de Lamé que dependen de las propiedades del material, en nuestro caso ambos valen <math> 1 </math>.
Definimos <math>\epsilon </math> y <math>\sigma</math> ( la divergencia ya la conocemos):
<math>\epsilon =\begin{pmatrix} \frac{4uv^3}{\sqrt{u^2+v^2}} & 10v^2\sqrt{u^2+v^2} \\ 10v^2\sqrt{u^2+v^2} & \frac {4uv(2u^2+3v^2}{\sqrt{u^2+v^2}} \end{pmatrix}</math>
Con la ayuda de MatLab obtenemos las tensiones normales en las direcciones <math>g_u</math> y <math>g_v </math>:
[[Archivo:gugv.jpg|500px|miniaturadeimagen|centro|Tensiones normales en las direcciones de <math>g_u</math> y <math>g_v </math> ]]

{{matlab|codigo=
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
sigmauu=12*uu.*vv.^2./sqrt(uu.^2+vv.^2)+2*4*vv.^3.*uu./sqrt(uu.^2+vv.^2);
sigmavv=12*uu.*vv.^2./sqrt(uu.^2+vv.^2)+2*4*uu.*vv.*(2*uu.^2+3*vv.^2)./sqrt(uu.^2+vv.^2);
sigmauv=12*uu.*vv.^2./sqrt(uu.^2+vv.^2)+2*10*vv.^2.*sqrt(uu.^2+vv.^2);
sigmavu=sigmauv;
subplot(1,2,1)
quiver(xx,yy,sigmauu,sigmauv)
axis([-1,1,-1,1])
subplot(1,2,2)
quiver(xx,yy,sigmauv,sigmavv)
axis([-1,1,-1,1])
}}
Podemos observar al comparar las tensiones y el rotacional que coinciden los puntos de mayor tensión con los de mayor rotacional.

==Masa de la placa==
Tengamos la siguiente función de densidad:
<math>d(x,y)=|x|e^{(-1/y^2)}</math>
Vamos a calcular la masa total del sólido integrando la función en el área del sólido de manera numérica dividiendo el sólido en pequeños rectángulos y sumando todo al final. Aquí vemos el codigo MatLab utilizado:
{{matlab|codigo=
h=1/100; %Tomamos un paso de 1/100 para conseguir un resultado mejor aproximado
u=1/3:h:1;
v=-1:h:1;
%mallado de la gráfica
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5.*((uu).^2-(vv).^2);
f=abs(xx).*exp(-1./yy.^2); %Definimos la función
a=h.^2*f; %Parte diferencial del sólido
masa=sum(sum(a)) %Suma de todas las partes del sólido
}}
Obtenemos una masa total de <math>Masa = 8.029 \cdot 10^-5</math>
[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]