MateWiki - Contribuciones del usuario [es]

Desintegración Radiactiva.Grupo 5

2015-03-01T16:42:38Z

Araceli Martin: /* Método Euler */

{{ 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.

Archivo:Eulerpaso0.0b1.png

2015-03-01T16:42:19Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-03-01T16:41:36Z

Araceli Martin: /* Método Euler */

{{ 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.01.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.

Archivo:Eulerpaso0.1b.png

2015-03-01T16:35:47Z

Araceli Martin:

Archivo:Euler paso 0.1 conc60.png

2015-03-01T16:34:27Z

Araceli Martin:

Archivo:Euler paso 0.1 conc20.png

2015-03-01T16:32:24Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-03-01T16:24:04Z

Araceli Martin: /* Método Euler */

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

 
[[Image:Eulerpaso0.1.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.01.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.

==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:Eulerpaso0.01.png

2015-03-01T16:22:31Z

Araceli Martin:

Archivo:Eulerpaso0.1.png

2015-03-01T16:11:04Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-26T10:23:02Z

Araceli Martin: /* Estudio de la evolución de concentraciones de compuestos interrelacionados */

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===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:Trapecioconstantescambiadas.png

2015-02-26T10:21:37Z

Araceli Martin:

Archivo:Eulerconstantesvariadas.png

2015-02-26T10:21:22Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-26T10:17:16Z

Araceli Martin:

Archivo:Metodo de eulerK.png

2015-02-26T10:14:53Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-26T10:13:14Z

Araceli Martin: /* Sistemas de ecuaciones lineales: Método del Trapecio. */

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===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)

M=[-5,0;5,-1];

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

===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. ]]

Desintegración Radiactiva.Grupo 5

2015-02-26T10:11:51Z

Araceli Martin:

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===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)

M=[-5,0;5,-1];

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

===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
%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: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. ]]

Archivo:Método del trapecioK.png

2015-02-26T10:10:11Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-26T10:06:35Z

Araceli Martin:

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===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)

M=[-5,0;5,-1];

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

===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
%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
}}

Desintegración Radiactiva.Grupo 5

2015-02-26T09:46:56Z

Araceli Martin: /* Sistemas de ecuaciones lineales: Método de Euler. */

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===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)

M=[-5,0;5,-1];

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

Desintegración Radiactiva.Grupo 5

2015-02-26T09:46:42Z

Araceli Martin: /* Sistemas de ecuaciones lineales: Método de Euler. */

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===Sistemas de ecuaciones lineales: Método de Euler.===

