• Author / Uploaded
• Frank Chac

Citation preview

METODOS NUMERICOS

UNIDAD 6 SOLUCION NUMÉRICA DE ECUACIONES DIFERENCIALES ORDINARIAS

Una ecuación diferencial es una ecuación en la que intervienen derivadas de una o más funciones. Dependiendo del número de variables independientes respecto de las que se deriva, las ecuaciones diferenciales se dividen en: • •

Ecuaciones diferenciales ordinarias: aquellas que contienen derivadas respecto a una sola variable independiente. Ecuaciones en derivadas parciales: aquellas que contienen derivadas respecto a dos o más variables.

Las ecuaciones diferenciales aparecen naturalmente al modelar situaciones físicas en las ciencias naturales, ingeniería, y otras disciplinas, donde hay envueltas razones de cambio de una ó varias funciones desconocidas con respecto a una ó varias variables independientes. Estos modelos varían entre los más sencillos que envuelven una sola ecuación diferencial para una función desconocida, hasta otros más complejos que envuelven sistemas de ecuaciones diferenciales acopladas para varias funciones desconocidas. Por ejemplo, la ley de enfriamiento de Newton y las leyes mecánicas que rigen el movimiento de los cuerpos, al ponerse en términos matemáticos dan lugar a ecuaciones diferenciales. Usualmente estas ecuaciones están acompañadas de una condición adicional que especifica el estado del sistema en un tiempo o posición inicial. Esto se conoce como la condición inicial y junto con la ecuación diferencial forman lo que se conoce como el problema de valor inicial. Por lo general, la solución exacta de un problema de valor inicial es imposible ó difícil de obtener en forma analítica. Por tal razón los métodos numéricos se utilizan para aproximar dichas soluciones. Las leyes que gobiernan los fenómenos de la naturaleza se expresan habitualmente en forma de ecuaciones diferenciales. Las ecuaciones del movimiento de los cuerpos (la segunda ley de Newton) es una ecuación diferencial de segundo orden, como lo es la ecuación que describe los sistemas oscilantes, la propagación de las ondas, la transmisión del calor, la difusión, el movimiento de partículas subatómicas, etc. Pocas ecuaciones diferenciales tienen una solución analítica sencilla, la mayor parte de las veces es necesario realizar aproximaciones, estudiar el 1   

comportamiento del sistema bajo ciertas condiciones. Así, en un sistema tan simple como un péndulo, la amplitud de la oscilación ha de ser pequeña y el rozamiento ha de ser despreciable, para obtener una solución sencilla que describa aproximadamente su movimiento periódico. Se estudia el procedimiento de Runge-Kutta que se aplica de forma directa a una ecuación diferencial de primer orden, pero veremos cómo se extiende a un sistema de ecuaciones de primer orden, a un ecuación diferencial de segundo orden y a un sistema de ecuaciones diferenciales de segundo orden. El procedimiento de Runge-Kutta se puede programar fácilmente en los ordenadores y además, se emplea mucho en la práctica, debido a la su exactitud relativamente elevada de la solución aproximada de la ecuación diferencial.

6.1 Metodos de Euler y Euler modificado.

Método de EULER En matemática y computación, el método de Euler, llamado así en honor de Leonhard Euler, es un procedimiento de integración numérica para resolver ecuaciones diferenciales ordinarias a partir de un valor inicial dado. El método de Euler es el más simple de los métodos numéricos, resuelve un problema del tipo siguiente:

Método de 4 pasos para resolver ecuaciones diferenciales de primer orden con valores iniciales mediante el método de EULER. Fórmulas usadas:

yi+1 = yi + h f(xi,yi)

xi+1 = xi + h

Donde: i=0, 1, 2, 3,… h = tamaño del incremento con x f(xi,yi) = segundo miembro de la ED de primer orden cuando tiene la forma: dy/dx = f(x,y) 2   

Esta fórmula es conocida como el método de Euler (punto medio). Se predice un nuevo valor de Y por medio de la pendiente (igual a la primera derivada en el valor original de X).

Pasos: i. ii. iii.

iv.

Escribimos la ED en la forma: dy/dx = f(x,y), para extraer su segundo miembro Definimos x0, y0 y h de acuerdo a los datos del problema. Plateamos la ecuación de Euler utilizando los datos iniciales, como sigue: y0+1=y0+h f(x0,y0) Y una vez obtenido este primer resultado repetimos el proceso iterativamente utilizando los nuevos datos: y1+1=y1+h f(x1,y1) Desarrollamos hasta el valor buscado en x.

Ejercicio: use el método de EULER para obtener una aproximación a cuatro decimales del valor indicado, ejecute a mano la ecuación de recursión yn+1=yn+h f(xn,yn), usando h=0.1 y′=2x–3y+1, y(1)=5, aproximar y(1.2) Pasos: i. Escribimos la ED en la forma: dy/dx=f(x,y), para extraer su segundo miembro dy/dx = f(x,y) = 2x–3y+1 ii. Definimos x0, y0 y h de acuerdo a los datos del problema x0=1, 3   

iii.

iv.

y0=5 h=0.1 Planteamos la ecuación de Euler utilizando los datos iniciales: y0+1 = y0 + h f(x0,y0) y1 = y0 + h (2x0–3y0+1) y1 = 5+(0.1)(2(1)–3(5)+1) Desarrollamos hasta el valor buscado en x, en este caso: x=1.2 (límite superior del intervalo). x1 = x0 + h = 1 + 0.1 = 1.1 y1 = 5+(0.1)(2–15+1) y1 = 5+(0.1)(–12) y1 = y(1.1) = 3.8000

