Diferenciasfinitas_13a20.pdf

Ecuaciones Diferenciales Ordinarias III ´ • Problemas de valores de contorno (P.V.C.). El metodo de ”shooting”. ´ • El metodo de diferencias finitas. ´ • El metodo de elementos finitos.

521230

-1-

´ DIM – Universidad de Concepcion

Problemas de valores de contorno (P.V.C.) ´ El metodo de “shooting” Consideremos el problema de valores de contorno (P.V.C.)

  y 00 (x) = f (x, y, y 0 ), x ∈ (a, b),  y(a) = α, y(b) = β. ´ ´ Para resolver este problema, el primer metodo que veremos es el metodo de “shooting” (o del disparo): Dado z

∈ R consideremos el P.V.I. asociado   y 00 (x) = f (x, y, y 0 ), x ∈ (a, b),  y(a) = α, y 0 (a) = z.

´ y denotemos por yz su solucion. 521230

-2-

´ DIM – Universidad de Concepcion

´ La idea del metodo de “shooting” es encontrar z¯ tal que yz¯(b)

= β. yz(b)

´ Para esto, si definimos la funcion

E(z) := yz (b) − β, z ∈ R,

yz (x)

E(z)=yz (b) − β

β

vemos que

E(z) = 0 ⇐⇒ yz (b) = β. Se trata entonces de encontrar que

α

z¯ tal

E(¯ z ) = 0. Es decir, encontrar

´ (en general una ra´ız para la ecuacion no lineal) E(z)

521230

y (x)

a

b

= 0.

-3-

´ DIM – Universidad de Concepcion

Ejemplo. Se lanza una piedra de masa 1 kg verticalmente hacia arriba. Sobre la piedra actuan la fuerza de gravedad y una fuerza de roce igual a 0.1 kg/s por la velocidad de la ´ ´ del tiempo, sabiendo que al cabo de piedra. Determinar la altura de la piedra como funcion 3 s la piedra esta´ exactamente a 50 m del punto de partida.

´ Sea h(t) la altura de la piedra en el instante t ≥ 0. El P.V.C. a resolver es: Solucion.   h ¨ = −g − 0.1h, ˙ t ∈ (0, 3),  h(0) = 0, h(3) = 50, donde g

521230

= 9.8 m/s2 .

-4-

´ DIM – Universidad de Concepcion

´ es la velocidad inicial z para que al cabo de 3 s la altura Para ello debemos averiguar cual de la piedra sea exactamente de 50 m. Esto conduce al P.V.I. asociado:

  h ¨ = −g − 0.1h, ˙ t ∈ [0, 3],  h(0) = 0, h(0) ˙ = z,

´ anal´ıtica puede calcularse en funcion ´ de z , a saber: cuya solucion

hz (t) = −(10z + 100g)e−0.1t + 10z + 100g − 10gt. ´ hz (3) Al resolver la ecuacion

= 50, se obtiene z¯ = 34.7 y sustituyendo este valor en la

´ del tiempo esta´ dada ´ de hz (t) se tiene que la altura de la piedra como funcion expresion por

hz¯(t) = −(347 + 100g)e

521230

−0.1t

+ 347 + 100g − 10gt = 1327 1 − e

-5-

−0.1t




− 98t.

´ DIM – Universidad de Concepcion

´ anal´ıtica de la funcion ´ E(z), pero, dado z En general no se tiene una expresion

∈ R,

´ ´ ´ yz puede usarse cualquiera de los metodos numericos estudiados para calcular la solucion del P.V.I. asociado

y con ella evaluar E(z)

  y 00 (x) = f (x, y, y 0 ), x ∈ (a, b),  y(a) = α, y 0 (a) = z, = yz (b) − β .

´ E(z) Para resolver la ecuacion

´ = 0 puede utilizarse cualquiera de los siguientes metodos:

´ ´ para lo cual deber´ıamos encontrar previamente dos valores • El metodo de biseccion, z0 y z1 tales que E(z0 ) < 0 < E(z1 ).

521230

-6-

´ DIM – Universidad de Concepcion

´ • El metodo de Newton-Raphson: ´ inicial, z0 : aproximacion E(zk ) , zk+1 = zk − 0 E (zk )

k = 1, 2, . . .

´ ´ complejo, pues requiere evaluar Usar el metodo de Newton-Raphson es mas