{{matlab|codigo=clear all
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)

M=[-5,0;5,-1];

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

Desintegración Radiactiva.Grupo 5

2015-02-26T09:46:16Z

Araceli Martin:

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

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

==Apartado 1==
{{matlab|codigo=clear all

%% TRABAJO 3 APARTADO 1
% M'(t)=-1,24*10^(-4)*M(t)
% M(0)=1

% M(t) representa la cantidad de C14 en un instante dado, como no nos
% dan una cantidad inicial M0 supondremos que esta es uno. Nuestro objetivo
% es determinar en que intanste queda un 8% de la cantidad inicial M0.7
% Como M0=1, es como si trabajaramos todo el rato en tanto por uno
% Después demostraremos que quedará un 8% de M0 pasado el tiempo que
% tratamos de determinar independientemente del M0 adoptado.

t0=0;
h=input('Inserte el valor del paso, por favor: ');
M0=input('Inserte la cantidad inicial de carbono 14, por favor: ');
t(1)=t0;
M(1)=M0;

% Euler explícito
i=1;
while M(i)>(0.08*M0)
M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
t(i+1)=t(i)+h;
i=i+1;
end

disp('Tiempo final:')
disp(t(end))
plot(t,M)
xlabel('Cantidad de Carbono 14');
ylabel('Tiempo (años)');
legend('Euler explícito','Location','best');
}}

==Resolución por el método del Trapecio==
{{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==
{{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==
===Sistemas de ecuaciones lineales: Método de Euler.===

{{matlab|codigo=clear all

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)

M=[-5,0;5,-1];

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

Desintegración Radiactiva.Grupo 5

2015-02-24T23:20:23Z

Araceli Martin:

Archivo:Rungekuttavidamedia.png

2015-02-24T23:17:57Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-24T23:17:36Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-24T23:02:17Z

Araceli Martin: /* Resolución por el método del Trapecio */

Desintegración Radiactiva.Grupo 5

2015-02-24T22:58:21Z

Araceli Martin: /* Resolución por el método del Trapecio */

Archivo:Resultado del tiempo del trapecio.png

2015-02-24T22:56:47Z

Araceli Martin:

Desintegración Radiactiva.Grupo 5

2015-02-24T22:55:27Z

Araceli Martin: /* Resolución por el método del Trapecio */

Desintegración Radiactiva.Grupo 5

2015-02-24T22:54:38Z

Araceli Martin: /* Apartado 1 */

Archivo:Trapecio3)cuandoy=0.08.png

2015-02-24T22:53:19Z

Araceli Martin:

Campos en Elasticidad

2014-12-06T12:33:45Z

Araceli Martin: /* Divergencia de un campo */

{{ TrabajoED |Campos escalares y vectoriales en elasticidad. Grupo 16-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Araceli Martín, Juan Carlos Durán, Francisco Javier Alcaraz, Álvaro Llera, Clara Callejo, Manuel Escudero }}
[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Para este análisis y representación de campos escalares en elasticidad nos situaremos en el contexto de una '''placa plana''' que ocupa la región comprendida entre las parábolas :

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

Para representar la placa se utilizará un sistema de coordenadas curvilíneas tal que:
* <math>x=uv</math>
* <math> y = \dfrac{1}{2}(u^2-v^2)</math>

En nuestro análisis, el dominio de '''<math>u</math>''' y '''<math>v</math>''' comprenderá:
<math>(u,v) \in [1/3,1]*[-1,1]</math>

==º Situación inicial de la placa==

Primero comenzaremos el estudio de la placa en cuestión planteando su representación gráfica mediante un mallado, utilizando, en el código, la conversión a coordenadas curvilíneas de '''<math>u</math>''' y '''<math>v</math>'''. El intervalo en el que representaremos comprende::

<math>(x,y) \in [-1,1]*[-1,1]</math>

Para discretizar los vectores <math>u</math> y <math>v</math> utilizaremos un paso <math> h = \dfrac{1}{20}</math>.

{{matlab|codigo=
h=1/20; % Paso de muestreo.
u=1/3:h:1; % Intervalo [1/3,1].
v=-1:h:1; % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v); % Malla de datos.
xx=uu.*vv; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
plot(xx,yy); % Imagen.
mesh(xx,yy,0*xx) % Mallado.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
}}

Los resultados de la representación se ven en la imagen:

[[Archivo:Pepemola.jpg|800px|thumb|centre|Mallado que representa la superficie de la placa comprendida entre las parábolas <math>P1</math> y <math>P2</math>.]]
 
=== Líneas Coordenadas===
Cuando se hace uso de cambio de coordenadas a ciertas coordenadas curvilíneas, es útil para el entendimiento de la transformación la representación de las '''Líneas Coordenadas'''. Las líneas coordenadas son aquellas que se obtienen variando una de las coordenadas de la transformación (<math>u</math> o <math>v</math>) y manteniendo fija la restante.
En este estudio hemos representado varias líneas coordenadas a base de dar un valor concreto a <math>u</math> o a <math>v</math>, dentro de sus respectivos intervalos:

{{matlab|codigo=
xx11=uu.*0.5 ;xx12=uu.*-0.5;xx13=uu.*1;xx14=uu.*-1;
xx15=uu.*0.75;xx16=uu.*-0.75;xx17=uu.*0; % Parametrización X fijando v (0.5/-0.5/1/-1/0.75/-0.75/0) y variando u.
yy11=(1/2).*((uu.^2)-(0.5.^2));yy12=(1/2).*((uu.^2)-((-0.5).^2));
yy13=(1/2).*((uu.^2)-(1.^2));yy14=(1/2).*((uu.^2)-((-1).^2));
yy15=(1/2).*((uu.^2)-(0.75.^2));yy16=(1/2).*((uu.^2)-((-0.75).^2));
yy17=(1/2).*((uu.^2)-(0.^2)); % Parametrización Y fijando v (0.5/-0.5/1/-1/0.75/-0.75/0) y variando u.
xx21=vv.*0.5;xx22=vv.*0.4;xx23=vv.*1;xx24=vv.*0.9;xx25=vv.*0.75;
xx26=vv.*0.65;xx27=vv.*(1/3); % Parametrización X fijando u (0.5/0.4/1/0.9/0.75/0.65/0.333) y variando v.
yy21=(1/2).*((0.5.^2)-(vv.^2));yy22=(1/2).*((0.4.^2)-(vv.^2));
yy23=(1/2).*((1.^2)-(vv.^2));yy24=(1/2).*((0.9.^2)-(vv.^2));
yy25=(1/2).*((0.75.^2)-(vv.^2));yy26=(1/2).*((0.65.^2)-(vv.^2));
yy27=(1/2).*(((1/3).^2)-(vv.^2)); % Parametrización X fijando u (0.5/0.4/1/0.9/0.75/0.65/0.333) y variando v.
subplot(1,2,1); % Muestra varias imágenes. 1ª Imagen.
hold on % Inicio superposición de gráficos.
mesh(xx11,yy11,0*xx);mesh(xx12,yy12,0*xx);mesh(xx13,yy13,0*xx);
mesh(xx14,yy14,0*xx);mesh(xx15,yy15,0*xx);mesh(xx16,yy16,0*xx);
mesh(xx17,yy17,0*xx); % Mallados.
axis([-1,1,-1,1]) % Selecciona la regíon a dibujar.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.
subplot(1,2,2); % Muestra varias imágenes. 2º Imagen.
hold on % Inicio superposición de gráficos.
mesh(xx21,yy21,0*xx);mesh(xx22,yy22,0*xx);mesh(xx23,yy23,0*xx);
mesh(xx24,yy24,0*xx);mesh(xx25,yy25,0*xx);mesh(xx26,yy26,0*xx);
mesh(xx27,yy27,0*xx); % Mallados.
axis([-1,1,-1,1]) % Selecciona la regíon a dibujar.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.
}}
A continuación se muestra las gráficas resultantes:
[[Archivo:Feofeo.jpg|1000px|thumb|centre|Líneas coordenadas fijando la variable <math>v</math> y <math>u</math>, respectivamente.]]
 

=== Base Natural===

Cuando se realiza una transformación a coordenadas curvilíneas (de <math>x</math> e <math>y</math> a <math>u</math> y <math>v</math>), el vector de posición <math> \vec{r_o}</math> se expresa como:

<math> \vec{r_o}=x\hat{e_1} +y \hat{e_2} =uv\hat{e_1} +\dfrac{1}{2}(u^2-v^2) \hat{e_2} </math>

La base natural <math>\vec{g_u}, \vec{g_v}</math> es la que tiene por vectores el resultado de derivar el vector posición <math> \vec{r_o}</math> según las nuevas coordenadas <math>u</math> y <math>v</math>. Dado que es un problema sobre una placa plana, estamos en una situación de dos dimensiones en la que para cualquier base, sólo se requieren dos vectores, (<math>\vec{g_u}, \vec{g_v}</math>). No obstante, para cuestiones que trataremos posteriormente será necesario considerar una tercera coordenada (por ejemplo, para el cálculo del rotacional), por ello, incluiremos en nuestra base natural el vector <math>\vec{g_w}</math>.

*<math> \vec{g_u}=v\hat{e_1} +u \hat{e_2}</math>
*<math> \vec{g_v}=u\hat{e_1} -v \hat{e_2}</math>

El código de Matlab utilizado para hallar los vectores <math>\vec{g_u}, \vec{g_v}</math> es :
{{matlab|codigo=
plot(xx,yy); % Imagen.
hold on % Inicio superposición de gráficos.
mesh(xx,yy,0*xx) % Mallado completo.
quiver(xx,yy,vv,uu); % Representación del primer vector de la base natural en cada punto.
quiver(xx,yy,uu,-vv); % Representación del segundo vector de la base natural en cada punto.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.

}}

La imagen obtenida es:
[[Archivo:Oleole.jpg|1000px|thumb|centre|Representación de los vectores de la base natural <math>\vec{g_u}</math> y <math>\vec{g_v}</math>]]
 

Resulta coherente porque los vectores dibujados son tangentes a las líneas coordenadas (Recordamos que <math>\vec{g_u}</math> y <math>\vec{g_v}</math> se obtienen al derivar el vector de posición <math> \vec{r_0}</math> con respecto a <math>u</math> y <math>v</math>).

==º Acción de la temperatura en la placa==
=== Influencia de un foco de calor ===
La temperatura proviene de un foco de calor dada por el campo escalar
<math>T(x,y)=(8-y^2+2y)e^{-x^2} </math>
{{matlab|codigo=
T=(8-yy.^2+yy.*2).*exp(-(xx.^2)); % Función Temperatura.
subplot(1,2,1); % Muestra varias imágenes. 1ª Imagen.
contour(xx,yy,T,20); % Define 20 líneas de nivel.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
subplot(1,2,2); % Muestra varias imágenes. 2ª Imagen.
surf(xx,yy,T); colorbar; % Visualización de superficie en 3D más leyenda en color.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
max(max(T)) % Valor máximo de la temperatura en toda la región.
}}

Al ejecutar el código se obtienen estas gráficas:
[[Archivo:Temperatura555.jpg|900px|thumb|centre|Gráficas de la variación de la temperatura en la placa.]]
 
El valor de temperatura máximo obtenido en la placa es '''8.7332''' y se da en el punto con el <math>|x|</math> mínimo y la <math>y</math> máxima.

===Gradiente de T y curvas de nivel===
El gradiente de una función escalar expresa la dirección en la cual el campo crece más rápido. Por otra parte, las curvas de nivel expresan los puntos que se encuentran a la misma altura, siendo en nuestro caso aquellos que tienen la misma temperatura. Por tanto, si tenemos dos curvas de nivel y estamos en un punto de la menor, el gradiente será el vector de mínima distancia a la otra curva de nivel, siendo en consecuencia perpendicular a ambas.
<math>\nabla T(x,y)= \frac{\partial T}{\partial x} \hat{e_1}+\frac{\partial T}{\partial y}\hat{e_2}=[-2xe^{-x^2}(8-y^2+2y)]\hat{e_1}+[(-2y+2)e^{-x^2}]\hat{e_2}</math>
{{matlab|codigo=
h=1/20; % Paso de muestreo.
u=1/3:h:1; % Intervalo [1/3,1].
v=-1:h:1; % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v); % Matrices.
xx=uu.*vv; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
f=(8-yy.^2+yy.*2).*exp(-(xx.^2)); % Función Temperatura.
fx =(-2.*xx).*(exp(-(xx.^2))).*(8-yy.^2+2.*yy); % Derivada con respecto a x de la función Temperatura.
fy=(exp(-(xx.^2))).*(-2.*yy+2); % Derivada con respecto a y de la función Temperatura.
hold on % Fin superposición de gráficos.
quiver(xx,yy,fx,fy) % Representación de los vectores gradiente.
contour(xx,yy,f,20);colorbar; % Visualización de superficie en 3D más leyenda en color.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.
}}

La imagen obtenida es la siguiente:
[[Archivo:Gradientetemperatura5.jpg|900px|thumb|centre|Gráfica que representa el gradiente superpuesto a las curvas de nivel.]]
 
Efectivamente en la imagen se aprecia la perpendicularidad de los vectores del gradiente respecto a las curvas de nivel.

==º Acción de una fuerza sobre el sólido==
=== Campo de desplazamientos ===
Una fuerza determinada aplicada sobre nuestro sólido ha provocado un desplazamiento del mismo que viene dado por <math> \vec{u}(x,y) </math> .
Este vector será <math>\vec u(u,v)=\vec a(\vec b\cdot\vec r_{o})</math> siendo ::
<math> \vec{a}= \frac{\vec{g}_u}{|\vec{g_u}|}=\frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}\qquad \vec{b}=-4 \frac{\vec{g}_u}{|\vec{g_u}|}=-4\frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}</math>.

Como hemos hallado anteriormente <math>\vec{g_u}=v \vec{e_1} +u \vec{e_2} </math>. Tomaremos:: <math>\vec{r_o}= x\vec {e_1}+y\vec {e_2}= uv\vec {e_1}+ \frac{1}{2}(u^2-v^2)\vec {e_2}</math>.

Con todo esto: <math>\vec{u}= \frac{\vec{g}_u}{|\vec{g_u}|}(-4 \frac{\vec{g}_u}{|\vec{g_u}|} \vec{r_o})= \frac{\vec{g_u}}{|\vec{g_u}|^2}(-4 \vec{g_u}\cdot\vec{r_o})=\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2} \vec{g_u}=\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v \vec{e_1} +u \vec{e_2})=-2u\vec{g_u}=-2uv \hat{e_1} -2u^2 \hat{e_2}</math>
La representación del campo de desplazamiento <math>\vec{u}</math> y la placa deformada será la siguiente:

{{matlab|codigo=
Ux=-2.*uu.*vv; % Componente x del campo de desplazamientos u.
Uy=-2.*uu.^2; % Componente y del campo de desplazamientos u.
subplot(1,3,1) % Muestra varias imágenes. 1ª Imagen.
quiver(xx,yy,Ux,Uy); % Representación del campo vectorial de desplazamientos u.
axis([-1,1,-1,1]); % Selecciona la región a dibujar.
xd=xx+Ux; % Componente x final del sólido deformado.
yd=yy+Uy; % Componente y final del sólido deformado.
subplot(1,3,2), mesh(xx,yy,0*xx); % Muestra varias imágenes (2ª Imagen) y mallado completo.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
subplot(1,3,3), mesh(xd,yd,0*xd); % Muestra varias imágenes (3ª Imagen) y mallado completo.
axis([-1,1,-2,1]); % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
}}

Las imágenes obtenidas fueron las siguientes:
[[Archivo:Todo0000.jpg|1000px|thumb|centre|Gráficas que muestran, de izquierda a derecha, el campo de vectores de desplazamiento, el mallado original de la placa, y el resultado final de la placa tras la aplicación de la fuerza.]]
 

===Divergencia de un campo===

La divergencia de un campo vectorial <math>\vec{u}</math> se halla por la expresión::

<math>\nabla\cdot\vec{u}= \frac{1}{ \sqrt{g} } \frac{\partial [\sqrt{g} u^{i}] }{\partial u^i}</math>

La divergencia controla la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial (sobre la superficie que rodea a un volumen de control), con lo cual, si la divergencia es positiva, el campo tiene "fuentes" , y si la divergencia es negativa se dice que tiene "sumideros" .

Por su parte, el campo <math>\vec{u}</math> , en este caso, recordamos que estaba definido por:

<math>\vec{u}=-2uv \hat{e_1} -2u^2 \hat{e_2}</math>

En el cálculo de la divergencia necesito las componentes contravariantes <math>u^u,u^v,u^w</math> del campo. Éstas se calculan por la expresión <math>u^i=\vec{u}{g^i}</math>. Recordando que <math>{g^i}</math> son las componentes contravariantes de la base natural <math>\{\vec{g_u},\vec{g_v}\}=\{ v\hat{e_1} +u \hat{e_2} , u\hat{e_1} -v \hat{e_2}\}</math>, éstas se hallan a partir de <math>{g^i}=G^{ij}g_j</math> siendo <math>G^{ij}</math> la matriz inversa de la matriz de Gram de la base natural <math>G_{ij}</math>:

<math>G_{ij}=\begin{bmatrix} \vec{g_u}\cdot\vec{g_u} & \vec{g_u}\cdot\vec{g_v} \\ \vec{g_v}\cdot\vec{g_u} & \vec{g_v}\cdot\vec{g_v}\end{bmatrix}=\begin{bmatrix} u^2+v^2 & 0 \\ 0 & u^2+v^2\end{bmatrix}\qquad G^{ij}=G^{-1}_{ij}</math>

Operando la anterior expresión se calculan las <math>g^{i}</math>:

*<math>g^{u}=\frac{1}{ u^2+v^2 }(v\hat{e_1} +u \hat{e_2}) </math>
*<math>g^{v}=\frac{1}{ u^2+v^2 }(u\hat{e_1} -v \hat{e_2}) </math>

De esta manera, las componentes contravariantes del campo quedan:
*<math>u^u=-4uv^2 -2u(u^2-v^2)</math>
*<math>u^v=0</math>

Y ya se pueden sustituir todos los términos en la expresión de la divergencia. El resultado final es <math>\nabla\cdot\vec u = -2\frac{3u^2+v^2}{u^2+v^2}</math>. Dado que el resultado es negativo se puede concluir que existen sumideros en el flujo que atraviesa el campo vectorial <math>\vec{u}</math>. Ésa es la expresión que se utilizará en el código de Matlab.

{{matlab|codigo=
Divergencia=(-2.*(3*uu.^2+vv.^2)./(uu.^2+vv.^2)); % Campo Divergencia.
surf(xx,yy,Divergencia);colorbar; % Visualización de superficie en 3D más leyenda en color.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
}}

