Ecuaciones Diferenciales Ordinarias

Ecuaciones Diferenciales Ordinarias

Una Ecuación Diferencial Ordinaria expresa la razón de cambio de la variable dependiente y con respecto al cambio de una variable independiente x.

La variable dependiente y puede ser la población de una bacteria, la temperatura de un cuerpo o la concentración de un producto en un reactor. La variable independiente x es aquella contra la que se mide la variable dependiente, por ejemplo, para los casos anteriores la variable independiente puede ser el tiempo. Entonces la ecuación diferencial expresa cómo cambia la población de una bacteria en el tiempo, cómo cambia la temperatura de un cuerpo en el tiempo, cómo cambia la concentración de un producto en el tiempo, etc.

Las ecuaciones diferenciales se resuelven eliminando las derivadas que contienen y obteniendo las funciones que cumplen con la relación expresada en la ecuación diferencial.

Por ejemplo, en la cinética de las reacciones, el problema habitual es investigar las reacciones químicas y la velocidad con la que un compuesto A (reactivo) se transforma en otro compuesto B (producto).

En los casos más sencillos ésta velocidad se puede expresar de la forma

Aquí se expresa el cambio de la concentración de A, la cual involucra una constante de proporcionalidad k, como la constante es negativa entonces la concentración de A disminuye.

Una ecuación diferencial tiene una infinidad de soluciones. Esto se entiende porque la ecuación diferencial solo expresa razón de cambio, si se conoce el valor inicial entonces se puede establecer una de las múltiples soluciones. En la gráfica siguiente se muestran las posibles soluciones de la ecuación diferencial y una solución particular para un valor inicial dado.

Si requerimos una solución en particular, se debe especificar la concentración inicial de A, condición inicial de la ecuación diferencial.

Los métodos aplicados a la solución de ecuaciones diferenciales en Ingeniería Química son muchos y dependen principalmente del tipo de ecuación diferencial generado en cada problema. De modo general podemos distinguir dos grandes grupos:

• Métodos Analíticos. Obtienen una solución exacta. Los métodos analíticos se ajustan a un solo tipo de ecuación diferencial, por lo que no hay un método general y en algunos casos no existen métodos analíticos aplicables, implican un conocimiento de cada método para aplicarlos.
• Métodos Numéricos. Obtienen una solución aproximada. Son métodos más genéricos y solo dependen de las condiciones iniciales, no es una solución continua y se requiere de hacer cálculos repetitivos

Problemas de Valor Inicial

Como hemos visto, las ecuaciones diferenciales expresan la razón de cambio de una variable dependiente con respecto a cada cambio de la variable independiente, la solución de la ecuación diferencial debe cumplir una condición inicial dada, ya que de lo contrario solo se tiene la solución general. Por ejemplo, para la ecuación diferencial

tiene como solución general

donde C es una constante arbitraria que puede ser cualquier valor real y cada uno es una solución de la ecuación diferencial. Como se requiere un valor inicial para obtener la solución particular entonces el problema se convierte en un problema de valor inicial. Se requieren de tres elementos para resolver un problema de valor inicial.

• Una ecuación diferencial de primer orden
• Un punto inicial conocido, y(x0)=y0
• El valor final al que se quiere llegar.

Euler

Para resolver una ecuación diferencial con un método numérico solo se cuenta con la ecuación diferencial f(x,y) que gráficamente es la pendiente de una recta tangente a la función en el punto inicial x0 y sabemos el valor de y0 para identificar la solución única. Con estos elementos lo que podemos construir es una linea recta con la pendiente, que nos da la ecuación diferencial y evaluar cuánto vale en el punto final ésto nos dará una aproximación al valor buscado.

De tal manera que el punto calculado en y(x0+h) se obtiene con:

El cual nos da una aproximación al valor de y(x0+h) donde h es el incremento de la variable independiente y donde se quiere conocer el valor de la variable dependiente. Es decir, si conocemos la concentración inicial en el tiempo 0 y deseamos calcular la concentración en el tiempo 0.6min, entonces h vale 0.6.

La diferencia entre la aproximación de Euler con el valor correcto se da porque se toma la derivada en el punto inicial x0=0 y se evalúa en x1=0.6, como se observa en la gráfica, el valor de la pendiente cambia en el intervalo x0 y x0+h, si tomamos cada cambio de la derivada se puede hacer una mejor aproximación, es decir, si en lugar de calcular y(0.6) en un solo paso con h=0.6, lo hacemos en 5 pasos entonces h=0.6/5=0.12 y tendremos que hacer 5 cálculos: y(0.12), y(0.24), y(0.36), y(0.48), y(0.60) y el resultado será mejor.