∂ ∂yz [yz (b) − β] = (b). E (z) = ∂z ∂z 0

´ anal´ıtica de la funcion ´ yz para calcular esa En general no se tiene una expresion derivada, pero derivando impl´ıcitamente respecto de z en el P.V.I. asociado se obtiene que uz

521230

:=

∂yz ´ del siguiente P.V.I.: es la solucion ∂z   u00 = ∂f (x, y, y 0 )u + ∂f (x, y, y 0 )u0 , x ∈ (a, b), ∂y ∂y 0  u(a) = 0, u0 (a) = 1. -7-

´ DIM – Universidad de Concepcion

´ Entonces, si se quiere aplicar el metodo de Newton, en cada paso debe resolverse ´ numericamente el sistema de E.D.O. de segundo orden

 00 0  (x) = f (x, y, y ), x ∈ (a, b), y    ∂f ∂f (x, y, y 0 )u + 0 (x, y, y 0 )u0 , x ∈ (a, b), u00 =  ∂y ∂y    y(a) = α, y 0 (a) = z, u(a) = 0, u0 (a) = 1,

para as´ı evaluar E(z)

´ calculada = yz (b) y E 0 (z) = uz (b), con (yz , uz ) la solucion

de este sistema.

521230

-8-

´ DIM – Universidad de Concepcion

´ • El metodo de la secante:

z0 , z1 : aproximaciones iniciales, E(zk ) zk+1 = zk − , E(zk ) − E(zk−1 ) zk − zk−1

k = 1, 2, . . .

´ ´ El metodo de “shooting” aplicando el metodo de la secante se puede resumir mediante el siguiente algoritmo:

521230

-9-

´ DIM – Universidad de Concepcion

1. Dar un par de estimaciones iniciales para z¯, z0 y z1 , y resolver los P.V.I. asociados con

z = z0 y z = z1 para calcular E(z0 ) y E(z1 ). 2. Para k (a)

zk+1

= 1, 2, . . . E(zk ) = zk − . E(zk ) − E(zk−1 ) zk − zk−1

´ (b) Resolver numericamente el P.V.I. asociado con z

= zk+1 :

  y 00 (x) = f (x, y, y 0 ), x ∈ (a, b),  y(a) = α, y 0 (a) = zk+1 .

´ numerica ´ (c) Utilizar la solucion obtenida yzk+1 para calcular

E(zk+1 ) = yzk+1 (b) − β . ˜ entonces terminar: la solucion ´ yzk+1 es (d) Si |E(zk+1 )| es suficientemente pequeno, la correcta. Si no, volver a (2).

521230

- 10 -

´ DIM – Universidad de Concepcion

Ejemplo. Consideremos un problema como el anterior pero en el que la fuerza de roce es no lineal, por ejemplo:

  h ¨ = −g − 0.01h˙ 2 , t ∈ [0, 3],  h(0) = 0, h(3) = 50. Aplicamos el algoritmo anterior con z0

= 10 y z1 = 20.

Utilizamos el comando ode45 para resolver el P.V.I.

Iteramos el procedimiento hasta que |E(zk+1 )|

521230

- 11 -

< 10−3 .

´ DIM – Universidad de Concepcion

60

Z0 Z1 Z2 Z3 Z4 Z5

50

k

zk

hzk (3) 40

0

10.0000

−16.4545

1

20.0000

10.0826

2

35.0421

40.1375

20

3

39.9782

48.3124

10

4

40.9972

49.9199

0

5

41.0480

49.9993

30

−10

−20

521230

0

- 12 -

0.5

1

1.5

2

2.5

3

´ DIM – Universidad de Concepcion

´ Metodo de Diferencias Finitas Supongamos que queremos resolver el P.V.C.

  −y 00 + ry 0 + ky = f (x),  y(a) = α, y(b) = β,

x ∈ (a, b),

donde k y r son constantes positivas. Supondremos que el problema tiene una unica ´ ´ la cual es al menos cuatro veces diferenciable. solucion,

´ El metodo de diferencias finitas se basa en las siguientes aproximaciones para las derivadas de y :

y(x + h) − y(x − h) , y (x) ≈ 2h 0

y(x + h) − 2y(x) + y(x − h) , y (x) ≈ h2 ´ ˜ las cuales son validas si h es suficientemente pequeno. 00

521230

- 13 -

´ DIM – Universidad de Concepcion

En efecto. Consideremos los siguientes desarrollos de Taylor:

h2 00 h3 000 y(x + h) = y(x) + y (x)h + y (x) + y (x) + O(h4 ), 2! 3! 0

2 3 h h y(x − h) = y(x) − y 0 (x)h + y 00 (x) − y 000 (x) + O(h4 ). 2! 3!

´ Sumandolos se tiene

y(x + h) + y(x − h) = 2y(x) + h2 y 00 (x) + O(h4 ), de donde

y 00 (x) =

521230

y(x + h) − 2y(x) + y(x − h) 2 + O(h ). 2 h

- 14 -

´ DIM – Universidad de Concepcion

Por otra parte, restando los desarrollos de Taylor anteriores,

h2 00 h3 000 y(x + h) = y(x) + y (x)h + y (x) + y (x) + O(h4 ), 2! 3! 0

h2 00 h3 000 y(x − h) = y(x) − y (x)h + y (x) − y (x) + O(h4 ), 2! 3! 0

se tiene

y(x + h) − y(x − h) = 2hy 0 (x) + O(h3 ), y de aqu´ı

y(x + h) − y(x − h) + O(h2 ). y (x) = 2h 0

´ Las aproximaciones de las derivadas que se obtuvieron reciben el nombre Observacion: de diferencias centradas.

521230

- 15 -

´ DIM – Universidad de Concepcion

´ uniforme del intervalo [a, b] en n subintervalos de igual Consideremos una particion longitud:

b−a , h= n

xi = a + ih, i = 0, . . . , n.

a x0

b x1

x i−1

xi

x i+1

x n−1

Entonces, las aproximaciones de y 0 e y 00 en los nodos interiores xi , i

xn

= 1, . . . , n − 1,

pueden escribirse como sigue:

521230

y 0 (xi ) =

y(xi+1 ) − y(xi−1 ) + O(h2 ), 2h

y 00 (xi ) =

y(xi+1 ) − 2y(xi ) + y(xi−1 ) 2 + O(h ). 2 h - 16 -

´ DIM – Universidad de Concepcion

Sustituyendo estas expresiones en el P.V.C.

  −y 00 + ry 0 + ky = f (x),  y(a) = α, y(b) = β,

para cada uno de los nodos interiores xi , i

x ∈ (a, b),

= 1, . . . , n − 1, se tiene:

 y(xi+1 ) − y(xi−1 ) y(xi+1 ) − 2y(xi ) + y(xi−1 )   +r   − h2 2h + ky(xi ) = f (xi ) + τi , i = 1, . . . , n − 1,     y(x ) = α, y(x ) = β, 0

´ donde los terminos τi

n

´ de = O(h2 ) provienen de los errores O(h2 ) de la aproximacion

las derivadas mediante las diferencias centradas.

521230

- 17 -

´ DIM – Universidad de Concepcion

´ ´ Al despreciar los terminos τi , que se denominan errores de truncamiento, el metodo de diferencias finitas permite aproximar los valores y(xi ) mediante valores yi que satisfacen

 y − 2yi + yi−1 yi+1 − yi−1   − i+1 + kyi = f (xi ), i = 1, . . . , n − 1, + r 2 h 2h   y = α, y = β. 0

n

´ numerica ´ De esta forma, la resolucion del P.V.C. por diferencias finitas se reduce a resolver un sistema de ecuaciones lineales.

Multiplicando por h2 , las ecuaciones pueden reescribirse como sigue:

    
rh rh  2 2  y y −1 − + 2 + kh + −1 + = h f (xi ), y  i−1 i i+1  2 2 i = 1, . . . , n − 1,     y = α, y = β. 0 n 521230

- 18 -

´ DIM – Universidad de Concepcion

As´ı, definiendo para i

ai := −1 −

= 1, . . . , n − 1, rh rh , bi := 2 + kh2 , ci := −1 + , fi = f (xi ), 2 2

llegamos a la forma matricial del sistema de ecuaciones:



b1

c1

0

···

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

···

0

0







´ de verificar que si h es suficientemente pequeno ˜ (h ´ Es facil Observacion:

< 2/r),

  a 2   0  .  ..  0

an−1







a1 α f1 y1       ..         f .   y2  0 2          ..    ..  2  ..    .  = h  .  −  . . 0                y f 0 n−2 n−2       cn−2  yn−1 fn−1 cn−1 β bn−1