x2 = x1 + h = 1.1 + 0.1 = 1.2 y1+1 = y1 + h f(x1,y1) y2 = y1 + h*(2x1–3y1+1) y2 = 3.8 + (0.1)*(2*1.1-3*3.8+1) y2 = y(1.2) = 2.9800 Nota . El hecho de que la x varíe en x0=1, x1=1.1, x2=1.2, etc, tras cada iteración es porque la x aumenta según la fórmula: xn+1=xn + h.

Ejercicio: Dada la siguiente ecuación diferencial con la condición inicial: y’=2xy, y(0) =1, aproximar y(0.5). NOTA: Primero observamos que esta ecuación sí puede resolverse por métodos tradicionales de ecuaciones diferenciales. Por ejemplo, podemos aplicar el método de separación de variables. Veamos las dos soluciones. Solución Analítica: dy/dx = 2x y

dy/y = 2 x dx

ᶴdy/y = ᶴ2 x dx

lnlyl = x2 + c

Sustituyendo la condición inicial: X=0, y=1 ln(1) = 02+c Por lo tanto, tenemos que la curva solución real está dada: y=exp(x2) y(0.5) = exp(0.52) = 1.28403 ln(y)=x2 Solución Numérica: Aplicamos el método de Euler y para ello, observamos que la distancia entre x0=0 y x1=0.5 no es lo suficientemente pequeña. Si dividimos esta distancia entre cinco 4   

obtenemos un valor de h=0.1 y por lo tanto, obtendremos la aproximación deseada en cinco pasos. De esta forma, tenemos los siguientes datos: x0=0 y0=1 h=0.1 f(x,y)=2xy Sustituyendo estos datos en la fórmula de Euler, tenemos, en un primer paso: x1 = x0+h = 0.1 y1 = y0+h f(x0,y0) = 1+0.1*(2*0*1) = 1 Aplicando nuevamente la fórmula de Euler, tenemos, en un segundo paso: x2 = x1+h = 0.2 y2 = y1+h f(x1,y1) = 1+0.1*(2*0.1*1) = 1.02 Y así sucesivamente hasta obtener y5. Resumimos los resultados en la siguiente tabla: n 0 1 2 3 4 5

xn 0 0.1 0.2 0.3 0.4 0.5

yn 1 1 1.02 1.0608 1.12445 1.2144

Concluimos que el valor aproximado, usando el método de Euler es: y(0.5) = 1.2144 Puesto que en este caso, conocemos el valor verdadero, podemos usarlo para calcular el error relativo porcentual que se cometió al aplicar la fórmula de Euler. Tenemos que: lϵvl = abs((1.28402 – 1.2144)/1.28402)*100% = 5.42%

Método de Euler mejorado Este método se basa en la misma idea del método anterior, pero hace un refinamiento en la aproximación, tomando un promedio entre ciertas pendientes. La fórmula es la siguiente: 5   

Donde

Para entender esta fórmula, analicemos el primer paso de la aproximación, con base en la siguiente gráfica:

En la gráfica, vemos que la pendiente promedio corresponde a la pendiente de la recta bisectriz de la recta tangente a la curva en el punto de la condición inicial y la "recta tangente" a la curva en el punto donde es la aproximación obtenida con la primera fórmula de Euler. Finalmente, esta recta bisectriz se traslada paralelamente hasta el punto de la condición inicial, y se considera el valor de esta recta en el punto

como la aproximación de Euler mejorada.

Ejercicio: Dada la siguiente ecuación diferencial con la condición inicial: y’=2xy, y(0) =1, aproximar.

Ejercicio: Aplicar el método de Euler mejorado, para aproximar y(0.5) si: y’=2xy, y(0) =1 Solución: Vemos que este es el mismo ejercicio del método anterior. Así que definimos h=0.1 y encontraremos la aproximación después de cinco iteraciones. A diferencia del método de Euler 1, en cada iteración requerimos de dos cálculos en vez de uno solo: el de yn* primero y posteriormente el de yn. Para aclarar el método veamos con detalle las primeras dos iteraciones. Primero que nada, aclaramos que tenemos los siguientes datos iniciales: 6   

x0=0 y0=1 h=0.1 f(x,y)=2xy En nuestra primera iteración tenemos: x1 = x0+h = 0.1 y1* = y0+h f(x0,y0) = 1+0.1*(2*0*1) = 1 y1 = y0+h ((f(x0,y0)+f(x1,y1*))/2) = 1.01 Nótese que el valor de y1* coincide con el y1 (Euler 1), y es el único valor que va a coincidir, pues para calcular y2* se usará y1 y no y1* Esto lo veremos claramente en la siguiente iteración: x2 = x1+h = 0.1+0.1 = 0.2 y2* = y1+h f(x1,y1) = = 1.0302 y2 = y1+h ((f(x1,y1)+f(x2,y2*))/2) = 1.040704 Nótese que ya no coinciden los valores de y2 (Euler 1) y el de y2*. El proceso debe seguirse hasta la quinta iteración. Resumimos los resultados en la tabla siguiente: n 0 1 2 3 4 5

xn 0 0.1 0.2 0.3 0.4 0.5

