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 con
Condiciones de Contorno o Frontera
José Luis de la Fuente O’Connor
[email protected] [email protected] Clase_integración_problemas_contorno_2016.pdf
1/34
2/34
Índice
Introducción
El método del disparo
El método de las diferencias finitas
Métodos de colocación y de elementos finitos
Introducción
3/34
En este tema analizamos los problemas que se pueden modelizar mediante una ecuación diferencial de segundo orden con condiciones en dos puntos dados, extremos de un intervalo de estudio Œa; b.
Son los denominados problemas de ecuaciones diferenciales con valores en la frontera, o problemas de contorno. Se formulan así: y 00 D f .t; y; y 0/;
y.a/ D ˛;
y.b/ D ˇ:
Para definir la solución se necesitan dos condiciones de contorno:
Si se dan para un mismo punto, se tiene un problema de valor inicial como los que hemos estudiado.
Si se dan para dos puntos, se tiene un problema de valores en la frontera.
El método del disparo
4/34
Se basa en encontrar la ecuación diferencial que tiene la misma solución del de contorno que se estudia. ChE-401: Computational Methods in Chem. Eng.
S
Shooting method
Con tal fin se genera una sucesión de problemas de valor inicial que converja al dado.
1
2
3
En cada uno se estima y 0.a/, la derivada –pendiente– de y en el extremo izquierdo del intervalo. Luego se resuelve (se dispara) hacia y.b/ y se obtiene el valor real yb en b. Chemical Engineering Dep.
Prof. Ibrahim S. Al-Mutaz
King Saud University
٢٥
Chemical
ChE-401: Computational Methods in Chem. Eng.
Se itera hasta que ySolution cerca deProblems y.b/- Shooting como method sea necesario. Boundary-Value b estéoftan
S
5/34
De forma más estructurada, si se define la función
F .s/ D
8 ˆ ˆ < Diferencia entre yb e y.b/; donde
y.t / es la solución del problema de valor inicial,
ˆ ˆ : con y.a/ D y y y 0 .a/ D s; a
el problema con valores de frontera se reduce a resolver la ecuación F .s/ D 0:
Para resolver esta ecuación se puede utilizar el método de la bisección, por ejemplo, informando al programa correspondiente de qué dos valores s0 y s1 hacen que F .s0/F .s1/ < 0.
6/34
Ejemplo Resolvamos mediante el método del disparo este problema: 8 00 ˆ ˆ < y D 4y
y.0/ D 1
ˆ ˆ : y.1/ D 3:
Escribamos la ecuación diferencial como un sistema de dos ecuaciones de orden uno: y0 D v v 0 D 4y:
En Matlab podría escribirse así:
function z=F(s) a=0; b=1; yb=3; ydot = @(t,y) [y(2); 4*y(1)]; [t,y]= ode45(ydot,[a,b],[1,s]); % Solver de Matlab z = y(end,1)-yb end
Comprobamos que para s D 1, F . 1/ 1;0512 y que F .0/ 0;7622 por lo que podemos invocar el programa Bisec_y.m en el intervalo Œ 1; 0: >> [ss,t,y]=Bisec_y(@F,[-1,0]); >> ss ss = -0.4203 >> plot(t,y) >> xlabel(’Tiempo’); ylabel(’y’) >> legend(’y(t)’, ’y’’(t)’,’Location’,’NorthWest’)
El resultado de la evolución en el tiempo de y e y 0 es este. 6 y(t) y'(t)
5
4
3
y
2
1
0
-1 0
0.1
0.2
0.3
0.4
0.5
Tiempo
0.6
0.7
0.8
0.9
1
7/34
Método de las diferencias finitas
Su intención es reemplazar las derivadas de la ecuación diferencial por su aproximación por diferencias fintas en una malla de puntos del intervalo de estudio.
Se consigue un sistema de ecuaciones más o menos grande según lo “tupida” que sea esa malla. Luego se resuelve.
Si y.t/ es una función con derivadas continuas al menos hasta cuarto orden, sabemos que 0
y .t / D
y.t C h/
y.t
h/
2h
h2 000 y ./ 6
y que 00
y .t / D
y.t C h/
2y.t/ C y.t h2
h/
h2 000 C y ./: 12
8/34
for approximations wi to the correct values yi , as shown in Figure 7.6. The boundary conditions are substituted in the system of equations where they are needed. y 1
w1 w2 ya
y1 y2
wn–1 wn yn–1
t0
t1
t 2 ...
yn
tn–1 tn
yb tn+1
t
Figure 7.6 The Finite Difference Method for BVPs. Approximations wi , i = 1, . . . , n for the correct values yi at discrete points ti are calculated by solving a linear system of equations. objetivo es reemplazar los y , y ; : : : ; y reales de la trayectoria por sus
El 1 2 n aproximaciones por diferencias finitas, wi . situations. If the original boundary value After the substitutions, there are two possible
problem was linear, then the resulting system of equations is linear and can be solved by Gaussian elimination or iterative methods. If the original problem was nonlinear, then Si el de contorno original es linealequations, el problema se convierte entonces theproblema algebraic system is a system of nonlinear requiring more sophisticated begin with linear example. en approaches. resolver unWesistema de aecuaciones lineales; si no es lineal, en lo propio.
E 7.8 Solve the BVP (7.7)
9/34
10/34
Ejemplo Resolvamos mediante diferencias finitas este problema: 8 00 ˆ ˆ < y D 4y
y.0/ D 1
ˆ ˆ : y.1/ D 3:
La aproximación de la ecuación diferencial y 00 D 4y por diferencias finitas en ti es wiC1 2wi C wi 1 4wi D 0; h2 2 o, wi 1 C 4h 2 wi C wi C1 D 0:
Si se divide el intervalo de estudio en tres tramos, n D 3, el parámetro h D 1=.n C 1/ D 1=4, por lo que habrá tres ecuaciones.
Introduciendo que w0 D 1 y que w4 D 3 se llega al sistema 4h 2 w1 C w2 D 0 4h2 2 w2 C w3 D 0 4h2 2 w3 C 3 D 0:
1C w1 C w2 C
Sustituyendo el valor de h a este sistema con matriz tridiagonal 2 9 3" # 2 3 1 0 1 4 w1 4 1 49 15 w2 D 4 05 : 0
11/34
2
1
9 4
w3
3
Resolviendo este sistema se obtiene la solución de esta tabla. i 0 1 2 3 4
ti 0;00 0;25 0;50 0;75 1;05
wi 1;0000 1;0249 1;3061 1;9138 3;0000
yi 1;0000 1;0181 1;2961 1;9049 3;0000
Las diferencias son O.10 2/. Para obtener más precisión habría que aumentar n.
12/34
En general, si h D .b a/=.n C 1/ D 1.n C 1/ el sistema de matriz tridiagonal que se obtiene en este ejemplo es de la forma 3 2 2 2 13 4h 2 1 0 0 0 0 2 3 4h2 2 1 0 0 0 7 w1 6 1 2 7 6 w2 7 6 007 6 0 1 4h 2 0 0 0 7 6 w3 7 6 :7 6 :: :: ::: ::: ::: 7 6 : 7 D 6 ::7 : 6 : 0 : 7 4 :: 5 6 7 6 ::: ::: 0 0 1 0 7 wn 1 6 0 4 05 5 4 0 ::: 0 0 0 4h2 2 1 0
0
0
0
1
4h
2
wn
2
3
Los errores en este método vienen de la aproximación de las derivadas y de la resolución del sistema de ecuaciones.
Para valores de h mayores que la raíz cuadrada de la precisión de la máquina, , prevalecen los de la aproximación, O.h2/.
Ejemplo Resolvamos ahora mediante diferencias finitas este otro problema: 8 00 ˆ ˆ sqrt(eps) % Iteraciones de Newton-Raphson w=w1; w1=w-jac(w,h,n)\f7(w,h,bv,n); end close all, plot ([a a+(1:n)*h b],[ya w1’ yb]); end function y=f7(w,h,bv,n) y=zeros(n,1); y(1)=bv(1)-(2+h^2)*w(1)+h^2*w(1)^2+w(2); y(n)=w(n-1)-(2+h^2)*w(n)+h^2*w(n)^2+bv(2); for i=2:n-1 y(i)=w(i-1)-(2+h^2)*w(i)+h^2*w(i)^2+w(i+1); end end
4.5
4
3.5
3
2.5
2
function a=jac(w,h,n) 1.5 a=zeros(n,n); for i=1:n a(i,i)=2*h^2*w(i)-2-h^2; 1 end 0 for i=1:n-1 a(i,i+1)=1; a(i+1,i)=1; end end
0.1
0.2
0.3
0.4
0.5
0.6
La respuesta se consigue con >> w=nlbvpfd([0 1],[1 4],40);
0.7
0.8
0.9
1
16/34
Métodos de colocación y de elementos finitos
La idea que comparten estos dos métodos es la de convertir el problema en el de resolver un conjunto de ecuaciones algebraicas mediante la aproximación de la ecuación diferencial por y.t / D c1 1.t / C C cn n.t /; donde las i son funciones base definidas en el intervalo Œa; b y los ci parámetros que hay determinar.
Las funciones base pueden ser polinomios, funciones trigonométricas, splines u otras funciones simples.
La determinación de una solución aproximada de la ecuación diferencial se reduce a encontrar los valores de los ci .
17/34
El enfoque de la colocación consiste en usar la aproximación de y.t / en un conjunto de n puntos a D t1 < < tn D b, denominados puntos de colocación, en los que se fuerza a que la solución aproximada y.t / satisfaga la ecuación diferencial.
Esto lleva a un sistema de ecuaciones lineales, o no lineales, en los que se calculan los ci . La idea es parecida a la de la interpolación.
El método de los elementos finitos aproxima mediante mínimos cuadrados los ci de tal forma que se minimice la suma al cuadrado de las desviaciones entre los valores exactos en unos puntos de la ecuación diferencial y los aproximados.
18/34
Método de colocación
Resolvamos la ecuación
8 00 ˆ ˆ < y D 4y
y.0/ D 1
ˆ ˆ : y.1/ D 3:
Se eligen como funciones base monomios j .t / D t j 1, para 1 j n. La solución será de la forma y.t / D
Pn
j D1 cj j .t /
D
Pn
j D1 cj t
j 1
.
Se formularán n ecuaciones con las n incógnitas c1; : : : ; cn. La primera y la última son las condiciones en la frontera: c1 D
n X
cj j .0/ D y.0/ D 1
j D1
c1 C C cn D
n X
cj j .1/ D y.1/ D 3:
j D1
Las otras n 2 ecuaciones resultan de evaluar la ecuación diferencial P y 00 D f .t; y; y 0/ en ti , 2 i n 1, siendo y.t / D jnD1 cj t j 1. Es decir 1 0 n n n X X X j 3 j 1 @ j 1 j 2 cj t D f t; cj t ; cj j 1 t j 2A : j D1
En el ejemplo, n h X j D1
j D1
j
1 j
j 2 ti
3
j 1 4ti
i
cj D 0;
j D1
i D 2; : : : ; n
1:
19/34
Las n ecuaciones forman un sistema Ac D b, donde la matriz de coeficientes A está definida así, 8 ˆ 1 0 0 0 ˆ ˆ ˆ < Aij D ˆ j 1 j 2 tij ˆ ˆ ˆ :1 1 1 1
fila i D 1 3
j 1
filas i D 2; : : : ; n
4ti
1
fila i D n
y b D Œ1 0 0 0 3T .
Los ti de la malla que se usan son los uniformemente espaciados ti D a C
i n
1 .b 1
a/ D (en el caso del ejemplo a)
i n
Para n D 2 el sistema queda
1 0 1 1
c1 c2
1
D 3 ;
y la solución es c D Œ1 2T . Es decir y.t / D c1 C c2t D 1 C 2t.
1 : 1
20/34
Para n D 4 el sistema es 2 4
1 4 4 1
0 4 1=3 2 4 2=3 2 1
21/34
32 3
0 0 c1 2 4 .1=3/ 6 1=3 4 .1=3/3 5 4c2 5 c3 4 .2=3/2 6 2=3 4 .2=3/3 c4 1 1
cuya solución es c D Œ1;0000 y.t / 1
23 D
1 405 ; 0 3
0;1886 1;0273 1;1613T . Es decir, 0;1886t C 1;0273t 2 C 1;161t 3 :
3 Exacta y(t)=1+2t 2
3
y(t)=1-0.1886t+1.0273t +1.1613t 2.5
El gráfico-resultado para el ejemplo:
y
2
1.5
1
0.5
0
0.1
0.2
0.3
0.4
0.5 Tiempo
0.6
0.7
0.8
0.9
1
22/34
El ejemplo que hemos utilizado para introducir el método era lineal, pues así lo era la ecuación diferencial. Si fuese no lineal, habría que utilizar el método de Newton-Raphson para resolver el sistema de ecuaciones no lineales que resultase.
Igual que hemos utilizado monomios como funciones base para aplicar el método, se puede recurrir a cualesquiera otras: polinomios de Chebyshev, funciones trigonométricas, etc.
23/34
Método de los elementos finitos de Galerkin
El método de los elementos finitos, formulado por Boris Grigoryevich Galerkin, Rusia 1871-1945, minimiza el cuadrado del error de la aproximación de la ecuación diferencial a lo largo de la solución. Aboca a un sistema de ecuaciones lineales diferente del anterior.
Si en el procedimiento general anterior se usan como funciones base unos splines nos referimos al método de los elementos finitos. Las splines son distintas de cero en un corto espacio de t .
En el caso de la colocación, los ci se obtenían forzando a que la solución aproximada satisficiese exactamente las condiciones de contorno y los puntos de colocación.
El método de los elementos finitos para resolver
24/34
8 00 0 ˆ ˆ < y D f .t; y; y /
y.a/ D y
a ˆ ˆ : y.b/ D y ; b
escoge una solución aproximada y de tal forma que el residuo r D y 00 lo más pequeño posible.
f sea
Análogamente a como vimos para mínimos cuadrados, esto se consigue eligiendo y de tal forma que el vector de residuos sea ortogonal al espacio vectorial del conjunto de las soluciones de la ecuación diferencial.
Para un intervalo Œa; b se define el espacio vectorial de las funciones integrables cuadradas ( ) ˇZ b ˇ L2Œa; bD funciones y.t / en Œa; bˇˇ y.t /2dt existe y es finita : a
25/34
En este espacio se define el producto interior Z b y1.t /y2.t /dt; hy1; y2i D a
con las propiedades habituales: 1. hy1; y2i 0; 2. h˛y1 C ˇy2; zi D ˛hy1; zi C ˇhy2; zi, siendo ˛ y ˇ escalares; 3. hy1; y2i D hy2; y1i. Dos funciones y1 e y2 son ortogonales en L2Œa; b si hy1; y2i D 0.
Como este espacio es de dimensión infinita, no podemos hacer r D y 00 f ortogonal al mismo en sentido estricto, pero si escoger una base que comprenda lo más posible del mismo en el sentido finito de lo cálculos.
Hay dos ideas principales detrás del método de Galerkin:
26/34
1º- Minimizar el vector residuo r forzando a que sea ortogonal, en el sentido del producto interior anterior, a un conjunto de base 0.t /, R b n C00 2 funciones 1.t /; : : :, nC1.t /. Es decir forzar a que a y f i dt D 0, o que b
Z
b
Z
00
f .t; y; y 0 / i .t/ dt
y .t/ i .t/ dt D
forma débil
a
a
para cada 0 i n C 1. 2º- Utilizar integración por partes para eliminar las segundas derivadas. Teniendo en cuenta que b
Z
ˇb y .t/ i .t/ dt D i .t/y 0 .t/ˇa 00
a
D i .b/y 0 .b/
Z
b
y 0 .t/ i0 .t/ dt Z b a i .a/y 0 .a/ y 0 .t/ i0 .t/ dt: a
Utilizando la forma débil y esta expresión resulta un conjunto de ecuaciones l
a
b
f .t; y; y 0 / i .t/ dt D i .b/y 0 .b/
i .a/y 0 .a/
l
a
b
y 0 .t/ i0 .t/ dt:
27/34
Resolviéndolas se determinan los ci de y.t / D
nC1 X
ci i .t /:
i D0
Esas dos ideas hacen conveniente usar funciones base i .t / —elementos finitos— muy simples, como lo son los B-splines.
Para ello, en los puntos t0 < t1 < < tn < tnC1 sobre el eje t, para i D 1; : : : ; n, se define
i .t / D
8 t ˆ ˆ ˆ t ˆ < i
ti 1 ti 1
ti C1 t tiC1 ti
ˆ ˆ ˆ ˆ :0
para ti
1
< t ti
para ti < t < ti C1 en otros casos.
28/34
Y, también, 0.t/ D
8 <
t1 t ti t0
:0
para t0 t < t1 en otros casos7.3
y nC1.t / D
8 <
t tn tnC1 tn
para tn < t tnC1
0 ElementenMethod otros casos. Collocation and the Finite | 369 :
y 1
0 1 2
t0
t1
t2
n–1 n n+1
t 3 ...
tn–1 tn
tn+1
t
Figure 7.10 Piecewise-linear B-splines used as finite elements. Each φi (t), for 1 ≤ i ≤ n, has support on the interval from ti−1 to ti+1 .
Estas funciones satisfacen que i .tj / D
1 si i D j 0 si i ¤ j:
for each i that can be solved for the ci in the functional form y(t) =
n+1
ci φi (t).
(7.23)
29/34
Para un conjunto de datos .ti ; ci /, se define la B-spline lineal por tramos S.t / D
nC1 X
ci i .t /:
i D0
Se puede comprobar que S.tj / D
PnC1 iD0
ci i .tj / D cj .
S.t/ es una función lineal por tramos que interpola los puntos .ti ; ci /. Además, los ci son las soluciones en los puntos ti de la aproximación de la ecuación diferencial.
30/34
Para calcular los ci : el primero y el último se obtienen por colocación: y.a/ D
nC1 X
ci i .a/ D c0 0.a/ D c0
i D0
y.b/ D
nC1 X
ci i .b/ D cnC1 nC1.b/ D cnC1:
i D0
Los demás, usando las ecuaciones de los elementos finitos: Z b Z b f .t; y; y 0/ i .t / dt C y 0.t / i0 .t / dt D 0; a
a
P
sustituyendo y.t / por ci i .t /, Z b Z X X i .t /f t; cj j .t /; cj j0 .t / dt C a
b a
i0 .t /
X
cj j0 .t / dt D 0:
31/34
Suponiendo que los puntos en t están uniformemente espaciados un paso h, hay que calcular todas estas integrales para cada i D 1; : : : ; n: b
Z
b
Z i .t/ i C1 .t/ dt D
a
a
D b
Z
t2 2h
Z b t t t 1 dt D h h h a ˇ 3 ˇh ˇt ˇ ˇ 3h2 ˇ D h6 h
t h
h
Z
2
a
a
0
b 0 .t/ dt i0 .t/ iC1
Z
b
0 2
i .t/ a
dt
0
i .t/ dt D 2 Z
t2 h2
0
1 h
Z
h
Z D
dt D 2 0
2
1 h
1 h 2
2 dt D h 3 dt D
1 h
2 dt D : h
Si la ecuación diferencial es lineal, las ecuaciones de los ci también.
32/34
Ejemplo Mediante elementos finitos, resolvamos
8 00 ˆ ˆ < y D 4y
y.0/ D 1
ˆ ˆ : y.1/ D 3:
Sustituyendo la ecuación diferencial en las expresiones de los ci se tiene para cada i l 1 ! nC1 nC1 0 D
4 i .t/
0
D
nC1 i j D0
X
j D0
cj j .t/ C
X
cj j0 .t/ i0 .t/ dt
j D0
1
Z cj 4
1
Z i .t/ j .t/ dt C
0
0
j0 .t/ i0 .t/ dt
:
Usando las expresiones de las B-splines para i D 1; : : : ; n y teniendo en cuenta que c0 D f .a/ y cnC1 D f .b/, las ecuaciones que resultan son estas:
33/34
2
1 h
c0 C
8
2 h
c1 C
2
1 h
c2 D 0
2
1 h
c1 C
8
2 h
c1 C
2
1 h
c1 D 0
h 3 h 3
hC 3 hC 3
h 3 h 3
:: : 0 1 h
2
h 3
cn
1
C
8
hC 3
2 h
cn C
2
h 3
1 h
cnC1 D 0:
Como c0 D ya D 1 y cnC1 D yb D 3 el sistema resultante es este 2 32 3 2 3 c1 ya ˇ 0 0 : : : : : : ::: 7 6 7 6 7 7 6 c2 7 6 0 7 6 7 6 7 ˇ : : : ˇ 07 7 6 ::: 7 D 6 ::: 7 ; 6 7 6 7 : : : : : : ˛ ˇ7 5 4cn 15 4 0 5
˛ ˇ 6ˇ ˛
6 60 6 6: 4 ::
0
donde
0
ˇ
˛
cn
2 2 8 y ˇD h ˛ D hC 3 h 3
yb ˇ
1 : h
Si recopilamos todo lo necesario en un fichero .m de Matlab para este caso:
34/34
function c=Galerkin_ef_1(inter,bv,n) % Elementos finitos de ejemplo clase % inter, intervalo de integración; bv, valores de contorno; % n, número de pasos. a=inter(1); b=inter(2); ya=bv(1); yb=bv(2); h=(b-a)/(n+1); alfa=(8/3)*h+2/h; beta=(2/3)*h-1/h; e=ones(n,1); M=spdiags([beta*e alfa*e beta*e],-1:1,n,n); d=zeros(n,1); d(1)=-ya*beta; d(n)=-yb*beta; c=M\d; end
Resolviendo con la orden c=Galerkin_ef_1([0 1],[1 3],3) se llega a la solución de la tabla que sigue. i 0 1 2 3 4
ti 0;00 0;25 0;50 0;75 1;05
ci 1;0000 1;0109 1;2855 1;8955 3;0000
yi 1;0000 1;0181 1;2961 1;9049 3;0000
Las diferencias son O.10 2/. Para obtener más precisión habría que aumentar n.