[[Archivo:Jauja.jpg|800px|thumb|centre|Imagen de la divergencia <math>\nabla\cdot\vec{u}</math>.]]
Podemos corroborar con la imagen de la placa deformada, que el mayor cambio de área se produce en la zona superior de la misma, la cual se deforma en sentido descendente; mientras que la menor variación se produce en los picos inferiores, los cuales se trasladan al lado opuesto, "doblándose" la placa sobre el eje de ordenadas.

===Cálculo del rotacional de un campo vectorial===

La expresión del rotacional de un campo vectorial <math>\vec{u}</math> se halla por la siguiente expresión:
<math>\nabla \times\vec{u}= \frac{1}{ \sqrt{g} } \begin{bmatrix} \vec{g_u} & \vec{g_v} & \vec{g_w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ u_{u} &u_{v}&u_{w}\end{bmatrix}</math>

De nuevo necesitaremos definir una tercera componente <math>\vec{g_w}</math> para el cálculo de ese determinante :
*<math> \vec{g_u}=v\hat{e_1} +u \hat{e_2}</math>
*<math> \vec{g_v}=u\hat{e_1} -v \hat{e_2}</math>
*<math> \vec{g_w}=\hat{e_3}</math> (suponiendo una tercera componente en la transformación <math> z=w</math>).

El término <math>g</math> es el determinante de la matriz de Gram <math>G</math>de la base natural:

<math>G=\begin{bmatrix} \vec{g_u}\cdot\vec{g_u} & \vec{g_u}\cdot\vec{g_v} & \vec{g_u}\cdot\vec{g_w} \\ \vec{g_v}\cdot\vec{g_u} & \vec{g_v}\cdot\vec{g_v} & \vec{g_v}\cdot\vec{g_w} \\ \vec{g_w}\cdot\vec{g_u} &\vec{g_w}\cdot\vec{g_v}&\vec{g_w}\cdot\vec{g_w}\end{bmatrix}=\begin{bmatrix} u^2+v^2 & 0 & 0 \\ 0 & u^2+v^2 & 0 \\ 0 &0&1\end{bmatrix}</math>

De modo que <math>g=(u^2+v^2)^2</math>

Por su parte, el campo <math>\vec{u}</math> había quedado definido por la expresión:

<math> \vec {u}(u,v)= \vec{a} (\vec{b}\cdot{r_{o}})=\frac{-4uv^2 -2u(u^2-v^2)}{ u^2+v^2 }(v\hat{e_1} +u \hat{e_2})</math>

Volviendo al Rotacional, necesitamos, por último hallar las componentes covariantes <math>u_{u},u_{v},u_{w}</math> del campo. Éstas se calculan por la expresión:

<math>u_{i}=\vec{u}{g_{i}}</math>

De esta manera:
*<math>u_{u}=-4uv^2 -2u(u^2-v^2)</math>
*<math>u_{v}=0</math>
*<math>u_{w}=0</math>

Y ahora se tienen todos los términos para sustituir en la expresión del rotacional:

<math>\nabla \times\vec{u}= \frac{1}{ \sqrt{g} } \begin{bmatrix} v\hat{e_1} +u \hat{e_2} & u\hat{e_1} -v \hat{e_2} & \hat{e_3} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ -4uv^2 -2u(u^2-v^2) &0&0\end{bmatrix}=\vec{0}</math>

El significado físico de esta situación se traduce en que el campo <math>\vec{u}</math> no tiene tendencia a rotación.

==º Tensiones sobre la placa==
===Tensor de tensiones===

En un medio elástico, lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones <math>\sigma ^ij=λ\nabla·\vec{u}1+2μ\epsilon</math> siendo:

*<math>\epsilon (\vec{u})</math> la parte simétrica del tensor gradiente de <math>\vec{u}</math>, <math>\nabla\vec{u}</math>::

<math>\epsilon (\vec{u})=\frac{\nabla\vec{u}+ \nabla\vec{u}^t}{2}</math>

*<math>λ</math> y <math>μ</math> los coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

Se utilizarán estas expresiones para dibujar las tensiones normales en la dirección que marca <math>\vec{g_u}</math> y la dirección que marca <math>\vec{g_v}</math>, es decir:

<math>\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}\qquad\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}</math>

Bajo estas instrucciones se empieza a definir <math>\sigma ^ij</math>. Para ello nos serviremos de la fórmula de las derivadas parciales covariantes <math>u^{i}, _j = \frac{\partial u^i}{\partial x^j}+\Gamma^i_{jk}u^k</math>

'''1. Tomaremos <math>λ=μ=1</math>.'''

'''2. Se calcula el tensor gradiente de <math>\vec{u}</math>, <math>\nabla\vec{u}</math>:'''

<math>\nabla\vec{u}= \begin{bmatrix} \frac{\partial u^u}{\partial u} + \Gamma ^1_{11} u^u & \frac{\partial u^u}{\partial v} + \Gamma ^1_{12} u^u & \frac{\partial u^u}{\partial w} + \Gamma ^1_{13} u^u \\ \frac{\partial u^v }{\partial u} + \Gamma ^{2}_{11}u^u & \frac{\partial u^v }{\partial v} + \Gamma ^{2}_{12}u^u & \frac{\partial u^v }{\partial w} + \Gamma ^{2}_{13}u^u \\ \frac{\partial u^w }{\partial u} + \Gamma ^{3}_{11}u^u &\frac{\partial u^w }{\partial v} + \Gamma ^{3}_{12}u^u&\frac{\partial u^w }{\partial w} + \Gamma ^{3}_{13}u^u\end{bmatrix}=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{-2uv}{u^2+v^2} & 0 \\\frac{2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix}\qquad\nabla \vec{u}^t=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{2uv}{u^2+v^2} & 0 \\\frac{-2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix}</math>:

'''3. Parte simétrica del tensor gradiente <math>\nabla\vec{u}</math>, <math>\epsilon (\vec{u})</math> :'''

Por lo tanto <math>\epsilon (\vec{u})=\begin{bmatrix}\frac{-4u^2-2v^2}{(u^2+v^2)} & 0 & 0 \\\ 0 & \frac{-2u^2}{(u^2+v^2)} & 0\\0&0&0 \end{bmatrix} </math>

'''4. Cálculo de <math>σ</math>:'''
Se recuerda que la divergencia de <math>\vec{u}</math> es
<math>\nabla\cdot\vec u = -2\frac{3u^2+v^2}{u^2+v^2}</math>

<math>\sigma ^ij=λ\nabla·\vec{u}1+2μЄ=\begin{bmatrix}\frac{-14u^2-6v^2}{u^2+v^2} & 0 & 0 \\\ 0 & \frac{-10u^2-2v^2}{u^2+v^2} & 0\\0 &0& \frac{-6u^2-2v^2}{u^2+v^2} \end{bmatrix}</math>