Se observa en las siguientes gráficas el método de Euler con 1, 5, 20 y 50 pasos, donde cada paso tiene una separación h cada vez más pequeña.

Definimos h como el tamaño de paso, y n como el número de pasos.

Cada paso será de tamaño h, de tal manera que

• x(1)=x(0)+h
• x(2)=x(1)+h,…,
• x(n)=x(n – 1)+h

y los valores de cada y(i+1) se obtienen con la siguiente ecuación:

Ejemplo. Calcular la concentración final de A, si la constante de velocidad de reacción k=2 y la concentración inicial es 1.5 en el tiempo 0, calcular la concentración de A en el tiempo t=0.6

function [xi,yi]=euler(f,x0 ,y0 ,x1 ,n)
     h=(x1 -x0)/n; % tamano de paso
     xi=zeros(n+1,1); % vector de x variable independiente
     yi=zeros(n+1,1); % vector de y variable dependiente
     xi(1)=x0; % tiempo iniicial
     yi(1)=y0; % concentracion inicial
     for i=1:n
         yi(i+1)= yi(i)+h*f(xi(i),yi(i));
         xi(i+1)= xi(i)+h;
     end
 end

f=@(x,y) -2*y; %Ecuación diferencial
x0=0; % tiempo inicial
y0=1.5; %valor inicial
xn=0.6; %valor final del tiempo
n=5; %numero de pasos
[tiempo,concentracion]=euler(f,x0,y0,xn,n); %llamada al metodo de Euler
concentracion(end) %valor final de la concentración
plot(tiempo,concentracion,'--o') %grafica de la solución
title('Reacción química')

Euler modificado

Como observamos en el método de Euler de un paso, el cálculo de y(i+1) se obtiene con la ecuación de la linea recta con pendiente f(xi,yi), la cual es una aproximación con un cierto error. Si calculamos la pendiente en el punto y(i+1) evaluando la ecuación diferencial como f(x(i+1), y(i+1)) y promediamos ambas pendientes, obtendremos un mejor resultado.

Observamos que aplicando el método de Euler simple con la pendiente m1 obtenemos el resultado de la gráfica m1, ahora evaluamos la ecuación diferencial en el punto final y obtenemos la gráfica m2, tomemos ahora el promedio de ambas pendientes, es decir, m3 y volvamos a calcular el resultado, observamos un mejor resultado que con m1 (método de Euler simple).

El método de Euler modificado implica calcular dos veces y(i+1), primero evaluando la ecuación diferencial en el punto anterior f(xi,yi), es decir, con la pendiente m1 y luego con el promedio de las dos pendientes m3. Al primero le llamamos y*(i+1) y al segundo y(i+1).

El método de Euler modificado también se puede mejorar dividiendo el intervalo en n pasos, entonces se deben hacer n veces con un tamaño de paso h cada vez.

Observe que se tiene un mejor resultado con 5 pasos, casi tan bueno como el método de Euler simple con 50 pasos.

function [xi,yi]=edo_euler_modificado(f,x0 ,y0 ,x1 ,n)
     h=(x1 -x0)/n; % tamano de paso
     xi=zeros(n+1,1); % vector de x variable independiente
     yi=zeros(n+1,1); % vector de y variable dependiente
     xi(1)=x0; % tiempo iniicial
     yi(1)=y0; % concentracion inicial

     for i=1:n
         xi(i+1)= xi(i)+h;
         yi(i+1)= yi(i)+h*f(xi(i),yi(i));
         yi(i+1)= yi(i)+h*(f(xi(i),yi(i))+f(xi(i+1), yi(i+1)) )/2;
     end
 end

Problemas de mezclas

Se requiere calcular la cantidad de una sustancia, C(t), que hay en un tanque en cada instante de tiempo t. Sabiendo que la derivada de C respecto a t expresa la razón de cambio de la sustancia presente en el tanque, se cumple la relación

Dada la velocidad a la que el fluido que contiene la sustancia entra en el tanque y la concentración de la sustancia, se cumple la relación:

Velocidad de entrada = Velocidad de flujo entrante x concentración

Suponiendo que la concentración de la sustancia es uniforme, para calcular la concentración se divide C(t) por el volumen de la mezcla que hay en el instante t. Así,

Velocidad de salida = Velocidad de flujo saliente x concentración.