yn 1 1.01 1.040704 1.093988 1.173192 1.28336

Se concluye entonces que la aproximación obtenida con el método de Euler mejorado es: y(0.5) = 1.28336 Con fines de comparación, calculamos el error relativo verdadero: lϵvl = abs((1.28402 – 1.28336)/1.28402)*100% = 0.05% Vemos que efectivamente se ha obtenido una mejor aproximación con este método, reduciendo el error relativo verdadero de un 5.4% hasta un 0.05%. En nuestro tercer método veremos cómo se reduce aún más este error prácticamente a un 0%!

7   

6.2 Método de Runge Kutta de cuarto orden

Los métodos de Taylor tienen la propiedad de un error local de truncamiento de orden superior, pero la desventaja de requerir el cálculo y la evaluación de las derivadas de f(t,y). Esto resulta algo lento y complicado, en la mayoría de los problemas, razón por la cual, en la práctica casi no se utilizan. El método de Euler, lamentablemente requiere de un paso muy pequeño para una precisión razonable. Los métodos de Runge kutta tienen el error local de truncamiento del mismo orden que los métodos de Taylor, pero prescinden del cálculo y evaluación de las derivadas de la función f(t,y). Se presenta de nuevo el problema de valor inicial cuya solución se intenta aproximar:

Como en los métodos anteriores, se determina primero la malla {t0,t1,...,tN} de paso h, donde t0 = a y tN = b. En estos puntos es donde se va a obtener la aproximación de la solución. En esencia, los métodos de Runge-Kutta son generalizaciones de la fórmula básica de Euler yi+1 = yi + h f(ti, yi) en los que el valor de la función f se reemplaza por un promedio ponderado de valores de f en el intervalo ti≤t≤ti+1, es decir,

En esta expresión las ponderaciones wi, i = 1,...,m son constantes para las que en general se pide que su suma sea igual a 1, es decir, w1+w2+...+wm = 1, y cada kj es la función f evaluada en un punto seleccionado (t,y) para el cual ti≤t≤ti+1. Se mostrará que los kj se definen en forma recursiva. Se define como orden del método al número m, es decir, a la cantidad de términos que se usan en el promedio ponderado. 8   

Runge-Kutta de primer orden. Si m = 1, entonces se toma w1 = 1 y la fórmula resulta

Igualando esta fórmula al desarrollo de Taylor de orden 1 de la función y(t), alrededor del punto ti, y calculado en el punto ti+1:

y teniendo en cuenta que yi=y(ti), resulta k1= f(ti,yi), obteniendo así la fórmula de Euler yi+1 = yi + h f(ti,yi). Por lo tanto, se dice también que el método de Euler es un método de Runge Kutta de primer orden. Runge-Kutta de segundo orden. Ahora se plantea, con m = 2, una fórmula del tipo: Donde

y las constantes a, b, se deben determinar, de manera que la expresión coincida con el desarrollo de Taylor de y de orden más alto posible. Para ello, utilizando un desarrollo de Taylor para funciones de dos variables, tenemos que:

donde el subíndice i indica que todas las derivadas están evaluadas en el punto (ti,yi). Reemplazando k1 y teniendo en cuenta la expresión de k2, tenemos que:

Agrupando los términos de por las potencias de h, y reemplazando el valor de k1 y k2, resulta

Reacomodando términos, resulta:

9   

Por otro lado, se hace un desarrollo de Taylor de orden 3 de la función y(t), calculado en el punto ti+1, obteniendo:

Aplicando regla de la cadena para las derivadas de f, se tiene:

Comparando las expresiones e igualando los coeficientes de h y h2, se tiene:

Sucede que se tienen cuatro incógnitas, pero tres ecuaciones, con lo que queda un grado de libertad en la solución del sistema. Se trata de usar este grado de libertad para hacer que los coeficientes de h3 en las expresiones coincidan. Esto obviamente no se logra para cualquier f. Hay muchas soluciones para el sistema, una de ellas es

obteniendo así la siguiente fórmula, del método de Runge Kutta de orden 2:

para i desde 0 hasta N-1, tomando un mallado {ti, i = 0,..,N}. Este método tiene un error local de O(h3), y global de O(h2). Mejora entonces el método de Euler, por lo que se espera poder usar con este método un paso mayor. El precio que debe pagarse en este caso, es el de evaluar dos veces la función en cada iteración. De la misma manera que se realizó arriba, se pueden derivar fórmulas de RungeKutta de cualquier orden, pero estas deducciones resultan excesivamente complicadas. Una de las más populares, y más utilizada por su alta precisión, es la de orden 4, que se presenta a continuación. 10   

Runge-Kutta de cuarto orden. Si ahora m = 4, se obtiene, con un desarrollo del tipo del anterior, la siguiente fórmula, para i desde 0 hasta N-1:

Si bien con facilidad se pueden deducir otras fórmulas, el algoritmo expresado arriba se denomina método de Runge-Kutta de cuarto orden, o método clásico de Runge-Kutta, abreviado como RK4. Este algoritmo es de uso extendido, y reconocido como una valiosa herramienta de cálculo, por la buena aproximación que produce. Esta fórmula tiene un error de truncamiento local de O(h5), y un error global de O(h4). De nuevo, el precio que se debe pagar por la mejora en el error, es una mayor cantidad de evaluaciones de la función, resultando en un mayor tiempo de cálculo si la función es complicada. Tiene la ventaja, sobre el método de Taylor de orden 4 (cuyo error global es también O(h4), que no requiere el cálculo de las derivadas de f. Implementación del método RK4. Se presenta a continuación el pseudocódigo del método RK4, para ser implementado en cualquier lenguaje de programación, o software simbólico.

11   

Ejercicio: Con el método RK4, obtener una aproximación del valor de y(1.5) para el siguiente problema de valor inicial, tomando un paso h = 0.1.

El primer paso para resolver este problema es determinar la malla de puntos en donde se va a obtener la solución. h=(xN-x0)/N Como en este caso h está dado, se tiene que N = (1.5 - 1)/0.1 = 5. Por lo tanto, los puntos en donde se va a determinar la solución, dados por la fórmula ti = 1 + 0.1*i, para i =1,2,3,4,5, son: t1 = 1.1, t2 = 1.2, t3 = 1.3, t4 = 1.4 y t5 = 1.5 Una vez establecida la malla del problema, tenemos, para i = 0:

Resulta entonces,

y aplicando sucesivamente la fórmula de RK4, para i desde 1 hasta 5, se obtienen los datos que se muestran en la siguiente tabla, donde además se muestra el valor de la solución exacta para cada punto de la malla.

12   

Al analizar la tabla anterior y comparar los resultados obtenidos con el método RK4 con los valores reales, se ve por qué es tan difundido este método. En la próxima tabla se comparan los métodos de Euler y Runge Kutta de orden 4 para el mismo problema.

6.3 Sistemas de dos ecuaciones y ecuaciones de orden superior.

Sistema de dos ecuaciones diferenciales de primer orden El procedimiento de Runge-Kutta es igualmente efectivo en la resolución de un sistema de dos ecuaciones diferenciales de primer orden. dx/dt = f(t,x,y) dy/dt = g(t,x,y) El procedimiento de aplicación del método de Runge-Kutta a cada una de las ecuaciones diferenciales, con la condición inicial siguiente, en el instante t0 el valor inicial de x es x0 el valor inicial de y es y0

13   

se esquematiza en la tabla adjunta. Como vemos además de los cuatro números k1, k2, k3, k4 para la primera ecuación diferencial precisamos otros cuatro números l1, l2, l3, l4 para la segunda ecuación diferencial. A partir del valor de x en el instante t, se determina el valor de x en el instante t+h, y a partir del valor de y en el instante t se determina el valor de y en el instante t+h mediante las fórmulas de la última fila de la tabla. dx/dt = f(t,x,y) k1 = h f(t,x,y) k2 = h f(t+½ h, x+½ k1, y+½l1,) k3 = h f(t+½ h, x+½k2, y+½l2) k4 = h f(t+h, x+k3, y+l3) x(t+h) = x(t)+1/6(k1+2k2+2k3+k4)

dy/dt = g(t,x,y) l1 = h g(t,x,y) l2 = h g(t+½h, x+½k1, y+½l1) l3 = h g(t+½h, x+½k2, y+½l2) l4 = h g(t+h, x+k3,y+l3) y(t+h) = y(t)+1/6(l1+2l2+2l3+l4)

Nos devuelve los vectores x e y para cada instante que se guarda en el vector t comprendido entre el instante inicial t0 y el final tf. function [t,x,y] =rk_2_1(f,g,t0,tf,x0,y0,n) h=(tf-t0)/n; t=t0:h:tf; x=zeros(n+1,1); %reserva memoria para n+1 elementos del vector x(i) y=zeros(n+1,1); x(1)=x0; y(1)=y0; for i=1:n k1=h*f(t(i),x(i),y(i)); l1=h*g(t(i),x(i),y(i)); k2=h*f(t(i)+h/2,x(i)+k1/2,y(i)+l1/2); l2=h*g(t(i)+h/2,x(i)+k1/2,y(i)+l1/2); k3=h*f(t(i)+h/2,x(i)+k2/2,y(i)+l2/2); l3=h*g(t(i)+h/2,x(i)+k2/2,y(i)+l2/2); k4=h*f(t(i)+h,x(i)+k3,y(i)+l3); l4=h*g(t(i)+h,x(i)+k3,y(i)+l3); x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6; y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6; end end Consideremos una serie radioactiva de tres elementos A-->B-->C en la que, una sustancia radiactiva A se desintegra y se transforma en otra sustancia radiactiva B, que a su vez se desintegra y se transforma en una sustancia C estable. Las ecuaciones diferenciales que gobiernan el proceso y sus soluciones analíticas son, respectivamente, 14   

dx/dt = −ax

x = x0 e−at

dy/dt = ax−by

y = a/(b−a) x0 (e−at−e−bt)

La solución analítica que aparece a la derecha, se ha obtenido con las condiciones iniciales t=0, x=x0 e y=0. La segunda solución se obtiene siempre que a sea distinto de b. En el caso de que a sea igual a b, la solución analítica para y es y=x0 a e−at La interpretación del sistema de ecuaciones diferenciales no es complicada. En la unidad de tiempo, desaparecen ax núcleos de la sustancia A al desintegrarse (primera ecuación). En la unidad de tiempo, se producen ax núcleos de la sustancia B y a su vez desaparecen bx núcleos de la sustancia B, que al desintegrarse se transforman en núcleos de la sustancia C estable (segunda ecuación). Escribimos el script radioactivo en el que definiremos las funciones f(t,x,y), g(t,x,y), las condiciones iniciales y llamaremos a la función rk_2_1 a=input('parámetro a: '); b=input('parámetro b: '); x0=input('valor inicial de x: '); y0=input('valor inicial de y: '); tf=input('tiempo final, tf: '); n=input('número de pasos, n: '); f=@(t,x,y) -a*x; g=@(t,x,y) a*x-b*y; %condiciones iniciales t0=0; [t,x,y]=rk_2_1(f,g,t0,tf,x0,y0,n); hold on plot(t,x,'b') plot(t,y,'r') xlabel('t') ylabel('x,y'); legend('x(t)','y(t)') title('dx/dt=-ax, dy/dt=ax-by') hold off En la ventana de comandos corremos el script radioactivo 15   

>> radioactivo parámetro a: 0.1 parámetro b: .2 valor inicial de x: 100 valor inicial de y: 0 tiempo final, tf: 10 número de pasos, n: 40

Ecuación diferencial de segundo orden Existen muchas situaciones en las que es necesario resolver una ecuación diferencial de segundo orden. d2x/dt2 = f(t,x,v) con las condiciones iniciales: x(t0) = x0

(dx/dt)t0 = v0

Una ecuación diferencial de segundo orden es equivalente a un sistema de dos ecuaciones diferenciales de primer orden, por lo que aplicaremos el mismo esquema. dx/dt =v k1 = h v k2 = h (v+½l1) k3 = h (v+½l2) k4 = h (v+l3) x(t+h) = x(t)+1/6 (k1+2k2+2k3+k4)

dv/dt = f(t,x,v) l1=h f(t,x,v) l2=h f(t+½h, x+½k1, v+½l1) l3=h f(t+½h, x+½k2, v+½l2) l4=h f(t+h, x+k3, v+l3) v(t+h) = v(t)+1/6 (l1+2l2+2l3+l4) 16 

 

Definimos la función rk_2 que resuelve la ecuación diferencial de segundo orden, cuando le pasamos: • la función f (t,x,v) • las condiciones iniciales: posición inicial x0 y velocidad inicial v0 en el instante t0 • el número n de pasos de integración entre t0 y el tiempo final tf Nos devuelve los vectores de las posiciones x y las velocidades v para cada instante que se guarda en el vector t comprendido entre el instante inicial t0 y el final tf. function [t,x,v] =rk_2(f,t0,tf,x0,v0,n) h=(tf-t0)/n; t=t0:h:tf; x=zeros(n+1,1); %reserva memoria para n+1 elementos del vector x(i) v=zeros(n+1,1); x(1)=x0; v(1)=v0; for i=1:n k1=h*v(i); l1=h*f(t(i),x(i),v(i)); k2=h*(v(i)+l1/2); l2=h*f(t(i)+h/2,x(i)+k1/2,v(i)+l1/2); k3=h*(v(i)+l2/2); l3=h*f(t(i)+h/2,x(i)+k2/2,v(i)+l2/2); k4=h*(v(i)+l3); l4=h*f(t(i)+h,x(i)+k3,v(i)+l3); x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6; v(i+1)=v(i)+(l1+2*l2+2*l3+l4)/6; end end La ecuación diferencial que describe un oscilador armónico amortiguado y su solución para unas condiciones iniciales fijadas es 2

d2x/dt2 + 2‫ ץ‬dx/dt + ω

0

x=0

2

ω = √ω − ‫ץ‬2 0

x=e−‫ץ‬t (A sin(ωt)+B cos(ωt)) 17   

v= dx/dt = − ‫ ץ‬e−‫ץ‬t (A sin(ωt)+B cos(ωt))+ω e−γt (A cos(ωt)−B sin(ωt)) Condiciones iniciales: en el instante t=0, la posición inicial es x0 y la velocidad inicial v0. x0 = B t=0 v0 = − ‫ץ‬B+ωA x=e−‫ץ‬t ((v0+‫ץ‬x0)/ ω sin(ωt)+x0 cos(ωt)) Escribimos el script oscilador en el que definiremos la función f(t,x,v), las condiciones iniciales y llamaremos a la función rk_2 w0=input('frecuencia angular w0: '); g=input('rozamiento, gamma: '); x0=input('posición inicial, x0: '); v0=input('velocidad inicial,v0: '); tf=input('tiempo final, tf: '); n=input('número de pasos, n: '); f=@(t,x,v) -2*g*v-w0*w0*x; %condiciones iniciales t0=0; hold on %solucion numérica [t,x,v]=rk_2(f,t0,tf,x0,v0,n); plot(t,x,'b') %solución analítica w=sqrt(w0*w0-g*g); x=((v0+g*x0)*sin(w*t)/w+x0*cos(w*t)).*exp(-g*t); plot(t,x,'r') grid on xlabel('t') ylabel('x'); legend('aproximado','exacto') title('oscilador amortiguado') hold off En la ventana de comandos corremos el script oscilador con distintas condiciones iniciales >> oscilador frecuencia angular, w0: 2 18   

rozamiento, gamma: 0.5 posición inicial, x0: 1.5 velocidad inicial, v0: 0 tiempo final, tf: 8 número de pasos, n: 100

No se aprecia tampoco diferencia entre la solución exacta y la numérica, aplicando el procedimiento de Runge_Kutta.

6.4 Aplicaciones de la solución numérica de ecuaciones diferenciales ordinarias.

Mecánica La física trata de la investigación de las leyes que gobiernan el comportamiento del universo físico. Por universo físico entendemos la totalidad de objetos a nuestro alrededor, no sólo las cosas que observamos sino también las que no observamos, tales como los átomos y moléculas. El estudio del movimiento de los objetos en nuestro universo es una rama de la mecánica llamada dinámica formulada mediante las leyes del movimiento de Newton. Para los objetos que se mueven muy rápido, cerca de la velocidad de la luz, no podemos usar las leyes de Newton. En su lugar debemos usar una versión revisada de estas leyes, desarrolladas por Einstein y conocidas como mecánica relativista, o mecánica de la relatividad. Para objetos de dimensiones atómicas las leyes de Newton tampoco son válidas. De hecho, para obtener descripciones precisas del movimiento de 19   

objetos de dimensiones atómicas, necesitamos establecer un conjunto de leyes denominadas mecánica cuántica Circuitos eléctricos Así como la mecánica tiene como base fundamental las leyes de Newton, la electricidad también tiene una ley que describe el comportamiento de los circuitos eléctricos, conocida como la ley de Kirchhoff. Realmente, la teoría de la electricidad está gobernada por un cierto conjunto de ecuaciones conocidas en la teoría electromagnética como las ecuaciones de Maxwell. La ley de Kirchhoff es adecuada para estudiar las propiedades simples de los circuitos eléctricos. El circuito eléctrico más simple es un circuito en serie, en el cual tenemos una fem (fuerza electromotriz), la cual actúa como una fuente de energía tal como una batería o generador, y una resistencia, la cual consume o usa energía, tal como una bombilla eléctrica, tostador, u otro electrodoméstico. En física elemental encontramos que la fem está relacionada con el flujo de corriente en el circuito. En forma simple, la ley dice que la corriente instantánea I (en un circuito que contiene sólo una fem E y una resistencia R) es directamente proporcional a la fem. Simbólicamente: I α E de donde, E = IR donde R es una constante de proporcionalidad llamada el coeficiente de resistencia o simplemente, resistencia. La ecuación anterior es conocida bajo el nombre de la ley de Ohm. Circuitos más complicados, pero para muchos casos más prácticos, son circuitos que contienen otros elementos distintos a resistencias. Dos elementos importantes son inductores y condensadores. Un inductor se opone a cambios en corriente. Tiene un efecto de inercia en electricidad de la misma manera que una masa tiene un efecto de inercia en mecánica. Un condensador es un elemento que almacena energía. En física hablamos de una caída de voltaje a través de un elemento. En la práctica podemos determinar esta caída de voltaje, o como se llama comúnmente, caída de potencial o diferencia de potencial, por medio de un instrumento llamado voltímetro. Aplicaciones a flujo de calor en estado estacionario. Considere una pieza de material de longitud indefinida acotada por dos planos paralelos A y B, según la figura, suponiendo que el material es uniforme en todas sus propiedades, por ejemplo, calor específico, densidad, etc. 20   

Considerando que los planos A y B se mantienen a 50º C y 100ºC, respectivamente, todo punto en la región entre A y B alcanza cierta temperatura que no cambia posteriormente. Así todos los puntos en el plano C entre A y B estarán a 75ºC,y en el plano E a 90ºC. Cuando la temperatura en cada punto de un cuerpo no varía con el tiempo decimos que prevalecen las condiciones de estado estacionario o que tenemos un flujo de calor en estado estacionario. Aplicaciones a problemas combinados de crecimiento y decrecimiento. La ecuación diferencial dt/dy = ay nos dice que la variación con el tiempo de una cantidad y es proporcional a y. Si la constante de proporcionalidad a es positiva siendo y positivo, entonces dy/dt es positivo e y aumenta. En este caso hablamos de que y crece, y el problema es de crecimiento. Por otro lado, si a es negativo siendo y positivo, entonces dy/dt es negativo e y decrece. Aquí el problema es uno que involucra decrecimiento. Puesto que la solución de dt/dy = ay identificada como una ecuación de variables separadas está dada por la función exponencial y = ceat , resolviendo mediante integración, definiéndose la ecuación diferencial dt/dy = ay como la ley de crecimiento exponencial si a > 0 y la ley de decrecimiento exponencial si a > carga_1 Resistencia R: 2 Capacidad C: 0.8 tiempo final, tf: 10 Sistema de dos ecuaciones diferenciales de primer orden 24   

Elaboramos el script titulado radiactivo_1 para integrar el sistema de dos ecuaciones diferenciales de primer orden que describe la serie de desintegración radioactiva. A-->B-->C donde C es un elemento estable. dx/dt = −ax dy/dt = ax−by En la matriz x que devuelve la función ode45, x(1) representará los sucesivos valores de la variable x y x(2) representará a la variable y. El mismo criterio se empleará para determinar el vector x0 de las condiciones iniciales. La definición de las funciones f(t,x,y) y g(t,x,y) aparecen en un vector columna, separadas por ; (punto y coma) fg=@(t,x) [-a*x(1);a*x(1)-b*x(2)]; % x(1) es x, x(2) es y El script radiactivo_1 será el siguiente: a=input('parámetro a: '); b=input('parámetro b: '); %condiciones iniciales en el vector x0 x0=zeros(1,2); x0(1)=input('valor inicial de x: '); x0(2)=input('valor inicial de y: '); tf=input('tiempo final, tf: '); tspan=[0 tf]; fg=@(t,x) [-a*x(1);a*x(1)-b*x(2)]; [t,x]=ode45(fg,tspan,x0); plot(t,x) xlabel('t') ylabel('x,y'); title('dx/dt=-ax, dy/dt=ax-by') En la ventana de comandos corremos el script radiactivo_1 >> radioactivo_1 parámetro a: 0.1 parámetro b: 0.2 valor inicial de x: 100 valor inicial de y: 0 tiempo final, tf: 20 Alternativamente, vamos a definir las funciones f(t,x,y) y g(t,x,y) en un fichero .m y le pasamos los valores de los parámetros a y b. 25   

function z=func_radioactivo(t,x,a,b) z=[-a*x(1);a*x(1)-b*x(2)]; % x(1) es x, x(2) es y end Elaboramos el script radioactivo_2 para establecer los valores de los parámetros a y b, las condiciones iniciales y llamar a la función que realiza la integración numérica ode45. El primer parámetro de ode45 es el handler (manejador de la función) a integrar que se especifica del siguiente modo @nombre_funcion. [t,x]=ode45(@func_radioactivo,tspan,x0); Ahora bien, func_radioactivo precisa de los valores de los parámetros a y b. Hay dos formas de hacerlo. La más sencilla es definir una función anónima fg en términos de func_radioactivo. fg=@(t,x) func_radioactivo(t,x,a,b); Véase la misma situación en la llamada a función fzero al final de la página Raíces de una ecuación a=input('parámetro a: '); b=input('parámetro b: '); %condiciones iniciales x0=zeros(1,2); x0(1)=input('valor inicial de x: '); x0(2)=input('valor inicial de y: '); tf=input('tiempo final, tf: '); tspan=[0 tf]; fg=@(t,x) func_radioactivo(t,x,a,b); [t,x]=ode45(fg,tspan,x0); plot(t,x) xlabel('t') ylabel('x,y'); title('dx/dt=-ax, dy/dt=ax-by') En la ventana de comandos corremos el script radioactivo_2 >> radioactivo_2 parámetro a: 0.1 parámetro b: 0.2 valor inicial de x: 100 valor inicial de y: 0 tiempo final, tf: 20 Opciones de ode45 Imaginemos que estudiamos el movimiento de caída de un cuerpo, no sabemos cuanto tiempo tardará en llegar al suelo, desconocemos el valor del elemento tf en el vector tspan. Sin embargo, queremos detener el proceso de integración 26   

numérica de la ecuación diferencial que describe el movimiento cuando la posición del móvil sea cero. La función MATLAB ode45 dispone de un parámetro adicional options donde podemos indicarlo, pero es bastante lioso e intentaremos explicarlo mediante ejemplos. Volvemos a resolver la ecuación diferencial que describe las oscilaciones amortiguadas y detendremos el proceso de integración cuando el móvil alcance la posición máxima, su velocidad es nula. Supongamos que el oscilador amortiguado estudiado anteriormente, de frecuencia natural ω0=2, constante de amortiguamiento γ=0.25, parte de la posición x0=2.5 con velocidad nula, queremos detener el proceso de integración cuando el móvil alcance la posición máxima, cuando su velocidad es nula, tal como se muestra en la figura, con lo que se completa un periodo.

Los pasos a seguir son los siguientes: 1.-Definimos la función cuyo nombre es opcion_ode45 function [detect,stopin,direction]=opcion_ode45(t,x) detect=[x(1) x(2)]; %[posición, velocidad] stopin=[0 1]; % 1 indice que detiene la integración cuando la velocidad se hace cero direction=[0 -1]; % 1 crece, -1 decrece, 0 no importa end 2.-Creamos la estructura opts con la llamada a la función odeset opts=odeset('events',@opcion_ode45); Cuando se utiliza options la función ode45 devuelve los tiempos te en los cuales ocurren los 'eventos' y los correspondientes valores de las variables dependientes 27   

(posición, velocidad) en el vector xe. Finalmente, ie es un índice que es útil cuando se vigilan varios eventos. 3.-Le pasamos opts a la función ode45 en su cuarto parámetro [t,x,te,xe,ie]=ode45(odefunc,tspan,x0,opts); Escribimos el script oscilador_2 para resolver la ecuación diferencial de segundo orden y detener el proceso de integración de acuerdo con lo estipulado en el parámetro opts. w0=2; g=0.25; x0=[2.5 0]; tf=10; f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; tspan=[0 10]; opts=odeset('events',@opcion_ode45); [t,x,te,xe,ie]=ode45(f,tspan,x0,opts); te,xe,ie plot(t,x(:,1),'r') grid on xlabel('t') ylabel('x'); title('oscilador amortiguado') Cuando corremos el script oscilador_2 en la ventana de comandos se imprime los siguientes datos relativos a los eventos. Tiempo, te Posición x(1) Velocidad x(2) Indice ie 0.0000 2.5000 -0.0000 2 2.4378 0.0000 2.7173 1 3.1662 1.1322 -0.0000 2 • Cuando parte de la posición inicial x(1)=2.5 la velocidad es cero x(2)=0, detecta velocidad (índice ie=2). • Cuando pasa por el origen x(1)=0 detecta posición (índice ie=1), pero no se detiene ya que en stopin se ha puesto un cero. • Cuando la posición es x(1)=1.1322 detecta velocidad nula x(2)=0, (índice ie=2) y la integración numérica se detiene ya que en stopin se ha puesto un uno y la velocidad decrece en direction se ha puesto un -1. La columna de tiempos nos proporcionan el periodo de la oscilación, te=3.1662. Vamos ahora a marcar en la representación gráfica de la oscilación amortiguada, las posiciones de máxima amplitud x(2)=0 y cuando pasa por el origen x(1)=0 sin detener el proceso de integración numérica. 28   

Definimos una nueva versión de la función opcion1_ode45 function [detect,stopin,direction]=opcion1_ode45(t,x) detect=[x(1) x(2)]; %[posición, velocidad] stopin=[0 0]; direction=[0 0]; end Creamos la estructura opts mediante odeset y se la pasamos al procedimiento de integración ode45. w0=2; g=0.25; x0=[2.5 0]; tf=10; f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; tspan=[0 7]; opts=odeset('events',@opcion1_ode45); [t,x,te,xe,ie]=ode45(f,tspan,x0,opts); hold on plot(t,x(:,1),'r') plot(te(ie==1),xe(ie==1),'o','markersize',6,'markerface color','k') plot(te(ie==2),xe(ie==2),'o','markersize',6,'markerface color','b') grid on xlabel('t') ylabel('x'); title('oscilador amortiguado') hold off 29   

Si solamente estamos interesados en los máximos definimos una nueva versión de la función opcion2_ode45 function [detect,stopin,direction]=opcion2_ode45(t,x) detect=x(2); stopin=0; direction=-1; end Modificamos el script oscilador_4 w0=2; g=0.25; x0=[2.5 0]; tf=10; f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; tspan=[0 7]; opts=odeset('events',@opcion1_ode45); [t,x,te,xe,ie]=ode45(f,tspan,x0,opts); te,xe,ie hold on plot(t,x(:,1),'r') plot(te(ie==1),xe(ie==1),'o','markersize',6,'markerfacecolor', 'b') grid on xlabel('t') ylabel('x'); title('oscilador amortiguado') hold off Corremos el script oscilador_4 en la ventana de comandos y observamos los resultados Paso de parámetros a la función Como hemos visto, a ode45 se le pasa la función (handle) a integrar en su primer argumento. Si la función contiene parámetros como la frecuencia angular ω0, no hay problema si la función se define como anónima, tal como se ha mostrado en los ejemplos previos. Si la función se define en un fichero entonces a la función se le pasan los valores de los parámetros en el quinto argumento params de ode45. Los pasos se explican en el siguiente ejemplo: El sistema de ecuaciones de Lorentz es un sistema de tres ecuaciones diferenciales de primer orden

30   

dx/dt = −σx+σy dy/dt = ρx−y−xz dz/dt = −βz+xy donde σ=10, β=8/3 y ρ=28 Vamos a resolver el sistema de tres ecuaciones diferenciales con las condiciones iniciales siguientes: en el instante t=0, x0=-8, y0=8 z0=27. 1. Escribir una función denominada lorentz(t,x,p) como fichero .m que contenga las tres ecuaciones y dependa de los tres parámetros σ=p(1), β=p(2) y ρ=p(3). Observar que la variable x se guarda en el primer elemento x(1), la variable y en el segundo x(2) y la variable z en el tercer x(3) elemento del vector x. 2. Escribir un script denominado lorentz_script que llame a la función MATLAB ode45, para resolver el sistema de ecuaciones diferenciales para las condiciones iniciales especificadas. 3. Pasar los parámetros σ, β y ρ como elementos de un vector p a la función ode45 para que a su vez se los pase a la función lorentz. 4. Dibujar el atractor de Lorentz de z en función de x hasta el instante tf=20 en una primera ventana gráfica. 5. Dibujar x, y y z en función del tiempo en tres zonas de una segunda ventana gráfica. 6. Examinar el comportamiento del sistema para otras condiciones iniciales, t=0, x0=1, y0=2 z0=3. Definimos la función lorentz como fichero .m function fg=lorentz(t,x,p) %x(1) es x, x(2) es y y x(3) es z % p(1) es sigma, p(2) es beta y p(3) es rho fg=[-p(1)*x(1)+p(1)*x(2); p(3)*x(1)-x(2)-x(1)*x(3); -p(2)*x(3)+x(1)*x(2)]; end

Elaboramos el script lorentz_script x0=[-8 8 27]; %valores iniciales tspan=[0 20]; p=[10 8/3 28]; %parámetros %no pasamos nada [] en el parámetro options de ode45 [t,x]=ode45(@lorentz,tspan,x0,[],p); figure plot(x(:,1),x(:,3),'r') xlabel('x') 31   

ylabel('z'); title('Atractor de Lorentz') figure subplot(3,1,1) plot(t,x(:,1)) ylabel('x'); subplot(3,1,2) plot(t,x(:,2)) ylabel('y'); subplot(3,1,3) plot(t,x(:,3)) ylabel('z'); xlabel('t') En la ventana de comandos corremos el script lorentz_script >> lorentz_script Aparecen dos ventanas gráficas, la primera con el atractor de Lorentz, la representación z(x) y la segunda ventana dividida en tres regiones que muestran x(t), y(t) y z(t)

32