Resolución Numérica de problemas de contorno no lineales con métodos de gradiente conjugado

Trabajo Fin de Grado Resolución Numérica de problemas de contorno no lineales con métodos de gradiente conjugado Julio Suárez Gómez 2020/2021 UNIVERSIDAD DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Resolución Numérica de problemas de contorno no lineales con métodos de gradiente conjugado Julio Suárez Gómez data presentación UNIVERSIDADE DE SANTIAGO DE COMPOSTELA ii iii iv Trabajo propuesto Área de Coñecemento: Matemática Aplicada Título: Resolución numérica de problemas de contorno no lineales con métodos de gradiente conjugado. Breve descrición do contido: Una de las estrategias para abordar la resolución numérica de los problemas de contorno de ecuaciones diferenciales elípticas no lineales es la formulación como un problema de mínimos cuadrados de su discretización con elementos finitos. Esta técnica está descrita en el artículo de Glowinski y Reinhart “Continuationconjugate gradient methods for the least squares solution of nonlinear boundary value problems”. El objetivo de este trabajo es el estudio de dicho artículo y la elaboración de un programa Matlab que resuelva problemas de contorno para ecuaciones diferenciales no lineales con un operador diferencial de tipo elíptico mediante la discretización con elementos finitos, la formulación del problema en mínimos cuadrados y el método de gradiente conjugado precondicionado. Recomendacións: Haber superado la materia Métodos Numéricos en Optimización y Ecuaciones Diferenciales y estar matriculado o haber superado Métodos Numéricos en EDP. Outras observacións: - - Índice general Resumen viii Introducción xi 1. Resultados y definiciones previas 1 2. Métodos de gradiente conjugado con Polak-Ribière 5 2.1. Introducción a la Búsqueda Lineal . . . . . . . . . . . . . . . . . . . . 7 2.2. Selección de la longitud del paso . . . . . . . . . . . . . . . . . . . . . 7 2.2.1. Búsqueda inexacta: Condiciones de Goldstein y Wolfe-Powell . 9 2.2.2. Regla de Wolfe-Powell . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3. Regla de Armijo-Goldstein . . . . . . . . . . . . . . . . . . . . 11 2.2.4. Método de interpolación cuadrática con tres puntos . . . . . . 12 2.3. Métodos de direcciones conjugadas . . . . . . . . . . . . . . . . . . . 14 2.4. Método de gradiente conjugado lineal . . . . . . . . . . . . . . . . . . 15 2.5. Método de gradiente conjugado no lineal . . . . . . . . . . . . . . . . 17 3. Resolución de problemas elípticos 1D mediante el método de elementos finitos 23 3.1. Introducción................................ 23 3.2. Descripción del MEF 1D . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2.1. Introducción al problema de Sturm-Liouville . . . . . . . . . . 24 3.2.2. Formulaciones diferencial y variacional del problema . . . . . . 24 3.2.3. Discretización del problema variacional . . . . . . . . . . . . . 25 3.2.4. Resolución del método de elementos finitos Langrange para K=1 27 3.2.5. Implementación de las condiciones de contorno . . . . . . . . . 36 v vi ÍNDICE GENERAL 4. Formulación del MEF para el caso no lineal 39 4.1. Introducción................................ 39 4.2. El problema de Bratu 1D . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.1. Ajuste del modelo no lineal . . . . . . . . . . . . . . . . . . . 40 4.2.2. Solución analítica del problema de Bratu . . . . . . . . . . . . 41 5. Resultados obtenidos 43 5.1. Resultados de los programas base . . . . . . . . . . . . . . . . . . . . 43 5.1.1. Gradiente conjugado . . . . . . . . . . . . . . . . . . . . . . . 43 5.1.2. Método de elementos finitos . . . . . . . . . . . . . . . . . . . 44 5.1.3. MEF con resolución mediante GC con búsqueda lineal inexacta 46 6. Conclusiones y trabajo futuro 49 7. ANEXOS 51 7.1. ANEXO I: Convergencia y existencia de los métodos de búsqueda inexacta .................................. 51 7.2. Anexo II: Existencia y unicidad de solución MEF . . . . . . . . . . . 56 7.3. Anexo III: Normas del error . . . . . . . . . . . . . . . . . . . . . . . 56 7.3.1. Normainfinito .......................... 56 7.3.2. Seminorma H1.......................... 56 7.3.3. Norma L2............................. 57 7.4. ANEXO IV: Métodos de cuadratura . . . . . . . . . . . . . . . . . . . 57 7.4.1. Cuadratura de Gauss con 2 puntos . . . . . . . . . . . . . . . 57 7.4.2. Cuadratura de Poncelet . . . . . . . . . . . . . . . . . . . . . 58 7.5. ANEXO V: Precondicionamiento en el método de gradiente conjugado 58 7.5.1. Factorización Incompleta . . . . . . . . . . . . . . . . . . . . . 59 7.5.2. Factorización incompleta de Cholesky . . . . . . . . . . . . . . 59 7.6. PROGRAMAS UTILIZADOS . . . . . . . . . . . . . . . . . . . . . . 59 Bibliografía 95 xiv INTRODUCCIÓN Capítulo 1 Resultados y definiciones previas En este capítulo inicial daremos las definiciones generales para el desarrollo posterior del trabajo. Para ellas, tomaremos como referencia [CCMT89],[SY06] y [BF84]. Definición 1.1. (Espacio de Hilbert) Dado un espacio dotado de un producto interior (V, h·,·i) decimos que este espacio es de Hilbert si el espacio métrico (V,d) es completo, donde V es un espacio normado (kvk=hv, vi2∀v∈V)y la distancia d se define de la siguiente manera: d(x, y) = kx−yk, con x, y ∈V(1.1) Definición 1.2. Definimos el espacio L2(a, b) := {u: (a.b)→Rtal que Rb au2(x)dx < ∞}. Con su norma correspondiente k·kL2:L2(a, b)→Rtal que k·kL2=Zb a u2(x)dx1 2 (1.2) y un producto escalar definido por h·,·iL2(a, b)×L2(a, b)→Rdonde hu, vi=Zb a u(x)v(x)dx (1.3) El espacio L2(a, b)es de Hilbert. Observación 1.3.Sean u, v dos funciones de L2(a, b)iguales en casi todo punto, o iguales salvo en un conjunto de medida 0, entonces u, v son consideradas funciones equivalentes en L2(a, b). 1 2CAPÍTULO 1. RESULTADOS Y DEFINICIONES PREVIAS Definición 1.4. (Espacio de Sobolev)H1(a, b) = {u(x)∈L2(a, b) : u0(x)∈ L2(a, b)}. Un espacio de Sobolev también será de Hilbert. Observación 1.5.Definimos también H1 0(a, b) = {u∈H1(a, b) : u(a) = u(b) = 0}. Definición 1.6. (Desigualdad de Cauchy-Schwarz). Sean u, v ∈L2(a, b)se cumple la siguiente desigualdad: Zb a|u(x)v(x)|dx ≤Zb a (u2(x)dx1 2Zb a v2(x)dx1 2 (1.4) Definición 1.7. Dada φ: (a, b)→R, el soporte de φ se define como supp(φ) = {x∈(a, b) : φ(x)6= 0} Definición 1.8. Se define D(a, b) := {φ: (a, b)→Rtal que φ ∈ C∞(a, b), donde supp(φ)⊂ (a, b)} Definición 1.9. Sea V⊂Rn, se dice que Ves convexo si ∀a, b ∈V y ∀t∈[0,1],(1− t)a+tb ∈V. Observación 1.10.Es sencillo demostrar que todo espacio de Hilbert es convexo [BF84]. Definición 1.11. (Convexidad) Dado un conjunto convexo V⊂Rncon V 6=∅,una función J:V→Ry un escalar α∈(0,1) se cumple que: 1. J convexa si J(αx + (1 −α)y)≤αJ(x) + (1 −α)J(y)∀x, y ∈V. 2. J estrictamente convexa si J(αx+(1−α)y)< αJ(x)+(1−α)J(y)∀x, y ∈V. 3. J uniformemente convexa si dado c > 0con c∈Rse cumple que J(αx + (1 − α)y)≤αJ(x) + (1 −α)J(y)−1 2cα(1 −α)kx−yk2∀x, y ∈V. Teorema 1.12. (Relación convexidad - gradiente). Dado un conjunto convexo V⊂Rncon V 6=∅y J :V→Runa función diferenciable se cumple que: 3 1. J convexa ⇐⇒ J(y) ≥J(x) + ∇J(x)T(y-x) ∀x, y ∈V. 2. J estrictamente convexa ⇐⇒ J(y) > J(x) + ∇J(x)T(y-x) ∀x, y ∈V, con x 6= y. 3. J uniformemente convexa ⇐⇒ dado c > 0, J(y) ≥J(x) + ∇J(x)T(yx)+1 2cky−xk2∀x, y ∈V y c ∈R. Definición 1.13. Una matriz A∈Rn×nse dice semi-definida positiva si ∀v∈Rn se cumple que vTAv ≥0. De la misma forma, diremos que es definida positiva si la desigualdad anterior es estricta vTAv > 0. Teorema 1.14. (Relación convexidad - hessiana). Dado un conjunto convexo V⊂Rncon V 6=∅y J :V→Runa función dos veces diferenciable se cumple que: 1. La matriz hessiana de J es semi-definida positiva ∀x∈V⇐⇒ J es convexa. 2. La matriz hessiana de J es definida positiva ∀x∈V⇒J es estrictamenteconvexa. 3. La matriz hessiana es uniformemente definida positiva (∃m > 0/mkuk2≤ uT∇2J(x)u, ∀x∈, u ∈Rn)∀x∈V⇐⇒ J uniformemente convexa . Observación 1.15.En el caso 2. la implicación hacia la izquierda no sería cierta, se puede tomar como contraejemplo la función J(x) = x4con x∈R, donde J00(0) = 0 y por tanto no sería estricta la desigualdad. Definición 1.16. Un punto y∈Rnes llamado mínimo global si f(y)≤f(x),∀x∈ Rn, con x 6=y. Este mínimo será estricto si la desigualdad también lo es. Definición 1.17. Dada J:Rn→Runa función continua y dos veces diferenciable en x∈Rndecimos que el vector d∈Rnes una dirección de descenso de J en x si verifica h∇J(x), di<0. Definición 1.18. Sea V un espacio normado se dice que J:V→Rnes una función coercitiva si l´ım kxkV→∞J(x) = +∞(1.5) 4CAPÍTULO 1. RESULTADOS Y DEFINICIONES PREVIAS Definición 1.19. Sea J:V→Rfunción sobre V un espacio de Hilbert, se dice que J es elíptica si es continuamente diferenciable en V y cumple que ∃α∈R+tal que h∇J(v)−∇J(u), v −ui ≥ αkv−uk2∀u, v ∈V(1.6) Teorema 1.20. (Condición para el mínimo relativo) Sea J:V⊂Rnfunción dos veces continuamente diferenciable en un abierto V. Un punto x∈Vserá un mínimo local si y solo sí ∇J(x)=0y∇2J(x)es semidefinida positiva. Teorema 1.21. Sea V un conjunto convexo sobre un espacio normado se cumple que: 1. Si J:V→Res una función convexa con un mínimo relativo en un punto x∈V, entonces tendrá un mínimo con respecto a todo el conjunto V. 2. Si J:V→Res estrictamente convexa con un mínimo relativo en x∈V, este conjunto tendrá como máximo un mínimo, además será estricto. 3. Dada J:U⊂V→Rsiendo V un conjunto abierto con U subconjunto abierto diferenciable en el punto x∈U, entonces J tendrá un mínimo en x con respecto a U si y solo si J0(x) = 0. Teorema 1.22. Sea J:Rn→Runa función diferenciable y convexa, entonces, x es un mínimo global si y solo si ∇J(x) = 0. A continuación, en un primer capítulo, utilizaremos algunos de los resultados anteriores para establecer un algoritmo de resolución del problema de optimización no lineal sin restricciones. Para ello desarrollaremos técnicas de búsqueda inexacta como complementario del método de gradiente conjugado no lineal con la variante de Polak-Ribiére. Capítulo 2 Métodos de gradiente conjugado con Polak-Ribière En este capítulo trataremos de dar solución al siguiente problema de optimización sin restricciones : Dada J :U→Rencontrar u ∈U tal que J(u) = ´ınf v∈UJ(v)(2.1) es decir, encontrar el mínimo de J(u) en U. Para ello trabajaremos durante el capítulo con el conjunto U como un espacio de Hilbert y una función elíptica J∈ C4(U). Comenzaremos enunciando dos teoremas que nos permitirán demostrar la existencia y unicidad de solución para el problema anterior. Teorema 2.1. (Existencia de solución). Sea U⊆Rnun conjunto no vacío y cerrado con J:Rn→Runa función continua, y coercitiva cuando U no está acotado. Entonces existe al menos un elemento u∈U que verifica el problema (P)u∈U y J(u) = ´ınf v∈UJ(v)(2.2) Demostración. Al ser conjunto no vacío, podemos tomar un punto cualquiera ˜u∈U, como J es una función coercitiva existirá r∈Rtal que kvk> r, de este modo J(˜u)< J(v). En consecuencia, el conjunto de soluciones del problema (P) coincide con el de las soluciones correspondientes al conjunto 5 6CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE ˜ U=U∩{v∈Rn:kvk ≤ r}(2.3) Por tanto tendremos ˜ Uun conjunto no vacío, cerrado y acotado. Uniendo esto a que J es una función continua garantizaremos la existencia de mínimo. Lema 2.2. Una función elíptica J:U→Res estrictamente convexa y coercitiva. Demostración. Tomamos u, v ∈U, dado que la función es elíptica es continuamente diferenciable, por tanto podremos escribir bajo la fórmula de Taylor (ver [CCMT89] pag.228-229) J(v)−J(u) = Z1 0h∇J(u+t(v−u)), v −uidt =h∇J(u), v −ui+ Z1 0h∇J(u+t(v−u)) −∇J(u), v −uidt > h∇J(u), v −ui. (2.4) De esta desigualdad obtenemos que la función es estrictamente convexa, pues J(v)> J(u) + h∇J(u), v −ui,∀u, v ∈U, con u 6=v. (2.5) Además, también vemos que es coercitiva, ya que h∇J(u), v −ui+Z1 0h∇J(u+t(v−u)) −∇J(u), v −uidt ≥ h∇J(u), v −ui+Z1 0 αtkv−uk2dt =h∇J(u), v −ui+α 2kv−uk2. (2.6) por tanto, J(v)≥J(0) + h∇J(0), viα 2kvk2≥J(0) −k∇J(0)kkvk+α 2kvk2.(2.7) Teorema 2.3. (Unicidad de solución). Sea V un espacio de Hilbert y J:U→Runa función elíptica, el problema (P)u∈U y J(u) = ´ınf v∈UJ(v)(2.8) tiene solución única. 2.1. INTRODUCCIÓN A LA BÚSQUEDA LINEAL 7 Demostración. La existencia de la solución en el problema (P) viene dada por el teorema anterior, ya que una función elíptica es coercitiva. Además, como una función elíptica es estrictamente convexa se obtiene la unicidad de la solución. A continuación comenzaremos describiendo las estrategias de búsqueda lineal. Estas serán necesarias para resolver el problema optimización en una dimensión, que está presente en casi cualquier método de resolución de problemas de optimización sin restricciones. Por tanto, será fundamental obtener un proceso eficiente de resolución de dichos problemas 1D para la resolución posterior mediante el método de gradiente conjugado. 2.1. Introducción a la Búsqueda Lineal En las estrategias de búsqueda lineal estableceremos una ecuación de búsqueda para aproximarnos a nuestro óptimo lo máximo posible y en el menor número de pasos posibles. Nuestra fórmula general para cada iteración vendrá dada por la siguiente expresión: xk+1 =xk+αk·dk(2.9) donde αkserá un escalar positivo que definirá el paso y dkla dirección de descenso en la iteración k. Por tanto, el método consistirá en escoger una dirección de descenso apropiada y movernos en esa misma dirección obteniendo valores xkcada vez más próximos a la solución. En cada iteración, será fundamental hallar el paso αminimizando la expresión anterior. Como podemos observar, el éxito de este algoritmo depende directamente de las selecciones de la dirección de descenso y del paso apropiados. Estudiaremos primero como seleccionar dicho paso. 2.2. Selección de la longitud del paso Para seleccionar el paso de forma exacta deberíamos escoger el mínimo de la función φ(α)con α > 0, φ(α) = J(xk+α·dk)(2.10) 8CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE pero este proceso podría ser excesivamente costoso, sobre todo si el punto inicial está lejos de la solución, dado que se requerirán de muchas evaluaciones de la función objetivo Jy de su gradiente. Por estas cuestiones estudiaremos un proceso de búsqueda inexacta, donde intentaremos aproximar el mínimo sin calcular ∇J. Para comenzar, discutiremos las condiciones que tendrá que cumplir la función objetivo para que cada etapa de iteración sea un avance hacia el mínimo. La primera y obvia de las condiciones es que la función objetivo Jse debe ir reduciendo, es decir J(xk+1)< J(xk)o lo que es lo mismo J(xk+1)< J(xk+αk·dk). Aun así, está condición no garantiza la convergencia al mínimo, pues se podría tener una secuencia de iterantes de la forma J(xk) = 1/k y que el mínimo fuera negativo, por lo que la reducción de esta función sería insuficiente. Figura 2.1: Descenso insuficiente de J A continuación describiremos algunas condiciones para garantizar un decrecimiento suficiente. 2.2. SELECCIÓN DE LA LONGITUD DEL PASO 9 2.2.1. Búsqueda inexacta: Condiciones de Goldstein y WolfePowell Como hemos mencionado anteriormente el proceso de búsqueda lineal puede ser a veces demasiado costoso, esto se debe mayoritariamente a tres factores: la propia función, la dirección de descenso o el iterante inicial seleccionado. Para ello, enunciaremos una condición que garantizará el decrecimiento suficiente de Jen el algoritmo, y así evitar los problemas relacionados con la reducción de la función objetivo. J(xk+αkdk)≤J(xk) + c1αk∇(Jk)Tdk, con c1∈(0,1) constante. (2.11) Nótese que para c1= 1 se tiene el lado derecho como una aproximación por Taylor de primer orden en el punto inicial. La condición anterior es conocida como Condición de Armijo. Esta condición que nos permite acotar el valor del nuevo punto J(xk+αkdk)por debajo del valor previo, ya que c1es positivo y por lo tanto el lado derecho de la desigualdad nos muestra un valor inferior a J(xk). Sin embargo, la condición de Armijo provoca que la reducción sea proporcional al tamaño del paso, lo cual funciona muy bien para pasos grandes, pero podría ser ineficiente para pasos pequeños no llegando a la convergencia. Por lo tanto, es necesario añadir una nueva restricción sobre el paso. Para solucionar este problema distinguiremos entre dos tipos de condiciones, dividiendo así en dos métodos alternativos: La Regla de Goldstein y el Método de Wolfe-Powell. Estos últimos se diferenciarán entre sí por el uso de la derivada, lo que hará a cada método útil en su contexto. 2.2.2. Regla de Wolfe-Powell Para obtener el método de Wolfe-Powell definiremos la Condición de curvatura, que exigirá a αkcumplir la siguiente desigualdad: ∇J(xk+αkdk)Tdk≥c2∇JT kdk, con c2∈(c1,1) constante. (2.12) Vemos que el lado izquierdo es φ0(αk), por tanto podemos ver que lo que se busca es que la pendiente sea c2veces mayor en el nuevo punto de la iteración. Esto tiene 16CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE Escogiendo los βidel mismo modo que (2.27) tendremos con j= 0,1, ..., k −1 βj=gT kGdj dT jGdj =gT k(gj+1 −gj) dT j(gj+1 −gj)(2.29) De esta forma βj= 0, j = 0,1, ..., k −2(2.30) βk−1=gT kgk gT k−1gk−1 (2.31) Estableciendo por tanto el método como xk+1 =xk+αkdk(2.32) dk=−gk+βk−1dk−1(2.33) Donde βk−1se calcula mediante la fórmula de Fletcher-Reeves 2.31 y el paso, en particular del caso cuadrático, será αk=−gT kdk dT kGdk (2.34) Como podemos observar, en los métodos de gradiente conjugado es fundamental seleccionar un buen βk−1, en nuestro caso introduciremos la fórmula de PolakRibière(PR) con la finalidad de optimizar dicho algoritmo para el caso lineal y no lineal. La fórmula de PR se considera mucho más eficiente ya que tiene la propiedad de reinicio automático, es decir, si el algoritmo va muy lento βk≈0, propiedad que en el caso de FR no tenemos. Notamos que nuestro αkvendrá delimitado por las técnicas de búsqueda lineal exacta e inexacta del los anteriores capítulos. Cabe mencionar a continuación la propiedad más importante que cumple el gradiente conjugado para el caso cuadrático, donde se garantiza una convergencia muy rápida. Teorema 2.8. (Propiedad del método de gradiente conjugado lineal) Para una función cuadrática definida positiva (2.23), el método de gradiente conjugado descrito por la fórmula de Fletcher-Reeves (2.31) y (2.32) con búsqueda lineal exacta converge en m≤npasos. Además, si el número de autovalores distintos de G es igual a m, entonces se cumplen las siguientes propiedades ∀i∈ {0, ..., m}: dT iGdj= 0, j = 0,1..., i −1(2.35) 2.5. MÉTODO DE GRADIENTE CONJUGADO NO LINEAL 17 gT igj= 0, j = 0,1..., i −1(2.36) dT igi=−gT igi, j = 0,1..., i −1(2.37) [g0, g1, ..., gi]=[g0, Gg0, ..., Gig0](2.38) [d0, d1, ..., di]=[g0, Gg0, ..., Gig0](2.39) Demostración. La demostración se puede consultar en [SY06]. 2.5. Método de gradiente conjugado no lineal Ahora adaptaremos el método para el caso no lineal siendo Juna función objetivo cualquiera. Para ello realizaremos algunos cambios sobre el algoritmo que describimos previamente para el caso lineal. Método de búsqueda inexacta no lineal de Polak-Ribiére Algoritmo 2.9. GC-PR Paso 1: Dado x0un punto inicial, evaluamos J(x0) = J0y∇J0=∇f(x0). Paso 2: Definimos de la misma manera que en el método lineal la dirección de descenso inicial d0← −∇J0y la etapa k←0. Paso 3: Comprobamos que ∇Jk6= 0. Paso 4: Realizamos la siguiente etapa computando αkyxk+1 =xk+αk·dk. Paso 5: Evaluamos de nuevo ∇Jk+1 y calculamos βk+1 mediante la fórmula de Polak-Ribiére para el caso no lineal. De esta forma escribimos βP R k+1 ←∇JT k+1(∇Jk+1 −∇Jk) ||∇Jk||2(2.40) dk+1 ← −∇Jk+1 +βk+1dk(2.41) 18CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE k←k+ 1 (2.42) Si escogemos Jcomo una función convexa cuadrática y αkcomo el mínimo exacto, el algoritmo previo se trasforma en el algoritmo del caso lineal. Si se quiere ver más detalles sobre el código consultar [Vog02]. Para completar el algoritmo anterior, debemos especificar cómo seleccionamos el paso αk. Para ello debemos recordar que para asegurar un descenso suficiente en el método estos términos αkdeben satisfacer ciertas condiciones. Gracias a la fórmula del producto escalar obtenemos: ∇JT kdk=−k∇Jkk2+βP R k∇JT kdk−1(2.43) Si la búsqueda lineal es exacta, como αk−1es un mínimo local de J en la dirección dk−1tenemos que ∇JT kdk−1=0, por tanto, por (2.43) tenemos que ∇JT kdk<0y entonces dkes una dirección de descenso. En el caso de que la búsqueda sea inexacta el término positivo del lado derecho de la ecuación (2.43) tiene más peso, por tanto ∇JT kdk<0. Para corregir esta situación debemos utilizar las anteriormente explicadas Condiciones fuertes de Wolfe-Powell. J(xk+αkpk)≤J(xk) + c1αk∇JT kdk(2.44) k∇J(xk+αkpk)Tpkk≤−c2∇JT kdk(2.45) Donde 0< c1< c2<1 2. Convergencia del método de Gradiente conjugado no lineal con PolakRibère Como mencionamos anteriormente, este algoritmo tiene una peculiaridad que lo diferencia de FR, el reinicio automático. Esto es, si el algoritmo es demasiado lento y ∇Jk+1 ≈ ∇Jk, entonces la fórmula provoca que βk≈0y por tanto dk+1 ≈ −∇Jk+1. Por tanto, esta propiedad es muy útil para evitar una tendencia lenta en el algoritmo. Se puede observar más concretamente en [Sha78]. 2.5. MÉTODO DE GRADIENTE CONJUGADO NO LINEAL 19 A continuación, enunciaremos los teoremas de convergencia global para el método de Polak-Ribère, al igual que en la propia tesis [Rei80] demostraremos los resultados utilizando búsqueda exacta. Lema 2.10. Dada una función J tal que ∇J(x)es uniformemente continua en nuestro subconjunto S⊆Uy considerando θkcomo el ángulo entre la dirección de descenso dky−∇J(xk). Si se cumple θk≤π 2−µ, con µ > 0,(2.46) entonces se tiene que ∇J(xk) = 0 para un cierto k, o J(xk)→ ∞ o∇J(xk)→0. Demostración. La prueba de este resultado se puede ver más detalladamente en [A+20]. Teorema 2.11. (Teorema del Valor Medio). Dada una función J:Rn→Rse tiene que para cualquier vector d∈Rn, J(x+αd) = J(x) + ∇J(x+αd)Td, α ∈(0,1).(2.47) Ahora estaremos en condiciones óptimas para demostrar el teorema de convergencia global. Teorema 2.12. Sea J:Rn→Runa función dos veces diferenciable en nuestro conjunto subconjunto S⊆Uacotado. Asumiendo que existe una constante positiva m > 0que verifica para todo x∈S, y ∈Rn mkyk2≤yT∇2J(x)y(2.48) Entonces la secuencia de {xk}generada por el método de Polak-Ribiére con búsqueda lineal exacta converge a un único mínimo x∗de la función J. Demostración. Por el lema anterior será suficiente probar (2.46), es decir, que existe una constante positiva ω > 0que verifique la siguiente expresión −gT k+1dk+1 ≥ωkgk+1kkdk+1k(2.49) esto es, cos θk≥ω≥0. Por tanto, con el lema mencionado se observa que gk→0yg(x∗)=0. De (2.48) se sigue que {xk} → x∗, mínimo, el cual es único 20CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE en la función J. A causa de que es búsqueda lineal exacta, se puede deducir de dk+1 =−gk+1 +βkdkygT k+1dk= 0 por tanto gT k+1dk+1 =−kgk+1k2. Entonces la expresión 2.49 será equivalente a kgk+1k kdk+1k≥ω(2.50) Recordamos que el valor de αkse obtiene de αk=−gT kdk dt kAkdk =kgkk2 dt kAkdk (2.51) donde Ak=Z1 0∇2J(xk+tαkdk)dt (2.52) Utilizando el Teorema del Valor Medio, se puede deducir de la fórmula anterior lo siguiente gk+1 −gk=∇J(xk+αkdk)−∇J(xk) = αkAkdk.(2.53) De esta forma la fórmula de βken el caso de PR se puede expresar βP R k=gT k+1(gk+1 −gk) gT kgk =αk gT k+1Adk kgkk2=gT k+1Akdk dT kAkdk (2.54) Por tanto, como el conjunto S es acotado, existe M > 0constante tal que para todo x∈Sy para cualquier y∈Rn yTA(x)y≤Mkyk2.(2.55) De este modo podemos establecer la siguiente cota |βP R k| ≤ kgk+1kkAkdkk mkdkk2≤M mkgk+1k kdkk(2.56) Además, kdk+1k ≤ kgk+1k+|βP R k|kdkk≤kdk+1k+M mkgk+1k= (1 + M m)kgk+1k,(2.57) esto es, 2.5. MÉTODO DE GRADIENTE CONJUGADO NO LINEAL 21 kgk+1k kdk+1k≥1 + M m−1 (2.58) De este modo, solo tendremos que tomar ω=m m+Mcumpliendo así 2.50. La experiencia práctica nos indica que PR tiende a confirmar que es un algoritmo más robusto y eficiente que Fletcher-Reeves u otras fórmulas para resolver problemas de optimización sin restricciones, de ahí su gran interés en la práctica. En lo relativo a la búsqueda inexacta, es importante mencionar que desde el punto de vista teórico este no garantiza que dksea dirección de descenso con las condiciones fuertes de Wolfe-Powell (2.11) y (2.13). Como en la programación del algoritmo utilizaremos PR con búsqueda lineal inexacta mencionaremos la siguiente formulación alternativa de βmencionada en varias de las fuentes consultadas para la realización de este trabajo. βk+1 = m´ax {βP R k+1,0}(2.59) Con esta sugerencia, propuesta ya por Powell en 1984, se desarrolla en [Noc92] para concluir la convergencia global de Polak-Ribière con las condiciones fuertes de Wolfe-Powell, es decir, con búsqueda inexacta. De esta forma se consigue alcanzar las propiedades de convergencia no lineal que cumplen otros métodos menos robustos como Fletcher-Reeves (estas se pueden ver en [SY06]). En el siguiente capítulo desarrollaremos el Método de Elementos Finitos, el cual complementaremos con los métodos de Gradiente Conjugado descritos en esta sección para obtener un algoritmo que resuelva eficientemente los problemas de optimización no lineales sin restricciones. 22CAPÍTULO 2. MÉTODOS DE GRADIENTE CONJUGADO CON POLAK-RIBIÈRE Capítulo 3 Resolución de problemas elípticos 1D mediante el método de elementos finitos 3.1. Introducción En esta sección estudiaremos uno de los métodos de discretizado más usados para resolver problemas de optimización en forma continua, usualmente asociados a campos de la física o la ingeniería, entre otros. El método que estudiaremos se conoce como el método de los elementos finitos (MEF), que desde la aparición de las computadoras ha gozado de un éxito más que notable para el estudio de problemas de optimización. Como mencionamos anteriormente, el MEF permite transformar un problema de optimización continuo en uno discreto, conociendo la forma aproximada de una función a través de un número finito de puntos. Notación 3.1.Durante este tema, para simplificar utilizaremos la notación u0(x)≡ du(x) dx . 23 24CAPÍTULO 3. RESOLUCIÓN DE PROBLEMAS ELÍPTICOS 1D MEDIANTE EL MÉTODO DE ELEMENTOS FINITOS 3.2. Descripción del MEF 1D 3.2.1. Introducción al problema de Sturm-Liouville En este apartado resolveremos con el MEF el problema de Sturm-Liouville de segundo orden. Comenzaremos resolviendo el problema para unas condiciones de contorno Dirichlet, es decir, considerando el siguiente problema (P) Dados α, β ∈Ryf∈L2(a, b) p∈H1(a, b), q ∈L∞(a, b), p(x)≥α > 0, q(x)≥0,∀x∈[a, b] determinar u∈H2(a, b)que verifique (P)(−(p(x)u0(x))0+q(x)u(x) = f(x)en (a, b) u(a) = α, u(b) = β 3.2.2. Formulaciones diferencial y variacional del problema Ahora consideraremos una función v∈D(a, b)que sea escalar, de clase infinito y con supp(v)⊂(a, b)un compacto, entonces v(a)=v(b)=0. Si multiplicamos la igualdad diferencial del problema (P)por vobtendremos la siguiente ecuación: −(p(x)u0(x))0v(x) + q(x)u(x)v(x) = f(x)v(x)(3.1) Realizamos las operaciones e integramos en (a,b) Zb a−(p(x)u0(x))0v(x) + q(x)u(x)v(x)dx =Zb a f(x)v(x)dx en (a, b)(3.2) Utilizando la linealidad de la integral −Zb a (p(x)u0(x))0v(x)dx +Zb a q(x)u(x)v(x)dx =Zb a f(x)v(x)dx en (a, b)(3.3) Ahora, podemos integrar por partes el primer sumando tomando u=v→du = v0dx y dv = (pu0)0→v=pu0 −Zb a (p(x)u0(x))0v(x)dx =−v(x)(p(x)u0(x))]b a−Zb a−p(x)u0(x)v0(x)dx (3.4) 3.2. DESCRIPCIÓN DEL MEF 1D 25 Como v(a)=v(b)=0 se obtiene =Zb a p(x)u0(x)v0(x)dx (3.5) Observación 3.2.Para un vcon las condiciones seleccionadas anteriormente y una función g∈ L∞se cumple la siguiente equivalencia g(x) = 0 en (a, b)⇐⇒ Rb ag(x)v(x)dx = 0 ∀v, entonces podremos formular el problema (P)de dos formas equivalentes. Para el resolver el problema (P) para condiciones tipo Dirichlet enunciaremos las formulaciones diferencial (o clásica)(D) y variacional (V) que enunciaremos a continuación, para posteriormente realizar el proceso de discretización. Definición 3.3. (Formulación diferencial o clásica[Dirichlet]) (D)(Determinar u ∈H1(a, b)tal que u(a) = α, u(b) = β, verificando Rb ap(x)u0(x)v0(x)dx +Rb aq(x)u(x)v(x)dx =Rb af(x)v(x)dx ∀v∈H1 0(a, b) Definición 3.4. (Formulación variacional[Dirichlet]) (V)     −(p(x)u0(x))0+q(x)u(x) = f(x)en (a, b) u(a) = α, u(b) = β Rb a(−(p(x)u0(x))0+q(x)u(x)−f(x))v(x)dx = 0 ∀v tal que v(a) = v(b)=0 3.2.3. Discretización del problema variacional Consideraremos una malla en [a, b]formada por Nelementos Ti= [ai, ai+1]que forman una partición del intervalo N ∪ i=1Ti= [a, b] a=a1< a2< ... < aN< aN+1 =b hi=ai+1 −ai=longitud de Ti,h= m´ax 1≤i≤Nhi i∈ {1, ..., N} Definición 3.5. Se define el espacio vectorial de polinomios de grado ≤kcomo Pk:= h1, x, x2, ..., xki, k ∈N.(dim Pk=k+ 1) 32CAPÍTULO 3. RESOLUCIÓN DE PROBLEMAS ELÍPTICOS 1D MEDIANTE EL MÉTODO DE ELEMENTOS FINITOS PASO 10: Cálculo de las integrales para la resolución del problema: Teniendo en cuenta lo desarrollado anteriormente podemos escribir las matrices elementales de Masa y Rigidez en función de dos integrales, por tanto, será preciso calcularlas para realizar el MEF. A(i) M=ZTi ~w(i)(x)q(x)(~w(i)(x))Tdx, matriz elemental de masa (3.7) A(i) R=ZTi D ~w(i)(x)p(x)(D ~w(i)(x))Tdx, matriz elemental de rigidez (3.8) A continuación, transformamos las expresiones anteriores a partir de las definiciones de las funciones vectoriales del Paso 3 A(i) M=ZTi w(i) j(x)q(x)(w(i) k(x))Tdx, matriz elemental de masa (3.9) A(i) R=ZTi (w(i) j)0(x)p(x)(w(i) k)0(x)dx, matriz elemental de rigidez (3.10) Ahora, para calcular estas matrices estableceremos un elemento referencia ˆ T= [0,1] y una transformación afín ϕi: [0,1] →[ai, ai+1] = Tital que ϕi(ˆx) = ˆ hi·ˆx+ai. Además, las funciones base para este elemento referencia se pueden expresar a través de las funciones definidas previamente Proposición 3.21. ˆwj(ˆx) = w(i) j(ϕi(ˆx)) (1 ≤j≤2,1≤i≤N).(3.11) Demostración. Para probar la igualdad anterior es suficiente ver que w(i) j◦ϕi es un polinomio de grado 1 que vale δjk en ˆxk(1≤k≤2). Como ˆwj=w(i) j◦ϕi(3.12) 3.2. DESCRIPCIÓN DEL MEF 1D 33 En consecuencia ( ˆwj)0(ˆx) = (w(i) j)0(ϕi(ˆx))hi(3.13) ( ˆwj)0(ϕi(ˆx)) = 1 hi ( ˆw(i) j)0(ˆx)(3.14) entonces, se tiene que para la matriz de masa ZTi (w(i) j(x)q(x)w(i) k(x))dx = |{z} x=ϕi(ˆx) hiZ1 0 (w(i) j(ϕi(ˆx))q(ϕi(ˆx))w(i) k(ϕi(ˆx)))dˆx= =hiZ1 0 ( ˆwj(ˆx) |{z} pol. de grado 1 q(ϕi(ˆx)) ˆwk(ˆx) |{z} pol. de grado 1 )dx (3.15) De equivalente forma, realizamos el mismo desarrollo para la matriz de rigidez ZTi ((w(i) j)0(x)p(x) (w(i) k)0(x)) dx = |{z} (3,13)(3,14) =1 hiZ1 0 ( ( ˆwj)0(ˆx) | {z } pol. de grado 0 p(ϕi(ˆx)) ( ˆwk)0(ˆx) | {z } pol. de grado 0 )dx (3.16) Observación 3.22.Si se desease calcular estas integrales de forma exacta será necesario que pyqcumplan condiciones muy restrictivas para los problemas a desarrollar, por ello, aproximaremos estas mediante fórmulas de cuadratura. Como en nuestro caso tratamos de resolver el Método de Elementos Finitos para Lagrange P1, sería suficiente usar fórmulas exactas en P0. En nuestro caso utilizaremos la fórmula de cuadratura de Poncelet para calcular estas integrales. Definición 3.23. (Fórmula del punto medio o de Poncelet.) Zb a f(x)dx ≈(b−a)f(x0), con x0=a+b 2(3.17) 34CAPÍTULO 3. RESOLUCIÓN DE PROBLEMAS ELÍPTICOS 1D MEDIANTE EL MÉTODO DE ELEMENTOS FINITOS Cálculo de A(i)=A(i) R+A(i) M Como hemos desarrollado anteriormente, se deduce que (A(i) R)ij =1 h1Z1 0 ( ˆwi)0(ˆx)p(ϕ1(ˆx))( ˆwj)0(ˆx)dˆx, donde h1=a2−a1, ϕ1(ˆx) = h1ˆx+a1 (3.18) Por tanto, (A(i) R)11 =1 h1Z1 0 ( ˆw1)0(ˆx) | {z } −1 p(ϕ1(ˆx)) ( ˆw1)0(ˆx) | {z } −1 dˆx= =1 h1Z1 0 p(ϕ1(ˆx))dˆx≈ |{z} (3,17) 1 h1 p(ϕ1(1 2)) = 1 h1 p(x(1) m) (3.19) Del mismo modo se puede desarrollar para el resto de puntos de la matriz elemental A(i) R: (A(i) R)12 =1 h1Z1 0 ( ˆw1)0(ˆx) | {z } −1 p(ϕ1(ˆx)) ( ˆw2)0(ˆx) | {z } +1 dˆx= =−1 h1Z1 0 p(ϕ1(ˆx))dˆx≈ |{z} (3,17) −1 h1 p(ϕ1(1 2)) = −1 h1 p(x(1) m) (3.20) (A(i) R)22 =1 h1Z1 0 ( ˆw2)0(ˆx) | {z } +1 p(ϕ1(ˆx)) ( ˆw2)0(ˆx) | {z } +1 dˆx= =1 h1Z1 0 p(ϕ1(ˆx))dˆx≈ |{z} (3,17) 1 h1 p(ϕ1(1 2)) = 1 h1 p(x(1) m) (3.21) En definitiva, tenemos la matriz elemental resultante A(i) R≈p(x(i) m) hi 1−1 −1 1 !(3.22) Notación 3.24.Escribimos el punto medio del intervalo Ticomo x(i) m=x(i) 1+x(i) 2 2. 3.2. DESCRIPCIÓN DEL MEF 1D 35 Hallaremos la matriz de masa elemental A(1) Mbajo el mismo método realizado para la matriz elemental de rigidez. (A(i) M)ij =h1Z1 0 ˆwi(ˆx)q(ϕ1(ˆx)) ˆwj(ˆx)dˆx≈ |{z} (3,17) h1ˆwi(1 2)q(ϕi(1 2)) ˆwj(1 2) = = |{z} ˆwi(1 2)=1 2,ˆwj(1 2)=1 2 h1 4q(x(1) m) (3.23) En definitiva, realizando los cálculos obtenemos la matriz elemental resultante A(1) M≈h1q(x(1) m) 4 1 1 1 1!(3.24) Análogamente, para cualquier A(i)se tiene que A(i)=A(i) M+A(i) R, con 2≤i≤N(3.25) A(i) R≈p(x(i) m) hi 1−1 −1 1 !(3.26) A(i) M≈hiq(x(i) m) 4 1 1 1 1!(3.27) Cálculo del vector ~ b(i)para 1≤i≤N Sabemos que ~ b(i)=ZTi ~w(i)(x)f(x)dx (3.28) Por tanto, dado j= 1,2podemos desarrollar por componentes del vector elemental al igual que en la matriz elementar A(i) b(i) j=ZTi w(i) j(x)f(x)dx = |{z} x=ϕi(ˆx) hiZ1 0 ˆwj(ˆx)f(ϕ(ˆx))dˆx≈ ≈ |{z} P oncelet hiˆwj |{z} 1 2 (1 2)f(x(i) m) = hif(x(i) m) 2 (3.29) 36CAPÍTULO 3. RESOLUCIÓN DE PROBLEMAS ELÍPTICOS 1D MEDIANTE EL MÉTODO DE ELEMENTOS FINITOS En consecuencia, ~ b(i)≈hif(x(i) m) 2 1 1!(3.30) 3.2.5. Implementación de las condiciones de contorno Resolvemos el problema enunciado mediante la formulación variacional (−(p(x)u0(x))0+q(x)u(x) = f(x)en (a, b) L1u(a) + L2u0(a) = α, L3u(b) + L4u0(b) = β Para abordar la implementación de las condiciones de contorno subdividiremos por casos los valores que podrán tomar los Lii= 1,2,3,4. Caso L1, L36= 0 yL2, L4= 0. (DIRICHLET-DIRICHLET) En este caso, el problema a resolver vendrá dado por (Hallar u ∈H1(a, b)tal que u(a) = α L1, u(b) = β L3 Rb ap(x)u0(x)v0(x)dx +Rb aq(x)u(x)v(x)dx =Rb af(x)v(x),∀v∈H1 0(a, b) Mediante el ensamblado se obtuvo A=                 a1,1a1,20. . . . 0 a2,1a2,2a2,30. . . 0 0a3,2a3,3a3,4. . . 0 . . . . . . . . . . . . . . . . . . . . . . . . 0. . . . aN,N−1aN,N aN,N+1 0. . . . 0aN+1,N aN+1,N+1                 ~ b=                 b1 b2 b3 . . . bN bN+1                 Por consiguiente los cambios que habrá que realizar sobre Ay~ bson a1,1= 1, a1,2= 0, aN+1,N = 0, aN+1,N+1 = 1 b1=α L1, bN+1 =β L3 3.2. DESCRIPCIÓN DEL MEF 1D 37 Con esto, obtenemos el sistema matricial resultante a través del método de elementos finitos con las condiciones de contorno indicadas, en nuestro caso serán Dirichlet. Para esta resolución matricial será fundamental el uso del método de gradiente conjugado explicado en el capitulo anterior. A continuación, detallaremos una serie de ejemplos que resolveremos con la utilización de Matlab a través de los métodos descritos en este trabajo. 38CAPÍTULO 3. RESOLUCIÓN DE PROBLEMAS ELÍPTICOS 1D MEDIANTE EL MÉTODO DE ELEMENTOS FINITOS Capítulo 4 Formulación del MEF para el caso no lineal 4.1. Introducción En este capítulo veremos como resolveremos algunos de los ejercicios propuestos en [Rei80], por tanto tendremos que adaptar el método de elementos finitos para el problema no lineal. En los ejercicios propuestos el coeficiente no lineal se encuentra en el vector ~ f, es decir, afectará a nuestra formulación del vector del segundo miembro en nuestro sistema. Finalmente, resolveremos el sistema utilizando el método de gradiente conjugado con la variante de Polak-Ribiére. Uno de los problemas no lineales estudiados en [Rei80] es problema de Bratu, que enunciaremos a continuación para el caso 1D. 4.2. El problema de Bratu 1D El problema de Bratu en dimensión uno se puede formular de la siguiente manera: (−u00(x) = λeu(x),∀x∈[0,1] y λ ∈R+ u(0) = 0, u(1) = 0 En consecuencia, el cambio se efectuará sobre f=λeu(x). Seguidamente, se desarrollará como modificar el método de elementos finitos para adaptarlo a un problema con término no lineal en f. 39 40 CAPÍTULO 4. FORMULACIÓN DEL MEF PARA EL CASO NO LINEAL 4.2.1. Ajuste del modelo no lineal A continuación describiremos como ajustar nuestro problema no lineal a la programación usual del MEF lineal. Por simplicidad tomaremos λ= 1, en el caso de tomar otro valor el procedimiento sería idéntico. En virtud de las condiciones de contorno, la formulación variacional del problema sería de la siguiente forma: Z1 0 u0(x)v0(x)dx =Z1 0 f(u(x))v(x)dx, ∀v∈V(4.1) siendo f(u(x)) = eu(x). Vemos que el término no lineal solo influye en el lado derecho de la ecuación, es decir, en el vector segundo miembro ~ b. Por ello idea será integrar ese mismo lado de la ecuación con una regla de integración que solo evalúe la función en los nodos, dejando exactamente igual los términos de la matriz A. Seguimos por tanto el mismo procedimiento que en el capítulo 3 del trabajo. Z1 0 f(uh(x))vh(x)dx = N+1 X i=1 ZTi f(uh(x))|Tivh(x)|Tidx = = N+1 X i=1 ZTi f(uh(x)|Ti)vh(x)|Tidx = (4.2) Sustituyendo ahora vh(x)|Ti= (~w(i)(x))T~v(i)yuh(x)|Ti= (~w(i)(x))T~u(i)se obtiene = N+1 X i=1 ZTi f((~w(i)(x))T~u(i))(~w(i)(x))T~v(i)dx = N+1 X i=1 (~v(i))TZTi (~w(i)(x))f((~w(i)(x))T~u(i))dx | {z } ~ b(u)(i) (4.3) Esta última integral la resolveremos aproximando con la ya utilizada Regla de Poncelet (o punto medio), ZTi (~w(i)(x))f((~w(i)(x))T~u(i))dx ≈hi[~w(i)(x(i) m))f((~w(i)(x(i) m))T~u(i)] = (4.4) 4.2. EL PROBLEMA DE BRATU 1D 41 Como ya vimos ~w(i)(x(i) m) = (0,5 0,5)Ty~u(i)= (uh(x(i) 1)uh(x(i) 2))Tcon lo cual se obtiene = hif (0,5 0,5) uh(x(i) 1) uh(x(i) 2))!! 2 1 1!(4.5) De esta forma modificamos el vector aproximando la función u por los valores en los nodos correspondientes. Esta técnica se puede consultar en más profundidad en la referencia [TM21]. 4.2.2. Solución analítica del problema de Bratu Para implementar el algoritmo será importante hallar las soluciones analíticas al problema de Bratu 1D. Para ello haremos un breve resumen de las posibilidades, estas se pueden consultar fácilmente en [Moh14]. Se define por tanto λc, el punto crítico de la ecuación diferencial a partir del cual varía la cantidad de soluciones del problema, que serán las siguientes: Si λ=λcentonces existirá una única solución exacta, que viene dada por la ecuación u(x) = −2 ln cosh ( θ 2(x−1 2)) cosh θ 4donde θ=√2λcosh θ 4. Si λ > λcentonces el problema no tendrá soluciones exactas. Si finalmente λ<λcel problema contará con dos soluciones. Estas se calculan a partir de la misma ecuación que en el primer caso. En nuestro caso, analizaremos en la práctica las soluciones del problema de solución única. Realizando los cálculos obtendremos los valores λc= 3,513830719 y θc= 4,79871456. Cabe recalcar que el enunciado del problema sería mucho más general que el presentado, pues se ha modificado las condiciones de λyxpara una mejor implementación del código. 48 CAPÍTULO 5. RESULTADOS OBTENIDOS Capítulo 6 Conclusiones y trabajo futuro Durante este trabajo se han revisado las técnicas de gradiente conjugado para resolver sistemas de ecuaciones tanto para el caso lineal y no lineal. Para ello se ha optado por la variante de Polak-Ribiére, recomendada en la literatura del artículo principal del trabajo y en el resto de referencias consultadas y ya mencionadas. Además, se han utilizado técnicas de búsqueda inexacta como la Regla de Wolfe-Powell o la Regla de Goldstein junto a la interpolación cuadrática. Esto ha sido de gran utilidad para profundizar y complementar el contenido impartido en las asignaturas de métodos numéricos que se han impartido en el grado. Además, se han mostrado algunos de los resultados del programa realizado con resultados satisfactorios. También se ha revisado la implementación del método de elementos finitos para la discretización de un problema 1D, así como modificado para la implementación del problema no lineal (Bratu) sugerido en el artículo [Rei80]. Para ello se han utilizado algunas de las técnicas sugeridas en la asignatura de Análisis Numérico de Ecuaciones en Derivadas Parciales del cuarto curso del grado. Aunque los resultados de la programación no han sido todo lo satisfactorio posible, esto nos puede servir para en el futuro realizar una programación mucho más eficiente y óptima del algoritmo basándonos en el estudio teórico mostrado durante el trabajo. Como complementario al trabajo, cabe resaltar que existen otras técnicas de discretización para los problemas vistos, como el método de diferencias finitas, en las que se puede profundizar para una comparación con los resultados obtenidos mediante el método de elementos finitos. Además, también existen otros métodos iterativos para la resolución del sistema matricial como por ejemplo el Método de 49 50 CAPÍTULO 6. CONCLUSIONES Y TRABAJO FUTURO Newton , revisado en varias de las asignaturas del grado. Estos métodos también son mencionados y estudiados en mucha de la bibliografía presentada en este trabajo. Finalmente, cabe destacar que se han resuelto problemas para el caso unidimensional, pero para muchos de las ecuaciones diferenciales utilizadas, como por ejemplo en el caso del problema de Bratu, existen enunciados alternativos en dimensión superior expuestos en [Rei80]. (∆u(x) = λeu(x)∀x∈Ω Con u(x) = 0 ∀x∈∂Ω donde ∆udenota el laplaciano de u. Por tanto, una ampliación al estudio realizado durante el trabajo sería desarrollar los métodos vistos para un caso de dimensión superior. Cabe mencionar que en el artículo se utilizan métodos de continuación para resolver los problemas en dimensión superior. Capítulo 7 ANEXOS 7.1. ANEXO I: Convergencia y existencia de los métodos de búsqueda inexacta A continuación comentaremos los teoremas más relevantes para la convergencia de los métodos de búsqueda inexacta con Wolfe-Powell o Regla de Goldstein. Lema 7.1. Supongamos J:Rn→Rcontinuamente diferenciable. Sea dkuna dirección de descenso en xky supongamos que f es acotada a lo largo del conjunto {xk+αdk:α > 0}. Entonces si 0< c1< c2<1, existen intervalos de longitudes de paso que cumplen las condiciones de Wolfe-Powell y las condiciones de Wolfe-Powell fuertes. Demostración. Como φ(α) = J(xk+αdk)es una función acotada para todo α > 0 con 0< c1<1, entonces si denotamos el lado derecho de la condición (2.11) por l(α) = J(xk) + αc1∇JT kdkeste debe intersecar al menos una vez con la gráfica de φ. Denotemos pues ¯α > 0el valor más pequeño de αpara el cual sucede dicha intersección. J(xk+ ¯αdk) = J(xk) + ¯αc1∇JT kdk(7.1) La condición de decrecimiento suficiente claramente es válida para longitudes de paso inferiores a ¯α. Utilizando el Teorema de Valor Medio, tenemos que existe ¯ ¯α∈(0,¯α)tal que 51 52 CAPÍTULO 7. ANEXOS J(xk+¯ ¯αdk)−J(xk) = ¯α∇J(xk+¯ ¯αdk)Tdk(7.2) Combinando las dos igualdades anteriores obtenemos la siguiente expresión: ∇J(xk+¯ ¯αdk)Tdk=c1∇JT kdk> c2∇JT kdk(7.3) ya que c1< c2y∇JT kdk<0. Por tanto, ¯ ¯αsatisface las condiciones de WolfePowell. Tomando nuestras hipótesis sobre f, sabemos que existe un entorno de ¯ ¯α para el cual se cumplen dichas condiciones de Wolfe-Powell. Además, dado que en la última igualdad (7.3) el lado izquierdo es negativo, las condiciones fuertes de Wolfe-Powell también se cumplen en el mismo intervalo. Además, las condiciones de Wolfe-Powell son escalarmente invariantes, es decir, si multiplicamos la función objetiva por una constante o realizamos una transformación afín no alteramos las condiciones. (Para más detalles mirar [NW06]). A continuación veremos la convergencia para los métodos de Wolfe-Powell y Goldstein. Convergencia del Algoritmo utilizado Regla de Goldstein o Wolfe-Powell Para demostrar el descenso de los métodos será fundamental evitar los casos en los que las direcciones de búsqueda (que definiremos sk=αkdk) estén próximas a la ortogonalidad con −gk, es decir, el ángulo θkentre sky−gkestá uniformemente delimitado θk≤π 2−µ, ∀k(7.4) donde µ >0, θ∈[0,π 2]está definido por cos(θk) = −gT ksk (kgkkkskk)(7.5) garantizando el descenso. Teorema 7.2. Dado un αkpor la Regla de Goldstein o Wolfe-Powell y un skque satisface las condiciones definidas anteriormente. Si ∇Jexiste y es uniformemente continuo en el conjunto {x:J(x)≤J(x0)}entonces ∇J(xk) = para algún k, o J(xk)→ −∞, o ∇J(xk)→0. Demostración. Ver [SY06]. 7.1. ANEXO I: CONVERGENCIA Y EXISTENCIA DE LOS MÉTODOS DE BÚSQUEDA INEXACTA53 Convergencia del algortimo de Interpolación cuadrática A continuación veremos el teorema que asegura la convergencia para el algritmo con Interpolaión cuadrática, para esto necesitaremos definir α∗como el punto que verifica φ0(α∗)=0yφ00(α∗)6= 0. Teorema 7.3. Dada J:R−→ Runa función continuamente diferenciable de orden 4 y el punto α∗con las condiciones anteriores. Entonces la sucesión {αk}generada por (2.18) y (2.19) es convergente con orden 1.32 Demostración. Definimos la fórmula de interpolación de Lagrange L(α)=φ1(α−α2)(α−α3) (α1−α2)(α1−α3) +φ2(α−α1)(α−α3) (α2−α1)(α2−α3)+φ3(α−α1)(α−α2) (α3−α1)(α3−α2),de esta forma, será equivalente calcular L0(α) = 0con el desarrollo de las fórmulas anteriores (2.18) y (2.18). Establecemos a partir de esto la ecuación φ(α) = L(α) + R(α)(7.6) Donde R(α) = 1 6φ00(ξ(α))(α−α1)(α−α2)(α−α3)(7.7) A partir de las hipótesis del teorema 0 = φ0(α∗) = L0(α∗)+R0(α∗), sustituyendo se obtiene: φ1 2α∗−(α2+α3) (α1−α2)(α1−α3)+φ2 2α∗−(α3+α1) (α2−α1)(α2−α3)+φ3 2α∗−(α1+α2) (α3−α1)(α3−α2)+R0(α∗) = 0 (7.8) Reescribimos (2.18) dividiendo por (α1−α2)(α2−α3)(α3−α1)en el denominador y numerador obteniendo el siguiente término de la sucesión: α4=1 2 φ1(α2+α3) (α1−α2)(α1−α3)+φ2(α3+α1) (α2−α1)(α2−α3)+φ3(α1+α2) (α3−α1)(α3−α2) J1 (α1−α2)(α1−α3)+φ2 (α2−α1)(α2−α3)+φ3 (α3−α1)(α3−α2) (7.9) Tomando (7.9) y (7.8) podemos ajustarlo tal que : α∗−α4=1 2 R0 3(α∗) φ1 (α1−α2)(α1−α3)+φ2 (α2−α1)(α2−α3)+φ3 (α3−α1)(α3−α2) (7.10) A continuación, establecemos ei=α∗−αi, i = 1,2,3,4.De esta forma, reescribimos la ecuación anterior por: 54 CAPÍTULO 7. ANEXOS e4[−φ1(e2−e3)−φ2(e3−e1)−φ3(e1−e2)] = −1 2R0(α∗)(e1−e2)(e2−e3)(e3−e1).(7.11) Sabiendo que φ0(α∗)=0se sigue desde la fórmula de expansión del polinomio de Taylor que φi=φ(α∗) + 1 2e2 iφ00(α∗) + O(e3 i)(7.12) Tomando las dos fórmulas anteriores (7.11 sin el término de tercer orden) y (7.12) llegamos a la expresión siguiente: e4=1 φ00(α∗R0(α∗).(7.13) Además, por la fórmula de Interpolación de Lagrange (7.7) tenemos la siguiente igualdad: R0(α∗) = 1 6J000(ξ(α∗))(e1e2+e2e3+e3e1) + 1 24J(4)(η)e1e2e3.(7.14) Ahora, eliminando el término de cuarto orden y sustituyendo (7.13) en (7.14) obtenemos la siguiente expresión e4=J000(ξ(α∗)) 6J00(α∗)(e1e2+e2e3+e3e1) = M(e1e2+e2e3+e3e1),(7.15) Donde M es una constante. Si generalizamos, tenemos la fórmula ek+2 =M(ek−1ek+ekek+1 +ek+1ek−1)(7.16) Ya que ek+1 =O(ek) = O(ek−1)cuando ek→0, entonces existe un ¯ M > 0 verificando que |ek+2| ≤ ¯ M|ek−1||ek|(7.17) como ¯ M > 0esto es lo mismo que escribir ¯ M|ek+2| ≤ ¯ M|ek−1|¯ M|ek|(7.18) Donde |ei|(i=1,2,3) es suficientemente pequeño para verificar la condición 7.1. ANEXO I: CONVERGENCIA Y EXISTENCIA DE LOS MÉTODOS DE BÚSQUEDA INEXACTA55 δ= m´ax{¯ M|e1|,¯ M|e2|,¯ M|e1|} <1,(7.19) Por un lado ¯ M|e1| ≤ ¯ M|e1|¯ M|e2| ≤ δ2.(7.20) Ajustamos por tanto ¯ M|ek| ≤ δqk(7.21) De esta forma ¯ M|ek+2| ≤ ¯ M|ek|¯ M|ek−1| ≤ δqkδqk−1∆ =δqk+2 (7.22) Por eso qk+2 =qk+qk−1,(k≥2) (7.23) donde q1=q2=q3, su ecuación característica será por tanto t3−t−1 = 0. Con la raíz real t1≈1,32 y otras dos raíces complejas |t2|=|t3|<1. La solución general de la ecuación (7.23) tiene la forma qk=A·t1 1+B·tk 2+C·tk 3,(7.24) donde A,B y C son coeficientes a determinar. Claramente cuando k→ ∞, qk+1 −t1qk=BtK 2(t2−t1) + CtK 3(t3−t1)→0 En este caso, cuando k es suficientemente grande qk+1 −t1qk≥ −0,1. Utilizando (7.21) tenemos |ek| ≤ (1 ¯ Mδqk∆ =Bk,(k≥1). Por tanto utilizando la condición para k suficientemente grande se obtiene Bk+1 Bk = δqk+1 ¯ M δt1qk (¯ M)t 1 =¯ Mt1−1δqk+1−t1qk≤δ−0,1¯ Mt1−1(7.25) Lo cual indica la convergencia con orden t1≈1,32. 56 CAPÍTULO 7. ANEXOS 7.2. Anexo II: Existencia y unicidad de solución MEF Teorema 7.4. (Lax-Milgram). Sea Hun espacio de Hilbert y sea A:H×H→R una forma bilineal, continua y coercitiva. Dada L:H→Rlineal y continua sobre H. Entonces existe un único u0∈Htal que A(u0, v) = L(v)∀v∈H A través del Teorema de Lax-Milgram garantizamos la existencia de una solución del problema discreto y variacional.(Para ver una prueba y desarrollo del teorema de Lax-Milgram se puede consultar [CM16],([LB10])) Además, podemos definir las hipótesis sobre las variables del problema de SturmLioville para la existencia y unicidad de la solución: p∈L∞(a, b)donde ∃c∈R:p(x)≥c > 0casi por doquier en (a, b) q∈L∞(a, b)donde p(x)≥0casi por doquier en (a, b) f∈L2(a, b) 7.3. Anexo III: Normas del error Desarrollaremos a continuación las normas utilizadas para la aproximación de soluciones de los métodos utilizados en este trabajo. Para ello hemos utilizado las normas infinito, seminorma H1y norma L2. 7.3.1. Norma infinito Dado un vector ~u = (x1, x2, x3, ..., xn)∈Rndefinimos la norma infinito como k~uk∞=max(|x1|,|x2|, ..., |xn|) = max i∈{1,...,n}|xi|(7.26) 7.3.2. Seminorma H1 Para cualquier solución exacta uy su aproximación uhen el dominio [a, b]se describe la norma de la energía (o semi-norma H1) como: ku−uhkH0 1= Zb a N+1 X i=1  du dx −duh dx  2 dx!1 2 (7.27) 7.4. ANEXO IV: MÉTODOS DE CUADRATURA 57 7.3.3. Norma L2 Para cualquier solución exacta uy su aproximación uhen el dominio [a, b]se describe la norma L2como: ku−uhkL2=Zb a|u−uh|2dx1 2 (7.28) A continuación de estas tres mediciones del error, producido a causa de la aproximación del método empleado, enunciaremos un teorema de vital importancia en referencia al error producido por el método de elementos finitos aplicado a las ecuaciones diferenciales elípticas. El siguiente Lema nos permitirá garantizar que uhes del orden de la mejor aproximación en el espacio que hemos diseñado, es decir, con la norma elegida es el elemento de Vhque menos dista de la solución exacta u. Lema 7.5. (Cea, cota del error). El lema de cea nos asegura que existe una constante C > 0, independiente de htal que ku(x)−uh(x)k ≤ C´ınfvh∈Vh0ku(x)−vh(x)k=Cd(u(x),Vh) Se puede consultar su demostración y desarrollo en el trabajo original de J.Cea ([Céa64]). 7.4. ANEXO IV: Métodos de cuadratura Durante este trabajo, sobre todo en la construcción del método de elementos finitos, se han resuelto diversas integrales mediante aproximación con métodos de cuadratura, fundamentales en la resolución de ecuaciones diferenciales. 7.4.1. Cuadratura de Gauss con 2 puntos La técnica de cuadratura Gaussiana será necesaria para el cálculo de las normas enunciadas en el apartado anterior, manteniendo así el orden buscado. Está regla será una buena aproximación de una integral de la función f(x)en el intervalo [a, b] mediante: Zb a f(x)dx ≈ n X i=1 wif(xi)(7.29) 64 CAPÍTULO 7. ANEXOS % % h es e l ve ct or que contiene l a s l o n gi tu de s de todos l o s elementos . % % vnp es un vecto r que contiene e l numero de puntos en cada s ubi nte rva lo % a bierto (a , vmesh (1)) (ONLY IF a < vmesh (1)) , (vmesh ( k ) , vmesh ( k+1)) % para todo k = 1 , . . . , length ( vmesh)−1, and ( vmesh ( end ) , b ) (ONLY IF % vmesh ( end ) < b ) . % len = length(vmesh ) ; %Numero de elementos que hay en e l vector de a c u m u l a c i n vnp = [ ] ; i f ( len==0) %Malla uniforme h ( 1 :N)=(b−a)/N; %Son N+1 elementos %Vector x con l o s nodos x=linspace(a , b ,N+1); else %caso no t r i v i a l %Construiremos una malla intermedia x , lo mas uniforme p o s i b l e %forzaremos a l o s puntos de vmesh ser nodos %la malla f i n a l debe tener N+1 nodos %Test de longitud de vmesh ( para que e ste contenida en x ) i f ( len>N+1) fprintf( of i , ’ \n\n∗␣La␣ l og it ud ␣de␣vmesh␣ supera ␣a␣ l o s ␣nodos : ␣ERROR’ ) ; error ( ’La␣ longitud ␣de␣vmesh␣ es ␣mayor␣que␣ e l ␣numero␣de␣puntos ’ ) end i f (vmesh(1)==a && vmesh( len)==b) %caso 1 NPL=N+1−len ; %Numero de puntos que quedan %Quiero meter en cada trozo , un numero de puntos %proporcional a la longitud de cada trozo vl=vmesh ( 2 : len)−vmesh ( 1 : len −1); %ve ct or de l o n gi tud es d i v i d i d a s por vmesh vnp=round(NPL.∗vl . / ( b−a ) ) ; % v ector puntos que debo meter en cada i n t e r v a l o 7.6. PROGRAMAS UTILIZADOS 65 %Realizamos una c o f r e c c i n para ver s i metimos todos lo s puntos %necesarios ( por u t i l i z a r round ) . suma=sum( vnp ) ; i f (suma<NPL) nee=NPL−suma ; % n de puntos que me f a l t a n ind=randperm( len −1); ind=ind ( 1 : nee ) ; vnp( ind)=vnp( ind )+1; elseif(suma>NPL) notnee=suma−NPL; % n de puntos que me sobran ind=randperm( len −1); ind=ind ( 1 : notnee ) ; vnp( ind)=vnp( ind ) −1; end %Construimos vector malla x x = [ ] ; for k=1: len −1 xk=linspace(vmesh (k ) , vmesh (k+1) ,vnp (k)+2); x=[x xk ] ; end x=unique (x ) ; %mismo v ecto r pero sin r e p e t i c i o n e s ( c o r r e c c i n de l i ns p a c e ) elseif(vmesh(1)>a && vmesh ( len)==b) NPL=N−len ; vl =[vmesh(1)−a , vmesh ( 2 : len )−vmesh ( 1 : len −1)]; vnp=round(NPL.∗vl . / ( b−a ) ) ; suma=sum( vnp ) ; i f (suma<NPL) nee=NPL−suma ; ind=randperm( len ) ; ind=ind ( 1 : nee ) ; vnp( ind)=vnp( ind )+1; 66 CAPÍTULO 7. ANEXOS elseif(suma>NPL) notnee=suma−NPL; ind=randperm( len ) ; ind=ind ( 1 : notnee ) ; vnp( ind)=vnp ( ind ) −1; end %Construimos vector x x=linspace(a , vmesh ( 1) , vnp (1)+2); for k=1: len −1 xk=linspace(vmesh (k ) , vmesh (k+1) ,vnp (k+1)+2); x=[x xk ] ; end x=unique (x ) ; elseif(vmesh(1)==a && vmesh( len)<b) NPL=N−len ; vl =[vmesh ( 2 : len)−vmesh ( 1 : len −1),b−vmesh( len ) ] ; vnp=round(NPL.∗vl . / ( b−a ) ) ; suma=sum( vnp ) ; i f (suma<NPL) nee=NPL−suma ; ind=randperm( len ) ; ind=ind ( 1 : nee ) ; vnp( ind)=vnp ( ind )+1; elseif(suma>NPL) notnee=suma−NPL; ind=randperm( len ) ; ind=ind ( 1 : notnee ) ; vnp( ind)=vnp ( ind ) −1; end %Construimos vector x x = [ ] ; 7.6. PROGRAMAS UTILIZADOS 67 for k=1: len −1 xk=linspace(vmesh (k ) , vmesh (k+1) ,vnp (k)+2); x=[x xk ] ; end xk=linspace(vmesh ( len ) ,b , vnp( len )+2); x=[x xk ] ; x=unique (x ) ; elseif(vmesh(1)>a && vmesh ( len)<b) NPL=N−len −1; vl =[vmesh(1)−a , vmesh ( 2 : len )−vmesh ( 1 : len −1) ,b−vmesh( len ) ] ; vnp=round(NPL.∗vl . / ( b−a ) ) ; suma=sum( vnp ) ; i f (suma<NPL) nee=NPL−suma ; ind=randperm( len +1); ind=ind ( 1 : nee ) ; vnp( ind)=vnp( ind )+1; elseif(suma>NPL) notnee=suma−NPL; ind=randperm( len +1); ind=ind ( 1 : notnee ) ; vnp( ind)=vnp( ind ) −1; end %Construimos vector x x=linspace(a , vmesh ( 1) , vnp (1)+2); for k=1: len −1 xk=linspace(vmesh (k ) , vmesh (k+1) ,vnp (k+1)+2); x=[x xk ] ; end xk=linspace(vmesh ( len ) ,b , vnp( len +1)+2); x=[x xk ] ; x=unique (x ) ; end 68 CAPÍTULO 7. ANEXOS % %Ahora esta creada l a malla intermedia % f p r i n t f ( ofi , ’\ n\n∗Malla intermedia \n ’ ) ; % f p r i n t f ( ofi , ’\ n %4.4 f ’ , x ) ; % % %Co e ficient e de uniformidad de la malla intermedia % h=x (2: length ( x))−x (1: l en g th ( x ) −1); % hmin=min(h ) ; hmax=max(h ) ; % f p r i n t f ( ofi , ’\ n\n∗Coefici ente de uniformidad de la malla intermedia ’ ) ; % f p r i n t f ( ofi , ’\ n\n∗( Cuanto mas proxima a 1 , mas uniforme ) = %−.6f . ’ , hmax/hmin ) ; %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− %Procederemos a crear l a malla f i n a l x f i n a l = [ ] ; i f ( vmesh(1)>a ) p1=a ; p2=vmesh ( 1 ) ; idx1=find (x>=p1 & x<=p2 ) ; z=affineA (x( idx1 ) , p1 , p2 ) ; %Nodos de [ p1 , p2 ] a [0 ,1 ] zz=fmesh1 ( z ) ; x f i n a l =[ x f i n a l affineAInv ( zz , p1 , p2 ) ] ; %Nodos de [ 0 , 1] a [ p1 , p2 ] end for k=1: len −1 p1=vmesh(k ) ; p2=vmesh(k+1); mk=(p1+p2 ) . / 2 ; idx1=find (x>=p1 & x<=mk) ; z=affineA (x( idx1 ) , p1 ,mk) ; zz=fmesh0 ( z ) ; x f i n a l =[ x f i n a l affineAInv ( zz , p1 ,mk) ] ; idx1=find (x>=mk & x<=p2 ) ; z=affineA (x( idx1 ) ,mk, p2 ) ; zz=fmesh1 ( z ) ; 7.6. PROGRAMAS UTILIZADOS 69 x f i n a l =[ x f i n a l affineAInv ( zz ,mk, p2 ) ] ; end i f (vmesh ( len)<b) p1=vmesh( len ) ; p2=b ; idx1=find (x>=p1 & x<=p2 ) ; z=affineA (x( idx1 ) , p1 , p2 ) ; zz=fmesh0 ( z ) ; x f i n a l =[ x f i n a l affineAInv ( zz , p1 , p2 ) ] ; end x=unique ( x f i n a l ) ; h=x ( 2 : length(x))−x ( 1 : length(x) −1); % f p r i n t f ( ofi , ’\ n\n∗Malla f i n a l x = ’ ); % f p r i n t f ( ofi , ’\ n %4.4 f ’ , x ) ; end fprintf( of i , ’ \n\n∗␣CREACION␣DE␣LA␣MALLA␣COMPLETADA␣\n ’ ) ; end function [ uh]=SolMatlab ( res ) global A bb %Resolvemos e l sistema con Matlab %se p o d r a hacer con matriz no s i m t r i c a , pero con bloqueo de 10^20 Matlab %detecta Matriz mal condicionada , pero s i u tilizamos un bloqueo de 10^10 no %lo detecta i f res ==0 uh = A\bb ; elseif res==1 %Cholesky ( la matriz debe ser definida p o s i t i v a ) L = chol(A) ; 70 CAPÍTULO 7. ANEXOS uh = L’\ bb ; uh = L\uh ; else fprintf( of i , ’ \n\n∗␣ p a r m e t r o ␣de␣ r e s o l u c i n ␣ matriz ␣ sin ␣cambios␣\n ’ ) ; end end function [ alpha , i0 , index ] = WolfeBiseccion ( f , g , x , d , . . . rho , sigma , . . . alphainit , gama , . . . f i d ) global A bb N R Rb % Wolfe Seleccion de paso para un metodo de descenso mediante la r egla % inexacta llamada r e gla de Wolfe . La generacion de un nuevo punto % alpha entre alpha_p y alpha_g se r e a l i z a por biseccion . % % Parametros de entrada : % f : Funcion coste ( psoblemente multi−dimensional ) . % g : Gradiente de la funcion coste . % x : Punto en e l que se encuentra e l algoritmo de descenso . % d : Direccion de descenso . % rho , sigma : Parametros que definen lo s c r i t e r i o s de la r e g l a de % Wolfe . Se supone que 0 < rho < sigma < 1. % Es a cons ejable que 0 < rho < 1/2. % a l p h a i n i t : Valor de alpha que se considera i nici alm ent e . % Deberia s a t i s f a c e r e l c r i t e r i o de paso grande . % gama : Parametro empleado para se l e c c i onar un paso i n i c i a l que % v e r i f i q u e e l c r i t e r i o de paso grande en caso de que % alphag no l o cumpla . 7.6. PROGRAMAS UTILIZADOS 71 % Se t r ata de un f a ct or de d i l a t a c i o n . En consecuencia debe % tomarse t a l que gama > 1. % f i d : I d e n t i f i c a d o r de un f i c h er o para la s a l i d a . % % Parametros de s a l i d a : % alpha : Paso se leccio nado por la reg la de Wolfe para e l metodo de % descenso. % i0 : Numero de evaluaciones de l funcinal coste . % index : Indicador de convergencia : 1 => exito , <= => fracaso . % Si index = 0 => No se ha podido i n i c i a l i z a r . % Si index = −1 => No se ha podido encontrar un alpha % admisible. % I n i c i a l i z a c i o n del algoritmo . % Seleccionamos un alpha que cumpla el c r i t e r i o de paso grande y otro que % cumpla e l c r i t e r i o de paso pequeno . index = 0; i0 = 0; % Numero de evaluaciones de l funcional coste . itmax1 = 10; % Numero maximo de i t e r a c i o n e s en l a i n i c i a l i z a c i o n . itmax2 = 50; % Numero maximo de i t e r a c i o n e s en l a busqueda . init = false ; alphap = 0; alpha = a l phai n i t ; j0 = f (x ) ; jp0 = g (x ) . ’ ∗d ; while ~ i n i t && ( i0 <= itmax1 ) xalpha = x + alpha ∗d ; jalpha = f ( xalpha ) ; jpalpha = g ( xalpha ) . ’ ∗d ; 72 CAPÍTULO 7. ANEXOS CPG = ( jalpha > j0 + rho ∗jp0 ∗alpha ) ; % Criterio de paso grande . CPP = ( jpalpha < sigma ∗jp0 ) && (~CPG) ; % Criterio de paso pequeno . i0 = i0 + 1; i f CPG % Alpha cumple CPG alphag = alpha ; i n i t = true ; elseif CPP % Alpha cumple CPP alphap = alpha ; alpha = alpha ∗gama ; else % Alpha cumple CPA index = 1; return;% Salimos con e x i t o de l a funcion end end i f ( i0 > itmax1 ) && ( i n i t == f a l s e ) % No se ha podido i n i c i a l i z a r Wolfe . % f p r i n t f ( fid , ’No se pudo i n i c i a l i z a r . Alphap= %e alphag= %e \n ’ , alphap , alphag ) ; return; end % Algoritmo i n i c i a l i z a d o . Disponemos de un alphap que cumple e l c r i t e r i o % de paso pequeno y un alphag que cumple e l c r i t e r i o de paso grande . % Buscamos un paso admisible entre ambos . fprintf( fid , ’ i n i c i a l i z a c i o n ␣ %e␣ %e␣\n ’ , alphap , alphag ) ; while ( i0 <= itmax2 ) % Seleccionamos e l nuevo alpha con e l metodo de b iseccion . alpha = ( alphag + alphap )/2; 7.6. PROGRAMAS UTILIZADOS 73 xalpha = x + alpha ∗d ; jalpha = f ( xalpha ) ; jpalpha = g ( xalpha ) . ’ ∗d ; CPG = ( jalpha > j0 + rho ∗jp0 ∗alpha ) ; % Criterio de paso grande . CPP = ( jpalpha < sigma ∗jp0 ) && (~CPG) ; % Criterio de paso pequeno . i0 = i0 + 1; i f CPG % Alpha cumple CPG alphag = alpha ; elseif CPP % Alpha cumple CPP alphap = alpha ; else % Alpha cumple CPA index = 1; fprintf( fi d , ’ s a l i d a ␣de␣ wolfe ␣ index= %i ␣\n ’ , index ) ; fprintf( fi d , ’ s a l i d a ␣de␣ wolfe ␣ alpha= %e␣\n ’ , alpha ) ; return; end end index = −1; return; end function [ x , i te r a c i o n ,K]= Interpolacion_cuadratica (x0 , d , c1 , c2 , f , g , itmax , tol , f i c 2 ) global A bb N R Rb %El m t o d o de i n t e r p o l quad busca aproximar phi ( alpha)=f ( x0+alpha ∗d) por un %polinomio q ( alpha)=a∗alpha^2 + b∗alpha+c , siendo d una direccion de iteracion=−1; %K=−10; %descenso 80 CAPÍTULO 7. ANEXOS end function z=H1semi (u , x , dfunf , iquad ) %ir e g = 1 malla regu lar %ir e g = 0 malla no regular num_pas=length(x)−1; z=0; %C l c u l o d e l cuadrado de la f u n c i n para obtener e l integrando i f iquad==1 %Regla del punto medio for i =2:num_pas+1 h=x( i )−x( i −1); aprox_derh=(u( i )−u( i −1))/h ; aprox_der=dfunf (( x( i )+x( i −1))∗0.5); z=z+h∗(aprox_derh−aprox_der )^2; end else i f iquad==2 %Cuadratura de Gauss con dos puntos for i =2:num_pas+1 h=x( i )−x( i −1); aprox_derh=(u( i )−u( i −1))/h ; xg1= (x( i )+x( i −1))∗0.5 + (x( i )−x( i −1))/(2∗sqrt (3)); xg2= (x( i )+x( i −1))∗0.5 −(x( i )−x( i −1))/(2∗sqrt (3)); aprox_der_1=dfunf ( xg1 ) ; aprox_der_2=dfunf ( xg2 ) ; z=z+h∗(( aprox_der_1−aprox_derh)^2+(aprox_der_2−aprox_derh )^2)∗0.5; end end 7.6. PROGRAMAS UTILIZADOS 81 end z=sqrt( z ) ; %calculamos por tanto l a norma L2 end function [ sol , t , index_conver ,A, bb ] = PNGC_parabola(y , par , pc , L , alph , bet , p , q , f , h) %y es la malla del FEM global N %Iniciamos e l gradi ente salida1=fopen( ’ gr ad iente ’ , ’w ’ ) ; f i c=fopen( ’ gradCuadWolfe ’ , ’w’ ) ; f i c 2=fopen( ’ gradInterp ’ , ’w ’ ) ; [ x0 , itmax , itmax2 , tol , tol2 , tol3 , rho , sigma , gama , a l phai n i t ]=datos_new ; %I ni c ializ a mos l o s p a r m e t r o s de s a l i d a i0 =1; index=−1; err =0; index_conver =0; s o l=x0 ; %punto i n i c i a l %−−−−−−−−−−−−−−−−−−−−−− %COMEZAMOS EL ENSAMBLADO MEF %−−−−−−−−−−−−−−−−−−−−−− %Ensamblado ( d ev uel ve t r e s d iago nale s de A y e l v ecto r segundo miembro ) % [ vd , vld , vud , bb]= assembly_Poncelet (N, x , h , p , q , f ) ; %Ensamblado de la matriz d e l problema A [ vd , vld , vud ,vmd, vmld ,vmud]=assembly_Poncelet_A (N, y , h , p , q ) ; %Creamos la Matriz A del m t o d o vd= vd+vmd; vld= vld+vmld ; 82 CAPÍTULO 7. ANEXOS vud=vud+vmud; %Imposicion de l a s condic iones de contorno para l a matriz A % y c o n s t r u c c i n de A con matriz sparse [A]=BC_A(N, vd , vld , vud ) ; for i =1:itmax %Ensamblado del vector de l segundo miembro b [ bb]=NL_assembly_Poncelet_b (N, sol , h , f ) ; %Imposicion de l a s condic iones de contorno para e l v ecto r b [ bb]=BC_b( alph , bet , bb ,L, vud , vld ) ; [ tf , tg , hg , hf ]= datos_extra (A, bb , pc ) ; i f i==1 %Datos i n i c i a l e s m t o d o de gradiente tgr ad ie nte = tg ( x0 ) ; td= −tg radi ente ; %con precondicionamiento hgradiente=hg ( x0 ) ; %No se cambia de signo l a direcci on i f pc==0 hgradiente=tg ( x0 ) ; %No se cambia de signo l a direccion end hd=−hgradiente ; % d i r e c c i n de descenso i n i c i a l nor_grad=sqrt( td ’∗hd ) ; % debemos minimizar l a f u n c i n g ( alpha)= f ( xk+alpha ∗dk ) ; fprintf(salida1 , ’∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗\n ’ ) ; fprintf( salida1 , ’ index=0␣␣norgrad=0␣ alpha=0␣\n ’ ) ; escribe_v ( ’ x=’ , sol , s a l ida1 ) ; escribe_v ( ’d=’ ,hd , s alida1 ) ; fprintf( salida1 , ’ e r r o r ␣ i n i c i a l ␣= %e␣ ’ , err ) ; escribe_v ( ’ norm_grad␣=’ , nor_grad , s a li da1 ) ; %escribe_cabecera(fic , sol ,d, err ); end 7.6. PROGRAMAS UTILIZADOS 83 %ESCOGEMOS UN ALPHA PTIMO i f par ==1 %Regla de Wolfe [ alpha , i0 , index ]= WolfeBiseccion ( tf , tg , sol , hd , rho , sigma , alphainit , gama , f i c ) ; elseif par==0 %I n t e r p o l a c i n c u a d r t i c a con r egla de Goldstein [ alpha , i0 , index ]= Interpolacion_cuadratica ( sol , hd , rho , sigma , hf , hg , itmax2 , tol2 , f i c 2 ) ; end %ACTUALIZAMOS EL PUNTO solold=sol ; s o l=s o l+alpha ∗hd ; %IMPLEMENTAMOS LA F RMULA DE POLAK−RIBI RE hgold=hg ( s o l o l d ) ; %guardamos e l gradient e ant er io r i f pc==0 hgold=tg ( s o l o l d ) ; end tgo ld=tg ( s o l o l d ) ; nor_grad=(tgold ’∗hgold ) ; %calculamos su norma hgradiente=hg ( s o l ) ; %gradiente nuevo con condicionante i f pc==0 hgradiente=tg ( s o l ) ; end tgr adi ent e=tg ( s o l ) ; %gradiente %Polak−R i b i r e beta=(tgradiente ’∗( hgradiente−hgold ))/ nor_grad ; i f beta<0 beta=0; end %NUEVA DIRECCI N DE DESCENSO hd=−hgradiente+beta∗hd ; %TEST DE PARADA MEDIANTE LA NORMA DEL GRADIENTE y DIFERENCIA CON EL ITERANTE ANTERIOR 84 CAPÍTULO 7. ANEXOS t=sqrt(nor_grad ); err=solold −s o l ; r=sqrt( err ’∗err ) ; fprintf(salida1 , ’∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗\n ’ ) ; fprintf( fi c 2 , ’ ∗∗∗∗∗∗∗∗∗ % i ∗∗∗∗∗∗∗∗∗\n ’ , i ) ; fprintf( fi c 2 , ’ alpha ␣= %e␣\n ’ , alpha ) ; fprintf( salida1 , ’ index= %i ␣␣ norgrad= %e␣ alpha= %e␣ i0= %i ␣ index= %i ␣ norerr= %e\n ’ , i , t , alpha , i0 , index , r ) ; escribe_v ( ’ x=’ , sol , s a l ida1 ) ; escribe_v ( ’d=’ ,hd , s alida1 ) ; escribe_v ( ’ x=’ , sol , f i c 2 ) ; escribe_v ( ’d=’ ,hd , f i c 2 ) ; %Primer t e s t de parada , norma gradiente i f t<t o l fprintf( ’ Se␣cumple␣ e l ␣ t e s t ␣de␣parada␣ del ␣ gr adi ent e \n ’ ) ; fprintf( ’ I t e r a c i o n ␣ %i \n ’ , i ) ; %escribe_vo ( ’ La s ol uc io n es : ’ , s o l ) ; %e s c r i b e por p a nt all a s o l fprintf( ’La␣norma␣ del ␣ gradiente ␣ %10.3e\n ’ , t ) ; fprintf( ’ D ife r en c ia ␣ i t e r a n t e ␣ a nt eri or ␣ %10.3e\n ’ , r ) ; index_conver=i ; %Escribimos en e l f ich e ro MEF e l vector s o l u c i n % f p r i n t f ( ofi , ’\ n\n∗El vector s o l u c i n es \n ’ ); % escribe_v ( ’ uh=’, sol , o f i ) ; residuo = A∗sol−bb ; residuo = residuo ’ ; normInf = norm(residuo ,Inf ) ; return; else %Segundo t e s t de parada , d ista nc i a entre i t e r a n t e s i f r<t o l 3 fprintf( ’ Se␣cumple␣ e l ␣ t e s t ␣de␣parada␣ del ␣ incremento \n ’ ) ; fprintf( ’ I t e r a c i o n ␣ %i \n ’ , i ) ; % escribe_vo ( ’ La s ol uc io n es : ’ , s o l ) ; %e s c r i b e por p a nta ll a s o l fprintf( ’ D ife r en c ia ␣ i t e r a n t e ␣ a nt eri or ␣ %10.3e\n ’ , r ) ; 7.6. PROGRAMAS UTILIZADOS 85 fprintf( ’La␣norma␣ del ␣ gradiente ␣ %10.3e\n ’ , t ) ; index_conver=i ; %Escribimos en e l f ich e ro MEF e l vector s o l u c i n % f p r i n t f ( ofi , ’\ n\n∗El vector s o l u c i n es \n ’ ); % escribe_v ( ’ uh= ’, sol , o f i ) ; residuo = A∗sol−bb ; residuo = residuo ’ ; normInf = norm(residuo ,Inf ) ; return; end end end i f index_conver==0 fprintf( ’No␣hubo␣ convergencia \n ’ ) ; end fclose( sali d a 1 ) ; fclose( f i c ) ; fclose( f i c 2 ) ; end function vaffineA=affineA (x , a , b) %Affine increasign map fomr [ a , b ] onto [ 0 , 1 ] vaffineA=(x−a ) . / ( b−a ) ; end function vaffineAInv=affineAInv (x , a , b) vaffineAInv=(b−a ) . ∗x+a ; end function [ b]=NL_assembly_Poncelet_b (N, uh , h , f ) %Realiza e l ensamblado d el vect or b %I ni c ializ a mos a cero todas l a s componentes de b b=zeros(1 ,N+1); 86 CAPÍTULO 7. ANEXOS for i = 1:N %Calculamos e l paso por Poncelet pto2=h( i ) / 2 ; ev=f (( uh( i )+uh( i +1))/2); %Calculamos vector b b( i )=b( i )+pto2∗ev ; b( i +1)=b( i +1)+pto2∗ev ; end % for i = 1:N % %Calculamos e l paso por Simpson % pto2=h( i )/6; % ev=f (( uh ( i )+uh( i +1))/2); % ev1=f (uh( i +1)); % ev2=f (uh( i +1)); % %Calculamos vector b % b ( i )=b ( i )+pto2 ∗( ev1+2∗ev ) ; % b ( i+1)=b ( i+1)+pto2 ∗(2∗ev+ev2 ) ; % end end function [ vd , vld , vud ,vmd, vmld ,vmud]=assembly_Poncelet_A (N, x , h , p , q) %Realiza e l c l c u l o y ensamblado de l a matriz de r i g i d e z (vd , vld , vud ) % y de masa (vmd , vmld , vmud) %Inicializamos a cero todas l a s diag onal es p r i n c i p a l e s vd=zeros(1 ,N+1); vld=zeros(1 ,N) ; vud=zeros(1 ,N) ; vmd=zeros(1 ,N+1); vmld=zeros(1 ,N) ; vmud=zeros(1 ,N) ; 7.6. PROGRAMAS UTILIZADOS 87 for i = 1:N %Calculamos punto medio xm=(x( i )+x( i +1))./2; pto=p(xm) . / h( i ) ; pto1=(h( i )∗q(xm) ) . / 4 ; %Calculamos diagonal de la matriz de r i g i d e z vd ( i )=vd ( i )+pto ; vd ( i +1)=vd( i+1)+pto ; %Calculamos subdiagonal de l a matriz de r i g i d e z vld ( i)=−pto ; %Calculamos superdiagonal de la matriz de r i g i d e z vud( i)=−pto ; %Calculamos diagonal de la matriz de masa vmd( i )=vmd( i )+pto1 ; vmd( i +1)=vmd( i +1)+pto1 ; %Calculamos subdiagonal de la matriz de masa vmld( i )=pto1 ; %Calculamos superdiagonal de la matriz de masa vmud( i )=pto1 ; end end function [A]=BC_A(N, vd , vld , vud) %Fijamos l a s cond iciones de contorno para sacar A, l a matriz d e l problema %Caso Dirichlet−Dirichlet vd (1)=1; vud (1)=0; vd (N+1)=1; vld (N)=0; vld(1)=0; vud(N)=0; 88 CAPÍTULO 7. ANEXOS %Construimos matriz A A = sparse ( 1 :N+1 ,1:N+1,vd ) ; A = A+sparse ( 1 :N, 2 :N+1,vud ,N+1,N+1)+sparse ( 2 :N+1 ,1:N, vld ,N+1,N+1); end function [ bb]=BC_b( alph , bet , bb , L , vud , vld ) global N a21 an_1n a21=vld ( 1 ) ; an_1n=vud (N) ; %Fijamos l a s cond iciones de contorno para e l vec tor b %Para e l vector de segundo miembro b solo hace f a l t a modificar %Los valores 1 y N+1 para l os nodos i n i c i a l y f i n a l . %Caso Dirichlet−Dirichlet bb(1)= alph/L ( 1 ) ; bb(N+1)=bet/L ( 3 ) ; bb (2) = bb(2)−alph∗a21/L ( 1 ) ; bb(N) = bb(N) −bet∗an_1n/L ( 3 ) ; bb=bb ’ ; end %ARCHIVO DE DATOS MEF global N % Extremo e s pac i al i z q u ierdo : fprintf( of i , ’ \n∗∗∗∗∗␣DATOS␣DEL␣PROBLEMA␣∗∗∗∗∗ ’ ) ; a = 0; fprintf( of i , ’ \n∗␣Extremo␣ e s p a c i a l ␣ i zq uie rdo : ␣a␣=␣ %−.15E. ’ , a ) ; 7.6. PROGRAMAS UTILIZADOS 89 % Extremo e s pa c i al derecho : b = 1; fprintf( of i , ’ \n\n∗␣Extremo␣ e s p a c i a l ␣ derecho : ␣b␣=␣ %−.15E. ’ , b ) ; %N mero de elementos N=100; fprintf( of i , ’ \n\n∗␣Numero␣de␣ elementos ␣N␣=␣ %d . ’ , N) ; %Funciones del problema funp= @p_NL_1; funq=@q_NL_1; funf=@f_NL_1; fexact=@exact_NL_1 ; dfunf=@df_NL_1; %Condiciones de contorno %L(1)∗u(a)+L(2)∗u ’ ( a ) = alph %L(3)∗u( b)+L(4)∗u ’ ( b ) = bet L=[1 0 1 0 ] ; fprintf( of i , ’ \n\n∗␣ vector ␣de␣ condiones ␣de␣ contorno ␣ ( Di rich let −Neumann−Dirichlet−Neumann) ’ ) ; fprintf( of i , ’ \n\n∗␣L=␣ %−.15E. ’ , L ) ; %Condiciones de contorno x=a & x=b alph = 0; fprintf( of i , ’ \n\n∗␣ alpha ␣=␣ %−.15E. ’ , alph ) ; bet =0; fprintf( of i , ’ \n\n∗␣ beta␣=␣ %−.15E. ’ , bet ) ; %Vector acumulador de puntos de la malla [ ] s i es uniforme vmesh = [ ] ; len=length(vmesh ) ; i f len == 0 fprintf( of i , ’ \n\n∗␣ e l ␣ vector ␣de␣puntos␣de␣acumulacion␣ es ␣ vacio ’ ) ; else fprintf( of i , ’ \n\n∗␣ vector ␣de␣puntos␣de␣acumulacion ’ ) ; fprintf( of i , ’ \n␣\n∗␣vmesh␣=␣ ’ ) ; 96 BIBLIOGRAFÍA [Rav83] Pierre-Arnaud Raviart. Introduction à l’analyse numérique des équations aux dérivées partielles. 1983. [Rei80] Laure Reinhart. Sur la résolution numérique de problemes aux limites non linéaires par des méthodes de continuation. PhD thesis, 1980. [Sha78] David F Shanno. Conjugate gradient methods with inexact searches. Mathematics of operations research, 3(3):244–256, 1978. [SY06] Wenyu Sun and Ya-Xiang Yuan. Optimization theory and methods: nonlinear programming, volume 1. Springer Science & Business Media, 2006. [TM21] Kevin Tolle and Nicole Marheineke. Extended group finite element method. Applied Numerical Mathematics, 162:1–19, 2021. [Vog02] Curtis R Vogel. Computational methods for inverse problems. SIAM, 2002.