Por ejemplo, en un tanque de 1000 litros de una solución que tiene 50kg de sal se hace entrar un flujo a una velocidad de 20 litros/min con una concentración de 0.08Kg/litro. La solución se mantiene bien agitada. Por otro lado sale un flujo de 20 litros/min, ¿cuál es la cantidad de sal pasados 30min?.

Condición inicial C0(0)=50 se requiere calcular yn(30)

f=@(x,y) 1.6-0.02*y; %ecuacion diferencial
x0=0; %tiempo inicial
y0=50; %concentracion inicial
xn=30; %tiempo final
n=30; %numero de pasos
[t,c]=edo_euler_modificado(f,x0,y0,xn,n); %llamada al metodo de Euler modificado
plot(t,c,'--.'),grid minor %grafica de la solucion
title('Mezclas')

Método de Runge-Kutta de 4to orden

Como observamos, el método de Euler mejorado tiene un resultado más cercano al analítico que el método de Euler simple, ¿por qué se mejora el resultado? la estrategia es hacer un promedio de las dos pendientes. Los métodos de Runge-Kutta básicamente aplican ese mismo principio, solo que hacen un promedio ponderado de las pendientes de puntos intermedios.

Como en los métodos de Euler y Euler mejorado, el resultado se aproxima más al correcto cuando se divide el intervalo [x0,xn] en n pasos. En el método de Runge-Kutta también se mejora dividiendo en n pasos.

function [xi,yi]=edo_rk4(f,x0,y0,x1,n)
     h=(x1-x0)/n; % tamano de paso
     xi=zeros(n+1,1); % vector de x variable independiente
     yi=zeros(n+1,1); % vector de y variable dependiente
     xi(1)=x0; % tiempo iniicial
     yi(1)=y0; % concentracion inicial

     for i=1:n
         k1=f(xi(i),yi(i));
         k2=f(xi(i)+h/2, yi(i)+k1*h/2);
         k3=f(xi(i)+h/2, yi(i)+k2*h/2);
         k4=f(xi(i)+h, yi(i)+k3*h);
         yi(i+1)=yi(i)+h*(k1/6+k2/3+k3/3+k4/6);
         xi(i+1)=xi(i)+h;
     end
end

Ejemplo. Un tanque esférico de radio R está inicialmente lleno de agua. En el fondo del tanque hay un agujero de radio r por el cual escapa el agua bajo la influencia de la gravedad. La ecuación diferencial que expresa la profundidad del agua como función del tiempo es:

donde g = 32.2ft / s^2, R = 12ft, r = 1/8ft. La condición inicial es que en t = 0, y = 22. Encontrar la altura del agua al minuto 1000.

La ecuación diferencial debe estar expresada en su forma básica

g =32.2;
R=12;
r=1/8;
f=@(x,y) -(r^2*sqrt(2*g))./(2*R*sqrt(y)-sqrt(y.^3) ); %ecuación diferencial en su forma básica
x0=0; %valor inicial del tiempo
y0=22; %valor inicial de la altura
xn=1000; %valor final del tiempo
n=10; %numero de pasos

[t,y]=edo_rk4(f,x0,y0,xn,n); %llamada al método de Runge-Kutta de 4to orden
plot(t,y,'--o') %gráfica de la solución
xlabel('tiempo');
ylabel('altura');

PROGRAM MetodoRungeKutta4
REAL :: a,b,y0,y
INTEGER :: n

PRINT *, "a : "
READ *, a
PRINT *, "b : "
READ *, b
PRINT *, "y0 : "
READ *, y0
PRINT *, "n : "
READ *, n

y=RungeKutta4(a,b,n,y0)

PRINT *, "y = ", y

CONTAINS

!Funcion
REAL FUNCTION f(x,y)
REAL :: x,y
f=2*x*y
END FUNCTION f

!Metodo RungeKutta4
REAL FUNCTION RungeKutta4(a,b,n,y0)
  REAL :: a,b,x,h,y0,y1
  REAL :: k1,k2,k3,k4
  INTEGER :: n
  h=(b-a)/n
  x=a
  DO i=1,n,1
    k1=f(x,y0)
    k2=f(x+h/2,y0+h/2*k1)
    k3=f(x+h/2,y0+h/2*k2)
    k4=f(x+h,y0+h*k3)
    y1=y0+h*(k1/6+k2/3+k3/3+k4/6)
    x=x+h
    y0=y1
  ENDDO
  RungeKutta4=y1
END FUNCTION RungeKutta4

END MetodoRungeKutta4