Solución Numérica de Problemas de Valor de Frontera para Ecuaciones Diferenciales Ordinarias 1

Universidad de Los Andes Facultad de Ciencias ´ticas Departamento de Matema ´n Grupo de Ciencias de la Computacio

Soluci´ on Num´ erica de Problemas de Valor de Frontera para Ecuaciones Diferenciales Ordinarias1

´rquez Br. Luis Jos´ e Berbes´ı Ma Trabajo Especial de Grado Para Optar al T´ıtulo de ´ticas. Licenciado en Matema ´ n. Tutor: Dr. Giovanni Caldero

M´ erida-Venezuela 2010 1

Este trabajo fue parcialmente financiado por Consejo de Desarrollo Cient´ıfico Human´ıstico y Tecnol´ ogico, CDCHT ULA, bajo el proyecto C - 1687 - 10 - 02 - ED

RESUMEN En muchas ´ areas de las Ciencias B´ asicas e Ingenier´ıa existen problemas donde es necesario encontrar la soluci´ on de una Ecuaci´ on Diferencial Ordinaria. Este tipo de problema, representa en la actualidad, uno de los temas b´ asicos del An´ alisis Num´erico. Un caso particular, el cual representar´ a la base de este trabajo, est´ a dado por encontrar la aproximaci´ on de la soluci´ on de un Problema de Valor de Frontera (PVF) de segundo orden, esto es, aproximar la soluci´ on de la ecuaci´ on diferencial ordinaria:   y ′′ = f (x, y, y ′ ), y(a) = α,  y(b) = β.

En los u ´ltimos a˜ nos, se han introducido en la literatura una variedad de m´etodos num´ericos para resolver este problema. En este trabajo se estudian cinco m´etodos para estimar la soluci´ on del problema planteado. Estos m´etodos son: el m´etodo del Disparo y Diferencias Finitas, considerados cl´ asicos en la literatura puesto que est´ an presentes en la mayor´ıa de los textos de An´ alisis Num´erico que abordan el problema propuesto, el m´etodo de los Elementos Finitos y Galerkin Discontinuo, los cuales parten de la teor´ıa de Espacios Normados y del An´ alisis Funcional, y el m´etodo Bvp4c, un m´etodo que usa la idea se superposici´on y el mismo se encuentra implementado en Matlab.

´Indice general

Introducci´ on

III

1. Preliminares

1 O(hn )

1.1. Orden de aproximaci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ´ 1.2. Teor´ıa de Algebra Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.3. Teor´ıa de Ecuaciones Diferenciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.4. Teor´ıa de An´ alisis Funcional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.4.1. Los espacios normados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.4.2. El espacio euclidiano y el espacio unitario . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.3. Funcionales lineales

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.4.4. Formas bilineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.4.5. El Teorema de Lax - Milgram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.5. M´etodo de Newton Raphson y Runge Kutta . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2

2. M´ etodos Num´ ericos para aproximar la soluci´ on de un Problema de Valor de Frontera de segundo orden 11 2.1. M´etodo del Disparo

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.1.1. Planteamiento para el caso

y ′′

2.1.2. Planteamiento para el caso

y ′′

=

p(x)y ′

=

f (x, y, y ′ )

+ q(x)y + r(x) . . . . . . . . . . . . . . . .

11 12

. . . . . . . . . . . . . . . . . . . . . . .

14

2.2. M´etodo de Diferencias Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2.1. Planteamiento para el caso

y ′′

2.2.2. Planteamiento para el caso

y ′′

=

p(x)y ′

=

f (x, y, y ′ )

+ q(x)y + r(x) . . . . . . . . . . . . . . . .

17

. . . . . . . . . . . . . . . . . . . . . . .

20

2.3. M´etodo de Elementos Finitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.3.1. Formulaci´ on d´ebil de un problema modelo unidimensional . . . . . . . . . . . . . .

23

2.3.2. El MEF para el problema modelo con funciones lineales a trozos . . . . . . . . . .

24

2.3.3. Integraci´ on num´erica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

2.4. M´etodo Galerkin Discontinuo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.4.1. Formulaci´ on d´ebil de un problema modelo unidimensional . . . . . . . . . . . . . .

28

i

´Indice General 2.4.2. Derivaci´ on del sistema lineal Aα = b mediante el uso de polinomios discontinuos a trozos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.3. Construcci´ on de la matriz A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31 32

2.4.4. Construcci´ on del vector b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5. M´etodo Bvp4c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34 36

3. Experimentaci´ on Num´ erica

39

A. Pruebas de algunas afirmaciones

65

A.1. Correspondientes al M´etodo del Disparo . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2. Correspondientes al M´etodo de Diferencias Finitas . . . . . . . . . . . . . . . . . . . . . . Referencias

65 66 68

ii

Introducci´ on

En muchas ´ areas de las Ciencias B´ asicas e Ingenier´ıa existen problemas donde es necesario encontrar la soluci´ on de una Ecuaci´ on Diferencial Ordinaria. Estas describen fen´ omenos que cambian frecuentemente. Com´ unmente, una soluci´ on de inter´es est´ a determinada especificando los valores de todas sus componentes en un punto x = a. Esto es un Problema de Valor Inicial. Sin embargo, en muchas ocasiones, una soluci´ on est´ a determinada en m´ as de un punto. Un problema de este tipo es denominado Problema de Valor de Frontera (PVF). Un PVF muy trabajado en la actualidad son los de segundo orden, es decir, los PVF que se especifican en dos puntos:   y ′′ = f (x, y, y ′ ), y(a) = α,  y(b) = β.

Los PVF de segundo orden suelen ser comunes en todas las ramas de las Ciencias Experimentales. Por ejemplo, en F´ısica, las leyes de Newton y muchas otras se expresan como un problema de este tipo, en Biolog´ıa aparecen en el modelado de din´ amica de poblaciones, en la Qu´ımica surgen en la evaluaci´ on de las concentraciones de diversos reactivos durante una reacci´ on, etc. A lo largo de los a˜ nos se han desarrollado t´ecnicas para encontrar la soluci´ on anal´ıtica a Ecuaciones Diferenciales Ordinarias, y por ende, la soluci´ on de los PVF que surgen de los modelos anteriormente mencionados. Sin embargo, resulta habitual que en la mayor´ıa de los casos no se conozca la soluci´ on anal´ıtica de la misma, y s´ olo estudios cualitativos de dichas ecuaciones son presentados en la literatura. Debido a esto, en la pr´ actica, resulta impetuoso usar m´etodos num´ericos para dar aproximaciones num´ericas de la soluci´ on (aproximaciones tan buenas como se quiera). Hoy en d´ıa, en la literatura, existen muchos m´etodos que ayudan a estimar la soluci´ on de un PVF de segundo orden, y es por esta raz´ on, que este trabajo se concentrar´ a en el estudio de cinco de ellos. Estos m´etodos son: el m´etodo del Disparo y Diferencias Finitas (ver [9, 18]), considerados cl´ asicos dentro de la literatura, el m´etodo de los Elementos Finitos y Galerkin Discontinuo (ver [3, 17] y [1], respectivamente), los cuales se fundamentan en resultados del An´ alisis Funcional, y el m´etodo Bvp4c (ver [14]), un m´etodo que usa la idea se superposici´ on y el mismo se encuentra implementado en Matlab. Las implementaciones num´ericas de los cinco m´etodos mencionados es llevada a cabo en Matlab. No obstante, por considerarlo irrelevante, los c´ odigos no son dados expl´ıcitamente. En general, el trabajo queda dividido de la siguiente manera: en el Cap´ıtulo 1 se desarrollan los preliminares te´oricos necesarios para abordar adecuadamente el estudio de los cinco m´etodos. Se introducen ´ resultados del Algebra Lineal (ver [2]), Ecuaciones Diferenciales (ver [7]) y An´ alisis Funcional (ver [3, 17]). As´ı mismo, se introducen el m´etodo de Newton Raphson, tanto para una ecuaci´ on como para un sistema de ecuaciones no lineales, y el m´etodo de Runge Kutta, para resolver Problemas de Valor Inicial de iii

Introducci´ on Ecuaciones Diferenciales Ordinarias, estudiados en un curso b´ asico de An´ alisis Num´erico, necesarios al momento de definir los m´etodos num´ericos. En el Cap´ıtulo 2 se expondr´ an cada uno de los cinco m´etodos que se emplear´ an para aproximar la soluci´ on de un PVF de segundo orden. Se presenta la construcci´ on de cada m´etodo, y en mucho de ellos, el an´ alisis del error. En el u ´ltimo cap´ıtulo se har´ a una breve discusi´ on, m´ as que comparaci´ on, sobre la eficiencia y propiedades de cada uno de estos m´etodos, esto en vista de que no se presentan las condiciones adecuadas para realizar una comparaci´ on detallada y ´ optima entre ellos. En el Ap´endice A se realizan las demostraciones de varios teoremas presentes en el Cap´ıtulo 2. Por u ´ltimo, la Bibliograf´ıa presenta las referencias literarias usadas.

iv

1 Preliminares

En este cap´ıtulo se expondr´ an los preliminares matem´ aticos necesarios para abordar adecuadamente el estudio de los m´etodos num´ericos usados para aproximar la soluci´ on de problema de valor de frontera de segundo orden. As´ı pues, el objetivo principal de este cap´ıtulo es el desarrollo de las herramientas matem´ aticas comunmente usadas en el planteamiento y an´ alisis de estos m´etodos num´ericos. Al mismo tiempo, se introducen el m´etodo de Newton Raphson, tanto para una ecuaci´ on como para un sistema de ecuaciones no lineales, y el m´etodo de Runge Kutta para resolver Problemas de Valor Inicial de Ecuaciones Diferenciales Ordinarias, estudiados en un curso b´ asico de An´ alisis Num´erico, necesarios al momento de definir los metodos usados en este trabajo. Es importante resaltar que los resultados que se desarrollan en el presente cap´ıtulo no ser´ an tratados de forma detallada por completo, pues el objetivo es s´ olamente dejar asentado un material de repaso, cuyo conocimiento es importante para entender el resto de este trabajo. Muchos de estos temas se tratan ´ de modo m´as amplio en algunos textos de An´ alisis Num´erico (ver [6, 12, 18]), Algebra Lineal (ver [2]), Ecuaciones Diferenciales (ver [7]) y An´ alisis Funcional (ver [3]).

1.1.

Orden de aproximaci´ on O(hn )

Definici´ on 1.1 Se dice que f (h) es de orden g(h) cuando h → 0 y h 6= 0, se denotar´ a por f (h) = O(g(h)), si existen n´ umeros reales M > 0 y k > 0 tales que, |f (h)| ≤ M |g(h)|, siempre que |h| < k. Ejemplo 1.1 Consideremos las funciones f (x) = x3 + 2x2 y g(x) = x2 . Puesto que, x3 < x2 para |x| ≤ 1, se obtiene, |x3 + 2x2 | < 3|x2 |. Por lo tanto, f (x) = O(g(x)). Definici´ on 1.2 Sean p y f funciones, se dice que p(h) aproxima a f (h) con un orden de aproximaci´ on O(hn ), lo que se denota por f (h) = p(h) + O(hn ), si existe un n´ umero real M > 0 y un n´ umero natural n tales que, |f (h) − p(h)| ≤ M, para h suficientemente peque˜ no. |hn | Al considerar el caso en que p(x) es la n-´esima aproximaci´ on por polinomios de Taylor a f (x) alrededor de x0 , es decir n X f (k) (x0 ) f n+1 (c) f (x) = (x − x0 )k + (x − x0 )n+1 , k! (n + 1)! k=0

1

Cap´ıtulo 1: Preliminares para alg´ un c entre x y x0 . Cuando x − x0 → 0, por la definici´ on 1.1 se tiene que O((x − x0 )n+1 ) =

f n+1 (c) (x − x0 )n+1 , (n + 1)!

as´ı, f (x) = p(x) + O((x − x0 )n+1 ), es decir, p(x), con p(x) el primer t´ermino del lado derecho de la igualdad, se aproxima a f (x) con un orden de aproximaci´ on O((x − x0 )n+1 ). Luego, el Teorema de Taylor se puede enunciar de la siguiente forma: Teorema 1.1 (Teorema de Taylor) Si f ∈ C n+1 [a, b] y x0 ∈ [a, b], entonces para cada x ∈ [a, b], f (x0 + h) =

n X f (k) (x0 ) k=0

1.2.

k!

hk + O(hn+1 ), donde h = x − x0 .

´ Teor´ıa de Algebra Lineal

Definici´ on 1.3 Se dice que una matriz A de n × n es invertible, si existe una matriz A−1 de n × n, con AA−1 = A−1 A = I, donde I es la matriz identidad. La matriz A−1 se llama inversa de A. Una matriz que no tiene inversa se le da el nombre de no invertible o singular. Definici´ on 1.4 (Matriz tridiagonal) siempre que |i − j| > 1. Esto es,  a11   a21    0 A=  ..  .   ..  . 0

Una matriz A = (aij ) de n × n se dice tridiagonal si aij = 0, a12

0

a22 a23

··· .. .

···

0 .. . .. .

.. . a32 a33 a34 .. .. .. .. . . . . 0 .. .. .. . . . an−1,1 ··· ··· 0 an,n−1 ann



     .     

Ahora bien, en muchas ocasiones estamos interesados en saber qu´e condiciones debe cumplir una matriz tridiagonal para que sea invertible. Dichas condiciones nos los da el siguiente teorema. Teorema 1.2 Sea A una matriz tridiagonal, con entradas [A]ij = ai,j . Supongamos que, para cada i = 2, 3, . . . , n − 1, se cumple ai,i−1 ai,i+1 6= 0. Si |a11 | > |a12 |, |aii | ≥ |ai,i−1 | + |ai,i+1 |, para i = 2, 3, . . . , n − 1 y |ann | > |an,n−1 |, entonces A es invertible. Demostraci´ on: Ver referencia [9], Cap´ıtulo 2, Secci´ on 3, p´ ag. 56.

1.3.

Teor´ıa de Ecuaciones Diferenciales

Es bien conocido en el quehacer cient´ıfico, bien sea en las ciencias puras o la ingenier´ıa, la utilidad de las ecuaciones diferenciales. El estudio de las ecuaciones diferenciales han influenciado el desarrollo de varias ´ areas de las matem´ aticas, como es el caso de la Topolog´ıa, el An´ alisis Real, y nuestra ´ area de estudio: el An´ alisis Num´erico. A lo largo de los a˜ nos se han desarrollado t´ecnicas para encontrar soluciones anal´ıticas a ecuaciones diferenciales ordinarias que modelan un sin fin de problemas que surgen de cualquiera de las ´ areas antes 2

Secci´ on 1.3: Teor´ıa de Ecuaciones Diferenciales mencionadas. Sin embargo, resulta habitual que en la mayor´ıa de los casos no se conozca la soluci´ on anal´ıtica de la misma, y s´ olo estudios cualitativos de dichas ecuaciones son presentados en la literatura. Debido a esto, en la pr´ actica, resulta impetuoso definir m´etodos num´ericos para dar aproximaciones num´ericas de la soluci´ on (aproximaciones tan buenas como se quiera). Vale se˜ nalar que muchos de estos estudios cualitativos se basan en aproximaciones num´ericas para iniciar la investigaci´ on del problema. Antes de continuar, es importante recordar algunas definiciones y conceptos b´ asicos de ecuaciones diferenciales. Una ecuaci´ on diferencial, el cual denotaremos mediante ED, es una ecuaci´ on que involucra una o m´ as variables dependientes, con respecto a una o m´ as variables independientes. Por ejemplo, en la ecuaci´ on: d2 x + 16x = 0, dt2 se tiene que el s´ımbolo x representa una variable dependiente mientras que la variable independiente es t. Si una ecuaci´ on contiene s´ olo las derivadas ordinarias de una o m´ as variables dependientes con respecto a una s´ ola variable independiente, entonces se dice que es una ecuaci´ on diferencial ordinaria (EDO). Las ecuaciones: dy + 5y = ex , dx d2 y dy + 6y = 0, − dx2 dx dx dy + = 2x + y, dt dt son ejemplos de ecuaciones diferenciales ordinarias. Es importante resaltar que en este trabajo las derivadas dy d2 y d3 y , ordinarias se escribir´ an con la notaci´ on de Leibniz: , ,. . . , o bien la notaci´ on del s´ımbolo dx dx2 dx3 ′ ′′ ′′′ de prima: y , y , y . En realidad, la notaci´ on de prima se emplea para denotar s´ olo las tres primeras derivadas, la cuarta derivada se escribir´ a y (4) en lugar de y ′′′′ ; m´ as general, la n-´esima derivada de y se dn y o y (n) . Por otro lado, el orden de la m´ as alta derivada de una ED se llama orden de la escribir´ a dxn ecuaci´ on. Por ejemplo: dy 3 d2 y + 5 − 4y = ex , dx2 dx es una EDO de segundo orden. En s´ımbolos, una EDO de n-´esimo orden de una variable dependiente, se puede expresar mediante la forma general F x, y, y ′ , . . . , y (n) = 0, (1.1)

donde F es una funci´ on de valores reales de n + 2 variables: x, y, y ′ ,. . . ,y (n) . Por razones pr´ acticas y te´ oricas supondremos de aqu´ı en adelante que es posible despejar en una EDO en forma u ´nica la derivada superior y (n) de una ED que est´e en la forma (1.1), en t´erminos de n + 1 variables restantes. La ecuaci´ on diferencial dn y = f x, y, y ′ , . . . , y (n−1) , n dx donde f es una funci´ on continua de valores reales, se denomina forma normal de (1.1). As´ı, cuando sea conveniente, se usan las formas normales dy = f (x, y) y dx

d2 y = f x, y, y ′ 2 dx

para representar las EDO de primer y segundo orden.

3

Cap´ıtulo 1: Preliminares Un concepto importante dentro del marco de las EDO es el concepto de linealidad. Se dice que una EDO de orden n es lineal si F es lineal en y, y ′ ,. . ., y (n) . Esto significa que una EDO de orden n es lineal cuando (1.1) se puede reescribir como an (x)y (n) + an−1 (x)y (n−1) + · · · + a1 (x)y ′ + a0 (x)y = g(x). Las ecuaciones diferenciales lineales se caracterizan por dos propiedades: 1. la variable independiente y junto con todas sus derivadas son de primer grado, es decir, la potencia de cada t´ermino en y es 1, 2. los coeficientes a0 , a1 ,. . . , an , de y, y ′ ,. . . ,y (n) , dependen s´ olo de la variable independiente x. Una ecuaci´ on que no es lineal se dice no lineal. Con frecuencia nos interesa problemas en los que se busca una soluci´ on y(x) de una ED de modo que y(x) satisfaga condiciones adicionales prescritas, es decir, condiciones impuestas en la y(x) desconocida o sus derivadas. Dos conceptos importantes dentro del marco de las EDO son los conceptos de problema de valor inicial y problema de valor de frontera. Para una ED, un Problema de Valor Inicial (el cual se denotar´ a como PVI) de n-´esimo orden es resolver y (n) = f x, y, y ′ , . . . , y (n−1) ,

sujeta a:

y(x0 ) = y0 , y ′ (x0 ) = y1 , . . . , y (n−1) (x0 ) = yn−1 .

Otro tipo de problemas consiste en resolver una ED de orden dos o mayor en el que la variable dependiente y o sus derivadas se especifican en diferentes puntos. Un problema que consiste en resolver y ′′ = f x, y, y ′

sujeta a:

y(a) = y0 , y(b) = y1 , se llama Problema de Valor de Frontera (el cual se denotar´ a mediante PVF). Los valores preescritos y(a) = y0 y y(b) = y1 se llaman condiciones de frontera. Ahora bien, estamos interesados en saber bajo qu´e condiciones un PVF de segundo orden tiene soluci´ on. El siguiente teorema establece condiciones generales que garantizan la existencia y unicidad de la soluci´ on a dicho problema de segundo orden. Teorema 1.3 Supongamos que la funci´ on f en el PVF y ′′ = f (x, y, y ′ ), a ≤ x ≤ b, y(a) = α, y(b) = β, es continua en el conjunto D = {(x, y, y ′ ) / a ≤ x ≤ b, −∞ < y < ∞, −∞ < y ′ < ∞}, y que fy y fy′ tambi´en son continuas en D. Si 1. fy (x, y, y ′ ) > 0, para toda (x, y, y ′ ) ∈ D, y 2. existe una constante M , con |fy′ (x, y, y ′ )| ≤ M, para todo (x, y, y ′ ) ∈ D, entonces el PVF tiene soluci´ on u ´nica. 4

Secci´ on 1.4: Teor´ıa de An´ alisis Funcional Demostraci´ on: Ver referencia [13].