Por último, para la representación de las tensiones normales en la dirección de <math>\vec{g_u}</math> y la dirección de <math>\vec{g_v}</math>, se necesita definir:<math>\frac{\vec{g}_u}{|\vec{g_u}|}= \frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}\qquad\frac{\vec{g}_v}{|\vec{g_v}|}=\frac{u\vec {e_1}-v\vec {e_2}}{\sqrt{v^2+u^2}}</math>

Sabiendo esto, se puede proceder en Matlab al cálculo de las matrices que permitirán posteriormente la representación de las tensiones normales en la dirección que marca <math>\vec{g_u}</math> y la dirección de <math>\vec{g_v}</math>:

*<math>\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}</math>
*<math>\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}</math>

===Tensión de Von Mises===
La tensión de Von Mises es un escalar que indica el comportamiento plástico del sólido.:
<math> \sigma _{vm}= \sqrt{ \frac{( \sigma _1- \sigma _2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2 }{2} }</math>
<math> \sigma_1,\sigma_2 </math> y <math> \sigma_3 </math> son los autovalores del tensor <math>σ</math> (tensiones principales).
{{matlab|codigo=
sigma1=(-14*uu.^2-6*vv.^2)./(vv.^2+uu.^2); % Autovalor 1.
sigma2=(-10*uu.^2-2*vv.^2)./(vv.^2+uu.^2); % Autovalor 2.
sigma3=(-6*uu.^2-2*vv.^2)./(vv.^2+uu.^2); % Autovalor 3.
Mises=sqrt(((sigma1-sigma2).^2+(sigma2-sigma3).^2+(sigma3-sigma1).^2)./2); % Fórmula de Von Misses.
surf(xx,yy,Mises); % Visualización de superficie en 3D.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
max(max(Mises)) % Valor máximo de la tensión de Von Mises en toda la región.
}}

Este fue el resultado:
[[Archivo:Elele.jpg|800px|thumb|centre|Visualización en 3D de la placa tras la tensión de Von Mises.]]
 

El valor máximo de la tensión de Von Mises es '''6.9282''' y se alcanza en los '''puntos de la placa pertenecientes al eje de ordenadas'''.

==º Masa de la placa==
En último lugar procedemos al cálculo de la masa de la placa. Para ello partimos de una función de densidad <math>xye^{-1/x^2}</math>.

Para hallar la masa de la placa, utilizaremos un proceso numérico basado en el '''Método del Trapecio para la aproximación de integrales'''. El cálculo se basa en aplicar la función en cada punto y obtener una matriz de valores de la densidad en cada punto evaluado de la malla. Una vez hallada la matriz, multiplicándola por un vector fila y columna y sumando todos los elementos de la matriz obtendremos la masa total de la placa.

Como <math>x</math> e <math>y</math> pueden ser valores negativos, puede resultar que la función densidad sea negativa en algunos puntos de la placa. Para evitar este absurdo, los resultados finales se obtendrán a partir de convertir todos los valores de la matriz de densidades a su valor absoluto, multiplicándolos posteriormente por los pasos y sumándolos todos entre sí.

{{matlab|codigo=
N1=200; N2=200; % Número de puntos.
a=1/3; b=1; c=-1; d=1; % Extremos de los intervalos.
h1=(b-a)/N1; h2=(d-c)/N2; % Pasos.
u=a:h1:b; v=c:h2:d; % Intervalos.
[uu,vv]=meshgrid(u,v); % Malla de datos.
xx=uu.*vv; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
d=(xx.*yy).*(exp(-1./(xx.^2))); % Función Densidad.
D=abs(d); % Función Densidad en valor absoluto.
w1=ones(N1+1,1); % Matriz de N1+1 elementos unidad.
w1(1)=1/2; w1(N1+1)=1/2; % Primer y último de la matriz w1.
w2=ones(N2+1,1); % Matriz de N2+1 elementos unidad.
w2(1)=1/2; w2(N1+1)=1/2; % Primer y último de la matriz w2.
result=h1*h2*w2'*D*w1 % Masa obtenida con la función densidad D.
}}

El valor final de la masa obtenido es '''0.0026'''

Campos en Elasticidad

2014-12-05T13:55:27Z

Araceli Martin: /* Tensor de tensiones */

{{ TrabajoED |Campos escalares y vectoriales en elasticidad. Grupo 16-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2014-15]] | Araceli Martín, Juan Carlos Durán, Francisco Javier Alcaraz, Álvaro Llera, Clara Callejo, Manuel Escudero }}
[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Para este análisis y representación de campos escalares en elasticidad nos situaremos en el contexto de una '''placa plana''' que ocupa la región comprendida entre las parábolas :

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

Para representar la placa se utilizará un sistema de coordenadas curvilíneas tal que:
* <math>x=uv</math>
* <math> y = \dfrac{1}{2}(u^2-v^2)</math>

En nuestro análisis, el dominio de '''<math>u</math>''' y '''<math>v</math>''' comprenderá:
<math>(u,v) \in [1/3,1]x[-1,1]</math>

==º Situación inicial de la placa==

Primero comenzaremos el estudio de la placa en cuestión planteando su representación gráfica mediante un mallado, utilizando, en el código, la conversión a coordenadas curvilíneas de '''<math>u</math>''' y '''<math>v</math>'''. El intervalo en el que representaremos comprende::

<math>(x,y) \in [-1,1]*[-1,1]</math>

Para discretizar los vectores <math>u</math> y <math>v</math> utilizaremos un paso <math> h = \dfrac{1}{20}</math>.

Código en Matlab:
{{matlab|codigo=
h=1/20; % Paso de muestreo.
u=1/3:h:1; % Intervalo [1/3,1].
v=-1:h:1; % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v); % Matrices de datos.
xx=uu.*vv ; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
plot(xx,yy); % Imagen.
mesh(xx,yy,0*xx) % Mallado.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
}}

Los resultados de la representación se ven en la imagen:

[[Archivo:Imagen_mallado.jpg|800px|thumb|centre|Mallado que representa la superficie de la placa comprendida entre las parábolas <math>P1</math> y <math>P2</math>.]]
 

Por otro lado, cuando se hace uso de cambio de coordenadas a ciertas coordenadas curvilíneas, es útil para el entendimiento de la transformación la representación de las '''Líneas Coordenadas'''. Las líneas coordenadas son aquellas que se obtienen variando una de las coordenadas de la transformación (<math>u</math> o <math>v</math>) y manteniendo fija la restante.
En este estudio hemos representado varias líneas coordenadas a base de dar un valor concreto a <math>u</math> o a <math>v</math>, dentro de sus respectivos intervalos:

{{matlab|codigo=
xx11=uu.*0.5 ; % ----------------------------------------------Parametrización X.
yy11=(1/2).*((uu.^2)-(0.5.^2)); % ----------------------------------------------Parametrización Y.
xx12=uu.*-0.5 ; % ----------------------------------------------Parametrización X.
yy12=(1/2).*((uu.^2)-((-0.5).^2)); % ----------------------------------------------Parametrización Y.
xx13=uu.*1 ; % ----------------------------------------------Parametrización X.
yy13=(1/2).*((uu.^2)-(1.^2)); % ----------------------------------------------Parametrización Y.
xx14=uu.*-1 ; % ----------------------------------------------Parametrización X.
yy14=(1/2).*((uu.^2)-((-1).^2)); % ----------------------------------------------Parametrización Y.
xx15=uu.*0.75 ; % ----------------------------------------------Parametrización X.
yy15=(1/2).*((uu.^2)-(0.75.^2)); % ----------------------------------------------Parametrización Y.
xx16=uu.*-0.75 ; % ----------------------------------------------Parametrización X.
yy16=(1/2).*((uu.^2)-((-0.75).^2)); % ----------------------------------------------Parametrización Y.
xx17=uu.*0 ; % ----------------------------------------------Parametrización X.
yy17=(1/2).*((uu.^2)-(0.^2)); % ----------------------------------------------Parametrización Y.}}