entonces la matriz del sistema es de diagonal dominante. Por lo tanto, puede utilizarse el algoritmo de Thomas para resolver el sistema. 521230

- 19 -

´ DIM – Universidad de Concepcion

Llamemos A a la matriz tridiagonal del sistema anterior, y

:= (y1 , . . . , yn−1 )t al vector

´ de incognitas y z al segundo miembro. Es decir,

Ay = z. Sea y ex

´ exacta. := (y(x1 ), . . . , y(xn−1 ))t el vector de valores de la solucion

Si no se desprecian los errores de truncamiento se tiene

Ay ex = z + τ ,

t 2

donde kτ k = (τ1 , . . . , τn−1 ) = O(h ).

Entonces el vector de errores y ex

− y satisface

A(y ex − y) = Ay ex − Ay = τ y, por lo tanto, y ex − y = A−1 τ .

−1 Se demuestra que A ≤ C , con C una constante independiente de h. Por lo tanto,

se tiene que

521230

−1 ky ex − yk ≤ A kτ k = O(h2 ). - 20 -

´ DIM – Universidad de Concepcion

´ Metodos de Elementos Finitos Supongamos que queremos resolver el P.V.C.

  −u00 + u = f, en (a, b),  u(a) = 0, u(b) = 0. ´ Para poder aplicar el metodo de elementos finitos (M.E.F.) primero debemos obtener la ´ v tal ´ debil ´ formulacion del problema. Para ello, multiplicamos la E.D.O. por una funcion que v(a)

= v(b) = 0, e integramos sobre el intervalo (a, b): Z b Z b Z b − u00 (x)v(x) dx + u(x)v(x) dx = f (x)v(x) dx. a

521230

a

- 21 -

a

´ DIM – Universidad de Concepcion

´ Integrando por partes el primer termino vemos que

−

Z

b 00

u (x)v(x) dx = a

Z

b a

b Z u0 (x)v 0 (x) dx − u0 (x)v(x) = a

b

u0 (x)v 0 (x) dx, a

´ anterior se tiene que la solucion ´ u del P.V.C. con lo que sustituyendo en la expresion satisface

Z

b

u0 (x)v 0 (x) dx + a

Z

b

u(x)v(x) dx = a

Z

b

f (x)v(x) dx. a

´ debil ´ Esto conduce a la formulacion del P.V.C.: Hallar

u : (a, b) → R tal que u(a) = u(b) = 0 y que satisfaga Z b Z b Z b u0 (x)v 0 (x) dx + u(x)v(x) dx = f (x)v(x) dx a

para toda v

521230

a

a

: (a, b) → R tal que v(a) = v(b) = 0. - 22 -

´ DIM – Universidad de Concepcion

´ ´ de la El objetivo del metodo de elementos finitos (M.E.F.) es buscar una aproximacion ´ finita Vh (llamado espacio discreto) que satisfaga ´ u en un espacio de dimension funcion ´ debil ´ con todas las funciones vh del mismo espacio Vh : la formulacion Hallar uh

Z

∈ Vh tal que

b a

u0h (x)vh0 (x) dx +

Z

b

uh (x)vh (x) dx = a

Z

b

f (x)vh (x) dx

∀vh ∈ Vh .

a

´ ´ debil ´ discreta del P.V.C. Esta se denomina la formulacion

521230

- 23 -

´ DIM – Universidad de Concepcion

´ sencillo de espacio discreto Vh (y el mas ´ usado) es el de las funciones El ejemplo mas ´ arbitraria del intervalo continuas y lineales a trozos en una particion

a = x0 < x1 < · · · < xn = b, que se anulan en los extremos a y b:  Vh := vh ∈ C(a, b) : vh |[xi−1 ,xi ] ∈ P1 , i = 1, . . . , n, y vh (a) = vh (b) = 0 .

´ h = max (xi − xi−1 ) denota la norma de la particion. 1≤i≤n

a=x 0 521230

x1

x2 - 24 -

x n−1

x n =b

x

´ DIM – Universidad de Concepcion

´ de Vh es n − 1. La dimension En efecto, la base natural de este espa-

Φ i(x)

cio esta´ formada por las llamadas fun- 1 ciones techo, φ1 , . . . , φn−1 , que satisfacen

  1, si i = j, φi (xj ) =  0, si i 6= j.

a=x 0

x i−1

xi

x i+1

xn =b

x