Como consecuencia del teorema anterior, tenemos el siguiente resultado. Corolario 1.1 Si el problema lineal con valor en la frontera y ′′ = p(x)y ′ + q(x)y + r(x), a ≤ x ≤ b, y(a) = α, y(b) = β, satisface: 1. p(x), q(x) y r(x) son continuas en [a, b], y 2. q(x) > 0 en [a, b], entonces el problema tiene soluci´ on u ´nica.

1.4.

Teor´ıa de An´ alisis Funcional

Para considerar los aspectos referentes al M´etodo de los Elementos Finitos y al M´etodo Galerkin Discontinuo, es importante conocer algunos resultados del An´ alisis Funcional.

1.4.1.

Los espacios normados

En un espacio lineal podemos asignar a cada elemento x la noci´ on de longitud por medio del n´ umero real ||x|| que es llamado norma de x. En particular, la norma es una funci´ on de valor real definida en un espacio lineal E que satisface las siguientes condiciones o axiomas: 1. ||x|| ≥ 0, con la particularidad de que ||x|| = 0 s´ olo si x = 0. 2. ||αx|| = |α|||x||, donde α es un n´ umero real. 3. ||x + y|| ≤ ||x|| + ||y||, x, y ∈ E. Ejemplo 1.2 Si E = R2 , la selecci´ on usual para la longitud o norma || · ||2 de un vector x = (x1 , x2 ) en el plano real es q ||x||2 = x21 + x22 , que es llamada norma 2 o norma euclidiana. Otro ejemplo de norma es la norma 1 que se define como ||x||1 = |x1 | + |x2 |, o la norma infinita que se define como ||x||∞ = m´ ax{|x1 |, |x2 |}. Ejemplo 1.3 Sea E = C[a, b]. La funci´ on de valor real definida por ||f || = m´ ax |f (t)| a≤t≤b

es una norma en E, conocida como norma uniforme o norma Chebyshev de f en el espacio C[a, b]. Una segunda norma en C[a, b], que es tambi´en muy importante, es la llamada norma L2 , que se denota Z b 1/2 2 || · ||L2 = [f (t)] dt . a

5

Cap´ıtulo 1: Preliminares Sea E un espacio normado. La sucesi´ on {un }∞ n=1 ⊂ E se dice acotada en E si existe un C > 0 tal que ||un || < C, para todo n. La sucesi´ on se dice convergente en E si existe un elemento u ∈ E tal que para todo 0 < ε ∈ R, existe un N0 ∈ N tal que ||u − un || < ε, para todo n > N0 . El elemento u es el l´ımite de la sucesi´ on {un }∞ a de forma indistinta n=1 . Usualmente se escibir´ l´ım un = u

o

n→∞

||u − un || → 0, para n → ∞.

La sucesi´ on {un }∞ a una sucesi´ on de Cauchy si para todo ε > 0, existe N0 ∈ N tal que n=1 ⊂ E ser´ ||un − um || ≤ ε, para todo n, m ≥ N0 . Un espacio lineal normado E se llama completo, cuando toda sucesi´ on de Cauchy de este espacio converge hacia un cierto elemento u ∈ E. Los espacios lineales normados completos se llaman espacios de Banach. Todo espacio lineal normado de dimensi´ on finita es completo.

1.4.2.

El espacio euclidiano y el espacio unitario

Una forma conocida de introducir una norma en un espacio lineal consiste en definir en ´este el producto escalar o producto interior. Se llama producto escalar en un espacio lineal real (o complejo) E a una funci´ on (·, ·) : E × E → R (o C) con las siguientes propiedades para todo par de elementos x, y ∈ E: 1. (x, y) = (y, x), 2. (x + y, z) = (x, z) + (y, z), 3. (λx, y) = λ(x, y), λ ∈ R (o C), 4. (x, x) ≥ 0, y (x, x) = 0 s´ olo si x = 0. En el caso de un espacio lineal real el primer axioma se reduce a la simetr´ıa (x, y) = (y, x). Ejemplo 1.4 Sean x(t) e y(t) dos funciones de C 2 [a, b]. Se pueden definir varios productos escalares: (x, y) =

Z

b

x(t)y(t)dt,

a

(x, y) =

Z bh a

(x, y) =

Z bh a

i x(t)y(t) + x′ (t)y ′ (t) dt,

i x(t)y(t) + x (t)y (t) + x (t)y (t) dt, ′

′

′′

′′

donde las integrales son integrales de Riemann. Lema 1.1 Sea E un espacio con producto interno. Entonces la funci´ on ||u|| = es una norma en E.

p

(u, u), u ∈ E

6

Secci´ on 1.4: Teor´ıa de An´ alisis Funcional Los espacios dotados de un producto interno son llamados espacio producto interno. Un espacio lineal real normado E, en el cual la norma est´ a generada por un producto escalar, es decir, p ||u|| = (u, u),

se llama espacio euclidiano; en el caso complejo recibe el nombre de espacio unitario. Un espacio unitario (euclidiano) completo es llamado espacio de Hilbert H. Ahora se introducen algunos espacios de Hilbert que resultan naturales para la formulaci´ on variacional de los problemas de valor de frontera que ser´ an considerados. Si Ω es un dominio acotado en R2 , se define el espacio de funciones de cuadrado integrables sobre Ω: Z L2 (Ω) = {v : Ω → R tal que v 2 dΩ < ∞}. Ω

El espacio L2 (Ω) es un espacio de Hilbert con el producto escalar Z (u, v)L2 (Ω) = vudΩ, Ω

y la correspondiente norma ||v||L2 (Ω) =

Z

Ω

[v(x)]2 dx

1/2

1/2

= (v, v)L2 (Ω) .

Se introduce adem´ as el espacio ∂v 2 2 H (Ω) = v ∈ L (Ω) : ∈ L (Ω), i = 1, 2 , ∂xi 1

dotado del producto escalar y la correspondiente norma Z (v, w)H1 (Ω) = (vw + ∇v∇w)dΩ, Ω

1/2

|| · ||H1 (Ω) = (·, ·)H1 (Ω) .

El espacio H1 (Ω) consiste de todas las funciones v definidas en Ω junto con sus primeras derivadas que son de cuadrado integrable, es decir, pertenecen a L2 (Ω). Adem´ as, se define el espacio H10 = v ∈ H1 (Ω) : v = 0 en ∂v ,

que se dota del mismo producto interno y norma que el espacio H1 (Ω). Para evitar sobrecargar la notaci´ on, siempre y cuando no se lleve a confusi´ on, se utilizar´ a la notaci´ on (·, ·) y ||·|| en lugar de (·, ·)L2 (Ω) y ||·||L2 (Ω) . De forma an´ aloga, (·, ·)1 y || · ||1 denotar´ an (·, ·)H1 (Ω) y || · ||H1 (Ω) , respectivamente. Los espacios L2 (Ω) y H1 (Ω) pertenecen a familias de espacios m´ as generales, conocidas como espacios con p ≥ 1,y espacios de Sobolev, respectivamente.

Lp (Ω),

1.4.3.

Funcionales lineales

Un funcional lineal es un caso particular de un operador lineal, es decir, es un operador lineal que transforma el espacio dado E en valores de R o C. En general, si E es un espacio lineal normado sobre un cuerpo K, entonces un funcional lineal f es una aplicaci´ on lineal de E → K tal que 1. f (u + w) = f (u) + f (v), ∀u, v ∈ E, 2. f (αu) = αf (u), ∀u ∈ E, α ∈ K. 7

Cap´ıtulo 1: Preliminares El espacio de todos los funcionales lineales de un espacio normado E sobre su cuerpo, L(E, K), se le llama espacio dual de E, y suele denotarse por E′ . Consideremos algunos ejemplos de funcionales. Ejemplo 1.5 Sea E = Rn . El operador F : E → R definido como el promedio de las componentes del vector n 1X f (v) = vi , para todo v ∈ E, n i=1

es un funcional lineal sobre E, es decir, f ∈ E′ . Ejemplo 1.6 El operador integral A : C[a, b] → R, definido por A(f ) =

Z

b

f (x)dx,

para todo v ∈ C[a, b],

a

es un funcional lineal sobre C[a, b]. Teorema 1.4 (Teorema de Representaci´ on de Riesz) Sea H un espacio de Hilbert y sea F un funcional lineal continuo de H. Entonces existe un u ´nico elemento u ∈ H tal que F (v) = (u, v),

∀v ∈ H.

Adem´ as ||F || = ||u||

1.4.4.

Formas bilineales

Otro tipo especial de operador que pueden ocurrir muy frecuentemente en el estudio de problemas de valor de frontera es uno que mapea un par de elementos a los n´ umeros reales, y que es lineal en cada uno de estos. Un operador B : U × V → R (U, V espacios lineales) se dice que es una forma bilineal si: 1. B(αu + βw, v) = αB(u, v) + βB(w, v),

u, w ∈ U, v ∈ V

2. B(u, αv + βw) = αB(u, v) + βB(u, w),

u ∈ U, v, w ∈ V

Formas bilineales continuas. Consideremos una forma bilineal B : U × V → R, donde U y V son espacios lineales normados. Si existe una constante k > 0 tal que |B(u, v)| ≤ k||u||||v||

para todo u ∈ U, v ∈ V,

entonces B es llamada una forma bilineal continua. Formas bilineales H - el´ıpticas. Dada una forma bilineal B : H × H → R, donde H es un espacio con producto interno, se dice que B es H - el´ıptica si existe una constante α > 0 tal que B(v, v) ≥ α||v||2 ,

∀v ∈ H.

As´ı una forma H - el´ıptica es acotada inferiormente. 8

Secci´ on 1.5: M´etodo de Newton Raphson y Runge Kutta

1.4.5.

El Teorema de Lax - Milgram

Teorema 1.5 (El Teorema de Lax - Milgram) Sea H un espacio de Hilbert y sea B : H × H → R una forma bilineal, H - el´ıptica, continua, definida en H. Entonces, dada cualquier funcional lineal F definida en H, existe un u ´nico elemento u ∈ H tal que B(u, v) = F (v),

∀v ∈ H,

adem´ as, ||u||1 ≤ α−1 ||f ||H′ .

1.5.

M´ etodo de Newton Raphson y Runge Kutta

En esta secci´ on se introducen dos m´etodos que ser´ an necesarios para el estudios de los m´etodos que se emplear´ an para aproximar la soluci´ on de un PVF de segundo orden. Estos m´etodos son el de Newton Raphson, tanto para una ecuaci´ on como para un sistema de ecuaciones no lineales, y el de Runge Kutta, para resolver PVI de ecuaciones diferenciales ordinarias. M´ etodo de Newton Raphson. Sea f una funci´ on dos veces diferenciable. Supongamos que α es una ra´ız simple de f (x) = 0 y x0 es una aproximaci´ on inicial a α. El m´etodo de Newton-Raphson consiste en conseguir la aproximaci´ on usando la recta tangente a la curva f , que pasa por el punto (x0 , f (x0 )). Esta recta intersecta al eje OX en el punto (x1 , 0), donde x1 ser´ a la aproximaci´ on a la ra´ız de f , luego x1 = x0 −

f (x0 ) . f ′ (x0 )

Para continuar, se renombra x0 = x1 y se repite el proceso tanto como desee. Teorema 1.6 (Convergencia del M´ etodo de Newton-Rapson) Sean f ∈ C 2 [a, b] y α ∈ [a, b] tal ′ on {xk }∞ que f (α) = 0. Si f (α) 6= 0, entonces existe δ > 0 tal que la sucesi´ k=0 definida por el proceso iterativo, f (xk−1 ) xk = g(xk−1 ) = xk−1 − ′ , para k = 1, 2, 3, . . . , f (xk−1 ) converge a α cualesquiera que sea la aproximaci´ on inicial x0 ∈ [α − δ, α + δ]. Demostraci´ on: Ver la referencia [18], Cap´ıtulo 2, Secci´ on 2.1, p´ ag. 69.

M´ etodo de Newton para sistemas de ecuaciones no lineales. En muchas ocasiones aparece la necesidad de resolver sistemas de ecuaciones algebraicas no lineales de la forma  f1 (x1 , x2 , . . . , xn ) = 0,    f2 (x1 , x2 , . . . , xn ) = 0,  ..  .    fn (x1 , x2 , . . . , xn ) = 0.

El m´etodo de Newton para sistema de ecuaciones no lineales es una extensi´ on directa del m´etodo del mismo nombre para buscar ceros de funciones de una variable. La idea es realizar el desarrollo de Taylor en n variables en el entorno de una raiz, ∂fi f(x) = f(x(0) ) + (x − x(0) ) + O(||x − x(0) ||2 ). ∂xj x(0) 9

Cap´ıtulo 1: Preliminares

Si x es una ra´ız aproximada, definiendo la matriz jacobiana como

J(x(0) )

=

∂fi ∂xj

encontramos que x(0)

x ≈ x(0) − J −1 (x(0) )f(x(0) ). En la pr´ actica, esta f´ ormula nos permite definir un proceso iterativo, x(k+1) = x(k) − J −1 (x(k) )f(x(k) ), que es el m´etodo de Newton en varias variables. A la hora de implementar el m´etodo, es costoso tener que calcular la inversa de J y por eso, normalmente, lo que se hace es resolver el sistema de ecuaciones lineales J(x(k) ) x(k) − x(k+1) = f(x(k) ),

on. del cual se puede despejar x(k+1) una vez conocida la soluci´ En [18], Cap´ıtulo 10, Secci´ on 10.2, se hace el estudio de la convergencia y la implementaci´on del m´etodo. M´ etodo de Runge Kutta. Es un m´etodo de resoluci´ on num´erica de ecuaciones diferenciales ordinarias, el cual fue inicialmente desarrollado alrededor del a˜ no 1900 por los matem´ aticos C. Runge 1 y M. W. 2 Kutta . El m´etodo de Runge Kutta no es s´ olo un m´etodo sino una importante familia de m´etodos iterativos tanto impl´ıcitos como expl´ıcitos. Un m´etodo cl´ asico es el m´etodo de RK expl´ıcito de cuarto orden. Para introducirlo definamos el PVI: y ′ = f (t, y), y(t0 ) = y0 . Entonces el m´etodo de RK de cuarto orden para este problema est´ a dado por la siguiente ecuaci´ on: yn+1 = yn + donde:

h k1 + 2k2 + 2k3 + k4 , 6

k1 = f (tn , yn ), h k2 = f tn + , yn + 2 h k3 = f tn + , yn + 2

k1 , 2 k2 , 2

k4 = f (tn + h, yn + k3 ). Por ser de cuarto orden, se tiene que el error de paso es de orden O(h5 ), mientras que el error total acumulado tiene orden O(h4 ). Para mayores detalles del m´etodo, ver [6], Cap´ıtulo 8, Secci´ on 8.3, p´ ag. 514.

1 C. Runge naci´ o el 30 de agosto de 1856 y muri´ o el 03 de enero de 1927. Fue un matem´ atico, f´ısico y espestroscopista alem´ an. Su primer nombre suele escribirse Carl. 2 Martin Wilhelm Kutta fue un f´ısico y matem´ atico alem´ an, naci´ o el 03 de noviembre de 1867 en Pitschen, Alta Silecia (en la actualidad pertenece a Polonia) y muri´ o el 25 de diciembre de 1944 en F¨ urstenfeldbruck, Alemania.

10

2 M´ etodos Num´ ericos para aproximar la soluci´ on de un Problema de Valor de Frontera de segundo orden

Las Ecuaciones Diferenciales Ordinarias describen fen´ omenos que cambian frecuentemente. Estas surgen en modelos totalmente matem´ aticos, f´ısicos o qu´ımicos. Un sistema de ecuaciones diferenciales ordinarias puede tener muchas soluciones. Com´ unmente, una soluci´ on de inter´es est´ a determinada especificando los valores de todas sus componentes en un punto x = a. Esto es un Problema de Valor Inicial. Sin embargo, en muchas aplicaciones, una soluci´on est´ a determinada en m´ as de un punto. Un problema de este tipo es denominado Problema de Valor de Frontera, el cual es nuestro problema a resolver, espec´ıficamente en dos puntos:   y ′′ = f (x, y, y ′ ), y(a) = α,  y(b) = β. Es importante resaltar que si el PVF anterior se expresa mediante

y ′′ = p(x)y ′ + q(x)y + r(x), a ≤ x ≤ b, y(a) = α, y(b) = β, entonces se dice que el PVF es del tipo lineal, en caso contrario, se dice que el PVF es del tipo no lineal. Problemas de este tipo suelen ser comunes en muchas ´ areas de la vida cotidiana, como, por ejemplo, en las Ciencias b´ asicas o en la Ingenier´ıa. Sin embargo, su principal dificultad radica en que la soluci´ on a dichos problemas no son f´ aciles de obtener anal´ıticamente. Debido a esto, resulta importante definir m´etodos num´ericos para dar estimaciones num´ericas de la soluci´ on. En este cap´ıtulo se expondr´ an cinco m´etodos num´ericos que nos ayudar´ an a encontrar estas estimaciones: el M´etodo del Disparo y el M´etodo de Diferencias Finitas, los cuales son presentados en las referencias [9, 18], el M´etodo de Elementos Finitos, planteado en [3, 17], el M´etodo Galerkin Discontinuo, planteado en la referencia [1], y el M´etodo Bvp4c, el cual se presenta en la referencia [14]. Es importante resaltar que, cuando el PVF es del tipo lineal, usaremos los cinco m´etodos mencionados, mientras que si el PVF es del tipo no lineal, entonces se emplear´ a s´ olo los m´etodos de Disparo, Diferencias Finitas y Bvp4c, en vista de que los m´etodos de Elementos Finitos y Galerkin Discontinuo ameritan un estudio mucho m´ as amplio que para el caso lineal, y esto conlleva a que la implementaci´ on de los c´ odigos requieran m´ as tiempo.

2.1.

M´ etodo del Disparo

Uno de los m´etodos num´ericos, considerado cl´ asico en la literatura, para aproximar la soluci´ on de un Problema de Valor de Frontera de segundo orden, es el m´etodo del Disparo. Este m´etodo, el cual se plantea en las referencias [9, 18], se reduce a resolver varios PVI para encontrar la soluci´ on al problema planteado. 11

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF

2.1.1.

Planteamiento para el caso y ′′ = p(x)y ′ + q(x)y + r(x)

Para aproximar la soluci´ on u ´nica del problema lineal y ′′ = p(x)y ′ + q(x)y + r(x), a ≤ x ≤ b, y(a) = α, y(b) = β,

(2.1)

garantizada por el cumplimiento de las hip´ otesis del Corolario 1.1, primero de considera los PVI: u′′ = p(x)u′ + q(x)u + r(x), a ≤ x ≤ b, u(a) = α, u′ (a) = 0

(2.2)

v ′′ = p(x)v ′ + q(x)v, a ≤ x ≤ b, v(a) = 0, v ′ (a) = 1.

(2.3)

y Seg´ un las hip´ otesis del Corolario 1.1, tanto (2.2) como (2.3) tienen soluci´ on u ´nica. La soluci´ on de (2.1) viene dada por β − u(b) y(x) = u(x) + v(x). (2.4) v(b) A partir de ´esta, β − u(b) ′ β − u(b) ′′ ′ ′ ′′ ′′ v (x), y (x) = u (x) + v (x). y (x) = u (x) + v(b) v(b) As´ı, y ′′ (x)

=

p(x)u′ (x)

+ q(x)u(x) + r(x) +

β − u(b) (p(x)v ′ (x) + q(x)v(x)) v(b)

β − u(b) ′ β − u(b) ′ = p(x) u (x) + v (x) + q(x) u(x) + v(x) + r(x) v(b) v(b) = p(x)y ′ (x) + q(x)y(x) + r(x). Como se puede apreciar, la soluci´ on (2.4) satisface las condiciones de frontera. En efecto, β − u(b) y(a) = u(a) + v(a) v(b) = α+

β − u(b) 0 v(b)

= α, y y(b) = u(b) +

β − u(b) v(b) v(b)

= u(b) + β − u(b) = β. Por tanto, y(x) es la soluci´ on u ´nica al problema de valor de frontera (sujeta a que v(b) 6= 0). El m´etodo del Disparo para las ecuaciones lineales se basa en la sustituci´ on del PVF por dos PVI. En el estudio de PVI para ecuaciones diferenciales ordinarias, se describen muchos m´etodos con los cuales se pueden aproximar las soluciones u(x) y v(x) (como ejemplo: el m´etodo de Euler, Trapecio, Punto Medio, Runge Kutta), y una vez que se cuenta con estas aproximaciones, la soluci´ on del PVF se aproxima por medio de la ecuaci´ on (2.4). Para efecto de nuestro estudio, el m´etodo que se usar´ a para aproximar (2.2) y (2.3) ser´ a el metodo de Runge Kutta. Desde el punto de vista gr´ afico, el m´etodo tiene el aspecto que se observa en la Figura 2.1. 12

Secci´ on 2.1: M´etodo del Disparo

b

β

b

v(x) y(x) b