A continuación se representa las gráficas resultantes:
[[Archivo:lineascoordenadas.jpg|800px|thumb|centre|Líneas coordenadas fijando la variable <math>v</math> y <math>u</math>, respectivamente.]]
 

=== Expresión en la Base Natural===

Cuando se realiza una transformación a coordenadas curvilíneas (de <math>x</math> e <math>y</math> a <math>u</math> y <math>v</math>), el vector de posición <math> \vec{r_0}</math> se expresa como:

<math> \vec{r_0}=x\hat{e_1} +y \hat{e_2} =uv\hat{e_1} +\dfrac{1}{2}(u^2-v^2) \hat{e_2} </math>

La base natural <math>\vec{g_u}, \vec{g_v}</math> es la que tiene por vectores el resultado de derivar el vector posición <math> \vec{r_0}</math> según las nuevas coordenadas <math>u</math> y <math>v</math>. Dado que es un problema sobre una placa plana, estamos en una situación de dos dimensiones en la que para cualquier base, sólo se requieren dos vectores, (<math>\vec{g_u}, \vec{g_v}</math>). No obstante, para cuestiones que trataremos posteriormente será necesario considerar una tercera coordenada (por ejemplo, para el cálculo del rotacional), por ello, incluiremos en nuestra base natural el vector <math>\vec{g_w}</math>.

*<math> \vec{g_u}=v\hat{e_1} +u \hat{e_2}</math>
*<math> \vec{g_v}=u\hat{e_1} -v \hat{e_2}</math>

El código de Matlab utilizado para hallar los vectores <math>\vec{g_u}, \vec{g_v}</math> es :
{{matlab|codigo=
subplot(3,3,4); % Muestra varias imágenes. 1º Imagen.
hold on % Inicio superposición de gráficos.
mesh(xx,yy,0*xx) % Mallado completo.
quiver(xx,yy,vv,uu); % Representación del primer vector de la base natural en cada punto.
quiver(xx,yy,uu,-vv); % Representación del segundo vector de la base natural en cada punto.
axis([-1,1,-1,1]) % Selecciona la regíon a dibujar.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.

}}

La imagen que se obtiene con este código es:
[[Archivo:basenatural.jpg|800px|thumb|centre|Representación de los vectores de la base natural <math>\vec{g_u}</math> y <math>\vec{g_v}</math>]]
 

Resulta coherente porque los vectores dibujados son tangentes a las líneas coordenadas (Recordamos que <math>\vec{g_u}</math> y <math>\vec{g_v}</math> se obtienen al derivar el vector de posición <math> \vec{r_0}</math> con respecto a <math>u</math> y <math>v</math>).

==º Acción de la temperatura en la placa==
=== Influencia de un foco de calor ===
La temperatura proviene de un foco de calor dada por el campo escalar
<math>T(x,y)=(8-y^2+2y)e^{-x^2} </math>
{{matlab|codigo=
f=(8-yy.^2+yy.*2).*exp(-(xx.^2)); % Función Temperatura.
subplot(1,2,1); % Muestra varias imágenes. 1ª Imagen.
contour(xx,yy,f,20); % Define 20 líneas de nivel.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
subplot(1,2,2); % Muestra varias imágenes. 2ª Imagen.
surf(xx,yy,f); colorbar; % Visualización de superficie en 3D más leyenda en color.
axis([-1,1,-1,1]) % Selecciona la región a dibujar.
view(2) % Ver imagen desde arriba.
max(max(f)) % Valor máximo de la temperatura en toda la región.
}}

Al ejecutar el código se obtienen estas gráficas:
[[Archivo:graficatemperatura.jpg|800px|thumb|centre|Gráficas de la temperatura en la placa.]]
 

===Gradiente de T y curvas de nivel===
El gradiente de una función escalar expresa la dirección en la cual el campo crece más rápido. Por otra parte,las curvas de nivel expresan los puntos que se encuentran a la misma altura, en nuestro caso, aquellos que tienen la misma temperatura. Por tanto, si tenemos dos curvas de nivel y estamos en un punto de la menor, el gradiente será el vector de mínima distancia a la otra curva de nivel, siendo por tanto perpendicular a ambas.
<math>\nabla T(x,y)= \frac{\partial T}{\partial x} \hat{e_1}+\frac{\partial T}{\partial y}\hat{e_2}=[-2xe^{-x^2}(8-y^2+2y)]\hat{e_1}+[(-2y+2)e^{-x^2}]\hat{e_2}</math>
{{matlab|codigo=
h=1/20; % Paso de muestreo.
u=1/3:h:1; % Intervalo [1/3,1].
v=-1:h:1; % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v); % Matrices.
xx=uu.*vv ; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
f=(8-yy.^2+yy.*2).*exp(-(xx.^2)); % Función Temperatura.
fx =(-2.*xx).*(exp(-(xx.^2))).*(8-yy.^2+2.*yy); % Derivada con respecto a x de la función Temperatura.
fy=(exp(-(xx.^2))).*(-2.*yy+2); % Derivada con respecto a y de la función Temperatura.
hold on % Fin superposición de gráficos.
quiver(xx,yy,fx,fy) % Representación de los vectores gradiente.
contour(xx,yy,f,20);colorbar; % Visualización de superficie en 3D más leyenda en color.
view(2) % Ver imagen desde arriba.
hold off % Fin superposición de gráficos.
}}

==Acción de una fuerza sobre sólido==
=== Campo de desplazamientos ===
Una fuerza determinada aplicada sobre nuestro sólido ha provocado un desplazamiento del mismo que viene dado por <math> \vec{u}(x,y) </math> .
Este vector será <math>\vec u(u,v)=\vec a(\vec b\cdot\vec r_{o})</math> siendo ::
<math> \vec{a}= \frac{\vec{g}_u}{|\vec{g_u}|}=\frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}\qquad \vec{b}=-4 \frac{\vec{g}_u}{|\vec{g_u}|}=-4\frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}</math>.

Como hemos hallado anteriormente <math>\vec{g_u}=v \vec{e_1} +u \vec{e_2} </math>. Tomaremos:: <math>\vec{r_o}= x\vec {e_1}+y\vec {e_2}= uv\vec {e_1}+ \frac{1}{2}(u^2-v^2)\vec {e_2}</math>.

