Universidad Politécnica de Madrid–Escuela Técnica Superior de Ingenieros Industriales Grado en Ingeniería en Tecnologías Industriales. Curso 2015-2016-3º Matemáticas de Especialidad–Ingeniería Eléctrica
Integración de
Ecuaciones Diferenciales Ordinarias
José Luis de la Fuente O’Connor
[email protected] [email protected] Clase_integración ODE_2016.pdf
1/57
2/57
Índice
Introducción
Integración numérica de ecuaciones diferenciales ordinarias (ODEs)
Método de Euler
Mejoras en el método de Euler
Método de Taylor
Métodos de Runge-Kutta
Métodos implícitos y problemas stiff
3/57
Introducción
Una ecuación diferencial en su forma más simple se define de forma explícita así: dy D y 0 D f .t; y/; dt Su solución
siendo y.0/ D y0: x
Z y.x/ D y0 C
f .t; y/ dt: 0
El orden de una ecuación diferencial es el de su máxima derivada.
Las ecuaciones diferenciales relacionan las variables de un problema a escala diferencial, en la cual son independientes de las condiciones de contorno (geometría, condiciones iniciales; : : :)
4/57
Son el resultado de modelizar una amplia familia de problemas de ingeniería en los que se estudia cómo y por qué varia en el tiempo el comportamiento de sistemas físicos, económicos, sociales, etc.
Ejemplo La ecuación logística, o función logística, y 0 D cy.1
y/
modeliza el crecimiento del tamaño de una población, que es proporcional al tamaño actual, y, y a la capacidad remanente de recursos disponibles.
Esta ecuación diferencial tiene infinitas soluciones. Si se le especifica una condición de partida, inmediatamente se identifica de qué solución se trata.
Para poder integrar una ecuación diferencial es necesario definir las condiciones de contorno.
5/57
Un problema de valor inicial es el que define la ecuación diferencial y una condición, o valor inicial, en un intervalo a t b. En conjunto, 8 0 ˆ < y D f .t; y/
y.a/ D y
a ˆ : t 2 Œa; b:
Resulta útil pensar en la evolución de una ecuación diferencial en el tiempo como la de un cuerpo dentro de un campo de fuerzas o de pendientes (como el de la figura, referida a la ecuación logística). y´ = 5*y*(1-y)
y´ = 5*y*(1-y)
1,2
1,2
1
1
0,8
0,8
En el lado derecho se han dibujado, en el mismo campo de fuerzas, las dos soluciones que corresponden a y.0/ D 0;2 y y.0/ D 1;26.
0,6
0,6 y
y 0,4
0,4
0,2
0,2
0
0
-0,2
-0,2
-0,4
-0,4 0
0,2
0,4
0,6
0,8
1
1,2 t
1,4
1,6
1,8
2
2,2
0
0,2
0,4
0,6
0,8
1
1,2 t
1,4
1,6
1,8
2
2,2
6/57
Tipos de ecuaciones diferenciales: Ecuaciones diferenciales ordinarias —EDO— De primer orden: una variable dependiente y una variable independiente. Problemas de valor inicial y problemas de contorno. Sistemas de ecuaciones diferenciales N variables dependientes relacionadas por N ecuaciones diferenciales. Ecuaciones en derivadas parciales Una variable dependiente y varias variables independientes. Problemas elípticos, parabólicos e hiperbólicos.
Integración numérica de EDO u ODE
7/57
Analíticamente sólo se pueden integrar algunas ecuaciones diferenciales.
La integración numérica no obtiene una función analítica y.t / que satisface la ecuación diferencial anterior, sino un conjunto de valores discretos yk que se corresponden con otros tantos, tk , de la variable independiente t.
Si se desean resultados para otros valores de t se pueden utilizar los métodos de interpolación vistos previamente.
Conocida la función yi en el instante ti los integradores numéricos calculan el valor yi C1 en ti C1 D ti C h mediante una expresión del tipo yi C1 D yi C h ; donde es un término que incorpora información de la derivada y 0.
8/57
Los métodos de integración numérica que veremos se clasifican en: De paso simple y paso múltiple Los de paso simple calculan la solución yi C1 en el instante ti C1 a partir del valor de la función yi en el instante ti . Por ejemplo, en el método de Euler yi C1 D yi C f .ti ; yi /h: Los de paso múltiple calculan la solución en tiC1 a partir del valor de la
función en ti , ti 1; : : : ; ti yiC1
nC1 .
Por ejemplo, en el método de Heun h D yi C f .ti ; yi / C f ti C1; yi C hf .ti ; yi / : 2
9/57
Métodos explícitos e implícitos Los métodos explícitos permiten hallar yi C1 directamente sin tener que resolver un sistema de ecuaciones no lineales. Por ejemplo, Runge-Kutta de orden dos yiC1 D yi C .a1k1 C a2k2/ h; donde k1 D f .ti ; yi / k2 D f .ti C p1h; yi C q11k1h/: Los métodos implícitos necesitan resolver un sistema de ecuaciones no
lineales pues yi C1 aparece a ambos lados de la ecuación. Por ejemplo, la regla trapezoidal yi C1
1 D yi C f .ti ; yi / C f .ti C1; yiC1/ h: 2
10/57
Método de Euler
Formulado por Leonhard Euler, Basilea 1707–San Petersburgo 1783.
Sigue la idea intuitiva de la geometría de los gráficos de los campos de fuerzas o pendientes. ¿Por qué no “resolver” una ecuación diferencial “siguiendo” las flechas, partiendo de un punto inicial?
Desde el punto de arranque hay que moverse una pequeña distancia, siguiendo una pendiente, llegándose a un nuevo punto .t1; y1/ en el que se revalúa la pendiente, luego hay que moverse con esa nueva pendiente otro poco, etc.
La aproximación será buena si la pendiente no cambia mucho súbitamente.
Poe ejemplo se trata de resolver
11/57
8 0 3 ˆ ˆ < y D ty C t
y.0/ D y
0 ˆ ˆ : t 2 Œ0; 1;
que discurre en el campo de “fuerzas” de la figura, partiendo de y.0/ D 1. y´ = t*y+t³ 2
1,8
1,6
1,4
1,2
y
1
0,8
0,6
0,4
0,2
0 0
0,1
0,2
0,3
0,4
0,5
0,6
t
0,7
0,8
0,9
1
12/57
La relación general de recurrencia del método de Euler la define yiC1 D yi C f .ti ; yi /h en la que f .ti ; yi / es la pendiente en cada punto –su derivada– y h el paso que se da con esa inclinación.
Si usamos Euler en el ejemplo, con h D 0;1, resulta la línea en rojo. y´ = t*y+t³ 2
1,8
1,6
1,4
1,2
y
1
0,8
0,6
0,4
0,2
0 0
0,1
0,2
0,3
0,4
0,5
0,6
t
0,7
0,8
0,9
1
Es muy fácil de aplicar aunque no muy preciso, pues la derivada en (ti, yi) no es suficientemente representativa del paso h. LaInterpretación interpretación geométrica de lo que sucede es esta: geométrica:
13/57
HUURU
YDORUH[DFWR YDORUDSUR[LPDGR
yi+1= yi+f(ti ,yi)h h
yi ti
ti+1
Se producen errores de truncamiento en cada punto, por retener sólo dos términos en el desarrollo de Taylor, y en el total, así como de redondeo en los cálculos del proceso. El global es O.h/.
14/57
Si ensayamos todo esto en un script de Matlab, puede resultar este: % Script del método de Euler con Matlab, Euler_2_S.m % y’(t)=t*y+t^3; solución analítica exacta y(t)=3*exp^(t^2/2)-t^2-2 h = 0.1; h1 = 0.01 t = 0:h:1; t1 = 0:h1:1; yx(1)=1.0; yx1(1)=1.0
% Pasos de integración h=0,1 y h=0,01
for i=1:length(t)-1 k1 = t(i)*yx(i)+t(i)^3; yx(i+1) = yx(i)+h*k1; end for i=1:length(t1)-1 k11 = t1(i)*yx1(i)+t1(i)^3; yx1(i+1) = yx1(i)+h1*k11; end
% Bucle del proceso de integración
% Punto de partida
y=3*exp(t.^2/2)-t.^2-2;
% Nuevo punto: valor de y % Bucle del proceso de integración % Nuevo punto: valor de y
% Solución exacta
% Gráficos plot(t,yx,’b--’,t,y,’r-’,t1,yx1,’k--’); legend(’Aproximación h=0,1’,’Exacto y’,’aproximación h=0,01’); title(’Paso de Euler, h=0,1 y h=0,01’); xlabel(’t’); ylabel(’y*(t), y(t)’); axis([0 1 0.95 2]); % Tabla for i=1:length(t) disp(sprintf(’t=%5.3f, end
y(t)=%6.4f,
y*(t)=%6.4f’,t(i),y(i),yx(i)));
15/57
Con pasos h D 0;1 y h D 0;01 se obtiene lo que sigue. Paso de Euler, h=0,1 y h=0,01 2 Aproximación h=0,1 Exacto y aproximación h=0,01
1.9 1.8 1.7
y*(t)=1.0000 y*(t)=1.0000 y*(t)=1.0101 y*(t)=1.0311 y*(t)=1.0647 y*(t)=1.1137 y*(t)=1.1819 y*(t)=1.2744 y*(t)=1.3979 y*(t)=1.5610 y*(t)=1.7744
1.6 y*(t), y(t)
>> Euler_2_S t=0.000, y(t)=1.0000, t=0.100, y(t)=1.0050, t=0.200, y(t)=1.0206, t=0.300, y(t)=1.0481, t=0.400, y(t)=1.0899, t=0.500, y(t)=1.1494, t=0.600, y(t)=1.2317, t=0.700, y(t)=1.3429, t=0.800, y(t)=1.4914, t=0.900, y(t)=1.6879, t=1.000, y(t)=1.9462,
1.5 1.4 1.3 1.2 1.1 1 0
0.1
0.2
0.3
0.4
0.5 t
0.6
0.7
Con el paso de integración h D 0;1 los errores son apreciables.
0.8
0.9
1
Mejoras del método de Euler
16/57
Método de Heun
Formulado por Karl Heun, Alemania, 1859-1929.
Su idea es utilizar en .ti ; yi / un valor promedio de la derivada en ese punto y en .ti C1; yiC1/.
Como no se conoce el valor yi C1, se calcula con el método de Euler, yiC1 D yi C hf .ti ; yi /.
El nuevo valor de yi C1 será yi C1 D yi C
h f .ti ; yi / C f .tiC1; yiC1/ : 2
17/57
El método se conoce como predictor-corrector. El error es O.h3/.
Una variante, denominada método del trapezoide, consiste en hacer iterar el yi C1, en cada nuevo paso, de acuerdo con la fórmula j
yi C1 D yi C
h 2
j 1
f .ti ; yi / C f ti C1 ; yiC1 ˇ ˇ ˇ j j 1ˇ hasta que se cumpla una tolerancia de ˇyi C1 yi C1 ˇ.
Ejemplo Integremos la ecuación y 0 D 4e 0;8t 0;5y en el intervalo t D Œ0; 4 de tal forma que en t D 0, y D 2. El paso ha de ser h D 1. Lo haremos con el método de Euler y con el de Heun, sin y con iteración (trapezoide). En Matlab:
18/57
function [t,yx,y,y_h,y_hi] = Heun_4(dydt,tspan,y0,h,es,maxit) % Método de Heun para integrar una ODE % dydt: función que evalúa la ODE; tspan: intervalos de cálculo % y0: punto de partida; h: paso; es: tolerancia de iteración corrector % maxit: máximo de iteraciones; yx: valor exacto; y: valor con Euler % y_h: valor con Heun sin itera.; y_hi: con iteración if nargin 0;
y0 D y.0/;
cuya solución analítica es y.t / D y0e
˛t
C t 2:
Si ˛ es un número real grande, la solución varía rápidamente hasta que el componente exponencial se desvanece o amortigua.
A partir de entonces prevalece el componente polinómico de variación lenta.
En la figura se observa lo apuntado. La línea continua azul es la solución de la ecuación, la discontinua verde –el componente exponencial– es la que se desvanece antes. La roja corresponde al componente polinómico de variación en el tiempo mucho más lenta.
40/57
1.4
1
0.9 1.2 0.8 1
0.7
0.6
y
y
0.8
0.5
0.6 0.4
0.3
0.4
0.2 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5 t
0.6
0.7
0.8
0.9
1
0
0
0.1
0.2
0.3
0.4
0.5 t
0.6
0.7
0.8
0.9
Figure 2.2.2: Euler solution (dash-dot) due to a perturbation (dashed) of = introduced at time = 0 4. y
t
1
2
t
(solid)
:
Stiness not only depends on the dierential equation under consideration, but the in-
Si se introduce una pequeña perturbación t 2of absolute en stability t Dof a numerical 0;4, la terval ofdel interest,componente the accuracy criteria, and the region method. For nonlinear systems solución exacta se restablece enseguida. En cambio, ely método de Euler = f ( y) produciría, con un h D 2, una perturbación se iría stiness will importante be related to the magnitudes que of the eigenvalues of the Jacobian f ( y). These eigenvalues may have vastly dierent sizes that vary as a function of . Thus, restableciendo después. detecting stiness can be a signicant problem in itself. 0
t
y
t
t
In order to solve sti problems we'll need a numerical method with a weaker stability
El método de Euler hacía atrás
41/57
Actúa de la siguiente manera: y0 D y.0/ yi C1 D yi C hf .ti C1; yiC1/ : La segunda expresión se obtiene de utilizar la aproximación atrasada de la derivada en ti C1, es decir y 0.tiC1/, en vez de la avanzada en ti como el método de Euler.
La información que se considera en cada paso tiene en cuenta los dos lados del intervalo Œti ; ti C h.
Euler hacia atrás es implícito, lo que significa que no da directamente una fórmula para la nueva aproximación yi C1 sino que ésta se obtiene mediante un proceso más elaborado.
42/57
Ese proceso, para la ecuación diferencial anterior, sería 2 yiC1 D yi C h ˛ yi C1 t C 2t D yi h˛yiC1 C h ˛t 2 C 2t : Al final
2
yi C1 D
yi C h ˛t C 2t : 1 C ˛h
Si a esta fórmula se le aplica un método iterativo de punto fijo –el de Newton por ejemplo– y C h ˛t 2 C 2t y ! g.y/ D ; 1 C ˛h como g 0.y/ D 1=.1 C ˛h/ < 1, para cualquier ˛ positivo, el proceso convergerá a un punto estacionario.
Apliquemos estas ideas en un código de Matlab. function [wi, win, y, ti] = back_euler_1(t0,x0,tf,N,TOL,Nmax) % % Backward Euler con y’=-alpha(y-t^2)+2t % Punto de arranque, t0, x0; final, tf; puntos N if nargin> Adams_Bas_Moul([0 2],0.5,20); es este. 1, Adams-Bash.-Moul. ode45
0,9 0,8 0,7 0,6
y
0,5 0,4 0,3 0,2 0,1 0, 0,
0,5
1,
1,5 Tiempo
2,
2,5
54/57
Sistemas de ecuaciones diferenciales ordinarias
En ingeniería, economía y ciencias sociales muchos sistemas se modelizan o simulan mediante sistemas de ecuaciones diferenciales. Basta para ello que el sistema tenga más de un grado de libertad o que haya que integrar una ecuación diferencial de orden superior a 1.
Ejemplos son los puntos moviéndose en el plano o en el espacio, el movimiento de un sólido rígido sin restricciones, mecanismos, vehículos, etc.
55/57
En general, un sistema de ecuaciones diferenciales ordinarias tiene la forma 9 0 y1 D f1 t; y1; y2; : : : ; yn > > > = 0 y2 D f2 t; y1; y2; : : : ; yn y 0 D f .t; y/: ::: > > > 0 yn D fn t; y1; y2; : : : ; yn ;
Al aplicar los métodos vistos hasta ahora hay que tener en cuenta que la variable dependiente y, la función f —y los coeficientes ki si se usa Runge-Kutta— son vectores.
El método de Euler se aplicaría en la forma vectorial y 0 D f .t; y/
!
y i C1 D y i C hf .ti ; y i /:
56/57
Ecuaciones diferenciales de orden superior
Una ecuación diferencial de orden n puede tener esta expresión d ny .n/ 0 00 000 .n D y D f t; y; y ; y ; y ; : : : y dt n
1/
y.t0/ D y0 y 0.t0/ D y00 y 00.t0/ D y000 .n y .n 1/.t0/ D y0
1/
Esta ecuación se puede transformar en un sistema de n ecuaciones diferenciales de orden uno.
57/57
Si se introducen las nuevas variables y1 D y; y2 D y 0; y3 D y 00; : : : ; yn D y .n
1/
se tendrá que y10 y20 y30
D D D ::: yn0 1 D yn0 D
y2 y3 y4
9 > > > > > > > > =
y 0 D f .t; y/;
> > > > > yn > > > f .x; y1; y2; : : : ; yn/ ;
obteniéndose así un sistema de n ecuaciones diferenciales de primer orden con n incógnitas, al que habría que añadirle las condiciones iniciales y.t0/ D y 0.