´ definidas por Estas funciones estan

 x−x i−1      xi − xi−1 xi+1 − x φi (x) :=  xi+1 − xi     0 521230

- 25 -

si

x ∈ [xi−1 , xi ],

si

x ∈ [xi , xi+1 ],

si

x 6∈ [xi−1 , xi+1 ].

´ DIM – Universidad de Concepcion

´ discreta uh Con respecto a esta base, la solucion

uh (x) =

n−1 X

∈ Vh puede representarse como

αj φj (x), con αj = uh (xj ).

j=1

´ de uh en la formulacion ´ debil ´ discreta Usando esta expresion

Z

b a

u0h (x)vh0 (x) dx +

Z

b

uh (x)vh (x) dx = a

Z

b

f (x)vh (x) dx

∀vh ∈ Vh

a

y tomando vh

= φi , i = 1, . . . , n − 1, por la linealidad del problema se llega a ! Z Z n−1 b b X φ0j (x)φ0i (x) dx + φj (x)φi (x) dx αi i=1

a

a

=

Z

b

f (x)φi (x) dx, i = 1, . . . , n − 1. a

´ ´ αi Este es un sistema de (n − 1) ecuaciones con (n − 1) incognitas 521230

- 26 -

= uh (xi ).

´ DIM – Universidad de Concepcion

1 Φ j (x)

a=x 0

Φ i (x)

xj

xi

x N=b

x

Como se ve en la figura, si j

∈ / {i − 1, i, i + 1}, entonces Z b Z b φ0j (x)φ0i (x) dx = 0 y φj (x)φi (x) dx = 0. a

a

´ debil ´ discreta se reduce a resolver un sistema de ecuaciones As´ı, vemos que la formulacion ´ ´ se demuestra es definida positiva (en ese con matriz tridiagonal y simetrica, que ademas ´ puede usarse el algoritmo de Thomas para resolver el sistema): caso tambien

521230

- 27 -

´ DIM – Universidad de Concepcion

 a1    b2   0  .  ..  0

donde

ai =

bi =

fi =

521230

Z Z

b a

b2

0

···

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

..

.

···

0

2

0

[φ0i (x)] dx +

b

Z





2

[φi (x)] dx, i = 1, . . . , n − 1, a

b

φ0i−1 (x)φ0i (x) dx +

a



f1 α1     ..       f .   α2  2        ..   ..   .  =  . , 0            α f n−2 n−2     bn−1  αn−1 fn−1 an−1

bn−1 Z



Z

b

φi−1 (x)φi (x) dx, i = 2, . . . , n − 1, a

b

f (x)φi (x) dx, i = 1, . . . , n − 1. a

- 28 -

´ DIM – Universidad de Concepcion

Para estimar el error, consideramos la siguiente norma:

kvk :=

(Z

b

2

|v(x)| dx a

)1/2

.

´ del P.V.C. y uh la solucion ´ calculada por el metodo ´ Teorema: Sea u la solucion de elementos finitos y supongamos que u tiene derivada segunda acotada en [a, b]. ´ del error: Entonces se tiene la siguiente estimacion

ku − uh k ≤ Ch2 max |u00 (x)|, 0≤x≤1

ku0 − u0h k ≤ Ch max |u00 (x)|, 0≤x≤1

donde C es una constante positiva independiente de h.

521230

- 29 -

´ DIM – Universidad de Concepcion

Observaciones: ´ • Si se aplica el metodo de elementos finitos en una particion uniforme y se calculan las ´ de los trapecios, se integrales que definen los coeficientes ai , bi y fi por el metodo ´ llega exactamente al mismo sistema de ecuaciones lineales que con el metodo de diferencias finitas. ´ • Ambos metodos, diferencias finitas y elementos finitos, se extienden a Ecuaciones en ´ Derivadas Parciales (E.D.P.), que es donde se manifiesta la mayor eficacia del metodo de elementos finitos. ´ • El metodo de “shooting”, en cambio, es muy eficiente para resolver P.V.C. para E.D.O., pero no se extiende a E.D.P. ´ acerca del metodo ´ ´ a la • Para mayor informacion de elementos finitos, y de su aplicacion ´ de E.D.P., se sugiere el siguiente texto: resolucion E.B. B ECKER , G.F. C AREY, J.T. O DEN . Finite Elements. An Introduction. Prentice Hall, 1981.

521230

- 30 -

´ DIM – Universidad de Concepcion