u(x)

α b

b

a

b

Figura 2.1: M´etodo del Disparo para el caso lineal.

Algoritmo del M´ etodo del Disparo para el caso lineal. Entrada: extremos a, b; condiciones de frontera α, β; n´ umero de subintervalos N . Salida: aproximaciones w1,i a y(xi ); w2,i a y ′ (xi ), para cada i = 0, . . . , N. 1. Tomar h = (b − a)/N ; u1,0 = α; u2,0 = 0; v1,0 = 0; v2,0 = 1. 2. Para i = 1, . . . , N − 1, hacer los pasos 3 y 4. 3. Tomar x = a + ih. 4. Aplicar Runge - Kutta para obtener u1,i+1 , u2,i+1 , v1,i+1 , v2,i+1 . 5. Tomar w1,0 = α; w2,0 = (β − u1,N )/v1,N . Mostrar (a, w1,0 , w2,0 ). 6. Para i = 1, . . . , N , tomar w1,i = u1,i + w2,0 v1,i , w2,i = u2,i + w2,0 v2,i , x = a + ih. Mostrar (x, w1,i , w2,i ). 7. Parar. Problemas por errores de redondeo. Desafortunadamente, esta t´ecnica, por errores de rodondeo, puede contener problemas ocultos. Si u(x) crece r´ apidamente cuando x recorre [a, b], entonces u1,N es grande; si adem´ as, β es peque˜ no en comparaci´ on con u1,N , entonces w2,0 =

β − u1,N v1,N

≈ −

u1,N . v1,N

Por lo tanto, w1,i = u1,i + w2,0 v1,i ≈ u1,i −

u1,N v1,i , v1,N 13

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF lo cual permite una posible p´erdida de d´ıgitos significativos debido a la cancelaci´ on. Pero como u1,i es una aproximaci´ on a u(xi ), se puede entonces vigilar el comportamiento de u(x), y si u1,i aumenta r´ apidamente de a a b, podemos aplicar hacia atr´ as el m´etodo del disparo, esto es, resolver en su lugar los problemas de valor inicial u′′ = p(x)u′ + q(x)u + r(x), a ≤ x ≤ b, u(b) = β, u′ (b) = 0, y

v ′′ = p(x)v ′ + q(x)v, a ≤ x ≤ b, v(b) = 0, v ′ (b) = 1.

Si este m´etodo de disparo inverso todav´ıa presenta la eliminaci´ on de los d´ıgitos significativos y si el aumento de precisi´ on no produce mayor exactitud, ser´ a necesario utilizar otros m´etodos num´ericos, como las que explicaremos m´ as adelante. An´ alisis del error. En muchas ocasiones estamos interesados en encontrar cotas para el error que se est´ a produciendo al momento de aproximar. El siguiente teorema nos proporciona dichas cotas. Para no recargar demasiado la exposici´ on, presentaremos la prueba en el Ap´endice A, al que referimos al lector interesado. Teorema 2.1 Sea h = (b − a)/N , xi = a + nh, con i = 0, 1, . . . , N. Sean u(x), v(x) y y(x) las soluciones respectivas de (2.2), (2.3) y (2.4). Si u ei y vei son aproximaciones de O(hn ) para ui = u(xi ) y vi = v(xi ), respectivamente, entonces wi ser´ a una aproximaci´ on de O(hn ) para yi = y(xi ). En particular, vi n |wi − yi | ≤ Kh 1 + , vN para alguna K apropiada.

Demostraci´ on: Ver Ap´endice A, secci´ on A.1.

2.1.2.

Planteamiento para el caso y ′′ = f (x, y, y ′)

El m´etodo del Disparo para el PVF no lineal y ′′ = f (x, y, y ′ ), a ≤ x ≤ b, y(a) = α, y(b) = β,

(2.5)

utiliza las soluciones de una sucesi´ on de PVI de la forma que contenga un par´ ametro t, para aproximar la soluci´ on al PVF y ′′ = f (x, y, y ′ ), a ≤ x ≤ b, y(a) = α, y ′ (a) = t. (2.6) Esto se hace escogiendo los par´ ametros t = tk de tal forma que

l´ım y(b, tk ) = y(b) = β,

k→∞

donde y(b, tk ) denotar´ a la soluci´ on al PVI (2.6) con t = tk y y(x) denota la soluci´ on al PVF (2.5). Esta t´ecnica se conoce con el nombre de m´etodo del Disparo, por la analog´ıa con el procedimiento de dispararle a objetos situados en un blanco fijo (ver Figura 2.2). Se comienza con un par´ ametro t0 que determina la elevaci´ on inicial a la cual se le dispara al objeto desde el punto (a, α) y a lo largo de la curva descrita por la soluci´ on al PVI: y ′′ = f (x, y, y ′ ), a ≤ x ≤ b, y(a) = α, y ′ (a) = t0 . Si y(b, t0 ) no est´ a suficientemente cerca de β, se cambia la aproximaci´ on seleccionando las elevaciones t1 , t2 , y as´ı sucesivamente, hasta que y(b, tk ) est´e bastante cerca de acertar en el blanco β (Ver Figura 2.3). Para determinar los par´ ametros tk , supongamos que (2.5) satisface las hip´ otesis del Teorema (1.3). Si y(x, t) denota la soluci´ on de (2.6), el problema consistir´ a en determinar t tal que y(b, t) − β = 0.

(2.7)

Esta es una ecuaci´ on no lineal, y para resolverlo se disponen de varios m´etodos. 14

Secci´ on 2.1: M´etodo del Disparo

β y(b, t0 )

(b, y(b, t0 )) b

y(x, t0 ) α

b

Pendiente t0

(a, α) a

b

Figura 2.2: M´etodo del Disparo para el caso no lineal.

M´ etodo de la Secante. Para emplear este m´etodo, se necesita elegir las aproximaciones iniciales t0 y t1 y luego generar los t´erminos restantes de la sucesi´ on mediante tk = tk−1 −

(y(b, tk−1 ) − β)(tk−1 − tk−2 ) , k = 2, 3, . . . y(b, tk−1 ) − y(b, tk−2 )

M´ etodo de Newton. Para generar la sucesi´ on {tk } con este m´etodo, que es m´ as poderoso, s´ olo se necesita una aproximaci´ on inicial t0 . Sin embargo, la iteraci´ on tiene la forma tk = tk−1 −

y(b, tk−1 ) − β , dy (b, tk−1 ) dt

(2.8)

dy (b, tk−1 ). Esto presenta un problema porque no se conoce una representaci´ on expl´ıcita dt de y(b, t); se conoce s´ olo los valores y(b, t0 ), y(b, t1 ), . . . , y(b, tk−1 ). Supongamos que se reescribe el problema de valor inicial (2.6), haciendo ´enfasis en que la soluci´ on se basa tanto en x como en el par´ ametro t: y requiere conocer

y ′′ (x, t) = f (x, y(x, t), y ′ (x, t)), a ≤ x ≤ b, y(a, t) = α, y ′ (a, t) = t.

(2.9)

Se conserva la notaci´ on prima para indicar la derivada respecto a x. Puesto que se necesita determinar dy (b, t) cuando t = tk−1 , primero se toma la derivada parcial de (2.9) respecto a t. Esto significa que dt ∂y ′′ (x, t) = ∂t =

∂f (x, y(x, t), y ′ (x, t)) ∂t ∂f ∂x ∂f ∂y (x, y(x, t), y ′ (x, t)) (x, t) + (x, y(x, t), y ′ (x, t)) (x, t) ∂x ∂t ∂y ∂t +