Con todo esto: <math>\vec{u}= \frac{\vec{g}_u}{|\vec{g_u}|}(-4 \frac{\vec{g}_u}{|\vec{g_u}|} \vec{r_o})= \frac{\vec{g_u}}{|\vec{g_u}|^2}(-4 \vec{g_u}\cdot\vec{r_o})=\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2} \vec{g_u}=\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v \vec{e_1} +u \vec{e_2})</math>
La representación del campo de desplazamiento <math>\vec{u}</math> será la siguiente:

Nuestra superficie sufre un desplazamiento en un instante <math>t_{0}</math> debido a una percusión. Dicho campo viene dado por el siguiente vector <math>\vec u</math>: <math>\vec u(u,v)=\vec a(\vec b\cdot\vec r_{o})</math>
===Divergencia de un campo===

La divergencia de un campo vectorial <math>\vec{u}</math> se halla por la expresión:

<math>\nabla\cdot\vec u = \frac{1}{\sqrt{g}} \frac{\part}{\part x^i}
\left(\sqrt{g} v^i \right)</math>

La divergencia controla la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial (sobre la superficie que rodea a un volumen de control), con lo cual, si la divergencia es positiva, el campo tiene "fuentes" , y si la divergencia es negativa se dice que tiene "sumideros" .

Por su parte, el campo <math>\vec{u}</math> , en este caso está definido por la expresión:

<math> \vec {u}(u,v)= \vec{a} (\vec{b} \vec{r_o})</math>

Donde

*<math> \vec{a}= \frac{\vec{g}_u}{|\vec{g_u}|}</math>
*<math> \vec{b}=-4 \frac{\vec{g}_u}{|\vec{g_u}|}</math>.
*<math> \vec{r_0}=x\hat{e_1} +y \hat{e_2} =uv\hat{e_1} +\dfrac{1}{2}(u^2-v^2) \hat{e_2} </math>

Dado que <math>{|\vec{g_u}|}=\sqrt{u^2+v^2}</math>,

*<math> \vec{a}= \frac{v\hat{e_1} +u \hat{e_2}}{\sqrt{u^2+v^2}}</math>
*<math> \vec{b}= -4\frac{v\hat{e_1} +u \hat{e_2}}{\sqrt{u^2+v^2}}</math>

Operando, la expresión del campo de desplazamientos <math>\vec{u}</math> queda:

<math>\vec{u}=\frac{-4uv^2 -2u(u^2-v^2)}{ u^2+v^2 }(v\hat{e_1} +u \hat{e_2})</math>

En el cálculo de la divergencia necesito las componentes contravariantes <math>u^u,u^v,u^w</math> del campo. Éstas se calculan por la expresión:

<math>u^i=\mathbf{\vec{u}}\cdot\mathbf{g^i}</math>

Recordando que <math>{g^i}</math> son las componentes contravariantes de la base natural <math>\{g_u , g_v\}=\{ v\hat{e_1} +u \hat{e_2} ; u\hat{e_1} -v \hat{e_2}\}</math>, éstas se hallan a partir de:

<math>{g^i}=G^{ij}g_j</math>

siendo <math>G^{ij}</math> la matriz inversa de la matriz de Gram de la base natural <math>G_{ij}</math>:

<math>G_{ij}=\begin{bmatrix} g_{u}g_{u} & g_{v}g_{u} \\ g_{v}g_{u} & g_{v}g_{v}\end{bmatrix}=\begin{bmatrix} u^2+v^2 & 0 \\ 0 & u^2+v^2\end{bmatrix}</math>

<math>G^{ij}=G^{-1}_{ij}</math>

Operando la anterior expresión se calculan las <math>g^{i}</math>:

*<math>g^{u}=\frac{1}{ u^2+v^2 }(v\hat{e_1} +u \hat{e_2}) </math>
*<math>g^{u}=\frac{1}{ u^2+v^2 }(u\hat{e_1} -v \hat{e_2}) </math>

De esta manera, las componentes contravariantes del campo quedan:
*<math>u^u=-4uv^2 -2u(u^2-v^2)</math>
*<math>u^v=0</math>

Y ya se pueden sustituir todos los términos en la expresión de la divergencia. El resultado final es:

<math>\nabla\cdot\vec u = -2\frac{3u^2+v^2}{u^2+v^2}</math>

Dado que el resultado es negativo se puede concluir que existen sumideros en el flujo que atraviesa el campo vectorial <math>\vec{u}</math>. Ésa es la expresión que se utilizará en el código de Matlab.

===Cálculo del rotacional de un campo vectorial===

La expresión del rotacional de un campo vectorial <math>\vec{u}</math> se halla por la siguiente expresión:
<math>\nabla \times\vec{u}= \frac{1}{ \sqrt{g} } \begin{bmatrix} g_{u} & g_{v} & g_{w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ u_{u} &u_{v}&u_{w}\end{bmatrix}</math>

De nuevo necesitaremos definir una tercera componente <math>\vec{g_w}</math> para el cálculo de ese determinante :
*<math> \vec{g_u}=v\hat{e_1} +u \hat{e_2}</math>
*<math> \vec{g_v}=u\hat{e_1} -v \hat{e_2}</math>
*<math> \vec{g_w}=\hat{e_3}</math> (suponiendo una tercera componente en la transformación: <math> z=w</math> ).

El término <math>g</math> es el determinante de la matriz de Gram <math>G</math>de la base natural:

<math>G=\begin{bmatrix} g_{u}g_{u} & g_{v}g_{u} & g_{w}g_{u} \\ g_{v}g_{u} & g_{v}g_{v} & g_{v}g_{w} \\ g_{w}g_{u} &g_{w}g_{v}&g_{w}g_{w}\end{bmatrix}=\begin{bmatrix} u^2+v^2 & 0 & 0 \\ 0 & u^2+v^2 & 0 \\ 0 &0&1\end{bmatrix}</math>

De modo que <math>g=(u^2+v^2)^2</math>

Por su parte, el campo <math>\vec{u}</math> había quedado definido por la expresión:

<math> \vec {u}(u,v)= \vec{a} (\mathbf{\vec{b}}\cdot\mathbf{r_{0}})=\frac{-4uv^2 -2u(u^2-v^2)}{ u^2+v^2 }(v\hat{e_1} +u \hat{e_2})</math>

Volviendo al Rotacional, necesitamos, por último hallar las componentes covariantes <math>u_{u},u_{v},u_{w}</math> del campo. Éstas se calculan por la expresión:

<math>u_{i}=\mathbf{\vec{u}}\cdot\mathbf{g_{i}}</math>

De esta manera:
*<math>u_{u}=-4uv^2 -2u(u^2-v^2)</math>
*<math>u_{v}=0</math>
*<math>u_{w}=0</math>

Y ahora se tienen todos los términos para sustituir en la expresión del rotacional:

<math>\nabla \times\vec{u}= \frac{1}{ \sqrt{g} } \begin{bmatrix} v\hat{e_1} +u \hat{e_2} & u\hat{e_1} -v \hat{e_2} & \hat{e_3} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ -4uv^2 -2u(u^2-v^2) &0&0\end{bmatrix}=0</math>

Dado que el el valor de <math>\nabla \times\vec{u}=0</math> no hay nada que representar ni hallar numéricamente. Él significado físico de esta situación se traduce en que el campo vectorial <math>\vec{u}</math> no tiene tendencia a rotar sobre ningún punto.

==Tensiones sobre la placa==
===Tensor de tensiones===

En un medio elástico, lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones <math>σij</math> a través de la fórmula:

<math>σ=λ\nabla·\vec{u}1+2μЄ</math>

siendo:

*<math>\epsilon (\vec{u})</math> la parte simétrica del tensor gradiente de <math>\vec{u}</math>, <math>\nabla\vec{u}</math>::

<math>\epsilon (\vec{u})=\frac{\nabla\vec{u}+ \nabla\vec{u}^t}{2}</math>

*<math>λ</math> y <math>μ</math> los coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

Se utilizarán estas expresiones para dibujar las tensiones normales en la dirección que marca <math>\vec{g_u}</math> y la dirección que marca <math>\vec{g_v}</math>,es decir:

*<math>\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}</math>
*<math>\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}</math>

