= qn q = 0
0 0 En este caso el dominio comprende = f0 x Lxg y el contorno = fx = 0;x = Lxg. Esta ecuación diferencial ordinaria puede resolverse analiticamente, y servirá para comparar con las soluciones aproximadas que obtendremos al aplicar un procedimiento discretio para resolverla.
2
X 1.2 Discretización Las soluciones numéricas a ecuaciones diferenciales ordinarias o a derivadas parciales por aplicación de procedimientos numéricos llevaran implícitos los conceptos de aproximación y de discretización. El concepto de aproximación estará involucrado porque se seleccionará una familia de funciones para aproximar a la función solución (analítica) con un procedimiento que guarda cierta similitud una combi- nación lineal. La selección del conjunto de funciones de aproximación es un punto clave para la solución de ecuaciones diferenciales utilizando métodos numéricos. Es el primer paso para replantear el problema mátemático a resolver, con el …n de obtener una forma puramente algebraica que solo involucre operaciones matemáticas básicas. Comenzaremos por elegir un conjunto de funciones de manera tal que la función de aproximación satisfaga las condiciones de borde exactamente. Si podemos encontrar una función que satisfaga las condiciones de borde en , y un conjunto de funciones de prueba Nm (m = 1::M, con M …nito) elegidas de manera tal que Nm = 0 en para todos los m. El conjunto de funciones de prueba Nm tiene la la función …nalidad de ajustar a la solución en el dominio. Entonces para todos los puntos del dominio ^ solución puede aproximarse por : ' ^ = + M amNm (11) m=1 am es un conjunto de parámetros que hay que calcular para obtener un buen ajuste. Las funciones de prueba también son conocidas como funciones de base o funciones de forma. La condición para elegir a este conjunto de ecuaciones es que si se aumenta el número M se tenga la garantía de obtener una mejor solución aproximada, con lo cual la función ^ convergerá a . Una condición para que se logre la convergencia es que el conjunto de funciones elegido tenga la posibilidad de representar cualquier variación de la función en el dominio . Esto es conocido como requerimiento de completitud. El concepto de discretización aparece cuando se reemplazan los in…nitos puntos donde se necesita conocer la función incógnita por un número …nito de parámetros desconocidos. En general este conjunto de parámetros estará asociado al conjunto de funciones que utilizaremos para aproximar a la función solución que estamos buscando, en el caso (11) serán los am. Solo resta establecer una metodología (o procedimiento) que nos permita determinar los valores de los parámetros. Con este …n se presenta, luego del ejemplo, el método de aproximación por residuos ponderados. EJEMPLO 1: Este ejemplo tiene la …nalidad de encontrar una función aproximante paramétrica como la (11) a una función dada. Supongamos que tenemos la función: = (0:1 + x=3) sin(1:7 x) Se propone la siguiente familia, o conjunto de funciones, para realizar el ajuste: = = 0:3505740306 x [ (0:1 + 1=3) sin(1:7 1)]x N1 N2 = = x(1 x2 (1 x) x) ^ = + a1 N1 + a2 N2
Como debemos calcular dos parámetros, pondremos las siguientes condiciones: ^ (1=3) = (1=3) y ^ (2=3) = (2=3), lo cual permitirá plantear un sistema de dos ecuaciones con dos incógnitas, resolverlo y obtener el valor de los parámetros a1 y a2. 3
X X X Z 2 9 2 9 2 27 4 27 a1 a1 = 0:2064978 +0:1310596 0:350574 1=3 0:350574 2=3 a1 a1 = 2:44824857 2:97944192 ^ = 0:3505740306 x + 2:44824857x(1 x) + 2:97944192×2 (1 x) 1.3 Aproximación por residuos ponderados Lo primero que se debe de…nir es el error o residuo de la aproximación (11), al que denominaremos R : R = ^ (12) Aqui debemos observar que R también dependerá de la posición x dentro del dominio , (dado que y ^ dependen de x). Lo que se quiere hacer es disminuir este residuo de una manera general en todo el dominio . Para ello vamos a requerir que un conjunto apropiado de integrales del error sobre el dominio pesandas (ponderadas) de diferentes maneras, sea cero. Matemáticamente: Z Wl R d = Z Wl ^ d =0 con l = 1::M (13) Wl es un conjunto de funciones de peso independientes. Para que se cumpla la condición general de convergencia ^ ! cuando M ! 1 , consiste en exigir que (13) se satisfaga para toda Wl cuando M ! 1. Esto será cierto solo si R ! 0 en todos los puntos del dominio (que es lo que se desea). Cuando se reemplaza (11) en (13) se obtiene un sistema de ecuaciones algebraicas lineal. Resolviendo el sistema de ecuaciones algebraicas obtenemos los valores de los parámetros am que se estaban buscando. Concretamente: Z Wl ^ d = Z Wl " M m=1 # amNm d (14) = Z Wl [ ]d Z Wl " M m=1 # amNm d = Z Wl [ ]d M m=1 am Z WlNmd = f Ka Donde: aT = (a1;a2;:::;aM) (15) Klm = Wl Nm d 1 l M y 1 m M fl = Z Wl ( ) d 1 l M (16) l sistema a resolver, escrito de la manera convencional es:
Ka = f
4
Donde K es la matriz de coe…cientes, a el vector de incógnitas y f el vector de términos independientes. Para poder resolver las integrales, previamente debemos elegir al conjunto de funciones de peso Wl. Varios conjuntos de funciones de peso pueden ser utilizadas, cada uno conduce a un método de apro- ximación por residuos ponderados diferente. A continuación se enumeran algunos de los más utilizados, para problemas unidimensionales. 1.3.1 Colocación por puntos Aqui el conjunto de funciones de peso Wl es:
Wl = (x xl) (17) donde (x xl) es la función Delta de Dirac, que tiene las siguientes propiedades: (18) (x (x xl) = 0 si x 6= xl xl) = 1 si x = xl x>xl G(x) (x xl) dx = G(xl) x< xl
Lo cual implica que R = 0 en un cierto número de puntos xl. Las componentes de la matriz K y del vector f que se obtienen al reemplazar (17) en (13) obteniendo las siguientes expresiones: Klm = Nm(xl) (19) fl = [ ]x=xl Nm(xl) :es la función de prueba valuada en xl [ ]x=xl :es la diferencia de las funciones valuada en xl 1.3.2 Colocación por subdominios En este caso, el conjunto de funciones de peso Wl se elige para que satisfaga las siguientes condiciones: Wl = 1 0 si xl < x < xl+1 si x < xl o x > xl+1 (20) Por propiedad de las integrales la (13) simplemente requiere que la suma de la integral del error sobre M subdominios que componen a sea cero. Reemplazando la función de peso en (13) se obtienen las Klm = l siguientes expresiones para las componentes de la matriz K y del vector f son en este caso:
xZ +1 Nm dx (21) fl = l xl xZ +1 xl ( ) dx 5
1.3.3 Mínimos cuadrados En general el método de mínimos cuadrados no es presentado como un método de residuos ponderados, pero puede mostrarse que pertenece a esta clase de métodos. La aproximación general por mínimos cuadrados es tratar de minimizar la suma de los cuadrados del residuo (o error) en cada punto del dominio . En este caso se requiere la minimización de: I (a1;a2;:::;aM) = Z ^ 2 d (22) tal que, se plantearemos que @I @al = 0 para l = 1;2;:::M (23) Llevando a cabo el proceso de derivación, y teniendo en cuenta (11) se obtiene: = Nl (24) @I @al Puede mostrarse que I es un mínimo cuando: Z ^ Nl d =0 (25) Klm = Esta última expresión es de la misma forma que (13) si se considera que Wl = Nl. Las expresiones para las componentes: Z Nl Nm d (26) fl = Z Nl ( ) d 1.3.4 Galerkin Klm = Este es sin dudas es el mas popular de los métodos de residuos ponderados. En este caso se eligen como funciones de peso a las mismas funciones que se utilizan como funciones de forma, es decir Wl = Nl. Haciendo el reemplazo en (13) las componentes de la matriz K y del vector f son: Z Nl Nm d (27) fl = Z Nl ( ) d Observaciones: | Este método tiene la ventaja computacional de que la matriz K resulta simétrica (generalmente). | Las expresiones (26) y (27) resultan ser equivalentes. | Si el conjunto de funciones de prueba Nm = sin(m x=Lx) para m = 1;2;:::M y = fxn 0 x Lxg se obtiene que: Klm = L Zx 0 sin l x Lx sin m x Lx dx fl = L Zx 0 ( ) sin m x Lx dx 6
= x (1 que al realizar las integrales da: Klm = Lx=2 si l = m 0 si l 6= m dando lugar a un sistema diagonal, resultando inmediatamente en la solución: am = 2 Lx L Zx ( ) sin m x Lx dx para m = 1;2;:::m 0 Que es idéntica a la representación de una serie seno de Fourier truncada, donde la particular simpli- cidad se debe a la propiedad de ortogonalidad de las funciones de forma. EJEMPLO 2: Encontraremos ahora las aproximaciones a la función del ejemplo 1 empleando el método de residuos = [0;1] ponderados utilizando las distintas funciones de peso presentadas. La función a aproximar en es: = (0:1 + x=3) sin(1:7 x) Se propone la siguiente aproximación a la función, que será utilizada en todos los casos: ^ = + a1 N1 + a2 N2 = 0:3505740306 x N1 N2 = x(1 2 x) x) Colocación por puntos Se adoptan los siguientes dos puntos de colocación: x1 = 0:25 y x2 = 0:75 con lo cual M = 2 Klm = Nm(xl) y fl = [ ]x=xl K11 = K12 = N1 (x1) = 0:1875 N2 (x1) = 0:046875 f1 = (x1) (x1) = 0:2659113 K21 = K22 = N1 (x2) = 0:1875 N2 (x2) = 0:140625 f2 = (x2) (x2) = 0:0032116 El sistema y su solución son: 0:1875 0:1875 0:046875 0:140625 a1 a2 = 0:2659113 0:0032116 ) a1 = 2:1358546 a2 = 2:87064426 ^ = 0:3505740306 x 2:13585478 x(1 x) + 2:87064417 x2 (1 x) Colocación por subdominios Se adoptan dos subdominios, el primer subdominio esta comprendido ente x1 = 0 y x2 = 0:5 y el segundo entre x2 = 0:5 y x3 = 1 con lo cual M = 2. Klm = l xZ +1 xl Nm dx y fl = l xZ +1 xl ( ) dx 7
= = = = 0: K11 =
K12 = 0: R5 0 R5 N1 (x)dx
N2 (x)dx = 0:083333
= 0:0260416 f1 = 0: R5 0 0 [ (x) (x)]dx = 0:11234019 1 K21 =
K22 = 1 R 0:5 R N1 (x)dx
N2 (x)dx = 0:0833333
= 0:0572916 f2 = 1 R 0:5 [ (x) (x)]dx = 0:02452481 0:5 El sistema algebraico y su solución son: 0:083333 0:0833333 0:0260416 0:0572916 a1 a2 = 0:11234019 0:02452481 ) a1 = 2:226236221 a2 = 2:8100923155 ^ = 0:3505740306 x 2:226236221 x(1 x) + 2:8100923155 x2 (1 x) Mínimos cuadrados En este caso se tiene que las expresiones para Klm y fl coinciden con las que corresponden a Galerkin: Klm = Z Nl Nm d y fl = Z Nl ( ) d 1 1 1 1 1 K11 =
K12 =
f1 =
K21 =
K22 =
f2 = 1 R 0
R R N1 (x)N1 (x)dx 0 R N1 (x)N2 (x)dx 0 N1 (x)[ (x) (x)]dx = R N2 (x)N1 (x)dx 0 R N2 (x)N2 (x)dx 0 N2 (x)[ (x) (x)]dx = 1 30
1 60
0:03053668742
1 60
1 105
0:01192414968 0 El sistema y su solución son: 1 30 1 60 1 60 1 105 a1 a2 = 0:03053668742 0:01192414968 ) a1 = 2:32066211519 a2 = 2:80912298519 ^ = 0:3505740306 x 2:32066211519 x(1 x) + 2:80912298519 x2 (1 x) En el siguiente gra…co se han superpuesto la función y las aproximaciones obtenidas. 8
k @n u ' u = X 1.4 Aproximación a la solución de ecuaciones diferenciales Consideremos que la ecuación dierencial puede ser escrita de la diguiente forma general: en (28) A(u) = L(u) + p = 0 Donde L() es un operador diferencial lineal p es una función independiente de u u variable independiente Las condiciones de borde pueden expresarse en forma general como: B(u) = M(u) + r en (29) Donde M(u) = u entonces r = u en u M(u) = @u entonces r = q en q 1.4.1 Condiciones de borde satisfechas por la elección de las funciones de prueba ^ expansión de la forma: ^ + M am Nm (30) En este caso la función m=1 será la encargada de cumplir las condiciones de borde y las funciones Nm las encargadas de ajustar a la función solución en el dominio. Estas funciones son elegidas de manera tal que cumplan con las siguientes condiciones: M( ) = r en (31) M(Nm) = 0 para m = 1;2;:::M en 9
' u = X @u X @2u @x @x X = A(u) = L(u) + p X X u ' u = X ^ los coe…cientes am. La aproximación (30) puede utilizarse para aproximar a las derivadas de u derivando y que todas sus derivadas, directamente, con tal que las funciones Nm sean continuas en el dominio que sean necesarias, existan: u ^ + M m=1 am Nm (32) @u @x ' ^
@x = @ @x + M
m=1 am @Nm @x @2u @x2 ' ^ 2 = @2 2 + M
m=1 am @2Nm @x2 Las aproximaciones a u, @u @x y @2u @x2 son reemplazadas en la ecuación diferencial A(u) para obtener la expresión del residuo. R (33) ^
^ = L( ) + M amL(Nm) + p m=1 Obtenido el residuo éste se minimizará utilizando el método de los residuos ponderados, para que R ' 0 en , lo cual lleva a plantear: Z Wl R d = Z " Wl L( ) + M # amL(Nm) + p = 0 (34) Klm = m=1 Agrupando adecuadamente los términos obtenemos:
Z Wl L(Nm) d 1 l M y 1 m M (35) fl = Z Wl p d + Z Wl L( ) d 1 l M Por lo cual obtenemos el siguiente sistema algebraico a resolver:
Ka + f = 0 1.4.2 Condiciones de borde no satisfechas por la elección de las funciones de prueba En muchos casos ocurre que encontrar una función que satisfaga las condiciones de borde, suele ser una tarea muy ardua y en ocasiones imposible, lo cual limita el rango de funciones de prueba admisibles. Para salvar este inconveniente se considerará un conjunto de funciones que no satisfagan las condiciones ^ ^ M am Nm (36) m=1 ^ en el dominio R . Tenemos ahora:
10
= A(u) = L(u) + p = B(u) = M(u) + r ' u = X @u X @2u X h i Wl L(u) + p d WlL(u) d + Wl L(u) d Wl E u C (Wl ) D u d (37) R R ^ ^
^ ^ en en Ahora deben ser minimizados ambos, para lo cual se propone plantear el método de los residuos ponderados de la siguiente manera: Z Wl R d + Z Wl R d = 0 (38) Las funciones de peso Wl (en el dominio) y las Wl (en el contorno) pueden ser elegidas en forma independiente. La aproximación (36) puede utilizarse para aproximar a las derivadas de u, con tal que las funciones Nm sean continuas en el dominio y que existan todas sus derivadas: u ^ M m=1 am Nm (39) @u @x ' ^
@x = M
m=1 am @Nm @x @2u @x2 ' ^ @x2 = M
m=1 am @2Nm @x2 El residuo en el dominio y en el borde se obtienen reemplazando las aproximaciones (39) en (28) y (29), para obtener: Klm = Z Wl L(Nm) d + Z Wl M(Nm) d 1 l M y 1 m M (40) fl = Z Wl p d Z Wl r d 1 l M 1.4.3 Forma débil del problema. Condiciones de borde naturales Se ha planteado la forma de encontrar el valor de los parámetros am aun cuando la solución aproximada ^ funciones a utilizar, pueden aparecer di…cultades cuando las condiciones de borde incluyen derivadas de ^ para evitar tales cálculos y proponer una forma mas general para tratar las condiciones de borde. Volviendo al planteo de residuos ponderados, se tenía que: Z Wl R d = Z ^ = Z ^ Z Wlp d El primer término pude ser integrado por partes (recordar que R udv = uv j R vdu) obteniendo: Z ^ = Z ^ d Z ^ (41) 11
Wl E u @x2 = 0 @x2 + + 1 = 0 para 0 Donde C , D y E son operadores diferenciales con un orden de diferenciación menor que L. A la expresión (41) se conoce como orma débil del problema. Si reemplazamos la expresión de la forma débil (41) en la expresión del residuo (38) y se eligen convenientemente las funciones de peso Wl el último término de (41) podrá cancelarse con parte del término que corresponde al residuo en el contorno. Este procedimiento es aplicable para tratar con las condiciones de borde naturales del problema, no asi para las condiciones esenciales. Z q Wl r d $ Z q ^ d Debido a que los nuevos operadores diferenciales son de menor orden, son menores también los re- querimientos que deben cumplir las funciones de prueba, lo cual es una ventaja adicional. Con el …n de aclarar los conceptos que han sido presentados de manera abstracta, plantearemos un par de problemas simples a los que aplicaremos los procedimientos recien vistos, estos se presentan en los siguientes ejemplos. EJEMPLO 3: Vamos a tratar con el siguiente problema de transferencia de calor unidimensional con una fuente distribuida gobernado por el siguiente sistema de ecuaciones: @2 @x2 + +1=0 = 0 en x = 0 = 1 en x = 1 Solución analítica La solución homogenea es: h(x) = Acos(x) + B sin(x) La solución particular: p(x) = @2 0+ +1=0 ) = 1 La solución general entonces (x) = h(x) + p(x) (x) = Acos(x) + B sin(x) 1 Utilizando las condiciones de contorno para encontrar las constantes A y B x =0 ; = 0 ) 0 = Acos(0) + B sin(0) 1 ) A =1 = 1 ) 1 = 1cos(1) + B sin(1) 1 ) B = 1:734697595 x = 1; La solución entonces es: 1 (x) = cos(x) + 1:734697595sin(x)
Soluciones aproximadas por residuos ponderados Con funciones que satisfacen las condiciones de borde El problema esta de…nido por: x 1 @2 = 0 en x = 0 = 1 en x = 1 12
| Página siguiente |