′ ∂f ′ (x, t)) ∂y (x, t). (x, y(x, t), y ∂y ′ ∂t

Dado que x y t son independientes,

∂x = 0, y ∂t

∂y ′′ ∂f ∂y ∂f ∂y ′ (x, t) = (x, y(x, t), y ′ (x, t)) (x, t) + ′ (x, y(x, t), y ′ (x, t)) (x, t), ∂t ∂y ∂t ∂y ∂t

(2.10) 15

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF y(b, t2 ) β y(b, t3 ) y(b, t1 ) y(b, t0 )

α

y(x, t2 ) b

b

b

b

y(x, t3 ) y(x, t1 ) y(x, t0 )

b

(a, α) a

b

Figura 2.3: Escogencia de los par´ ametros tk .

con a ≤ x ≤ b. Las condiciones iniciales dan ∂y (a, t) = 0, ∂t ∂y ′ (a, t) = 1. ∂t ∂y (x, t) y si se supone que el orden de derivaci´ on ∂t de x y t puede invertirse, con las condiciones iniciales, la ecuaci´ on (2.10) se convierte en el PVI: Si se simplifica la notaci´ on usando z(x, t) para denotar

z ′′ (x, t) =

∂f ∂f (x, y, y ′ )z ′ (x, t) + (x, y, y ′ )z(x, t), a ≤ x ≤ b, z(a, t) = 0, z ′ (a, t) = 1. ′ ∂y ∂y

(2.11)

As´ı pues, el m´etodo de Newton requiere que los PVI (2.9) y (2.11), sean resueltos en cada iteraci´ on. Entonces, conforme a la ecuaci´ on (2.8), tk = tk−1 −

y(b, tk−1 ) − β . z(b, tk−1 )

(2.12)

Para efecto de nuestro estudio, se emplear´ a el m´etodo de Runge - Kutta de cuarto orden para aproximar las soluciones de (2.9) y (2.11). Algoritmo del M´ etodo del Disparo para el caso no lineal. Entrada: extremos a, b; condiciones de frontera α, β; n´ umero de subintervalos N ≥ 2; tolerancia Tol, n´ umero m´ aximo de iteraciones M . Salida: aproximaciones w1,i a y(xi ); w2,i a y ′ (xi ), para cada i = 0, . . . , N , o bien un mensaje de que se excedi´ o el n´ umero m´ aximo de iteraciones. 1. Tomar h = (b − a)/N ; u1,0 = α; k = 1, T K = (β − α)/(b − a). 2. Mientras k ≤ M , hace los pasos 3 - 10. 3. Tomar w1,0 = α, w2,0 = T K, u1 = 0, u2 = 1. 4. Para i = 1, . . . , N , hacer los pasos 5 - 6. 16

Secci´ on 2.2: M´etodo de Diferencias Finitas 5. Tomar x = a + (i − 1)h.

6. Aplicar Runge - Kutta para obtener w1,i , w2,i , u1 , u2 . 7. Si |w1,N − β| < Tol, entonces hacer los pasos 8 y 9.

8. Para i = 0, 1, . . . N , tomar x = a + ih, mostrar (x, w1,1 , w2,i ). 9. Procedimiento terminado, PARAR.

10. Tomar k = k + 1 y T K = T K − (w1,N − β)/u1 . El m´etodo de Newton se emplea para calcular T K. 11. N´ umero m´ aximo de iteraciones excedido, procedimiento terminado sin ´exito, PARAR. El valor t0 = T K escogido en el paso 1 es la pendiente de la recta que pasa por (a, α) y por (b, β). Como el m´etodo del Disparo emplea el m´etodo de Runge Kutta de cuarto orden para aproximar la soluci´ on de dos problemas de valor inicial, entonces la soluci´ on encontrada aproxima la soluci´ on exacta con un orden de aproximaci´ on O(h4 ) (ver [9], Cap´ıtulo 8, Secci´ on 7.2, p´ ag. 426).

2.2.

M´ etodo de Diferencias Finitas

El m´etodo de Diferencias Finitas, al igual que el Disparo, es considerado cl´ asico dentro de la literatura. Este m´etodo, presentado en las referencias [9, 18], se basa en la idea de aproximar las derivadas que aparecen en un problema de contorno o valor de frontera mediante una f´ ormula de diferencias sobre una malla de puntos. Luego se sustituye dichas aproximaciones en la ecuaci´ on diferencial para posteriormente resolver un sistema de ecuaciones lineales o no lineales, dependiendo si el PVF es lineal o no lineal.

2.2.1.

Planteamiento para el caso y ′′ = p(x)y ′ + q(x)y + r(x)

Para el problema de contorno y ′′ = p(x)y ′ + q(x)y + r(x), a ≤ x ≤ b, y(a) = α, y(b) = β, se procede a discretizar el intervalo [a, b] en una malla de puntos a = x0 < x1 < · · · < xn+1 = b, el cual se tomar´ a uniforme, esto es, xi+1 − xi = h, aunque, en general, los puntos no necesitan estar distribuidos de esta manera. En los puntos interiores de la malla (xi = a + ih, i = 1, 2, . . . , n), la ecuaci´ on diferencial a aproximar es y ′′ (xi ) = p(xi )y ′ (xi ) + q(xi )y(xi ) + r(xi ). (2.13) Si se usa el polinomio de Taylor de orden 3 para y(x), alrededor de xi , y se eval´ ua en xi+1 y xi−1 , entonces, suponiendo que y ∈ C 4 [xi−1 , xi+1 ], se tiene y(xi+1 ) = y(xi + h) = y(xi ) + hy ′ (xi ) +

h2 ′′ h3 h4 y (xi ) + y (3) (xi ) + y (4) (ξi+ ), 2 6 24

para alg´ un ξi+ en (xi , xi+1 ), y y(xi−1 ) = y(xi − h) = y(xi ) − hy ′ (xi ) +

h2 ′′ h3 h4 y (xi ) − y (3) (xi ) + y (4) (ξi− ), 2 6 24

para alg´ un ξi− en (xi−1 , xi ). 17

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF Si sumamos adecuadamente estas dos u ´ltimas igualdades, h2 (4) + (4) − y (ξi ) + y (ξi ) , y(xi+1 ) + y(xi−1 ) = 2y(xi ) + h y (xi ) + 24 2 ′′

y al despejar y ′′ (xi ), 1 h2 (4) + (4) − y (ξi ) + y (ξi ) . y (xi ) = 2 y(xi+1 ) − 2y(xi ) + y(xi−1 ) − h 24 ′′

Ahora, aplicando el teorema de valor intermedio, tenemos que 1 h2 ′′ y (xi ) = 2 y(xi+1 ) − 2y(xi ) + y(xi−1 ) − y (4) (ξi ), h 12

(2.14)

para alg´ un ξi en (xi−1 , xi+1 ). A este esquema se le llama f´ ormula de diferencias centradas para y ′′ (xi ). De manera semejante se obtiene una f´ ormula de este tipo para y ′ (xi ), el cual viene dada por 1 h2 ′ y(xi+1 ) − y(xi−1 ) − y (3) (ηi ), (2.15) y (xi ) = 2h 6 para alg´ un ηi en (xi−1 , xi+1 ). Formaci´ on del sistema de ecuaciones. la igualdad

Sustituyendo (2.14) y (2.15) en (2.13), entonces se obtiene

y(xi+1 ) − y(xi−1 ) = p(xi ) + q(xi )y(xi ) + r(xi ) 2h

y(xi+1 ) − 2y(xi ) + y(xi−1 ) h2

h2 (3) (4) − 2p(xi )y (ηi ) − y (ξi ) . 12 El m´etodo de Diferencias Finitas con error de truncamiento de orden O(h2 ) se obtiene empleando esta ecuaci´ on junto con las condiciones de frontera y(a) = α y y(b) = β para definir w0 = α,

wn+1 = β

y

wi+1 − 2wi + wi−1 wi+1 − wi−1 − p(xi ) − q(xi )wi = r(xi ), h2 2h para cada i = 1, 2, . . . , n. Reescribiendo la ecuaci´ on anterior, se tiene que h h 2 − 1 + p(xi ) wi−1 + 2 + h q(xi ) wi − 1 − p(xi ) wi+1 = −h2 r(xi ). 2 2

(2.16)

El sistema de ecuaciones resultante se expresa en la forma Aw = z, donde,



     A=     

γ1

δ1

0

ǫ2

γ2

δ2

0 .. . .. . 0

ǫ3 .. .

γ3 .. . .. . ···

···

··· .. . δ3 .. . .. . 0

··· ..

.

0 .. . .. .

..

.

0

..

.

δn−1 γn

ǫn



     ,     

(2.17)

γi = 2 + h2 q(xi ),

1 ≤ i ≤ n,

δi = −1 +

h p(xi ), 1 ≤ i ≤ n − 1, 2

ǫi = −1 −

h p(xi ), 2 ≤ i ≤ n, 2 18

Secci´ on 2.2: M´etodo de Diferencias Finitas



w1 w2 w3 .. .

    w=    wn−1 wn

        



y

    z=   

−h2 r(x1 ) + 1 + h2 p(x1 ) w0 −h2 r(x2 ) −h2 r(x3 ) .. . −h2 r(xn−1 ) 1 + h2 p(xn ) wn+1

−h2 r(xn ) +



    .   

El siguiente teorema establece bajo cu´ ales condiciones el sistema lineal (2.17) tiene soluci´on u ´nica. Su demostraci´ on es consecuencia del Teorema 1.2. Teorema 2.2 Supongamos que p(x), q(x) y r(x) son continuas en [a, b]. Si q(x) ≥ 0 en [a, b], en2 tonces el sistema lineal tridiagonal (2.17) tiene una soluci´ on u ´nica siempre y cuando h < , donde L L = m´ ax p(x) . a≤x≤b

Demostraci´ on: Ver secci´ on A.2 del Ap´endice A.

Algoritmo del M´ etodo de Diferencias Finitas para el caso lineal. Entrada: extremos a, b; condiciones de frontera α, β; N ≥ 2. Salida: aproximaciones wi a y(xi ), para cada i = 0, . . . , N + 1. 1. Tomar

2. Para i = 2, . . . , N − 1, tomar

3. Tomar

4. Tomar

h x a1 b1 d1

= = = = =

(b − a)/(N + 1), a + h, 2 + h2 q(x), −1 + (h/2)p(x), −h2 r(x) + (1 + (h/2)p(x))α. x ai bi ci di

x = aN = cN = dN =

= = = = =

a + ih, 2 + h2 q(x), −1 + (h/2)p(x), −1 − (h/2)p(x), −h2 r(x).

b − h, 2 + h2 q(x), −1 − (h/2)p(x), −h2 r(x) + (1 − (h/2)p(x))β. l1 = a1 , u1 = b1 /a1 , z1 = d1 /l1 .

(Los pasos 4 - 8 resuelven un sistema lineal tridiagonal utilizando el algoritmo de Crout). 19

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF 5. Para i = 2, . . . , N − 1, tomar

6. Tomar

li = ai − ci ui−1 , ui = bi /li zi = (di − ci zi−1 )/li . lN zN

= aN − cN uN −1 , = (dN − cN zN −1 )/lN .

7. Tomar

w0 = α, wN +1 = β, wN = zN .

8. Para N − 1, . . . , 1, tomar wi = zi − ui wi+1 . 9. Para i = 0, . . . , N + 1, tomar x = a + ih. Mostrar (x, wi ). 10. Parar. An´ alisis del error El siguiente teorema nos ayuda a deducir cotas para el error de la aproximaci´ on de las ecuaciones diferenciales mediante f´ ormula de diferencias finitas de segundo orden. La prueba est´ a detallada en la Secci´ on A.2 del Ap´endice A. Teorema 2.3 Sean Q∗ tal que q(x) ≥ Q∗ > 0, para a ≤ x ≤ b, L = m´ ax p(x) , M3 = m´ ax y (3) (x) y a≤x≤b a≤x≤b M4 = m´ ax y (4) (x) . Entonces, para i = 0, 1, . . . , n + 1, se cumple: a≤x≤b

! M + 2LM 4 3 2 wi − y(xi ) ≤ h . 12Q∗

Demostraci´ on: Ver Ap´endice A, Secci´ on A.2.

El teorema anterior nos dice que la soluci´ on de diferencias finitas converge a la soluci´ on exacta cuando h → 0, y el error es O(h2 ). Observaci´ on. Cuando, en la ecuaci´ on diferencial lineal, p(x) = 0, entonces se puede demostrar que el error cometido es del orden O(h4 ).

2.2.2.

Planteamiento para el caso y ′′ = f (x, y, y ′)

Para el caso del problema no lineal general con valor de frontera: y ′′ = f (x, y, y ′ ), a ≤ x ≤ b, y(a) = α, y(b) = β, el m´etodo de Diferencias Finitas se parece al que se aplic´ o en la subsecci´ on anterior para el caso lineal. Sin embargo, el sistema de ecuaciones ser´ a no lineal y, por tanto, se requiere un proceso iterativo para resolverlo. Nuevamente, se procede a discretizar el intervalo [a, b] en una malla de puntos a = x0 < x1 < · · · < xn+1 = b, que tomaremos equiespaciada o uniforme, es decir, xi+1 − xi = h. En los puntos interiores de la malla (xi = a + ih, i = 1, 2, . . . , n), la ecuaci´ on diferencial a aproximar

es: y ′′ (xi ) = f (xi , y(xi ), y ′ (xi )).

(2.18) 20

Secci´ on 2.2: M´etodo de Diferencias Finitas Formaci´ on del sistema de ecuaciones. Sustituyendo (2.14) y (2.15) en (2.18), obtendremos: y(xi+1 ) − 2y(xi ) + y(xi−1 ) y(xi+1 ) − y(xi−1 ) h2 (3) h2 (4) − y y (ξi ), = f x , y(x ), (η ) + i i i h2 2h 6 12 para alg´ un ξi y ηi en el intervalo (xi−1 , xi+1 ), con i = 1, 2, . . . , n. El m´etodo de Diferencias Finitas se obtiene empleando esta ecuaci´ on junto con las condiciones de frontera y(a) = α y y(b) = β para definir w0 = α,

wn+1 = β,

y

wi+1 − 2wi + wi−1 wi+1 − wi−1 = 0, − f xi , wi , h2 2h para cada i = 1, 2, . . . , n. El sistema no lineal de n × n obtenido con este m´etodo Fi (w1 , . . . , wn ) = 0,

(2.19)

con

wi+1 − wi−1 , i = 1, . . . , n, Fi (w1 , . . . , wn ) = wi+1 − 2wi + wi−1 − h f xi , wi , 2h ∂f 2 tiene soluci´ on u ´nica si h < , donde L es tal que ′ (x, y, y ′ ) ≤ L (ver [13], p´ ag. 86). L ∂y 2

Resoluci´ on del sistema no lineal.

Para resolver el sistema no lineal (2.19), se usa el m´etodo de (0) (0) T Newton para sistemas de ecuaciones no lineales. Si tomamos inicialmente w(0) = w1 , . . . , wn , la iteraci´ on mediante este m´etodo es: w(j+1) = w(j) − J −1 w(j) F w(j) , donde J w(j) es la matriz jacobiana. En la pr´ actica es m´ as conveniente a efectos computacionales iterar resolviendo las ecuaciones: J w(j) w(j+1) − w(j) = −F w(j) . h i ∂Fp Calculemos el jacobiano J(w(j) ) = , ∂wq p,q ∂f wk+1 − wk−1 2 [J]k,k = −2 − h xk , wk , , ∂y 2h [J]k,k+1

[J]k+1,k

h ∂f wk+1 − wk−1 = 1− xk , wk , , 2 ∂y ′ 2h h ∂f wk+1 − wk−1 , = 1+ xk , wk , 2 ∂y ′ 2h

con k = 1, . . . , n. El resto de las entradas de la matriz J son nulos. En particular, como todos los elementos con |p − q| > 1 son nulos, entonces J es tridiagonal. En cada iteraci´ on del sistema lineal: J w(j) v = −F w(j) , con v = (v1 , . . . , vn )T , se requiere obtener v, porque

w(j+1) = w(j) + v. Puesto que J es tridiagonal, se puede utilizar un m´etodo para resolver sistemas de ecuaciones lineales. Para efecto de la implementaci´ on del c´ odigo, se usar´ a el algoritmo de reducci´ on de Crout para sistemas tridiagonales (ver referencia [18], Cap´ıtulo 6, Secci´ on 6.6, p´ ag. 408). 21

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF Algoritmo del M´ etodo de Diferencias Finitas para el caso no lineal. Entrada: extremos a, b; condiciones de frontera α, β; entero N ≥ 2; tolerancia Tol, n´ umero m´ aximo de iteraciones M . Salida: aproximaciones wi a y(xi ), para cada i = 0, . . . , N +1, o un mensaje de que se excedi´ o el n´ umero m´ aximo de iteraciones. 1. Tomar h = (b − a)/(N + 1); w0 = α; wN +1 = β. β−α 2. Para i = 1, . . . , N , tomar wi = α + i h. b−a

3. Tomar k = 1.

4. Mientras k ≤ M , hacer los pasos 5 - 16. 5. Tomar

x t a1 b1 d1

6. Para i = 2, . . . , N − 1, tomar x t ai bi ci di

7. Tomar

= = = = = =

= = = = =

a + h, (w2 − α)/2h, 2 + h2 fy (x, w1 , t), −1 + (h/2)fy′ (x, w1 , t), −(2w1 − w2 − α + h2 f (x, w1 , t)).

a + ih, (wi+1 − wi−1 )/2h, 2 + h2 fy (x, wi , t), −1 + (h/2)fy′ (x, wi , t), −1 − (h/2)fy′ (x, wi , t), −(2wi − wi+1 − wi−1 + h2 f (x, wi , t)).

x = b − h, t = (β − wN −1 )/2h, aN = 2 + h2 fy (x, wN , t), cN = −1 − (h/2)fy′ (x, wN , t), dN = −(2wN − wN −1 − β + h2 f (x, wN , t)).

8. Tomar l1 = a1 ; u1 = b1 /a1 ; z1 = d1 /l1 . (Los pasos 8-12 resuelven un sistema lineal tridiagonal utilizando el algoritmo de Crout). 9. Para i = 2, . . . , N − 1, tomar li = ai − ci ui−1 ; ui = bi /li ; zi = (di − ci zi−1 )/li . 10. Tomar lN = aN − cN uN −1 ; zN = (dN − cN zN −1 )/lN . 11. Tomar vN = zN ; wN = wN + vN . 12. Para i = N − 1, . . . , 1, tomar vi = zi − ui vi+1 ; wi = wi + vi . 13. Si ||v|| ≤ Tol, hacer los pasos 14-15. 14. Para i = 0, . . . , N + 1, tomar x = a + ih. Mostrar (x, wi ). 15. Parar. Procedimiento terminado con ´exito. 16. Tomar k = k + 1. 17. M´ aximo de iteraciones, procedimiento terminado sin ´exito. Parar. Al igual que el caso lineal, la soluci´ on encontrada por el m´etodo de Diferencias Finitas aproxima a la soluci´ on exacta con un orden de aproximaci´ on O(h2 ) (ver [9], Cap´ıtulo 8, Secci´ on 7.2, p´ ag. 432 - 433). 22

Secci´ on 2.3: M´etodo de los Elementos Finitos

0

1

Figura 2.4: La funci´ on y(x) = −3x2 + 3x es un ejemplo de una funci´ on en V.

2.3.

M´ etodo de Elementos Finitos

El m´etodo de Elementos Finitos (MEF), el cual se plantea en las referencias [3, 17], constituye hoy en d´ıa el procedimiento habitual para la aproximaci´ on num´erica de la soluci´ on de problemas pr´ acticos que surgen en Ingenier´ıa o las Ciencias. Se trata de un m´etodo general para la soluci´ on de problemas de contorno de ecuaciones diferenciales ordinarias o parciales. En esencia se trata de una t´ecnica que sustituye el problema diferencial por otro algebraico, aproximadamente equivalente, para el cual se conocen t´ecnicas generales de resoluci´ on. Para ello hace uso de la discretizaci´ on o subdivisi´ on de una regi´ on sobre la cual est´ an definidas las ecuaciones en formas geom´etricas simples denominadas elementos finitos. Una de las ventajas de este m´etodo es su facilidad de implementaci´ on en un programa computacional, que a su vez es una condici´ on b´ asica para su utilizaci´ on ya que el tratamiento de un problema, en particular, debe efectuarse en un n´ umero muy elevado de operaciones para resolver sistemas algebraicos del orden de cientos o miles de ecuaciones. No obstante, esta cantidad no es una limitaci´ on con las computadoras est´ andar de hoy. El MEF fue al principio desarrollado en 1943 por R. Courant, qui´en utiliz´ o el m´etodo de Ritz de an´ alisis num´erico y minimizaci´ on de las variables de c´ alculo para obtener soluciones aproximadas a un sistema de vibraci´ on. En la d´ecada de 1960 el m´etodo fue generalizado para la soluci´ on aproximada de problemas de an´ alisis de tensi´ on, flujo de fluidos y transferencia de calor. El primer libro sobre Elementos Finitos fue publicado en 1967 por Zienkiewicz y Cheung ([19]). En la d´ecada de 1970 el m´etodo fue extendido al an´ alisis de problemas no lineales de la mec´ anica del continuo. Hoy el m´etodo permite resolver pr´ acticamente cualquier situaci´ on f´ısica que pueda formularse mediante un sistema de ecuaciones diferenciales.

2.3.1.

Formulaci´ on d´ ebil de un problema modelo unidimensional

Consideremos el problema de valor de frontera siguiente: −u′′ (x) = f (x),

(2.20)

u(0) = u(1) = 0, dv donde v ′ = y f es una funci´ on continua. Se introduce el siguiente espacio vectorial de funciones: dx V = v : v ∈ C[0, 1], v ′ es continua a trozos y acotada en [0, 1], v(0) = v(1) = 0 .

Si se multiplica la ecuaci´ on (2.20) por una funci´ on cualquiera de V (ver Figura 2.4) e integramos en el intervalo [0, 1] usando la f´ ormula de integraci´ on por partes, entonces se tiene: Z 1 Z 1 Z 1 ′′ ′ ′ ′ ′ − u (x)v(x)dx = −u (1)v(1) + u (0)v(0) + u (x)v (x)dx = u′ (x)v ′ (x)dx, 0

0

0

23

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF b

b

b

b

b

0

1

Figura 2.5: Ejemplo de una funci´ on vh ∈ Vh , 5 subintervalos.

y concluimos que para toda funci´ on v ∈ V, la soluci´ on u del problema modelo (2.20) verifica: B(u, v) = F(v), donde B(u, v) =

Z

0

1

′

′

u (x)v (x)dx es una forma bilineal y F(v) =

(2.21) Z

1

f (x)v(x)dx es un funcional lineal.

0

La formulaci´ on anterior (2.21) se llama formulaci´ on d´ebil o variacional del problema de partida (2.20). Las dos formulaciones no son estrictamente equivalentes. En efecto, toda soluci´ on de (2.20) es soluci´ on de (2.21). Sin embargo, para que toda soluci´ on de (2.21) sea soluci´ on de (2.20), se debe exigir que dicha soluci´ on sea dos veces continuamente diferenciable. La existencia y unicidad de la soluci´ on de (2.21) queda garantizada en virtud del Teorema de Lax Milgram y el Teorema de Representaci´ on de Riesz (Ver Preliminares, Teor´ıa de An´ alisis Funcional).

2.3.2.

El MEF para el problema modelo con funciones lineales a trozos

El espacio V es un espacio de funciones de dimensi´ on infinita. La idea del m´etodo de los elementos finitos es buscar soluciones de (2.21) en un subespacio de V m´ as sencillo, en particular, en un espacio vectorial de funciones que sea de dimensi´ on finita. Esto nos permitir´ a representar cualquier funci´ on de este subespacio como una combinaci´ on lineal de elementos de una base. En primer lugar, construiremos un subespacio Vh de V. Para ello sea 0 = x0 < x1 < · · · < xM < xM +1 = 1, una partici´ on del intervalo (0, 1) en subintervalos Ij = (xj−1 , xj ) de longitud hj = xj − xj−1 , j = 1, 2, . . . , M + 1, y sea h = m´ ax hj . La cantidad h es una medida de lo fina que es la partici´ on. Consideremos ahora Vh = {v : v es lineal en cada subintervalo Ij , v ∈ C[0, 1] y v(0) = v(1) = 0}. Un ejemplo de una funci´ on en Vh se aprecia en la Figura 2.5. Una base del espacio Vh est´ a constituido por el siguiente conjunto de funciones ϕj ∈ Vh , j = 1, . . . , M , definidas por: 1 si i = j, ϕj (xi ) = 0 si i 6= j. es decir, ϕj es continua y lineal a trozos, y toma el valor 1 en el nodo xj y el valor 0 en los otros nodos. 24

Secci´ on 2.3: M´etodo de los Elementos Finitos 1

b

b

b

b

b

b

b

b

b

b

0

b

1

Figura 2.6: Ejemplo de una funci´ on de la base de Vh , 10 subintervalos.

De forma m´ as precisa:

ϕj (x) =

  0        x − xj−1      xj − xj−1   xj+1 − x     xj+1 − xj        0

si x ≤ xj−1 , si xj−1 < x ≤ xj , si xj < x < xj+1 , si xj+1 ≤ x.

Una funci´ on de esta base se aprecia en la Figura 2.6.

Una funci´ on vh ∈ Vh se escribe de forma u ´nica como combinaci´ on lineal de funciones de la base M {ϕj }j=1 , M X vh (x) = vh (xi )ϕj (x), x ∈ [0, 1]. j=1

El m´etodo de los elementos finitos para el problema (2.20) se formula de la siguiente manera: encontrar uh ∈ Vh tal que para todo vh ∈ Vh verifique B(uh , vh ) = F(vh ).

(2.22)

Este problema, como vamos a ver, es equivalente a resolver un sistema algebraico lineal. En efecto, por una parte la soluci´ on uh (x) se expresa en funci´ on de los elementos de la base {ϕi }M i=1 , esto es, uh (x) =

M X

ui ϕi (x),

i=1

donde ui = uh (xi ). Por otro lado, para que se verifique (2.22) con cualquier vh (x), es necesario y suficiente que se verifique para cualquier funci´ on de la base. Teniendo en cuenta la linealidad de la integral, resulta que el problema a resolver es: M X i=1

B(ϕi (x), ϕj (x))ui = F(ϕj ),

para j = 1, . . . , M .

En forma matricial, el sistema se puede escribir como: Au = b,

(2.23) 25

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF donde [A]ij = B(ϕi (x), ϕj (x)) y bj = F(ϕj ).

Los elementos [A]ij se calculan f´ acilmente. Primero notemos que si |i − j| > 1, entonces [A]ij = 0, pues, en este caso, para x ∈ [0, 1], ϕi (x) o ϕj (x) es igual a cero. Por lo tanto, la matriz A es tridiagonal, esto es, u ´nicamente los elementos de la diagonal principal y los elementos de las dos diagonales adyacentes pueden ser diferentes de cero. As´ı, para j = 1, . . . , M , Z 1 ϕ′j (x)ϕ′j (x)dx B(ϕj (x), ϕj (x)) = 0

=

Z

xj

xj−1

=

ϕ′j (x)ϕ′j (x)dx +

Z

xj+1

xj

ϕ′j (x)ϕ′j (x)dx

1 1 + , hj hj+1

y para j = 2, . . . , M , B(ϕj (x), ϕj−1 (x)) = =

Z

1

0

Z

ϕ′j (x)ϕ′j−1 (x)dx

xj

xj−1

= −

1 dx h2j

1 . hj

La matriz A es sim´etrica y definida positiva, esto es, xAxT > 0, para todo vector fila x de dimensi´ on M no nulo. En efecto, notemos que [A]ij = B(ϕi (x), ϕj (x)) = B(ϕj (x), ϕi (x)) = [A]ji . Por otro lado, para un vector fila x arbitrario no nulo, tenemos ! M M X X vp B(ϕp (x), ϕk (x)) vk xAxT = p=1

k=1

=

M X k=1

= B

B M X

M X

!

vp ϕp (x), ϕk (x) vk

p=1

vp ϕp (x),

p=1

M X k=1

!

vk ϕk (x) ,

= B(e v (x), ve(x)) > 0.

As´ı, como A es definida positiva, entonces A es invertible, por tanto, el sistema (2.23) tiene soluci´on. En el caso especial en que los subintervalos Ij tengan todos la misma longitud h = M1+1 , el sistema tiene la forma:      2 −1 0 . . . 0 u1 b1  .      −1 2 −1 . . . ..        1     . .  . . . . . =  . . .. .. .. .       h       −1   uM bM 0 0 . . . −1 2 26

Secci´ on 2.3: M´etodo de los Elementos Finitos Observaci´ on: Si consideramos la ecuaci´ on diferencial u′′ (x) = −f (x), con u(0) = u(1) = 0, entonces la matriz A, del m´etodo de Diferencias Finitas para el caso lineal, coincide con la matriz A del m´etodo de los Elementos Finitos.

2.3.3.

Integraci´ on num´ erica

A la hora de calcular el segundo miembro del sistema (2.23), aparecen integrales de la forma Z

xj

f (x)ϕj (x)dx + xj−1

Z

xj+1

f (x)ϕj (x)dx,

(2.24)

xj

con j = 1, . . . , M . Salvo en el caso en que la funci´ on f sea una funci´ on constante a trozos, la forma pr´ actica de calcular estas integrales ser´ a mediante la integraci´ on num´erica. Para efecto de este trabajo, se usar´ a la f´ ormula de cuadratura Gaussiana, espec´ıficamente la de dos puntos (Ver [12], Cap´ıtulo 5, Secci´ on 5.3, p´ ag. 270). Dicha f´ ormula dice que Z Por otro lado, una integral

Z

1 −1

f (ξ)dξ ≈ f

1 −√ 3

1 √ +f . 3

b

g(x)dx en un intervalo arbitrario [a, b] se puede transformar en otra en

a

[−1, 1], usando el cambio de variables t= el cual equivale a decir que

2x − a − b , b−a

1 x= (b − a)t + a + b . 2

Esto permite aplicar la f´ ormula de cuadratura Gaussiana a cualquier intervalo [a, b], ya que Z

b a

b−a g(x)dx = 2

Z

(b − a)t + b + a g dt. 2 −1 1

Por comodidad pongamos ψj (x) = f (x)ϕj (x), j = 1, . . . , M . Para (2.24), se cumple bj =

Z

xj

ψj (x)dx +

xj−1

=

≈

1 donde pe = √ . 3

hj 2

Z

1

−1

ψj

Z

xj+1

ψj (x)dx xj

Z hj t + xj + xj−1 hj+1 1 hj+1 t + xj+1 + xj dt + ψj dt 2 2 2 −1

hj −hj pe + xj + xj−1 hj pe + xj + xj−1 ψj + ψj 2 2 2

hj+1 −hj+1 pe + xj+1 + xj hj+1 pe + xj+1 + xj + ψj + ψj , 2 2 2

27

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF Descripci´ on de un algoritmo del M´ etodo de Elementos Finitos. del m´etodo de Elementos Finitos requiere, en general, de cuatro etapas:

El desarrollo de un algoritmo

1. El problema debe reformularse en forma variacional. 2. El dominio de la variable independiente debe dividirse mediante una partici´ on en subdominios, llamados elementos finitos. Asociada a la partici´ on anterior se construye un espacio vectorial de dimensi´ on finita, llamado espacio de elementos finitos, siendo la soluci´ on num´erica una combinaci´ on lineal en dicho espacio vectorial. 3. Se obtiene la proyecci´ on del problema variacional original sobre el espacio de elementos finitos obtenido de la partici´ on. Esto da lugar a un finito sistema de ecuaciones con un n´ umero elevado de inc´ ognitas. El n´ umero de ecuaciones y de inc´ ognitas ser´ a igual a la dimensi´ on del espacio vectorial. En general, cuanto mayor sea dicha dimensi´ on tanto mejor ser´ a la aproximaci´ on num´erica obtenida. 4. El u ´ltimo paso es el c´ alculo num´erico de la soluci´ on del sistema de ecuaciones. Es importante resaltar que la soluci´ on estimada por el m´etodo de Elementos Finitos aproxima la soluci´ on exacta con un orden de aproximaci´ on O(h2 ) (ver [17]).

2.4.

M´ etodo Galerkin Discontinuo

El m´etodo Galerkin Discontinuo, presentado en la referencia [1], forma parte de una clase de m´etodos num´ericos para solucionar ecuaciones diferenciales parciales y ordinarias. Combina muchas caracter´ısticas del m´etodo de Elementos Finitos y se ha aplicado con ´exito a problemas hiperb´ olicos, el´ıpticos y parab´ olicos. Este m´etodo fue propuesto y analizado en los a˜ nos 70 como una t´ecnica para solucionar num´ericamente ecuaciones diferenciales parciales que surgen en muchas aplicaciones. Este m´etodo num´erico sigue la misma estructura del m´etodo de Elementos Finitos. La idea consiste en reformular el problema dado en forma d´ebil o forma variacional. Se procede a dividir el dominio de la variable independiente en subdominios. Asociada a la partici´ on anterior se construye un espacio vectorial de dimensi´ on finita, siendo la soluci´ on num´erica una combinaci´ on lineal en dicho espacio vectorial, lo que conllevar´ a a un sistema de ecuaciones lineales, es decir, encontrar la soluci´ on de una ecuaci´ on diferencial se limita a encontrar la soluci´ on de un sistema de ecuaciones lineales.

2.4.1.

Formulaci´ on d´ ebil de un problema modelo unidimensional

Consideremos el problema siguiente en el intervalo (0, 1) −(K(x)p′ (x))′ = f (x), p(0) = 1, p(1) = 0,

(2.25)

donde K ∈ C 1 (0, 1), es positiva y acotada, y f ∈ C 0 (0, 1). Vamos a reformular el problema modelo, pero antes se introducir´ a algunos conceptos que ser´ an u ´tiles para todo lo que sigue. Sea 0 = x0 < x1 < · · · < xN = 1 una partici´ on εh de (0, 1), pongamos In = (xn , xn+1 ), y definamos hn = xn+1 − xn , hn−1,n = m´ ax(hn−1 , hn ), h =

m´ ax

0≤n≤N −1

hn .

Denotemos por Dk (εh ) el espacio de los polinomios discontinuos a trozos de grado k: Dk (εh ) = v : v|In ∈ Pk (In ), ∀n = 0, . . . , N − 1 ,

28

Secci´ on 2.4: M´etodo Galerkin Discontinuo 1

0

0.25

0.5

0.75

1

Figura 2.7: Ejemplo de una funci´ on v en D2 (εh ), 4 subintervalos de igual longitud. En I0 e I2 v es un abola cuadr´ atica. segmento de recta, mientras que en I1 e I3 v es una par´

donde Pk (In ) es el espacio de los polinomios de grado k en el intervalo In . Sean v(x+ ım v(xn + ǫ) n ) = l´ ǫ→0

y v(x− ım v(xn − ǫ), donde ǫ > 0. Definiremos el salto y promedio de v en los extremos de In , n ) = l´ ǫ→0 respectivamente, como 1 + − + [v(xn )] = v(x− n ) − v(xn ), {v(xn )} = (v(xn ) + v(xn )), n = 1, . . . , N − 1, 2 mientras que el salto y promedio en los extremos del intervalo (0, 1) se define como: + − − [v(x0 )] = −v(x+ 0 ), {v(x0 )} = v(x0 ), [v(xN )] = v(xN ), {v(xN )} = v(xN ).

Ahora, se introducen los t´erminos de salto de la soluci´ on y su derivada. Estos t´erminos vienen dados por J0 (v, w) =

N X

σ0

h n=0 n−1,n

[v(xn )][w(xn )],

J1 (v, w) =

N −1 X n=1

σ1 hn−1,n

[v ′ (xn )][w′ (xn )],

donde σ 0 y σ 1 son dos n´ umeros reales no negativos. Ahora bien, ya estamos en condici´ on de reformular el problema modelo. Sea v una funci´on en Dk (εh ) (Ver Figura 2.7). Si se multiplica la ecuaci´ on (2.25) por v y se integra en In , con n = 0, . . . , N − 1, usando la f´ ormula de integraci´ on por partes, entonces Z xn+1 Z xn+1 − ′ ′ ′ ′ + K(x)p (x)v (x)dx − K(xn+1 )p (xn+1 )v(xn+1 ) + K(xn )p (xn )v(xn ) = f (x)v(x)dx. xn

xn

Sumando las N ecuaciones anteriores y usando la definici´ on de salto para v en cada nodo de la partici´ on, se tiene Z 1 N −1 Z xn+1 N X X K(x)p′ (x)v ′ (x)dx − [K(xn )p′ (xn )v(xn )] = f (x)v(x)dx. n=0

xn

0

n=0

Notemos que, para n = 1, . . . , N − 1, se cumple

[K(xn )p′ (xn )v(xn )] = {K(xn )p′ (xn )}[v(xn )] + {v(xn )}[K(xn )p′ (xn )].

(2.26)

Si se usa (2.26) y sabiendo que la soluci´ on exacta p satisface [K(xn )p′ (xn )] = 0 en los nodos interiores de la partici´ on εh (esto gracias a la continuidad de p′ ), se obtiene N −1 Z xn+1 X n=0

xn

′

′

K(x)p (x)v (x)dx −

N X

′

{K(xn )p (xn )}[v(xn )] =

n=0

Z

1

f (x)v(x)dx.

0

29

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF Como la soluci´ on exacta p es tambi´en continua, entonces [p(xn )] = 0. Notemos que p(x0 ) = p(0) = 1 y p(xN ) = p(1) = 0, con lo que N −1 Z xn+1 X n=0

xn

K(x)p′ (x)v ′ (x)dx − =

=

Z Z

1 0 1 0

N X

{K(xn )p′ (xn )}[v(xn )] + ǫ

n=0

N X

{K(xn )v ′ (xn )}[p(xn )]

n=0

f (x)v(x)dx − ǫK(x0 )v ′ (x0 )p(x0 ) + ǫK(xN )v ′ (xN )p(xN ) f (x)v(x)dx − ǫK(x0 )v ′ (x0 ).

El tercer t´ermino del primer miembro de la igualdad anterior es casi siempre cero. Adem´ as, ǫ puede ser cualquier n´ umero real, sin embargo, nos restringiremos al caso en que ǫ = −1, ǫ = 0 o ǫ = 1. As´ı, definimos la forma bilineal del m´etodo Galerkin Discontinuo Bǫ : Dk (εh ) × Dk (εh ) → R, como N −1 Z xn+1 X

Bǫ (w, v) =

n=0

xn

+ǫ

N X

′

′

K(x)w (x)v (x)dx −

N X

{K(xn )w′ (xn )}[v(xn )]

n=0

{K(xn )v ′ (xn )}[w(xn )] + J0 (w, v) + J1 (w, v).

n=0

La forma bilineal Bǫ tiene las siguientes propiedades: 1. Para ǫ = −1, la forma es sim´etrica, y adem´ as B−1 (v, v) =

N −1 Z xn+1 X n=0

xn

K(x)(v ′ (x))2 dx − 2

N X

n=0

{K(xn )v ′ (xn )}[v(xn )] + J0 (v, v) + J1 (v, v).

2. Para ǫ = 0 o ǫ = 1, la forma es no sim´etrica, y cumple B1 (v, v) =

N −1 Z xn+1 X

B0 (v, v) =

N −1 Z xn+1 X

n=0

n=0

xn

xn

K(x)(v ′ (x))2 dx + J0 (v, v) + J1 (v, v) ≥ 0,

′

2

K(x)(v (x)) dx −

N X

{K(xn )v ′ (xn )}[v(xn )] + J0 (v, v) + J1 (v, v).

n=0

La idea de este m´etodo consiste en encontrar Pe ∈ Dk (εh ) tal que Bǫ (Pe, v) = L(v),

(2.27)

para cada v ∈ Dk (εh ), donde L : Dk (εh ) → R es el funcional lineal definido por Z 1 σ0 L(v) = f (x)v(x)dx − ǫK(x0 )v ′ (x0 ) + v(x0 ). h0,1 0 La existencia y unicidad de la soluci´ on de (2.27) queda garantizada en virtud del Teorema de Lax Milgram y el Teorema de Representaci´ on de Riesz (Ver Preliminares, Teor´ıa de An´ alisis Funcional). 0 1 Dependiendo de la elecci´ on de los par´ ametros ǫ, σ y σ , el m´etodo Galerkin Discontinuo recibe varios nombres. 30

Secci´ on 2.4: M´etodo Galerkin Discontinuo Si ǫ = −1, σ 1 = 0 y σ 0 est´ a acotado inferiormente por un n´ umero de gran magnitud, el m´etodo resultante es llamado symmetric interior penalty Galerkin (SIPG), introducido despu´es de 1970 por Wheeler (ver [16]) y Arnold (ver [8]). Si ǫ = −1 y σ 0 = σ 1 = 0, el m´etodo resultante es llamado global element, introducido en 1979 por Delves y Hall (ver [15]). Sin embargo, la matriz asociada a la forma bilineal es indefinida, la parte real de los autovalores no son todos positivos y, por tanto, el m´etodo es inestable. Si ǫ = 1, σ 1 = 0 y σ 0 = 1, el m´etodo resultante es llamado nonsymmetric interior penalty Galerkin (NIPG), introducido en 1999 por Rivi`ere, Wheeler and Girault (ver [4]). Si ǫ = +1 y σ 0 = σ 1 = 0, el m´etodo resultante fue introducido por Oden, Babu˘ska y Baumann en 1998 (ver [10]). Haremos referencia a este m´etodo como NIPG 0, ya que corresponde al caso particular de NIPG, con σ 0 = 0. Si ǫ = 0, obtenemos el m´etodo incomplete interior penalty Galerkin (IIPG), el cual fue introducido por Dawson, Sun and Wheeler en 2004 (ver [5]). Cabe preguntarnos, ¿qu´e pasa si ǫ = 0 y σ 0 = σ 1 = 0?. Resulta que el m´etodo no es convergente y es inestable, m´ as a´ un, no se puede probar la existencia y unicidad de la soluci´ on.

2.4.2.

Derivaci´ on del sistema lineal Aα = b mediante el uso de polinomios discontinuos a trozos

En este apartado se deriva el sistema lineal que se obtiene al emplear el esquema Galerkin Discontinuo para el problema modelo, tomando el caso simple cuando K(x) es id´enticamente igual a 1 y σ 1 = 0. Al mismo tiempo, y por efectos de c´ alculo, se consideran los polinomios cuadr´ aticos discontinuos a trozos, es decir, cuando se trabaja en D2 (εh ). Una base del espacio P2 (In ) est´ a contituido por el siguiente conjunto de funciones φnj ∈ P2 (In ), j = 0, 1, 2, definidas por φn0 (x) = 1,

φn1 (x) = 2

x−x en , xn+1 − xn

φn2 (x) = 4

(x − x en )2 , (xn+1 − xn )2

1 xn +xn+1 es el punto medio del intervalo In . Notemos que φn0 , φn1 y φn2 son las traslaciones 2 respectivas de las funciones g0 (x) = 1, g1 (x) = x y g2 (x) = x2 , desde el intervalo (−1, 1) al intervalo In (Ver Figura 2.8). Con efecto de simplificar un poco los c´ alculos, supongamos que εh es una partici´ on uniforme de (0, 1). Por tanto, las funciones de base local y sus derivadas se pueden simplificar 2 2 1 4 1 n n n φ0 (x) = 1, φ1 (x) = x− n+ h , φ2 (x) = 2 x − n + h , (2.28) h 2 h 2 2 8 1 n ′ n ′ n ′ h . (2.29) (φ0 ) (x) = 0, (φ1 ) (x) = , (φ2 ) (x) = 2 x − n + h h 2 donde x en =

Las funciones de base global {Φni } para el espacio D2 (εh ) son obtenidas de las funciones de base local extendi´endolas a cero: n φi (x), x ∈ In , n Φi (x) = 0, x 6∈ In .

Ahora, para cada x ∈ (0, 1), notemos que la soluci´ on aproximada Pe viene dada por Pe (x) =

N −1 X 2 X

m αm i Φi (x),

(2.30)

m=0 i=0

31

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF φn0 (x)

1

φn2 (x) 0

xn+1

xn

1

φn1 (x) −1

Figura 2.8: Funciones de base local para el espacio P2 (In ).

donde los coeficientes αm i son desconocidos y requieren ser calculados. Si se sustituye (2.30) en (2.27), entonces, gracias a la linealidad de Bǫ y L, se tiene N −1 X 2 X m=0 i=0

m n n αm i Bǫ (Φi (x), Φj (x)) = L(Φj (x)),

para n = 0, . . . , N − 1 y j = 0, 1, 2. La igualdad anterior no es m´ as que un sistema lineal de la forma m n Aα = b, donde α es el vector con componentes αj , A es la matriz con entradas Bǫ (Φm i (x), Φj (x)), y b es n el vector con componentes L(Φj (x)).

2.4.3.

Construcci´ on de la matriz A

Las entradas de la matriz global A pueden ser obtenidas calculando y ensamblando matrices locales, por tanto, vamos a describir c´ omo calcular dichas matrices. Se reagrupar´ an los t´erminos de la definici´ on Bǫ en tres grupos: los t´erminos que involucran las integrales sobre In , los t´erminos que involucran los nodos interiores xn , y los relacionados con los nodos de la frontera x0 y xN . Primero se considera los t´erminos correspondientes a las integrales sobre los intervalos In . En cada elemento In , la soluci´ on Pe es un polinomio de grado 2, y se puede expresar como Pe(x) = αn0 φn0 (x) + αn1 φn1 (x) + αn2 φn2 (x),

(2.31)

para cada x ∈ In . Usando (2.31) y elegiendo v = φnj (x), para cada j = 0, 1, 2, se obtiene Z

In

Pe′ (x)(φnj )′ (x)dx =

2 X

αni

i=0

Z

In

(φni )′ (x)(φnj )′ (x)dx,

Este sistema lineal puede ser reescrito como An αn , donde  n  Z α0 αn =  αn1  , [An ]ij = (φni )′ (x)(φnj )′ (x)dx. In αn2 La matriz An es f´ acil de determinar,



0 0 1 0 4 An =  h 0 0

 0 0  . 16 3 32

Secci´ on 2.4: M´etodo Galerkin Discontinuo En segundo lugar, se considera los t´erminos que involucran los nodos interiores xn . Usando la definici´ on promedio y salto de la funci´ on v para dichos nodos, se tiene

donde

σ0 − Pe ′ (xn ) v(xn ) + ǫ v ′ (xn ) Pe(xn ) + Pe(xn ) v(xn ) = bn + cn + dn + en , h bn =

1 e′ + ǫe + ′ + σ0 e + P (xn )v(x+ P (x P (xn )v(x+ ) − )v (x ) + n n n n ), 2 2 h

1 ǫe − ′ − σ0 e − − cn = − Pe′ (x− P (xn )v(x− n )v(xn ) + P (xn )v (xn ) + n ), 2 2 h

1 ǫe + ′ − σ0 e + − dn = − Pe′ (x+ P (xn )v(x− n )v(xn ) − P (xn )v (xn ) − n ), 2 2 h en =

1 e′ − ǫe − ′ + σ0 e − P (xn )v(x+ ) + P (x )v (x ) − P (xn )v(x+ n n n n ). 2 2 h

an las matrices Usando (2.31) y con el cambio v = φnj (x), los cuatros t´erminos anteriores determinar´ locales Bn , Cn , Dn y En , respectivamente. Las entradas de dichas matrices son 1 n ′ + n + ǫ σ0 n + n + n ′ + (φi ) (xn )φj (xn ) − φni (x+ φ (x )φ (x ), n )(φj ) (xn ) + 2 2 h i n j n

[Bn ]ij

=

[Cn ]ij

= −

0 1 n ′ − n − ǫ n )′ (x− ) + σ φn (x− )φn (x− ), (φi ) (xn )φj (xn ) + φni (x− )(φ n n j 2 2 h i n j n

[Dn ]ij

= −

1 n ′ + n − ǫ σ0 n + n − n ′ − (φi ) (xn )φj (xn ) − φni (x+ φ (x )φ (x ), n )(φj ) (xn ) − 2 2 h i n j n

[En ]ij

=

0 1 n ′ − n + ǫ n )′ (x+ ) − σ φn (x− )φn (x+ ). (φi ) (xn )φj (xn ) + φni (x− )(φ n n j 2 2 h i n j n

Las entradas ij de cada matriz local 3 × 3 se puede calcular f´ acilmente usando (2.28) y (2.29): Bn =

Cn =

Dn =

En =

  σ0 1 − σ0 −2 + σ 0 1 −ǫ − σ 0 −1 + ǫ + σ 0 2 − ǫ − σ0  , h 0 0 2ǫ + σ 1 − 2ǫ − σ −2 + 2ǫ − σ 0 

 σ0 −1 + σ 0 −2 + σ 0 1 ǫ + σ 0 −1 + ǫ + σ 0 −2 + ǫ + σ 0  , h 2ǫ + σ 0 −1 + 2ǫ + σ 0 −2 + 2ǫ + σ 0   −σ 0 −1 + σ 0 2 − σ0 1 −ǫ − σ 0 −1 + ǫ + σ 0 2 − ǫ − σ 0  , h −2ǫ − σ 0 −1 + 2ǫ + σ 0 2 − 2ǫ − σ 0 

 −σ 0 1 − σ0 2 − σ0 1 ǫ + σ0 −1 + ǫ + σ 0 −2 + ǫ + σ 0  . h 0 −2ǫ − σ 1 − 2ǫ − σ 0 2 − 2ǫ − σ 0

33

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF Para finalizar, se considera los t´erminos relacionados con los nodos de frontera x0 y xN : σ0 f0 = Pe′ (x0 )v(x0 ) − ǫPe (x0 )v ′ (x0 ) + Pe(x0 )v(x0 ), h

fN

σ0 = −Pe ′ (xN )v(xN ) + ǫPe(xN )v ′ (xN ) + Pe(xN )v(xN ). h

Los dos t´erminos anteriores determinar´ an, respectivamente, dos matrices locales, F0 y FN , las cuales vienen dadas por:   σ0 2 − σ0 −4 + σ 0 1 F0 = −2ǫ − σ 0 −2 + 2ǫ + σ 0 4 − 2ǫ − σ 0  , h 4ǫ + σ 0 2 − 4ǫ − σ 0 −4 + 4ǫ + σ 0 FN

=



 σ0 −2 + σ 0 −4 + σ 0 1 2ǫ + σ 0 −2 + 2ǫ + σ 0 −4 + 2ǫ + σ 0  . h 4ǫ + σ 0 −2 + 4ǫ + σ 0 −4 + 4ǫ + σ 0

Estas matrices locales son independientes del intervalo In . Una vez que las matrices han sido calculadas, ellas son ensambladas dentro de la matriz global. El ensamble depende del orden de las incognitas αni . Si se supone que las incognitas est´ an ordenadas de la siguiente manera −1 N −1 N −1 α00 , α01 , α02 , α10 , α11 , α12 , α20 , α21 , α22 , . . . , αN , , α , α 0 1 2 entonces la matriz global viene dada por:  M0 D1  E1 M D2   ··· ··· ···   ··· ···   EN −2

 ··· M DN −1 EN −1 MN

donde

M = An + Bn + Cn+1 ,

M0 = A0 + F0 + C1 ,

   ,   

MN = AN −1 + FN + BN −1 .

Observaci´ on: Ya que el par´ ametro de penalizaci´ on es constante, las matrices locales son independientes de los subintervalos. Por tanto, ellas pueden ser definidas antes de ensamblar la matriz global.

2.4.4.

Construcci´ on del vector b

Recordemos que el vector b es un vector con componentes L(Φnj (x)), donde Z 1 σ0 n L(Φj (x)) = f (x)Φnj (x)dx − ǫ(Φnj )′ (x0 ) + Φnj (x0 ), h 0 para n = 0, . . . , N − 1 y j = 0, 1, 2. En virtud de la definici´ on de Φnj (x), se tiene Z 1 Z xn+1 n f (x)Φj (x)dx = f (x)φnj (x)dx. 0

xn

Despu´es de un cambio de variable (ver [18], Cap´ıtulo 4, Secci´ on 4.7, p´ ag. 224), se obtiene Z 1 Z h 1 h 1 f (x)Φnj (x)dx = f t+ n+ h tj dt. 2 2 2 0 −1 34

Secci´ on 2.5: M´etodo Galerkin Discontinuo Salvo en el caso en que la funci´ on f sea una funci´ on constante a trozos, la forma pr´ actica de calcular la integral anterior ser´ a mediante la integraci´ on num´erica, y, al igual que en el MEF, usaremos la f´ ormula de cuadratura Gaussiana de dos puntos (Ver [12], Cap´ıtulo 5, Secci´ on 5.3, p´ ag. 270). Por lo tanto Z

1 0

f (x)Φnj (x)dx

2 hX h 1 ≈ f sp + n + h (sp )j , 2 2 2 p=1

donde s1 = − √13 y s2 = √13 . Escribiremos las componentes del vector b en un orden consecuente con el orden de las incognitas αni −1 −1 −1 b00 , b01 , b02 , b10 , b11 , b12 , b20 , b21 , b22 , . . . , bN , bN , bN , 0 1 2 donde las tres primeras componentes vienen dadas por b00

=

2 hX h h σ0 sp + f + , 2 2 2 h p=1

b01

=

2 hX h h 2 σ0 f sp + sp − ǫ − , 2 2 2 h h p=1

b02 =

2 hX h h 4 σ0 f sp + (sp )2 + ǫ + , 2 p=1 2 2 h h

mientras que el resto de las 3(N − 1) componentes son bni

2 hX h 1 = f sp + n + h (sp )i . 2 2 2 p=1

A continuaci´ on se presenta un breve algoritmo de este m´etodo, el cual emplea las ideas b´ asicas anteriormente descritas. Algoritmo del M´ etodo Galerkin Discontinuo. Entrada: N´ umero de subintervalos N ≥ 2; par´ ametros ǫ y σ 0 . Salida: Vector de inc´ ognitas αni , i = 0, 1, 2; n = 0, . . . , N . 1. Tomar h = 1/N . 2. Construir las matrices locales An , Bn , Cn , Dn , En , F0 , FN . 3. Ensamblar las matrices locales dentro de la matriz global A. 4. Aplicar cuadratura Gaussiana de dos puntos para calcular y construir el vector b. 5. Tomar A−1 b. Para efecto de la implementaci´ on de los c´ odigos, se trabajar´ a con SIPG, esto es, cuando ǫ = −1 y σ 0 = 106 . A ra´ız de esto, la soluci´ on encontrada por este m´etodo aproxima la soluci´ on exacta con un 3 orden de aproximaci´ on O(h ) (ver [1], p´ ag. 13). 35

Cap´ıtulo 2: M´etodos Num´ericos para aproximar la soluci´ on de un PVF

2.5.

M´ etodo Bvp4c

Hasta los momentos se han planteado cuatro m´etodos para aproximar la soluci´ on de un Problema de Valor de Frontera de segundo orden: los m´etodos de Disparo y Diferencias Finitas son los m´etodos com´ unmente m´ as usados, con frecuencia aparecen en todas las referencias b´ asicas de An´ alisis Num´erico, mientras que los m´etodos de los Elementos Finitos y Galerkin Discontinuo, un poco m´ as complejos que los anteriores, parten de la teor´ıa de espacios normados y se fundamentan en resultados del An´ alisis Funcional. En esta secci´ on se plantear´ a el u ´ltimo de los m´etodos num´ericos que se emplear´ a para la resoluci´on del problema planteado: el M´etodo Bvp4c, presentado en [14], el cual se encuentra implementado en Matlab. El M´etodo Bvp4c pone en pr´ actica el principio de superposici´ on para aproximar la soluci´ on del Problema de Valor de Frontera de segundo orden:   y ′′ = f (x, y, y ′ ), y(a) = α, (2.32)  y(b) = β.

Para conseguir esto, se reescribe (2.32) como un sistema de ecuaciones diferenciales ordinarias de primer orden: (2.33) y′ = f(x, y), a ≤ x ≤ b, sujeta a condiciones generales de frontera no lineales:

g(y(a), y(b)) = 0. La soluci´ on aproximada del sistema (2.33), S(x), es una funci´ on continua. M´ as en general, la soluci´ on S(x) es un polinomio c´ ubico en cada subintervalo In = [xn , xn+1 ], de una malla de puntos a = x0 < x1 < · · · < xN = b. Esta soluci´ on aproximada satisface las condiciones de frontera que cumple y, es decir: g(S(a), S(b)) = 0. As´ı mismo, la soluci´ on S(x) satisface la ecuaci´ on diferencial en los extremos y punto medio de cada subintervalo In , esto es, S′ (xn ) = f(xn , S(xn )), ′

S

xn + xn+1 2

xn + xn+1 xn + xn+1 = f ,S , 2 2

S′ (xn+1 ) = f(xn+1 , S(xn+1 )). Estas ecuaciones conllevan a un sistema de ecuaciones algebraicas no lineales para los coeficientes de S(x). A diferencia del m´etodo del Disparo para el caso no lineal, la soluci´ on y(x) es aproximada sobre todo el intervalo [a, b] y las condiciones de frontera son tomadas en cuenta todo el tiempo. Las ecuaciones algebraicas no lineales son resueltas iterativamente por linealizaci´ on. La soluci´ on S(x) es una aproximaci´ on de cuarto orden para y(x), es decir, ||y(x) − S(x)|| ≤ Ch4 (ver [11]). Aqu´ı, h es el m´ aximo tama˜ no de paso hn = xn+1 − xn y C una constante.

Este m´etodo adapta la malla para obtener as´ı una soluci´ on num´erica acertada con un modesto n´ umero de puntos de dicha malla. Poder hacer una aproximaci´ on suficientemente buena es la parte m´ as dif´ıcil al momento de resolver un PVF. La continuidad de S(x) en [a, b] y la superposici´ on en los extremos y punto medio de cada intervalo implica que la derivada de S(x) tambi´en es continua en [a, b]. El m´etodo Bvp4c controla el error cometido al momento de aproximar y(x) mediante S(x). Para el control de dicho error, se introduce el residuo o resto de la ecuaci´ on diferencial. Este viene dado por r(x) = S′ (x) − f(x, S(x)). Como se puede apreciar, en los extremos y punto medio de cada intervalo In , el resto es cero. 36

Secci´ on 2.5: M´etodo Bvp4c Sintaxis. Como se dijo, el m´etodo Bvp4c es un m´etodo que usa la idea se superposici´ on y el mismo se encuentra implementado en Matlab, el cual permite resolver num´ericamente un PFV del tipo y ′′ (x) = f (x, y(x), y ′ (x)), x ∈ [a, b], sujeta a g(y(a), y(b)) = 0. La sintaxis b´ asica, empleando Matlab, es la siguiente: sol = bvp4c(@odef un, @bcf un, solinit), donde, 1. @odef un es el nombre del fichero de la ecuaci´ on diferencial ordinaria que se quiere resolver, expresado como un sistema de ecuaciones diferenciales ordinarias de primer orden, dv = h(x, v). Aqu´ı x es un escalar, y adem´ as, v

= [v1 v2 ]T = [y(x) y ′ (x)]T ,

dv = [v1′ v2′ ]T = [y ′ (x), f (x, y(x), y ′ (x))]T . 2. @bcf un es el nombre del fichero que contiene las condiciones de contorno. 3. solinit es el dato inicial de la ecuaci´ on diferencial en un conjunto de nodos ordenados a = x0 < x1 < · · · < xN = b. Para ello se usa la funci´ on bvpinit: solinit = bvpinit(xmesh, yinit), para un mallado inicial. La salida sol es una estructura con sol.x, que es el mallado seleccionado por Bvp4c, sol.y, aproximaci´ on a y(x) en los puntos del mallado sol.x, sol.yp, aproximaci´ on a y ′ (x) en los puntos del mallado sol.x.

37

Cap´ıtulo 2: Experimentaci´ on Num´erica

38

3 Experimentaci´ on Num´ erica

Como ya se ha dicho, los problemas de valor de frontera de segundo orden suelen ser comunes en todas las ramas de las Ciencias Experimentales e Ingenier´ıa. Por ejemplo, en F´ısica, las leyes de Newton y muchas otras se expresan como un problema de este tipo, en Biolog´ıa aparecen en el modelado de din´ amica de poblaciones, en la Qu´ımica surgen en la evaluaci´ on de las concentraciones de diversos reactivos durante una reacci´ on, etc. Sin embargo, la dificultad radica en que las ecuaciones diferenciales que surgen de estos modelos no pueden resolverse anal´ıticamente excepto en algunos casos especiales. A ra´ız de esto, resulta necesario usar m´etodos num´ericos para dar aproximaciones (suficientemente buenas) de la soluci´ on de dichos modelos. En ese sentido, se han introducido, en este trabajo, cinco m´etodos num´ericos con la finalidad de estimar esas aproximaciones, los cuales son: el M´etodo del Disparo y Diferencias Finitas (ver [9, 18]), el M´etodo de los Elementos Finitos (ver [3, 17]), el M´etodo Galerkin Discontinuo (ver [1]) y el M´etodo Bvp4c, un m´etodo que usa la idea se superposici´ on y el mismo se encuentra implementado en Matlab (ver[14]). En muchas ocasiones estamos interesados en tener conocimento sobre costo de error y costo computacional (tiempo de CPU usado) que ocurre al momento de usar uno de los m´etodos anteriormente mencionados, es decir, nos interesa saber qu´e tan grande es la magnitud del error que se comete cuando se aproxima la soluci´ on y el tiempo de m´ aquina que se requiere para realizar tal aproximaci´ on. Es por esto que en este cap´ıtulo se har´ a una breve discusi´ on, m´ as que comparaci´ on, sobre la eficiencia y propiedades de cada uno de estos m´etodos, esto en vista de que no se presentan las condiciones adecuadas para realizar una comparaci´ on detallada y ´ optima entre ellos. Por ejemplo, el m´etodo Bvp4c realiza proceso adaptativo mientras que los otros m´etodos no realizan. Otro ejemplo es que el m´etodo de los Elementos Finitos y Galerkin Discontinuo s´ olo son aplicables al caso en que la ecuaci´ on diferencial es del tipo lineal, mientras que los m´etodos del Disparo, Diferencias Finitas y Bvp4c se pueden emplear en un PVF ya sea lineal como no lineal. Esta discusi´ on se har´ a usando los c´ odigos de estos m´etodos sobre un conjunto de ejemplos. Las implementaciones num´ericas de estos m´etodos es llevada a cabo en Matlab. No obstante, por considerarlo irrelevante, los c´ odigos no son dados expl´ıcitamente en este trabajo. Igualmente, en este cap´ıtulo, se estimar´ a el error cometido al momento de aproximar la soluci´ on con los m´etodos que se han introducido. Para ello se emplear´ a la norma del m´ aximo o norma infinita, esto es, si x es un vector de n componentes, entonces ||x||∞ = m´ ax{|x1 |, |x2 |, . . . , |xn |}. De la misma maner´ a, se usar´ a la norma L2 , dada por Z b 1/2 2 ||eh (x)||L2 = (eh (x)) dx , a

con eh (x) = u(x) − u e(x), donde u(x) es la soluci´ on exacta del PVF planteado y u e(x) su aproximaci´ on. En el caso en que se emplee el m´etodo Galerkin Discontinuo, u e(x) es un polinomio de grado 2 en cada 39

Cap´ıtulo 3: Experimentaci´ on Num´erica subintervalo de la partici´ on, mientras que si se emplea alguno de los otros m´etodos, entonces u e(x) es un polinomio lineal en cada subintervalo de dicha partici´ on.

En lo que sigue, se ponen a prueba los c´ odigos de los m´etodos num´ericos sobre un conjunto de seis ejemplos. Los tres primeros son del tipo lineal, lo que permite emplear los cinco m´etodos, mientras que los ejemplos 4, 5 y 6 son del tipo no lineal, lo que conlleva a usar s´ olo el m´etodo del Disparo, Diferencias Finitas y Bvp4c. En vista de que el m´etodo del Disparo y Diferencias Finitas, para el caso no lineal, usan un m´etodo iterativo para encontrar el cero de una funci´ on no lineal y resolver un sistema de ecuaciones no lineales, respectivamente, estos se ejecutar´ an usando una tolerancia de 10−3 . De la misma manera, como el m´etodo Bvp4c realiza proceso adaptativo, este se ejecutar´ a usando una tolerancia tambi´en de 10−3 . Para el m´etodo Galerkin Discontinuo, se trabajar´ a en el caso particular en que ε = −1 y σ 0 = 106 , es decir, se trabajar´a con el m´etodo SIPG. Se estimar´ a el error mediante la norma del m´ aximo y la norma L2 , los cuales se denotar´ an mediante ||ε||∞ y ||ε||L2 , respectivamente, as´ı como tambi´en el tiempo de m´ aquina, en segundos, requerido para realizar las distintas operaciones, y estos datos se representar´ an en cuadros. Para finalizar, se presentan un conjunto de gr´ aficas que revelan el comportamiento de los errores producidos para una partici´ on particular del dominio donde est´ a definida la soluci´ on de la ecuaci´ on diferencial planteada.

Ejemplo 1 Consideremos el siguiente problema de valor de frontera lineal: y ′′ = −27x +

y(0) = 1, y(1) = 0,

28 , 3

el cual tiene por soluci´ on exacta 14 7 9 y(x) = − x3 + x2 − x + 1. 2 3 6 La gr´ afica de la soluci´ on exacta y(x) se refleja en la Figura 3.1. Si se particiona el intervalo [0, 1] en 10, 50 y 100 subintervalos de igual longitud, entonces, en virtud del Cuadro 3.1, se logra apreciar que los errores, empleando la norma del m´ aximo, son relativamente buenos. Notemos que el m´etodo del Disparo logra obtener una mejor aproximaci´ on que los restantes m´etodos, a pesar que emplea casi el doble del tiempo que el m´etodo de Diferencias Finitas, que es m´ as r´ apido. Esto es previsible debido a que Disparo emplea el m´etodo de Runge - Kutta de cuarto orden, el cual ofrece una exactitud O(h4 ) a las soluciones de los problemas de valor inicial. Para un tama˜ no de paso de h = 0.02, el comportamiento del error cometido por cada uno de los m´etodos se ilustra en la Figura 3.1. En las gr´ aficas se reflejan como, para Diferencias Finitas, Elementos Finitos y Galerkin Discontinuo, al aumentar x, el error va aumentando, hasta llegar a un cierto punto a partir del cual el error empieza a disminuir. Para Disparo y Bvp4c, el error empieza a aumentar y disminuir en peque˜ na magnitud, pero a partir de un cierto valor el error presenta oscilaciones considerables. Para el error empleando la norma L2 , el m´etodo Galerkin Discontinuo logra obtener una mejor aproximaci´ on sobre los otros m´etodos, en vista de que este emplea polinomios de grado 2, a diferencias de los otros m´etodos que emplean polinomios de grado 1. Es importante resaltar que el m´etodo Bvp4c no realiza proceso de adaptatividad en las tres mallas iniciales dadas, ya que en cada subintervalo de cada una de las tres particiones, la norma del residuo, que resulta de aproximar la soluci´ on mediante polinomios c´ ubicos, est´ a por debajo de la tolerancia dada. Si se toma como referencia la malla de 10 subintervalos de igual longitud, entonces, en vista del Cuadro 3.1, 40

Ejemplo 1

8

x 10

−16

7 6 5 y 4 3 2 1 0 0

0.1

(a) Gr´ afica de la soluci´ on exacta 8

x 10

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0.9

1

(b) Error absoluto producido con Disparo

−15

1.2

x 10

−14

7

1 6

0.8

5

y

y 4

0.6

3

0.4 2

0.2

1 0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

0 0

1

(c) Error absoluto producido con Diferencias. Finitas 8

x 10

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(d) Error absoluto producido con Elementos Finitos

−9

1.8

7

1.6

6

1.4

x 10

−15

1.2

5 y

y

4

1 0.8

3

0.6 2

0.4 1 0 0

0.2 0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

(e) Error absoluto producido con Galerkin Discontinuo

1

0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(f) Error absoluto producido con Bvp4c

Figura 3.1: Ejemplo 1. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 41

Cap´ıtulo 3: Experimentaci´ on Num´erica

8

x 10

−15

7 6 5 y 4 3 2 1 0 0

0.2

0.4

0.6

0.8

1

x

(a) Gr´ afica de la soluci´ on exacta 1

x 10

(b) Error relativo producido con Disparo

−14

1.4

x 10

−14

1.2 0.8 1 0.6

0.8

y

y 0.6

0.4

0.4 0.2 0.2 0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x

(c) Error relativo producido con Diferencias. Finitas 8

x 10

0.6

0.8

1

x

(d) Error relativo producido con Elementos Finitos

−9

6

x 10

−15

7 5 6 4 5 y

y

3

4 3

2 2 1 1 0 0

0.2

0.4

0.6

0.8

x

(e) Error relativo producido con Galerkin Discontinuo

1

0 0

0.2

0.4

0.6

0.8

1

x

(f) Error relativo producido con Bvp4c

Figura 3.2: Ejemplo 1. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 42

Ejemplo 2 N = 10

M´etodo Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

Tiempo 0.043 0.021 0.062 0.115 0.979

N = 50

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.046 0.024 0.080 0.127 1.028

N = 100

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.049 0.026 0.104 0.148 1.064

||ε||∞ 4.441 · 10−16 5.551 · 10−16 8.882 · 10−16 1.788 · 10−9 7.772 · 10−16

||ε||L2 7.340 · 10−3 7.340 · 10−3 7.340 · 10−3 2.165 · 10−4 7.340 · 10−3

6.661 · 10−16 1.088 · 10−14 1.410 · 10−14 2.195 · 10−8 2.442 · 10−15

7.365 · 10−5 7.365 · 10−5 7.365 · 10−5 2.170 · 10−7 7.365 · 10−5

7.772 · 10−16 7.994 · 10−15 1.177 · 10−14 7.703 · 10−9 1.665 · 10−15

2.946 · 10−4 2.946 · 10−4 2.946 · 10−4 1.732 · 10−6 2.946 · 10−4

Cuadro 3.1: Ejemplo 1. Tiempo empleado y errores obtenidos para un mallado de 10, 50 y 100 subintervalos. Bvp4c emplea 0.979 segundos para obtener la aproximaci´ on con un error m´ aximo de 7.772 · 10−16 . Para que los restantes cuatros m´etodos requieran un tiempo relativamente igual, necesitan emplear el n´ umero de intervalos que vienen dados en el Cuadro 3.2. M´etodo

Intervalos

Disparo

3610

Dif. Finitas

5635

Elementos Finitos G. Discontinuo

1100 524

||ε||∞

6.151 · 10−14 5.358 · 10−12 7.480 · 10−13 1.570 · 10−7

Cuadro 3.2: Ejemplo 1. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 0.979 segundos. Como se puede apreciar en el Cuadro 3.2, el m´etodo de Elementos Finitos logra obtener una aproximaci´ on casi tan buena como la aportada por Disparo, en una cantidad de intervalos mucho menor, mientras que, si se hace la comparaci´ on con Diferencias Finitas, la aproximaci´ on es mejor y emplea s´ olo casi una quinta parte de los intervalos que amerita este u ´ltimo m´etodo. El m´etodo Galerkin Discontinuo requiere menos intervalos pero el error obtenido aumenta significativamente en relaci´ on cuando se usan menos intervalos.

Ejemplo 2 Trabajemos ahora con el siguiente problema de valor de frontera: 2

y ′′ = −(4x3 − 4x2 − 6x + 2)e−x , y(0) = 1, y(1) = 0, 43

Cap´ıtulo 3: Experimentaci´ on Num´erica N = 10

M´etodo Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

Tiempo 0.041 0.023 0.061 0.115 0.987

N = 50

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.045 0.024 0.084 0.141 1.018

N = 100

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.049 0.027 0.103 0.151 1.075

||ε||∞ 1.915 · 10−6 1.203 · 10−3 1.277 · 10−6 1.277 · 10−6 1.915 · 10−6

||ε||L2 1.115 · 10−3 8.653 · 10−4 1.114 · 10−3 3.796 · 10−5 1.115 · 10−3

1.902 · 10−10 1.211 · 10−5 1.268 · 10−10 9.361 · 10−9 1.902 · 10−10

1.118 · 10−5 8.637 · 10−6 1.118 · 10−5 3.846 · 10−8 1.118 · 10−5

3.042 · 10−9 4.845 · 10−5 2.028 · 10−9 1.755 · 10−9 3.042 · 10−9

4.472 · 10−5 3.455 · 10−5 4.472 · 10−5 3.041 · 10−7 4.472 · 10−5

Cuadro 3.3: Ejemplo 2. Tiempo empleado y errores obtenidos para un mallado de 10, 50 y 100 subintervalos. con soluci´ on

2

y(x) = (1 − x)e−x . La gr´ afica de la soluci´ on exacta de este problema lineal se refleja en la Figura 3.3. Particionando uniformente el intervalo [0, 1] en 10, 50 y 100 subintervalos, se aprecia, en virtud del Cuadro 3.3, que las aproximaciones, empleando la norma del m´ aximo, son ´ optimas. Para el caso del m´etodo de Diferencias Finitas, los resultados son menos exactos que los restantes m´etodos, esto se debe a que presenta un error del orden O(h2 ), cosa que no sucede, por citar, con el m´etodo del Disparo o Bvp4c, que presentan un error del orden O(h4 ).

El comportamiento del error al dividir el intervalo [0, 1] en 50 subintervalos de igual longitud se ilustra en la Figura 3.3. Se puede apreciar en las gr´ aficas como, para Disparo, Diferencias finitas, Elementos Finitos y Bvp4c, el error empieza a aumentar hasta alrededor de x = 0.5, punto en el cual empieza a disminuir en forma considerable. Para el m´etodo Galerkin Discontinuo, el error se comporta de manera oscilatoria hasta un entorno de x = 0.5, para luego empezar a disminuir de manera estable. Al igual que en el ejemplo anterior, el m´etodo Galerkin Discontinuo logra dar una mejor aproximaci´ on a la soluci´ on del PVF empleando la norma L2 , por encima de los otros cuatros m´etodos.

Para cada una de las tres mallas que se usan en el Cuadro 3.3, el m´etodo Bvp4c no realiza proceso de adaptatividad, en vista de que en cada subintervalo de esas tres mallas, la norma del residuo, que resulta de aproximar la soluci´ on mediante polinomios c´ ubicos, est´ a por debajo de la tolerancia dada. Tomando como referencia la malla de 10 subintervalos de igual longitud, entonces, en vista del Cuadro 3.3, el m´etodo Bvp4c emplea 0.987 segundos para obtener la aproximaci´ on con un error m´ aximo de 1.915 · 10−6 . Para que los restantes cuatros m´etodos requieran un tiempo relativamente igual, necesitan emplear el n´ umero de intervalos que se expresan en el Cuadro 3.4. Como resultado se tiene, empleando el Cuadro 3.3 y el Cuadro 3.4, que tanto el m´etodo del Disparo como Diferencias Finitas requieren una cantidad considerable de intervalos para usar el tiempo empleado por Bvp4c para 10 subintervalos. Sin embargo, esto implica que los errores, empleando la norma del 44

Ejemplo 2

3.5

x 10

−9

3 2.5 2 y 1.5 1 0.5 0 0

0.1

(a) Gr´ afica de la soluci´ on exacta 5

x 10

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0.9

1

(b) Error absoluto producido con Disparo

−5

2.5

4

2

3

1.5

y

x 10

−9

y 2

1

1

0.5

0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

0 0

1

(c) Error absoluto producido con Diferencias. Finitas 1.8

x 10

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(d) Error absoluto producido con Elementos Finitos

−9

3.5

x 10

−9

1.6 3

1.4 2.5

1.2 y

1

y

0.8

2 1.5

0.6 1

0.4 0.5

0.2 0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

(e) Error absoluto producido con Galerkin Discontinuo

1

0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(f) Error absoluto producido con Bvp4c

Figura 3.3: Ejemplo 2. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 45

Cap´ıtulo 3: Experimentaci´ on Num´erica

2.5

x 10

−8

2

1.5 y 1

0.5

0 0

0.1

(a) Gr´ afica de la soluci´ on exacta 6

x 10

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(b) Error relativo producido con Disparo

−4

x 10

1.6

−8

1.4 5 1.2 4 1 y

y 0.8

3

0.6 2 0.4 1 0.2 0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x

(c) Error relativo producido con Diferencias. Finitas 2

x 10

0.6

0.8

1

x

(d) Error relativo producido con Elementos Finitos

−9

2.5

x 10

−8

2 1.5 1.5 y

y

1

1 0.5 0.5

0 0

0.2

0.4

0.6

0.8

x

(e) Error relativo producido con Galerkin Discontinuo

1

0 0

0.2

0.4

0.6

0.8

1

x

(f) Error relativo producido con Bvp4c

Figura 3.4: Ejemplo 2. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 46

Ejemplo 3 M´etodo

Intervalos

Disparo Dif. Finitas

3554 5559

Elementos Finitos

1108

G. Discontinuo

514

||ε||∞

4.796 · 10−14 3.918 · 10−9 3.221 · 10−12 3.160 · 10−7

Cuadro 3.4: Ejemplo 2. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 0.987 segundos. m´ aximo, aunmenten considerablemente, en relaci´ on a cuando se usan menos intervalos. Por otro lado, y a diferencia del m´etodo del Disparo y Diferencias Finitas, tanto el m´etodo de Elementos Finitos como el m´etodo Galerkin Discontinuo requieren menos intervalos para obtener el tiempo de prueba, aunque el error aumenta levemente. Como es de apreciar, el m´etodo de Elementos Finitos aproxima mejor la soluci´ on en menos intervalos en comparaci´ on con Diferencias Finitas. Para finalizar, Galerkin Discontinuo emplea menos intervalos, m´ as el error obtenido no disminuye, por el contrario, aumenta debido al gran volumen de operaciones que se ven involucrados en la implementaci´ on del c´ odigo.

Ejemplo 3 Consideremos el problema de valor de frontera lineal: y ′′ = −4

sin(2x) + 4 cos2 (x) − 4x cos2 (x) − 5 + 5x , cos6 (x)

y(0) = 1, y(1) = 0, con soluci´ on exacta y(x) =

1−x . cos4 (x)

La gr´ afica de la soluci´ on exacta y(x) viene dada en la Figura 3.5. En el Cuadro 3.5, se reflejan los datos obtenidos para una partici´ on uniforme de [0, 1] en 10, 50 y 100 subintervalos. Para el mallado inicial de 10 subintervalos, el m´etodo Bvp4c adapta la malla con la finalidad de que el residuo satisfaga la tolerancia dada, obteniendo as´ı un nuevo mallado de 15 subintervalos (no todos de la misma longitud). Al igual que los dos ejemplos anteriores, las aproximaciones arrojadas por cada m´etodo son relativamente optimas, haciendo ´enfasis, nuevamente, en que el m´etodo de Diferencias Finitas proporciona errores un ´ poco mayor que los dados por los otros m´etodos, esto, en parte, al bajo orden de aproximaci´ on (O(h2 )) que presenta. Es importante resaltar que el proceso de adaptatividad por parte de Bvp4c conlleva a que se realicen realizen m´ as operaciones con el fin de mejorar los resultados. Es por eso, por ejemplo, que el tiempo empleado al usar la malla inicial de 10 subintervalos es mucho mayor que el tiempo usado en la malla de 50 o 100 subintervalos. Para el caso particular en que el mallado consta de 51 puntos distribuidos uniformemente, el comportamiento del error cometido se refleja en la Figura 3.5. Aqu´ı se puede notar como los errores, para cada uno de los cinco m´etodos, se comportan de manera similar. Los errores empiezan a aumentar de forma lineal hasta un entorno de x = 0.8, punto donde aproximadamente alcanzan el m´ aximo, para luego disminuir de manera r´ apida. En esta malla, el m´etodo Bvp4c usa un tiempo de 1.045 segundos para realizar las operaciones involucradas en la aproximaci´ on, cometiendo un error m´ aximo de 3.073 · 10−6 . Para tener un tiempo aproximadamente igual, los cuatros m´etodos necesitar´ an aproximar la soluci´ on sobre el n´ umero de intervalos que se reflejan en el Cuadro 3.6. 47

Cap´ıtulo 3: Experimentaci´ on Num´erica

3.5

x 10

−6

3 2.5

y

2 1.5 1 0.5 0 0

0.1

(a) Gr´ afica de la soluci´ on exacta 3.5

x 10

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0.9

1

(b) Error absoluto producido con Disparo

−3

2.5

x 10

−6

3

2 2.5

y

1.5

2

y 1.5

1

1

0.5 0.5 0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

0 0

1

(c) Error absoluto producido con Diferencias. Finitas 2.5

x 10

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(d) Error absoluto producido con Elementos Finitos

−6

3.5

x 10

−6

3

2 2.5

1.5 y

y

2 1.5

1

1

0.5 0.5

0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

(e) Error absoluto producido con Galerkin Discontinuo

1

0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(f) Error absoluto producido con Bvp4c

Figura 3.5: Ejemplo 3. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 48

Ejemplo 3

5

x 10

−6

4

3 y 2

1

0 0

0.2

0.4

0.6

0.8

1

x

(a) Gr´ afica de la soluci´ on exacta 5

x 10

(b) Error relativo producido con Disparo

−3

3.5

x 10

−6

3 4 2.5 3

2 y

y

1.5

2

1 1 0.5 0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x

(c) Error relativo producido con Diferencias. Finitas 3.5

x 10

0.6

0.8

1

x

(d) Error relativo producido con Elementos Finitos

−6

5

x 10

−6

3 4 2.5 3

2 y

y 1.5

2

1 1 0.5 0 0

0.2

0.4

0.6

0.8

x

(e) Error relativo producido con Galerkin Discontinuo

1

0 0

0.2

0.4

0.6

0.8

1

x

(f) Error relativo producido con Bvp4c

Figura 3.6: Ejemplo 3. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 49

Cap´ıtulo 3: Experimentaci´ on Num´erica N = 10

M´etodo Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

Tiempo 0.043 0.021 0.059 0.116 1.118

N = 50

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.045 0.024 0.077 0.129 1.045

N = 100

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.049 0.027 0.102 0.148 1.070

N = 15

||ε||∞ 1.706 · 10−3 7.850 · 10−2 1.128 · 10−3 1.128 · 10−3 1.120 · 10−4

||ε||L2 2.340 · 10−2 6.295 · 10−2 2.274 · 10−2 2.877 · 10−3 6.232 · 10−3

1.928 · 10−7 8.543 · 10−4 1.285 · 10−7 1.159 · 10−7 1.928 · 10−7

2.466 · 10−4 6.933 · 10−5 2.465 · 10−4 2.901 · 10−6 2.466 · 10−4

3.073 · 10−6 3.407 · 10−3 2.048 · 10−6 2.052 · 10−6 3.073 · 10−6

9.846 · 10−4 2.765 · 10−3 9.833 · 10−4 2.321 · 10−5 9.846 · 10−4

Cuadro 3.5: Ejemplo 3. Tiempo empleado y errores obtenidos para un mallado de 10, 50 y 100 subintervalos. Para la malla inicial de 10 subintervalos, el m´etodo Bvp4c adapta la malla con la finalidad de mejorar la estimaci´ on de la soluci´ on, obteniendo as´ı una malla final de 15 subintervalos. M´etodo

Intervalos

Disparo

3654

Dif. Finitas

5801

Elementos Finitos

1109

G. Discontinuo

533

||ε||∞

1.020 · 10−13 2.542 · 10−7

8.607 · 10−12 2.698 · 10−7

Cuadro 3.6: Ejemplo 3. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 1.045 segundos. Nuevamente, al igual que los Ejemplos 1 y 2, se obtiene, en virtud del Cuadro 3.5 y del Cuadro 3.6, que los errores se incrementan a medida que aumentan los intervalos de la partici´ on. Tanto el m´etodo del Disparo como Elementos Finitos aportan errores relativamente iguales, con la diferencia de que Elementos Finitos requiere menos intervalos para aproximar la soluci´ on. De la misma manera, tanto Diferencias Finitas como Galerkin Discontinuo obtienen errores similares, pero se diferencian en el hecho de que Galerkin Discontinuo necesita menos de la d´ecima parte del total de intervalos que necesita Diferencias Finitas.

Ejemplo 4 Trabajemos ahora con la siguiente ecuaci´ on diferencial: −12(2x − 1)p , (p + (2x − 1)2 )5/2 y(0) = 1, y(1) = 0, y ′′ =

50

Ejemplo 4 el cual tiene por soluci´ on exacta √ 1 (2 p + 1 + p + 1)(2x − 1) 1 y(x) = p − + , p+1 2 p + (2x − 1)2 2 2x − 1

donde p = 10−3 . En la Figura 3.7, se ilustra la gr´ afica de la soluci´ on y(x). Para una partici´ on uniforme del intervalo [0, 1] en 100, 200 y 500 subintervalos, los errores y tiempo requerido por cada m´etodo se expresan en el Cuadro 3.7. Es importante mencionar que para el mallado inicial que consta de 101 puntos distribuidos uniformemente, el m´etodo Bvp4c realiza proceso adaptativo, para obtener as´ı una malla final que consta de 65 puntos, mientras que para el mallado de 201 puntos, este m´etodo realiza nuevamente proceso adaptativo para obtener una mallado constitu´ıdo por 139 puntos, no necesariamente distribu´ıdos uniformente.

N = 100

N = 64 N = 200

N = 138 N = 500

M´etodo Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

Tiempo 0.066 0.032 0.120 0.199 1.582

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.091 0.047 0.164 0.376 1.669

Disparo Dif. Finitas Elementos Finitos G. Discontinuo Bvp4c

0.171 0.064 0.417 1.231 1.744

||ε||∞ 1.230 · 10−3 2.522 · 10−2 8.797 · 10−4 8.796 · 10−4 1.115 · 10−3

||ε||L2 4.154 · 10−3 2.438 · 10−3 4.479 · 10−3 4.479 · 10−3 1.408 · 10−3

2.144 · 10−6 1.150 · 10−3 1.431 · 10−6 1.435 · 10−6 2.144 · 10−6

1.761 · 10−4 9.045 · 10−6 1.758 · 10−4 1.758 · 10−4 1.761 · 10−4

9.264 · 10−5 6.931 · 10−3 6.216 · 10−5 6.223 · 10−5 9.849 · 10−4

1.096 · 10−3 1.479 · 10−4 1.087 · 10−3 1.087 · 10−3 1.404 · 10−3

Cuadro 3.7: Ejemplo 4. Tiempo empleado y errores obtenidos para un mallado de 100, 200 y 500 subintervalos. Para las mallas iniciales de 100 y 200 subintervalos, el m´etodo Bvp4c adapta dichas mallas con la finalidad de mejorar la estimaci´ on de la soluci´ on, obteniendo as´ı dos mallas finales de 64 y 138 subintervalos, respectivamente

En la Figura 3.7 se muestran las gr´ aficas de los distintos errores absolutos al emplear cada uno de los cinco m´etodos num´ericos, trabajando con un mallado de 200 subintervalos. Como se dijo anteriormente, para esta malla en m´etodo Bvp4c realiza proceso adaptativo con la finalidad de mejorar la estimaci´ on, para ello crea una malla final que va a constar de 139 nodos. Para realizar dicha adaptatividad, este m´etodo emplea 1.669 segundos, cometiendo un error m´ aximo de 9.849 · 10−4 . Para que los restantes cuatros m´etodo usen un tiempo similar, necesitar´ an usar aproximadamente el n´ umero de subintervalos que se ilustran en el Cuadro 3.8 51

Cap´ıtulo 3: Experimentaci´ on Num´erica

1

x 10

−4

0.8

0.6 y 0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

x

(a) Gr´ afica de la soluci´ on exacta 7

x 10

(b) Error absoluto producido con Disparo

−3

7

6

6

5

5

4

x 10

−5

4

y

y 3

3

2

2

1

1

0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x

(c) Error absoluto producido con Diferencias. Finitas 7

x 10

0.6

0.8

1

x

(d) Error absoluto producido con Elementos Finitos

−5

1

x 10

−3

6 0.8 5 0.6

4 y

y 3

0.4

2 0.2 1 0 0

0.2

0.4

0.6

0.8

x

(e) Error absoluto producido con Galerkin Discontinuo

1

0 0

0.2

0.4

0.6

0.8

1

x

(f) Error absoluto producido con Bvp4c

Figura 3.7: Ejemplo 4. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 200 subintervalos, empleando la norma del m´ aximo. 52

Ejemplo 4

7

x 10

−3

6 5 4 y 3 2 1 0 0

0.1

(a) Gr´ afica de la soluci´ on exacta

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0.9

1

(b) Error relativo producido con Disparo

1.6 4.5

1.4

y

0.2

x 10

−3

4

1.2

3.5

1

3 2.5

0.8

y 2

0.6 1.5

0.4 1

0.2

0.5

0 0

0.5 0.6 0.7 0.8 0.9 x (c) Error relativo producido con Diferencias. Finitas

4.5

x 10

0.1

0.2

0.3

0.4

0 0

1

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(d) Error relativo producido con Elementos Finitos

−3

0.25

4

0.2 3.5 3

0.15

2.5

y

y 2

0.1

1.5 1

0.05

0.5 0 0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

(e) Error relativo producido con Galerkin Discontinuo

1

0 0

0.1

0.2

0.3

0.4

0.5 0.6 0.7 0.8 x (f) Error relativo producido con Bvp4c

0.9

1

Figura 3.8: Ejemplo 4. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas, Elementos Finitos, Galerkin Discontinuo y Bvp4c, para un mallado de 50 subintervalos, empleando la norma del m´ aximo. 53

Cap´ıtulo 3: Experimentaci´ on Num´erica M´etodo

Intervalos

Disparo Dif. Finitas

3841 6930

Elementos Finitos

1301

G. Discontinuo

631

||ε||∞

6.173 · 10−10 5.960 · 10−6 3.123 · 10−8 9.876 · 10−7

Cuadro 3.8: Ejemplo 4. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 1.669 segundos. N = 20

M´etodo Disparo Dif. Finitas Bvp4c

Tiempo 0.045 0.042 0.886

N = 40

Disparo Dif. Finitas Bvp4c

0.052 0.047 1.046

N = 80

Disparo Dif. Finitas Bvp4c

0.059 0.052 1.087

||ε||∞ 1.672 · 10−6 3.643 · 10−5 3.087 · 10−7

||ε||L2 7.431 · 10−5 5.530 · 10−5 7.384 · 10−5

4.807 · 10−7 2.276 · 10−6 3.079 · 10−7

4.785 · 10−6 3.454 · 10−6 4.759 · 10−6

8.711 · 10−7 9.106 · 10−6 3.077 · 10−7

1.875 · 10−5 1.382 · 10−5 1.857 · 10−5

Cuadro 3.9: Ejemplo 5. Tiempo empleado y errores obtenidos para un mallado de 20, 40 y 80 subintervalos.

Ejemplo 5 Consideremos la siguiente ecuaci´ on diferencial: y ′′ = −y 3 + y(0.2) =

0.2 , 1.04

y(0.8) =

0.8 , 1.64

el cual tiene por soluci´ on exacta y(x) =

3x3 − 6x , (1 + x2 )3

x . 1 + x2

Como se puede notar, la ecuaci´ on planteada es un problema de valor de frontera del tipo no lineal. En la Figura 3.9, se ilustra la gr´ afica de la soluci´ on y(x). Particionando uniformemente el intervalo [0.2, 0.8] en 20, 40 y 80 subintervalos, entonces, en virtud del Cuadro 3.9, se aprecia que los errores, empleando las dos normas, son relativamente buenos. Se puede notar que las estimaciones aportadas tanto por Disparo como por Bvp4c son mejores que las aportadas por Diferencias Finitas, en vista de que el orden de aproximaci´ on 4 2 de los dos primeros es el doble que el orden de aproximaci´ on de Diferencias Finitas (O(h ) contra O(h )).

Para un tama˜ no de paso h = 0.0075, el comportamiento del error que se comete con los tres m´etodos se refleja en la Figura 3.9. Como es de apreciar, los errores, tanto para Diferencias Finitas (que amerita dos iteraciones por parte del m´etodo de Newton para sistemas de ecuaciones no lineales) como Bvp4c, empiezan a aumentar considerablemente, alcanzando el m´ aximo alrededor de x = 0.4 y x = 0.5, respectivamente, 54

Ejemplo 5

5

x 10

−7

4

3 y 2

1

0 0.2

0.3

0.5 0.6 0.7 x (b) Error absoluto producido con Disparo

(a) Gr´ afica de la soluci´ on exacta 2.5

x 10

0.4

0.8

−7

−6

3.5

x 10

3

2 2.5

1.5 y

y

2 1.5

1 1

0.5 0.5

0 0.2

0.3

0.4

0.5 0.6 0.7 0.8 x (c) Error absoluto producido con Diferencias. Finitas

0 0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(d) Error absoluto producido con Bvp4c

Figura 3.9: Ejemplo 5. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 80 subintervalos, usando la norma del m´ aximo.

55

Cap´ıtulo 3: Experimentaci´ on Num´erica

1

x 10

−6

0.8

0.6 y 0.4

0.2

0 0.2

0.4

0.6

0.8

x

(a) Gr´ afica de la soluci´ on exacta 7

x 10

(b) Error relativo producido con Disparo

−6

8

x 10

−7

7

6

6

5 5

4

y

y

4

3 3

2

2

1 0 0.2

1

0.4

0.6

0.8

x

(c) Error relativo producido con Diferencias. Finitas

0 0.2

0.4

0.6

0.8

x

(d) Error relativo producido con Bvp4c

Figura 3.10: Ejemplo 5. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenido con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 80 subintervalos, usando la norma del m´ aximo.

56

Ejemplo 6 para luego empezar a disminuir. En el caso del m´etodo del Disparo (que requiere tres iteraciones del m´etodo Newton Raphson) el error crece de forma lineal. Es importante resaltar que tanto en el m´etodo de Diferencias Finitas como en Bvp4c, la soluci´ on es aproximada sobre todo el intervalo [0.2, 0.8] y las condiciones de frontera son tomadas todo el tiempo, cosa que no ocurre con Disparo. Es por esta raz´ on, para este u ´ltimo m´etodo, que el error producido para x = 0.8 no necesariamente es cero, como se podr´ıa pensar. Para las tres mallas iniciales dadas, el m´etodo Bvp4c no realiza proceso de adaptatividad, lo que hace pensar que las aproximaciones usando los polinomios c´ ubicos son buenas. Particularmente, cuando se trabaja en una malla uniforme que consta de 21 nodos, el m´etodo Bvp4c usa un tiempo de 0.886 segundos para realizar las operaciones involucradas en la aproximaci´ on, cometiendo un error m´ aximo de 3.087 · 10−7 . Para que Disparo y Diferencias Finitas requieran un tiempo similar, necesitan emplear el n´ umero de intervalos que vienen dados en el Cuadro 3.10. M´etodo

Intervalos

Disparo

3601

Dif. Finitas

3391

||ε||∞

1.050 · 10−7 1.267 · 10−9

Cuadro 3.10: Ejemplo 5. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 0.886 segundos. Empleando el Cuadro 3.9 y el Cuadro 3.10, se logra apreciar que el error obtenido por el m´etodo del Disparo no mejora, sino que se mantiene aproximadamente igual pese a aumentar la partici´ on del dominio de la ecuaci´ on diferencial. Por el contrario, las estimaciones aportadas por Diferencias Finitas son mejores a medida que se incrementa la partici´ on del dominio. El error que comete Diferencias Finitas es mucho menor en comparaci´ on con el m´etodo del Disparo, con la ventaja de que el mallado requiere menos intervalos.

Ejemplo 6 Consideremos el siguiente PVF no lineal: y ′′ = −yy ′ − 2 sin(x) − x cos(x) + x cos2 (x) −

y(−1) = − cos(−1), y(2) = 2 cos(2),

x2 sin(2x) , 2

con soluci´ on exacta y(x) = x cos(x). La gr´ afica de la soluci´ on del PVF planteado se muestra en la Figura 3.11. Nuevamente, al igual que el ejemplo anterior, se realiza una partici´ on uniforme del intervalo [−1, 2] en 20, 40 y 80 subintervalos. En virtud del Cuadro 3.11, se logra notar que los errores usando cada una de las dos normas introducidas son aceptables. Las estimaciones aportadas por Disparo y Bvp4c siguen siendo mejores que las aportadas por Diferencias Finitas, debido a que este u ´ltimo m´etodo posee un orden de aproximaci´ on bajo, en comparaci´ on con los otros dos. Para un tama˜ no de paso h = 0.0375, es decir, dividiendo uniformemente el intervalo [−1, 2] en 80 subintervalos, el comportamiento del error cometido viene dada en la Figura 3.11. Para el m´etodo del Disparo, que amerita cuatro iteraciones por parte del m´etodo de Newton Raphson, el error aumenta considerablemente a lo largo de todo el intervalo [−1, 2], sucediendo lo mismo que en el ejemplo anterior, en el cual el error en la segunda condici´ on de frontera no necesariamente es cero. Para Diferencias Finitas, que requiere cuatro iteraciones por parte del m´etodo de Newton para sistemas de ecuaciones no lineales, 57

Cap´ıtulo 3: Experimentaci´ on Num´erica

6

x 10

−5

5 4 y 3 2 1 0 −1

0.5 1 1.5 x (b) Error absoluto producido con Disparo

(a) Gr´ afica de la soluci´ on exacta x 10

−4

1.8

x 10

−0.5

0

2

−0.5

0

2

−6

1.6 1.4 2

1.2 y

y

1 0.8

1

0.6 0.4 0.2

0 −1

−0.5

0

0.5 1 1.5 2 x (c) Error absoluto producido con Diferencias. Finitas

0 −1

0.5 1 1.5 x (d) Error absoluto producido con Bvp4c

Figura 3.11: Ejemplo 6. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 80 subintervalos, usando la norma del m´ aximo.

58

Ejemplo 6

2.5

x 10

−3

2

1.5 y 1

0.5

0 −1

(a) Gr´ afica de la soluci´ on exacta

−0.5

0

0.5 x

1

1.5

2

(b) Error relativo producido con Disparo

0.012 2

x 10

−5

0.01 1.5

0.008 y 0.006

y

1

0.004 0.5

0.002 0 −1

−0.5

0

0.5 1 1.5 2 x (c) Error relativo producido con Diferencias. Finitas

0 −1

−0.5

0

0.5 x

1

1.5

2

(d) Error relativo producido con Bvp4c

Figura 3.12: Ejemplo 6. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 80 subintervalos, usando la norma del m´ aximo.

59

Cap´ıtulo 3: Experimentaci´ on Num´erica N = 20

M´etodo Disparo Dif. Finitas Bvp4c

Tiempo 0.049 0.046 1.070

N = 40

Disparo Dif. Finitas Bvp4c

0.053 0.049 1.136

N = 80

Disparo Dif. Finitas Bvp4c

0.062 0.059 1.230

||ε||∞ 6.304 · 10−4 3.210 · 10−3 3.095 · 10−6

||ε||L2 5.778 · 10−3 4.085 · 10−3 5.356 · 10−3

5.792 · 10−5 2.010 · 10−4 1.769 · 10−6

2.995 · 10−4 2.551 · 10−4 3.359 · 10−4

8.948 · 10−5 8.038 · 10−4 1.901 · 10−6

1.282 · 10−3 1.020 · 10−3 1.340 · 10−3

Cuadro 3.11: Ejemplo 6. Tiempo empleado y errores obtenidos para un mallado de 20, 40 y 80 subintervalos.

el error aumenta hasta la proximidades de x = 0.75, punto en el cual aproximadamente alcanza su m´ aximo, para posteriormente empezar a disminuir. Para el m´etodo Bvp4c, el error se comporta de manera oscilatoria hasta una vecindad de x = 0.5, punto donde alcanza su m´ aximo, y luego empieza a disminuir de forma estable. Es importante acotar, que para las tres mallas iniciales dadas, el m´etodo Bvp4c no realiza proceso de adaptatividad, con lo que se concluye que las estimaciones usando polinomios c´ ubicos son buenas. En particular, trabajando en una malla uniforme de 21 puntos, el m´etodo Bvp4c usa un tiempo de 1.070 segundos para realizar las operaciones relacionadas con el c´ alculo de la aproximaci´ on de −6 la soluci´ on, cometiendo un error m´ aximo de 3.095 · 10 . Para requerir un tiempo aproximadamente igual, tanto Disparo como Diferencias Finitas necesitar´ an aproximar la soluci´ on sobre el n´ umero de intervalos que se muestran en el Cuadro 3.12.

M´etodo

Intervalos

Disparo

3841

Dif. Finitas

3092

||ε||∞

1.021 · 10−5 1.347 · 10−7

Cuadro 3.12: Ejemplo 6. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 1.070 segundos.

En virtud del Cuadro 3.11 y del Cuadro 3.12, se obtiene, al igual que el Ejemplo 5, que el error que proporciona el m´etodo del Disparo se sigue manteniendo igual, pese a aumentar considerablemente el total de subintervalos de la partici´ on del dominio donde est´ a definido el PVF. Todo lo contrario ocurre con Diferencias Finitas, ya que su error disminuy´ o relativamente a medida que se consider´ o m´ as intervalos en la partici´ on del dominio. 60

Ejemplo 7

Ejemplo 7 Para finalizar, consideremos ahora el PVF no lineal: p p 2 − 3p + x3 p + x2 ) x(px p + x p y ′′ = −y 2 + , (p + x2 )5 y(−0.1) = − √ y(0.1) =

√

0.1 , p + 0.01

0.1 , p + 0.01

con soluci´ on y(x) = p

x p + x2

,

donde p = 10−5 . La soluci´ on gr´ afica del PVF planteado se ilustra en la Figura 3.13. Nuevamente, al igual que el ejemplo anterior, se realiza una partici´ on uniforme del intervalo [−1, 2], pero esta vez en 200, 300 y 500 subintervalos. En virtud del Cuadro 3.13, se aprecia que los errores, empleando cada una de las dos N = 200

M´etodo Disparo Dif. Finitas Bvp4c

Tiempo 0.157 0.079 1.249

N = 300

Disparo Dif. Finitas Bvp4c

0.214 0.097 1.383

N = 500

Disparo Dif. Finitas Bvp4c

0.333 0.134 1.618

||ε||∞ 9.359 · 10−5 6.934 · 10−3 9.266 · 10−5

||ε||L2 4.901 · 10−4 6.614 · 10−5 4.901 · 10−4

2.709 · 10−6 1.151 · 10−3 2.153 · 10−6

7.874 · 10−5 4.047 · 10−6 7.874 · 10−5

1.685 · 10−5 3.171 · 10−3 1.613 · 10−5

2.184 · 10−4 1.901 · 10−5 2.184 · 10−4

Cuadro 3.13: Ejemplo 7. Tiempo empleado y errores obtenidos para un mallado de 200, 300 y 500 subintervalos. normas introducidas, no son los esperados para la gran cantidad de subintervalos de las tres particiones. De nuevo, las estimaciones aportadas por Disparo y Bvp4c siguen siendo mejores que las aportadas por Diferencias Finitas, debido a que ambos presentan un orden de aproximaci´ on de O(h4 ), contra el orden 2 de aproximaci´ on O(h ) que presenta Diferencias Finitas. Usando un tama˜ no de paso h = 0.0004, es decir, dividiendo uniformemente el intervalo [−0.1, 0.1] en 500 subintervalos, el comportamiento del error que se produce se refleja en la Figura 3.13. Para el m´etodo del Disparo, que amerita tres iteraciones por parte del m´etodo de Newton Raphson, el error tiene un comportamiento lineal, pero alrededor de x = 0 presenta oscilaciones, esto debido a la alta pendiente que presenta la soluci´ on en un entorno de ese punto. Para el m´etodo de Diferencias Finitas (que requiere tres iteraciones por parte del m´etodo de Newton para sistemas de ecuaciones no lineales), y el m´etodo Bvp4c, el error tiene un comportamiento estable en casi todo en intervalo [−0.1, 0.1], presentando una mayor magnitud en un entorno de x = 0, donde presenta oscilaciones considerablemente altas. Es importante resaltar que en las tres mallas iniciales dadas, el m´etodo Bvp4c no realiza proceso de adaptatividad. En particular, cuando se trabaja en una malla uniforme que consta de 201 nodos, el 61

Cap´ıtulo 3: Experimentaci´ on Num´erica

3

x 10

−6

2.5 2 y 1.5 1 0.5 0 −0.1 −0.08 −0.06 −0.04 −0.02

0 0.02 0.04 0.06 0.08 x (b) Error absoluto producido con Disparo

(a) Gr´ afica de la soluci´ on exacta 1.2

x 10

−3

2.5

1

x 10

0.1

−6

2

0.8 1.5

y 0.6

y 1

0.4 0.2 0 −0.1 −0.08 −0.06 −0.04 −0.02

0 0.02 0.04 0.06 0.08 0.1 x (c) Error absoluto producido con Diferencias. Finitas

0.5

0 −0.1 −0.08 −0.06 −0.04 −0.02

0 0.02 0.04 0.06 0.08 x (d) Error absoluto producido con Bvp4c

0.1

Figura 3.13: Ejemplo 7. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error absoluto obtenidos con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 500 subintervalos usando la norma infinita.

62

Ejemplo 7

1.6

x 10

−5

1.4 1.2 1 y 0.8 0.6 0.4 0.2 0 −0.1

(a) Gr´ afica de la soluci´ on exacta 4

x 10

−0.05

0 x

0.05

0.1

(b) Error relativo producido con Disparo

−3

1.2

x 10

−5

3.5 1 3 0.8 2.5 y

y

2

0.6

1.5 0.4 1 0.2 0.5 0 −0.1

−0.05

0 x

0.05

(c) Error relativo producido con Diferencias. Finitas

0.1

0 −0.1

−0.05

0 x

0.05

0.1

(d) Error relativo producido con Bvp4c

Figura 3.14: Ejemplo 7. De arriba hacia abajo y de izquierda a derecha se muestran la soluci´ on exacta y el error relativo obtenidos con Disparo, Diferencias Finitas y Bvp4c, para un mallado de 500 subintervalos, usando la norma del m´ aximo.

63

Cap´ıtulo 3: Experimentaci´ on Num´erica m´etodo Bvp4c usa un tiempo de 1.249 segundos para realizar las operaciones involucradas, cometiendo un error m´ aximo de 9.266 · 10−5 . Para que Disparo y Diferencias Finitas requieran un tiempo similar, necesitan emplear el n´ umero de intervalos que vienen dados en el Cuadro 3.14. M´etodo

Intervalos

Disparo

1889

Dif. Finitas

4250

||ε||∞

7.399 · 10−7 1.585 · 10−5

Cuadro 3.14: Ejemplo 7. N´ umeros de intervalos necesarios y error obtenido por los m´etodos para requerir un tiempo pr´ oximo a 1.249 segundos. Como se logra apreciar, pese a aumentar ampliamente la cantidad de intervalos de la partici´on del dominio donde est´ a definida la soluci´ on del PVF, los errores s´ olo disminuyen levemente. El m´etodo del Disparo logra aportar mejores estimaciones en menos subintervalos que el m´etodo de Diferencias Finitas.

64

A Pruebas de algunas afirmaciones

A.1.

Correspondientes al M´ etodo del Disparo

Teorema A.1 Sea h = (b − a)/N , xi = a + nh, con i = 0, 1, . . . , N. Sean u(x), v(x) y y(x) las soluciones respectivas de (2.2), (2.3) y (2.4). Si u ei y vei son aproximaciones de O(hn ) para ui = u(xi ) y vi = v(xi ), respectivamente, entonces wi ser´ a una aproximaci´ on de O(hn ) para yi = y(xi ). En particular, vi n |wi − yi | ≤ Kh 1 + , vN para alguna K apropiada.

Demostraci´ on: Definamos

ei = wi − yi , (u)

ei

(v)

ei

= u ei − ui , = vei − vi ,

para i = 0, 1, . . . , N. En virtud de (2.4), en cada nodo xi , se tiene, yi = ui + Cvi , donde C=

(A.1)

β − u(b) . v(b)

Por otra parte, donde

Restando (A.1) a (A.2),

evi , wi = u ei + Ce eN e = β−u C . veN

evi − Cvi wi − yi = u ei − ui + Ce

evi − Cv e i + Cv e i − Cvi = u ei − ui + Ce

e vi − vi ) + (C e − C)vi . = u ei − ui + C(e 65

(A.2)

Ap´endice Por tanto, (u)

ei = ei Para j = N , notemos que eN = 0, as´ı,

(v)

e + Ce i

e − C)vi . + (C

(u) e (v) eN + Ce N e . C −C = − vN

Sustituyendo (A.4) en (A.3),

(A.3)

(A.4)

vi (v) (u) (v) e e ei = + Cei − eN + CeN . vN (u) (v) Tomando valor absoluto, y sabiendo que ei ≤ M1 hn y ei ≤ M2 hn , para ciertas constantes positivas M1 y M2 , tenemos vi n n n n e e |ei | ≤ M1 h + C M2 h + M1 h + C M2 h vN (u) ei

=

e M2 hn 1 + vi . M1 + C vN

e M2 , concluimos que wi aproxima a yi con un orden O(hn ). As´ı, eligiendo K = M1 + C

A.2.

Correspondientes al M´ etodo de Diferencias Finitas

Teorema A.2 Supongamos que p(x), q(x) y r(x) son continuas en [a, b]. Si q(x) ≥ 0 en [a, b], entonces el sistema lineal tridiagonal (2.17) tiene una soluci´ on u ´nica siempre y cuando h < 2/L, donde L = m´ ax |p(x)|. a≤x≤b

Demostraci´ on: Basta verificar que se cumplen las hip´ otesis del teorema (1.2). Como L < 2/h, entonces, para cada i = 1, 2, . . . , n, se tiene que |p(xi )| < 2/h, con lo que, −1 < De la u ´ltima desigualdad se obtiene que 0 < 1 −

h p(xi ) < 1 2

h h p(xi ) = |δi | y 0 < 1 + p(xi ) = |ǫi |. Por tanto, 2 2

1. Para i = 2, 3, . . . , n − 1, se verifica que ǫi δi 6= 0. 2. Usando (A.5) y sabiendo que q(x) ≥ 0, en [a, b], se tiene, |δ1 | = 1 −

h p(x1 ) 2

< 2 ≤ 2 + h2 q(x1 ) = |γ1 |.

3. Para i = 2, . . . , n − 1,

|ǫi | + |δi | = 1 +

h h p(xi ) + 1 − p(xi ) 2 2

= 2 ≤ 2 + h2 q(xi ) = |γi |.

(A.5)

Secci´ on A.2: Correspondientes al M´etodo de Diferencias Finitas 4. Nuevamente, usando (A.5), se obtine, |ǫn | = 1 +

h p(xn ) 2

< 2 ≤ 2 + h2 q(xn ) = |γn |.

As´ı, el sistema (2.17) tiene soluci´ on u ´nica.

Teorema A.3 Sean Q∗ tal que q(x) ≥ Q∗ > 0, para a ≤ x ≤ b, L = m´ ax p(x) , M3 = m´ ax y (3) (x) y a≤x≤b a≤x≤b ax y (4) (x) . Entonces, para i = 0, 1, . . . , n + 1, se cumple: M4 = m´ a≤x≤b

! M + 2LM 4 3 2 wi − y(xi ) ≤ h . 12Q∗

Demostraci´ on: Para i = 0, 1, . . . , n + 1, pongamos ei = wi − y(xi ), y sea e = Existen ξi , ηi ∈ (xi−1 , xi+1 ) tales que, y(xi+1 ) − 2y(xi ) + y(xi−1 ) h2

m´ ax |ei |.

0≤i≤n+1

y(xi+1 ) − y(xi−1 ) = p(xi ) + q(xi )y(xi ) + r(xi ) 2h h2 (3) (4) − 2p(xi )y (ηi ) − y (ξi ) , 12

Se sabe que,

wi+1 − 2wi + wi−1 wi+1 − wi−1 = p(xi ) + q(xi )wi + r(xi ). h2 2h

(A.6)

(A.7)

Restando (A.6) de (A.7) y usando la definici´ on de ei , ei+1 − 2ei + ei−1 ei+1 − ei−1 h2 (3) (4) = p(xi ) + q(xi )ei + 2p(xi )y (ηi ) − y (ξi ) . h2 2h 12 Multiplicando la igualdad anterior por h2 y manipulando algebr´ aicamente se obtiene, h h h4 (4) 2 (3) (2 + h q(xi ))ei = 1 − p(xi ) ei+1 + 1 + p(xi ) ei−1 − y (ξi ) − 2p(xi )y (ηi ) , 2 2 12 para i = 1, 2, . . . , n. Por tanto, h h h4 2 (4) (3) (2 + h q(xi ))|ei | ≤ 1 − p(xi ) e + 1 + p(xi ) e + |y (ξi )| + 2|p(xi )||y (ηi )| 2 2 12 ≤ 2e +

h4 M4 + 2LM3 , 12

Notemos que e0 = en+1 = 0. As´ı, la desigualdad anterior es v´ alida para todo i = 0, 1, . . . , n + 1. Adem´ as h2 q(xi ) ≥ h2 Q∗ . Por lo tanto, h4 2 ∗ M4 + 2LM3 , h Q e≤ 12 de donde, h4 e≤ M4 + 2LM3 . 12Q∗

Referencias

[1] B´eatrice Rivi`ere, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equation, SIAM (2008) 3-17. [2] Charles W. Curtis, Linear Algebra: An Introductory Approach, Editorial Board, 1974. [3] B. Dayanad Reddy, Functional Analysis and Boundary-Value Problems: an Introductory treatment. John Wiley & Sons, Inc. 1986. [4] B. Rivi`ere, M. Wheeler, and V. Girault, Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems. Part I, Computational Geosciences, 3 (1999), pp. 337 - 360. [5] C. Dawson, S. Sun, and M. Wheeler, Compatible algorithms for coupled flow and transport, Computer Methods in Applied Mechanics and Engineering, 193 (2004), pp. 2565 - 2580. [6] David Kincaid y Ward Cheney, An´ alisis Num´erico: Las matem´ aticas del c´ alculo cient´ıfico. Addison Wesley Iberoamericana, 1994. [7] Dennis G. Zill, Ecuaciones diferenciales con aplicaciones de modelado, 7a. ed., Thomson Learning, 2002. [8] D. N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM Journal on Numerical Analysis, 19 (1982), pp. 742 - 760. [9] Eugene Isaacson y Helbert Bishop Keller, Analysis of Numerical Methods, John Wiley & Sons, Inc., 1966. [10] J. Oden, I. Babu˘ska, and C. Baumann, A discontinuous hp finite element method for diffusion problems, Journal of Computational Physics, 146 (1998), pp. 491 - 519. [11] J. Kierzenka, Studies in the Numerical Solution of Ordinary Differential Equations, PhD thesis, Southern Methodist University, Dallas, TX, 1998. [12] Kendall E. Atkinson, An introduction to NUMERICAL ANALYSIS, Second Edition, John Wiley & Sons, Inc., 1989. [13] Keller, H. B., Numerical methods for two-point boundary-value problems, Blaisdell, Waltham, MA, 1968. 69

Referencias [14] Lawrence F. Shampine, Jacek Kierzenka and Mark W. Reichelt, Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c. El tutorial y programas est´ an disponibles en http : //www.mathworks.com/bvp tutorial. [15] L. Delves and C. Hall, An implicit matching principle for global element calculations, Journal of the Institute of Mathematics and its Applications, 23 (1979), pp. 223 - 234. [16] M. F. Wheeler, An elliptic collocation - finite element method with interior penalties, SIAM Journal on Numerical Analysis, 15 (1978), pp. 152 - 161. ˇ ın, Partial Differential Equations and the Finite Element Method. John Wiley & Sons, Ltd. [17] Pavel Sol´ 2006 [18] Richard L. Burden y J. Douglas Faires, An´ alisis num´erico, 7a. ed., Thomson Learning, 2002. [19] Zienkiewicz and Cheung, The Finite Element Method in Structural and Continuum Mechanics, 1967.

70

Solución Numérica de Problemas de Valor de Frontera para Ecuaciones Diferenciales Ordinarias 1

Recommend Stories

Story Transcript

Get in touch

Social