Bajo estas instrucciones se empieza a definir <math>σij</math>:

'''1.Tomaremos <math>λ=μ=1</math>.'''

'''2.Se calcula el tensor gradiente de <math>\vec{u}</math>, <math>\nabla\vec{u}</math>:'''

<math>\nabla\vec{u}= \begin{bmatrix} \frac{\partial u_{u}}{\partial u} & \frac{\partial u_{u}}{\partial v} & \frac{\partial u_{u}}{\partial w} \\ \frac{\partial u_{v}}{\partial u} & \frac{\partial u_{v}}{\partial v} & \frac{\partial u_{v}}{\partial w} \\ \frac{\partial u_{w}}{\partial u} &\frac{\partial u_{w}}{\partial v}&\frac{\partial u_{w}}{\partial w}\end{bmatrix}</math>
Reutilizando las componentes covariantes de <math>\vec{u}</math>:

*<math>u_{u}=-4uv^2 -2u(u^2-v^2)</math>
*<math>u_{v}=0</math>
*<math>u_{w}=0</math>

<math>\nabla\vec{u}= \begin{bmatrix} -6u^2-2v^2 & -4uv & 0 \\ 0 & 0 & 0 \\ 0 &0&0\end{bmatrix}</math>

'''3.Parte simétrica del tensor gradiente <math>\nabla\vec{u}</math>, <math>\epsilon (\vec{u})</math> :'''

<math>\epsilon (\vec{u})=\frac{\nabla\vec{u}+ \nabla\vec{u}^t}{2}=\begin{bmatrix} -6u^2-2v^2 & -2uv & 0 \\ -2uv & 0 & 0 \\ 0 &0&0\end{bmatrix} </math>

'''4. Cálculo de <math>σ</math>:'''
Se recuerda que la divergencia de <math>\vec{u}</math> es
<math>\nabla\cdot\vec u = -2\frac{3u^2+v^2}{u^2+v^2}</math>

<math>σ=λ\nabla·\vec{u}1+2μЄ=\begin{bmatrix} \frac{-6u^2-2v^2}{u^2+v^2}-12u^2-4v^2 & -4uv & 0 \\ -4uv & \frac{-6u^2-2v^2}{u^2+v^2} & 0 \\ 0 &0&\frac{-6u^2-2v^2}{u^2+v^2}\end{bmatrix}</math>

Por último,para la respresentación de las tensiones normales en la dirección de <math>\vec{g_u}</math> y la dirección de <math>\vec{g_v}</math>,se necesita definir

*<math>\frac{\vec{g}_u}{|\vec{g_u}|}= \frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}</math>
*<math>\frac{\vec{g}_v}{|\vec{g_v}|}=\frac{u\vec {e_1}-v\vec {e_2}}{\sqrt{v^2+u^2}}</math>

Sabiendo esto,se puede prodecer en el código de Matlab al cálculo de las matrices que permitirán posteriormente la respresentación de las tensiones normales en la dirección que marca <math>\vec{g_u}</math> y la dirección de <math>\vec{g_v}</math>:

*<math>\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}</math>
*<math>\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}</math>

===Tensión de Von Mises===
==Masa de la placa==
{{matlab|codigo=
N1=200; N2=200; % Número de puntos.
a=1/3; b=1; c=-1; d=1; % Extremos de los intervalos.
h1=(b-a)/N1; h2=(d-c)/N2; % Pasos.
u=a:h1:b; v=c:h2:d; % Intervalos.
[uu,vv]=meshgrid(u,v); % Matrices.
xx=uu.*vv; % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2)); % Parametrización Y.
d=(xx.*yy).*(exp(-1./(xx.^2))); % Función Densidad.
D=abs(d); % Función Densidad en valor absoluto.
w1=ones(N1+1,1); % Matriz de N1+1 elementos unidad.
w(1)=1/2; w(N1+1)=1/2; % Primer y último de la matriz w.
w2=ones(N2+1,1); % Matriz de N2+1 elementos unidad.
w(1)=1/2; w(N2+1)=1/2; % Primer y último de la matriz w.
result=h1*h2*w2'*d*w1 % Masa obtenida con la función densidad d.
result=h1*h2*w2'*D*w1 % Masa obtenida con la función densidad D.
}}

Campos en Elasticidad

2014-12-04T23:50:25Z

Araceli Martin: /* Cálculo del rotacional de un campo vectorial */

Campos en Elasticidad

2014-12-04T23:49:44Z

Araceli Martin: /* Divergencia de un campo */

Campos en Elasticidad

2014-12-04T23:40:20Z

Araceli Martin: /* Cálculo del rotacional de un campo vectorial */

Campos en Elasticidad

2014-12-04T23:37:05Z

Araceli Martin: /* Divergencia de un campo */

Campos en Elasticidad

2014-12-04T22:52:15Z

Araceli Martin: /* Cálculo del rotacional de un campo vectorial */

Campos en Elasticidad

2014-12-03T18:40:47Z

Araceli Martin:

Campos en Elasticidad

2014-12-03T16:51:27Z

Araceli Martin: /* Campo de desplazamientos */

Campos en Elasticidad

2014-12-03T16:10:29Z

Araceli Martin: /* Base Natural */

Campos en Elasticidad

2014-12-03T16:09:52Z

Araceli Martin: /* Líneas coordenadas */

Campos en Elasticidad

2014-12-03T16:09:32Z

Araceli Martin: /* Representación de la placa */

Campos en Elasticidad

2014-12-03T16:06:22Z

Araceli Martin:

Campos en Elasticidad

2014-12-03T16:02:28Z

Araceli Martin: /* Base Natural */

Campos en Elasticidad

2014-12-03T16:01:09Z

Araceli Martin: /* Base Natural */

Campos en Elasticidad

2014-12-03T15:54:05Z

Araceli Martin: /* Base Natural */

Campos en Elasticidad

2014-12-03T15:52:22Z

Araceli Martin: /* Líneas coordenadas y vectores de la base natural */

Campos en Elasticidad

2014-12-03T15:51:55Z

Araceli Martin: /* Base Natural */

Campos en Elasticidad

2014-12-03T15:49:49Z

Araceli Martin: /* º Líneas coordenadas */

Campos en Elasticidad

2014-12-03T15:47:21Z

Araceli Martin: /* Representación del mallado del sólido */