Aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa

Aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Vicente García Marín Aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Vicente García Marín Memoria presentada como parte de los requisitos para la obtención del título de Grado en Matemáticas por la Universidad de Sevilla. Tutorizada por Prof. Tomás Chacón Rebollo Prof. Faustino Maestre Caballero Dpto. Ecuaciones Diferenciales y Análisis Numérico. Índice general English Abstract 1 Resumen 3 1. Introducción al problema de optimización 5 1.1. Resolución de una EDP elíptica no lineal . . . . . . . . . . . . . . . . 5 1.2. Caracterización del problema de contorno como problema del Cálculo deVariaciones............................... 8 1.2.1. Discretización........................... 9 2. Repaso de resultados del análisis funcional 11 2.1. Espaciosdefunciones........................... 11 2.2. Distribuciones............................... 14 2.3. Espacios de Sobolev. Propiedades . . . . . . . . . . . . . . . . . . . . 16 2.4. Algunos teoremas de interés . . . . . . . . . . . . . . . . . . . . . . . 18 3. F-diferenciabilidad y G-diferenciabilidad 19 4. Minimización de funcionales convexos 31 4.1. Convexidadfuerte............................. 33 ii aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa 4.2. Existencia y unicidad de óptimos . . . . . . . . . . . . . . . . . . . . . 35 5. Condiciones de optimalidad del problema discreto 39 5.1. Análisis de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.2. Estimación del error a priori del problema elíptico . . . . . . . . . . . 43 6. Métodos numéricos para la resolución del problema de optimización 47 6.1. Métododelgradiente ........................... 47 6.2. Método del gradiente conjugado . . . . . . . . . . . . . . . . . . . . . 51 6.2.1. Método del gradiente conjugado para funcionales cuadráticos 51 6.2.2. Método del gradiente conjugado no lineal . . . . . . . . . . . 54 6.3. MétododeNewton ............................ 57 7. Aplicación numérica: Resolución de una EDP elíptica no lineal 63 7.1. Resolución con FreeFem++ . . . . . . . . . . . . . . . . . . . . . . . . 64 A. Código FreeFem++ 69 English Abstract Generally, elliptic equations correspond to stationary models, i.e , models that do not depend on time. In this memory, we will solve PDEs that formulate a physical principle of minimum energy in mathematical terms. Specifically, the PDE is the optimality condition for the minimum of an energy functional. This arises in many different contexts, often describing the total energy of the state of a system, whether it be kinetic energy, potential energy, or other relevant forms. Despite its importance, we ocasionally find ourselves with the problem of solving this type of PDEs from an analytical and exact perspective, so we must resort to other methods that give us an approximate solution to the problem. The Calculus of Variations and Numerical Analysis will be used to solve this issue. Once it is formulated as an optimization problem and its characteristics are studied, conditions will be given to allow us to ensure the existence and uniqueness of a solution. Finally, some numerical algorithms will be developed in order to solve the variational problem associated with our nonlinear elliptic PDE. Resumen En general, las ecuaciones elípticas corresponden a modelos estacionarios, es decir, modelos que no dependen del tiempo. En este trabajo se resolverán EDPs que formulan en términos matemáticos un principio físico de energía mínimo. En concreto, la EDP es la condición de optimalidad del mínimo de un funcional de energía. Este surge en muchos contextos diferentes, pero a menudo describen la energía total del estado de un sistema, tratándose de energía cinética, potencial u otras formas relevantes para el mismo. Pese a su importancia, con frecuencia no es posible resolver este tipo de EDPs desde el punto de vista analítico y exacto, por lo que hay que recurrir a otros métodos que nos den una solución aproximada del problema. Se recurrirá al Cálculo de Variaciones y al Cálculo Numérico para resolver dicho problema. Una vez se formule como un problema de optimización y se estudien las características del mismo, se darán condiciones que nos permitan asegurar la existencia y unicidad de solución. Por último, se desarrollarán diversos algoritmos numéricos para la resolución del problema variacional asociado a nuestra EDP elíptica no lineal. 2 Repaso de resultados del análisis funcional En este capítulo se van a definir y recordar algunos conceptos básicos del análisis funcional que van a ser de suma importancia durante la elaboración del trabajo. Estas primeras definiciones se pueden encontrar en [9]y[15]. Definición 2.1. Sean XeYdos espacios normados sobre R. Diremos que la aplicación T:X−→ Yes lineal si para cualesquiera x, y de Xy cualesquiera α, β de R, se tiene: T(αx +βy) = αT(x) + βT (y). Definición 2.2. Sean XeYdos espacios normados sobre R. Diremos que la aplicación T:X−→ Yes continua si para toda sucesión {xn}n≥1⊂Xtal que xn→xen X, entonces se tiene que T(xn)→T(x). Definición 2.3. Se define L(X, Y )como el espacio vectorial de los operadores lineales y continuos del espacio Xen el espacio Y. Cuando X=Y, se denota simplemente L(X). Se define la norma del operador Tcomo: ∥T∥L(X,Y )= sup u∈X u=0 ∥Tu∥Y ∥u∥X = sup ∥u∥X=1 ∥Tu∥Y. Observación 2.1.Cuando Y=R, notamos L(X, R) = X′, el espacio dual de X. 2.1 Espacios de funciones Se definen en esta sección los principales espacios y sus propiedades, de las que vamos a hacer uso. Para ello, se han seguido generalmente [9]y[7] En adelante, Ω 12 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa denota un conjunto abierto no vacío de Rn. Definición 2.4. Sea φ∈C0(Ω). Se llama soporte de φal conjunto sop φ={x∈Ω|φ(x)= 0}. A partir de este espacio, es necesario considerar los siguientes conjuntos: C0 c(Ω) = {φ∈C0(Ω) |sop φ⊂Ωes compacto} Ck c(Ω) = C0 c(Ω) ∩Ck(Ω) D(Ω) = C0 c(Ω) ∩C∞(Ω) AD(Ω) se le conoce con el nombre de espacio de las funciones "test". Definición 2.5. Sea p∈[1,+∞]y supongamos que nuestro abierto Ωes medible Lebesgue. Se define el espacio Lp(Ω) como: Lp(Ω) = f: Ω −→ Rfes medible ,ZΩ |f(x)|p≤ ∞,si p < +∞. L∞(Ω) = f: Ω −→ Rfes medible y existe C > 0 : |f(x)|< C e.c.t Ω. Una vez definidos estos espacios, se definen las normas asociadas a ellos como sigue: ∥f∥Lp(Ω) =ZΩ |f(x)|pdx1/p ,∀f∈LP(Ω),si p < +∞. ∥f∥L∞(Ω) = ´ınf C > 0|f(x)|< C, e.c.t Ω. El siguiente resultado de densidad se puede encontrar probado en [7], y es de gran utilidad por los resultados que se obtienen de él. Teorema 2.1. D(Ω) es denso en Lp(Ω), para p∈[1,+∞). Como consecuencia se obtiene: Observación 2.2.Lp(Ω) es un espacio de Banach separable para p∈[1,+∞) Corolario 2.1.Sea f∈L2(Ω). Si ZΩ f(x)φ(x)dx = 0 ,∀φ∈D(Ω), entonces u= 0 e.c.t Ω. 2. repaso de resultados del análisis funcional 13 Definición 2.6. Sea funa función medible-Lebesgue en Ω. Se dice que fes localmente integrable en Ωsi se tiene que la restricción f|K∈L1(Ω), para cada K⊂Ω compacto. Se denomina L1 loc(Ω) al conjunto de las funciones localmente integrables en Ω. Se presentan a continuación una serie de importantes desigualdades que serán de gran utilidad, en especial, para demostrar acotaciones y continuidades de aplicaciones, como veremos más adelante. Definición 2.7. Sea p > 1. Se dice que p′>1es el exponente conjugado de psi 1 p+1 p′= 1 En particular, el exponente conjugado de p= 1 es p′= +∞. Desigualdad de Young: Sea p > 1yp′>1su exponente conjugado. Entonces: αβ ≤1 pα+1 p′β . Desigualdad de Hölder: Sea p > 1yp′>1su exponente conjugado. Entonces, para cualquier f∈Lp(Ω) y cualquier g∈Lp′(Ω) se verifica: ZΩ |f(x)g(x)|dx ≤ ∥f∥Lp(Ω)∥g∥Lp′(Ω) . En particular, esta desigualdad se utiliza para demostrar la desigualdad de Minkowski, que se presenta a continuación. Desigualdad de Minkowski: Sea p > 1y sean f, g ∈Lp(Ω). Entonces, se verifica: ZΩ |f(x) + g(x)|pdx1/p ≤ZΩ |f(x)|pdx1/p +ZΩ |g(x)|pdx1/p 14 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa 2.2 Distribuciones Para elaborar esta sección, se han seguido principalmente [9]y[7]. Para cada k≥0, tenemos las siguientes contenciones: D(Ω) ⊂Ck c(Ω) ⊂Ck(Ω) ⊂Ck(Ω) ⊂C0(Ω) ⊂L1 loc(Ω) , y para cada p∈[1,+∞], D(Ω) ⊂Lp(Ω) ⊂L1 loc(Ω) . Todos estos espacios se pueden dotar de estructuras topológicas adecuadas. En el caso de los espacio Lp(Ω), basta considerar la norma que hemos definido anteriormente. Igualmente, para cada k≥0, se puede dotar Ck(Ω) de la siguiente norma de forma que se convierte en un espacio de Banach separable: ∥φ∥Ck:= m´ax |α|≤km´ax x∈Ω∂αφ(x),∀φ∈Ck(Ω) . Para D(Ω), Ck(Ω) yL1 loc(Ω) en particular no se puede conseguir que sean completos con la norma. No obstante, en estos casos, es posible definir una topología natural compatible con la estructura de espacio vectorial, de los que estudiamos su convergencia. Definición 2.8. Sea {vn}n≥1una sucesión de L1 loc(Ω) y sea v∈L1 loc(Ω). Se dice que vn→ven L1 loc(Ω) si: ZK |vn−n|dx →0,∀K⊂Ωcompacto. Definición 2.9. Sea {φn}n≥1una sucesión de Ck(Ω) y sea φ∈Ck(Ω). Se dice que φn→φen CK(Ω) si: φn→φen Ck(K),∀K⊂Ωcompacto, es decir, si para todo multi-índice α= (α1, ..., αn)tal que |α|=α1+... +αn≤k, se tiene que ∂αφn→∂αφuniformemente en K. Definición 2.10. Sea {φn}n≥1una sucesión de Ck c(Ω) y sea φ∈CK c(Ω). Se dice que φn→φen CK(Ω) si existe K⊂Ωcompacto que contiene los soportes de todas las φny el soporte de φy, además se tiene: ∂αφn→∂αφuniformemente en K, para todo αmulti-índice tal que |α| ≤ k. 2. repaso de resultados del análisis funcional 15 Definición 2.11. Sea {φn}n≥1una sucesión de D(Ω) y sea φ∈D(Ω). Se dice que φn→φen D(Ω) si la propiedad anterior se cumple para todo valor de k≥0. Observación 2.3.Por las definiciones anteriores, tenemos que si {φn}n≥1⊂C0(Ω) es una sucesión tal que φn→φen C0(Ω), entonces φn→φen L1 loc(Ω). Diremos así que C0(Ω) se inyecta en L1 loc(Ω), con inyeción secuencialmente continua y se denotará como C0(Ω) ,→L1 loc(Ω). De igual forma, esto puede aplicarse a D(Ω) en Ck c(Ω),Ck c(Ω) en Ck(Ω), etc y por tan- to de nuevo diremos que estas inyeciones son continuas y se denotará como D(Ω) ,→ Ck c(Ω) ,→Ck(Ω), etc. Definición 2.12. Una distribución en Ωes cualquier aplicación T:D(Ω) −→ R lineal y secuencialmente continua. Observación 2.4.El conjunto de las distribuciones en Ωtiene estructura de espacio vecorial y se denota D′(Ω). Dada una distribución T, el valor asociado a φ∈D(Ω) por Tse denotará ⟨T, φ⟩. Definición 2.13. Se dice que Tn→Ten D′(Ω) si para todo φ∈D(Ω) se tiene ⟨Tn, φ⟩→⟨T, φ⟩. Definición 2.14. Sea T∈D′(Ω) y sea α= (α1, ..., αn)un multi-índice. Se denomina derivada distribucional o derivada en el sentido de las distribuciones de orden αde Ta la distribución ∂αTde D′(Ω) dada por ⟨∂αT, φ⟩= (−1)|α|⟨T, ∂αφ⟩,∀φ∈D(Ω). Para cada α, la aplicación T7→ ∂αTes lineal y secuencialmente continua. Este concepto de derivada distribucional extiende la definición de derivada parcial conocida. Proposición 2.1.Sea u∈C0(Ω) y supongamos que para todo xde Ω ∃∂iu(x),con ∂iu∈C0(Ω). Entonces, ∂iTu=T∂iu. 16 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa 2.3 Espacios de Sobolev. Propiedades En esta sección se definen y trabajan las principales propiedades de las que vamos a hacer uso sobre los espacios de Sobolev H1(Ω), H1 0(Ω) yH−1(Ω). Estos resultados se detallan en profunidad en [7]y[2]. Para más, ver [1]. Como anteriormente, Ωdenota un abierto no vacío de RN Definición 2.15. Se define el espacio de Sobolev H1(Ω) como: H1(Ω) = u∈L2(Ω) : ∃∂iu∈L2(Ω),∀i: 1 ≤i≤N, donde las derivadas se entienden como las derivadas distribucionales. Proposición 2.2.El espacio H1(Ω) dotado del siguiente producto escalar (u, v)H1(Ω) = (u, v)L2(Ω) + (∇u, ∇v)L2(Ω) ,∀u, v ∈H1(Ω),(2.1) y la norma inducida por dicho producto escalar ∥u∥H1(Ω) =∥v∥2 L2(Ω) +∥∇v∥2 L2(Ω)1/2 ,∀u, v ∈H1(Ω),(2.2) es un espacio de Hilbert. Observación 2.5.Por lo general, D(Ω) no es denso en H1(Ω). Sin embargo, cuando Ω = RNsi lo es. Por otro lado, si Ωes un abierto conexo de frontera regular, se tiene que D(Ω) es denso en H1(Ω). Definición 2.16. Se define el espacio de Sobolev H1 0(Ω) como la adherencia de D(Ω) en H1(Ω), es decir, H1 0(Ω) = D(Ω)∥·∥H1(Ω) . Proposición 2.3.H1 0(Ω) es un subespacio cerrado de H1(Ω) que, dotado de (2.1) se convierte en un espacio de Hilbert separable. En la práctica, la forma más habitual de caracterizar H1 0(Ω) es a través de la aplicación traza. Teorema 2.2. Sea Ω⊂RN, o bien Ω = RN +, o bien Ωes un abierto no vacío y acotado de frontera regular. Entonces existe una única aplicación γ0∈ L(H1(Ω); L2(∂Ω)) que llamaremos aplicación traza tal que γ0(φ) = φ|∂Ω,∀φ∈C1 c(Ω) El operador traza se puede extender por continuidad a H1(Ω) −→ L2(∂Ω) y de denotará ˜γ0, aunque en la práctica por abuso de notación se escribirá simplemente γ0. 2. repaso de resultados del análisis funcional 17 Teorema 2.3. Sea Ω⊂RNun abierto no vacío y acotado de frontera regular. Entonces, H1 0(Ω) = ker(γ0) = v∈H1(Ω) : γ0(v) = 0 Teorema 2.4 (Desigualdad de Poincaré). Sea Ω⊂RNun abierto no vacío que está acotado (al menos) en una dirección. Entonces existe una constante C > 0que solo depende de Ωtal que: ZΩ |v(x)|2dx ≤CZΩ |∇v(x)|2dx , ∀v∈H1 0(Ω).(2.3) Observación 2.6.Por lo general, la Desigualdad de Poincaré no es cierta para las funciones del espacio H1(Ω). Gracias a (2.3) se puede definir un nuevo producto escalar y por tanto una nueva norma equivalente en H1 0(Ω), (u, v)H1 0(Ω) = (∇u, ∇v)L2(Ω) ,∀u, v ∈H1 0(Ω). Corolario 2.2.La seminorma ∥v∥H1 0(Ω) =ZΩ |∇v(x)|2dx1/2 ,∀v∈H1 0(Ω) es una norma en H1 0(Ω) equivalente a la norma ∥·∥H1(Ω). Por tanto, (H1 0(Ω),∥·∥H1 0(Ω)) es un espacio de Hilbert. Definición 2.17. Se denomina H−1(Ω) al dual topológico de H1 0(Ω). Observación 2.7.Se puede ver la norma de H−1(Ω) como ∥T∥H−1(Ω) = sup v∈H1 0(Ω) v=0 ⟨T, v⟩ ∥v∥H1 0(Ω) ,∀T∈H−1(Ω). Observación 2.8.H−1(Ω) es un nuevo espacio de Hilbert, isométrico e isomorfo a H1 0(Ω). Con el siguiente resultado podemos ver una caracterización de H−1(Ω). Teorema 2.5. Sean φ0, ..., φn∈L2(Ω) y sea L:H1 0(Ω) −→ Runa forma lineal definida de la siguiente forma: L(v) = ZΩ φ0v+ n X i=1 φi ∂v ∂xi!dx , ∀v∈H1 0(Ω). Entonces L∈H−1(Ω). 18 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Teorema 2.6. Sea Ω⊂RNun abierto no vacío y acotado de frontera regular. Sean u, v ∈H1(Ω), entonces se satisface ZΩ u(x)∂v ∂xi (x)dx =−ZΩ v(x)∂u ∂xi (x)dx +Z∂Ω γ0u(x)γ0v(x)ni(x)dσ donde nicon i= 1, ..., N denota la i-ésima componente del vector normal exterior n= (ni)N i=1. La expresión anterior se conoce como primera identidad de Green sobre Ω, o simplemente fórmula de Green. 2.4 Algunos teoremas de interés Se enuncian a continuación algunos teoremas que se han utilizado a lo largo del trabajo y que son de gran interés. Este primer resultado se ha extraído de [15]. Teorema 2.7 (Hahn-Banach). Sea Mun subespacio de un espacio vectorial normado Xy sea ρ∈M′un operador lineal. Entonces, existe otra forma lineal Λ∈X′de forma que: Λx=ρ(x),∀x∈M ∥Λ∥X′=∥ρ∥M′ El siguiente resultado, así como su demostración, ha sido extraído de [11]. Teorema 2.8 (Convergencia dominada de Lebesegue). Sea {fn}una sucesión de funciones integrables. Supongamos que existe otra función integrable g≥0de forma que |fn| ≤ gpara todo n. Supongamos además que la sucesión {fn}converge casi por doquier hacia otra función f. Entonces, fes integrable y además l´ım n→∞ Zfndµ =Zfdµ. 3 F-diferenciabilidad y G-diferenciabilidad Las nociones de diferenciabilidad juegan un papel muy importante en el Cálculo de Variaciones pues son la base que permiten establecer criterios para afrontar y resolver los problemas de extremos funcionales. Para introducir en este capítulo las nociones del Cálculo diferencial en dimensión infinita, se han seguido principalmente [8]y[2]. A partir de ahora, VyZdenotan espacios de Banach. Definición 3.1. Dados dos espacios, VyZ, y una aplicación J:V−→ Z, diremos que Jes G-diferenciable en v∈Ven la dirección w∈Vsi existe d dt J(v+tw)t=0 = l´ım t→0 J(v+tw)−J(v) t. En tal caso, se denomina derivada direccional de Jen ven la dirección de wal valor d dt J(v+tw)t=0, y se denota por J′(v)w. Si Z=R, también se denota J′(v)w=⟨J′(v), w⟩. Definición 3.2. Si para todo w∈Vexiste J′(v)wy además la aplicación J′(v) : w∈V−→ J′(v)wes lineal y continua, se dice entonces que Jes G-derivable en v. Lema 3.1.El funcional de energía que se definió en el Capítulo 1es G-derivable. Demostración. Sea J:H1 0(Ω) −→ Rdado por J(u) = 1 2ZΩ |∇u|2+ZΩ B(u)−ZΩ fu. Entonces, J(u+tv) = 1 2ZΩ |∇(u+tv)|2+ZΩ B(u+tv)−ZΩ f(u+tv). 26 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Luego σ(ξ+λ) = σ(ξ) + λσ′(ξ) + λ2 2σ′′(ξ+θλλ),con 0< θλ<1 De aquí, dada h∈Lp(Ω), J(v+h) = ZΩ [σ(v(x)) + h(x)σ′(v(x)) + h2(x) 2σ′′(v(x) + θλh(x))] dx = =J(v) + ZΩ p|v(x)|p−1signo(v(x))h(x)dx +ZΩ α(x)dx = =J(v) + DJ(v)h+ZΩ α(x)dx con α(x) = h2(x) 2σ′′(v(x) + θxh(x)). Sea ahora γ(x) = p|v(x)|p−1signo(v(x)). Entonces, |γ(x)|p′=pp′|v(x)|(p−1)p′, pero (p−1)(p′−1) = 1 ⇒(p−1)p′= 1+(p−1) = p. Por tanto, |γ(x)|p′=pp′|v(x)p|. De aquí se se deduce que RΩγ(x)h(x)dx =DJ(v)hy además DJ(v)∈ L(Lp(Ω),R). Por último, queda probar que RΩα(x)dx =∥h∥Lp(Ω)ε(h), donde l´ım h→0ε(h)=0: Escribimos ahora α=h(hδ), donde δva a ser δ=1 2σ′′(v+θxh). Así, integrando y aplicando la desigualdad de Holder, ZΩ α dx ≤ZΩ |h|pdx 1/p ZΩ |hδ|p′dx 1/p′ . Ahora bien, δp′= [ p(p−1) |v+θxh|p−2]p′≤[cp(|v|p−2+|h|p−2)]p′≤ ≤cp′|v|(p−2)p′+|h|(p−2)p′. Y como (p−2)p′=p−p′,entonces: |hδ|p′≤cp′|v|p−p′|h|p′+|h|p. De nuevo, integrando y aplicando la desigualdad de Holder, ZΩ |hδ|p′dx ≤ZΩ |v|α(p−p′)dx 1/ρ ZΩ |h|α′p′dx 1/ρ′′ +cp′∥h∥p p.(3.3) 3. f-diferenciabilidad y g-diferenciabilidad 27 Se impone ρ′p′=p;ρ′=p p′. Tiene que darse que ρ′sea mayor o igual que 1, pero esto se cumple puesto que anteriormente se ha impuesto p≥2. Con esta condición, se llega a que α(p−p′) = p. Sustituyendo así en (3.3): ZΩ |hδ|p′dx 1/p′ ≤cp′∥v∥p−p′ Lp(Ω)∥h∥p′ Lp(Ω) +∥h∥p Lp(Ω)1/p′ . Así, ∥h∥Lp(Ω)ε(h)≤1 2∥h∥Lp(Ω)cp′∥v∥p−p′ Lp(Ω)∥h∥p′ Lp(Ω) +∥h∥p Lp(Ω)1/p′ . Y se concluye que l´ım ∥h∥Lp(Ω)→0ε(h) = 0. Por último, volvamos al caso de estudio del tema y veamos la F-derivabilidad del funcional (1.2). En este caso, se consideran Ωun abierto acotado contenido en R3, los espacios V=H1 0(Ω) ,Z=Ry el funcional J:H1 0(Ω) −→ Rtal que J(v) = 1 2ZΩ |∇v|2+ZΩ B(v)−ZΩ fv. Se prueba que Jes F-derivable en todo v∈H1 0(Ω) usando la definición 3.4 J(v+h) = 1 2ZΩ |∇(v+h)|2+ZΩ B(v+h)−ZΩ f(v+h). Desarrollando en serie de Taylor el integrando B(v+h), se llega a J(v+h) = J(v) + ZΩ ∇v∇h+ZΩ B′(v)h−ZΩ fh +1 2ZΩ |∇h|2+ +ZΩ 1 2B′′(v(x) + θxh(x))h2(x)dx, donde ZΩ ∇v∇h+ZΩ B′(v)h−ZΩ fh =⟨J′(v), h⟩, y por tanto ε(h)∥h∥H1 0(Ω) =1 2ZΩ |∇h|2+ZΩ 1 2B′′(v(x) + θxh(x))h2(x)dx . Es decir, para probar la F-derivabilidad, ha de verificarse: J′(v)∈ L(H1 0(Ω),R) 28 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa l´ım ∥h∥H1 0(Ω)→0ε(h)=0 En primer lugar, para probar que J′(v)∈ L(H1 0(Ω),R)veamos que se trata de una transformación lineal y continua. Fijamos J′(v)y tomamos α1, α2∈Ryh1, h2∈ H1 0(Ω). ⟨J′(v), α1h1+α2h2⟩=ZΩ ∇v∇(α1h1+α2h2) + ZΩ B′(v)(α1h1+α2h2)− −ZΩ f(α1h1+α2h2) = α1ZΩ ∇u∇h1+α2ZΩ ∇u∇h2+ +α1ZΩ B′(v)h1+α2ZΩ B′(v)h2−α1ZΩ f h1−α2ZΩ f h2= =α1⟨J′(v), h1⟩+α2⟨J′(v), h2⟩. (3.4) Por otro lado, |⟨J′(v), h⟩| ≤ ZΩ ∇v∇h+ZΩ B′(v)h+ZΩ fh ≤ ≤ZΩ |∇v| |∇h|+ZΩ |B′(v)| |h|+ZΩ |f| |h| ≤ ≤ZΩ |∇v|21/2ZΩ |∇h|21/2 +ZΩ |B′(v)|6/55/6ZΩ |h|61/6 + +ZΩ |f|21/2ZΩ |h|21/2 . Por construcción y por (1.1), |B′(t)|=|β(t)| ≤ c1+c3|t|5. Sustituyendo este nuevo valor en la integral ZΩ |B′(v)|6/55/6 ≤ZΩ |c1+c3|v|5|6/55/6 ≤c4+c5∥v∥L6(Ω)5/6≤c4+c5∥v∥H1(Ω)5/6. 3. f-diferenciabilidad y g-diferenciabilidad 29 Por último, basta aplicar la desigualdad de Poincaré: |⟨J′(v), h⟩| ≤ ∥v∥H1 0(Ω) ∥h∥H1 0(Ω) +c4+c5∥v∥H1(Ω)5/6∥h∥H1(Ω) +∥f∥L2(Ω)∥h∥H1(Ω) ≤ ≤ ∥v∥H1 0(Ω) ∥h∥H1 0(Ω) +c4+c5∥v∥H1(Ω)5/6∥h∥2 L2(Ω) +∥h∥2 H1 0(Ω)1/2+ +∥f∥L2(Ω) ∥h∥2 L2(Ω) +∥h∥2 H1 0(Ω)1/2≤ ∥v∥H1 0(Ω) ∥h∥H1 0(Ω)+ +c4+c5∥v∥H1(Ω)5/6C1∥h∥H1 0(Ω) +∥f∥L2(Ω) C2∥h∥H1 0(Ω) = =∥v∥H1 0(Ω) +c4+c5∥v∥H1(Ω)5/6C1+∥f∥L2(Ω) C2∥h∥H1 0(Ω). (3.5) De (3.4) se deduce la linealidad de J′(v)y de (3.5) la continuidad del mismo. En consecuencia, J′(v)∈ L(H1 0(Ω),R) Una vez probado esto, faltaría ver que l´ım ∥h∥H1 0(Ω)→0ε(h) = 0. Por un lado, B′′ =β′y |β′(t)| ≤ c1+c2|t|4. Por tanto, |B′′(v+θh)| ≤ c1+c2|v+θh|4≤c1+c3(|v|4+|h|4) Luego, sustituyendo en la integral y tomando valores absolutos: ZΩ B′′(v+θh)h2≤ZΩ |B′′(v+θh)h2| ≤ c1ZΩ h2+c3ZΩ |v|4|h|2+c3ZΩ |h|6≤ ≤c(c1+c3)∥h∥2 H1 0(Ω) +∥h∥6 H1 0(Ω)+c3ZΩ |v|4α1/α ZΩ |h|2α′1/α′ (3.6) donde en el último paso se tienen en cuenta las inyecciones de Sobolev H1 0(Ω) ,→ L2(Ω) yH1 0(Ω) ,→L6(Ω) y se aplica la desigualdad de Hölder. Así, tomando 4α= 6 se tiene α=6 4=3 2. Por tanto, α′= 3. Sustituyendo y aplicando normas en (3.6), ZΩ B′′(v+θh)h2≤c4∥h∥2 H1 0(Ω) +∥h∥6 H1 0(Ω)+c5∥v∥4 H1 0(Ω)∥h∥2 H1 0(Ω) ⇒ | ε(h)∥h∥H1 0(Ω) | ≤ 1 2∥h∥2 H1 0(Ω) +C∥h∥2 H1 0(Ω) +∥h∥6 H1 0(Ω) +∥v∥4 H1 0(Ω)∥h∥2 H1 0(Ω) Basta ahora tomar límites en la expresión anterior y l´ım ∥h∥H1 0(Ω)→0ε(h)=0 Con lo que se concluye que Jes por tanto F-derivable en todo v∈H1 0(Ω). 4 Minimización de funcionales convexos En este capítulo, se introducen los conceptos de convexidad y convexidad fuerte de un funcional, que serán de gran importancia a la hora de probar la existencia y unicidad del problema sobre el que se está trabajando. Para la elaboración, se han seguido principalmente [2], [5]y[8]. Definición 4.1. Sea un espacio de Hilbert V,Ω⊂Vun abierto y U⊂Ωun convexo no vacío. Dado un funcional J: Ω →R, diremos que Jes convexo en Usi J(θu + (1 −θ)v)≤θJ(u) + (1 −θ)J(v),∀u, v ∈U, ∀θ∈(0,1) Se dirá que Jes estrictamente convexo si la desigualdad es estricta cuando u=v. Lema 4.1.Supongamos que Jes G-derivable en Ω. Entonces: 1. Jes convexo en Usi y sólo si para todas u, v ∈Use verifica J(v)−J(u)≥ ⟨J′(u), v −u⟩(4.1) 2. Jes estrictamente convexo si la desigualdad anterior es estricta cuando u=v. Demostración. Supongamos en primer lugar que Jes convexo. Se puede ver que para todo k≥1 J1−1 2ku+1 2kv≤1−1 2kJ(u) + 1 2kJ(v), de donde 2kJu+1 2k(v−u)−J(u)≤J(v)−J(u). 32 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Tomando ahora límites con k tendiendo a +∞, se llega a (4.1). De manera recíproca, supongamos que (4.1) se cumple y tratemos de ver que Jes convexo en U. Sean u, v ∈Uy tomamos θ∈(0,1). Entonces: J(u)≥J(v+θ(u−v)) + ⟨J′(v+θ(u−v),(1 −θ)(u−v))⟩ J(v)≥J(v+θ(u−v)) + ⟨J′(v+θ(u−v),−θ(u−v))⟩ Si ahora multiplicamos las desigualdades anteriores por θy1−θrespectivamente, y teniendo en cuenta el carácter homogéneo de la G-derivada respecto de la dirección, se llega a J(θu + (1 −θ)v)≤θJ(u) + (1 −θ)J(v), lo que prueba la convexidad de Jen U. Esta misma prueba sirve para demostrar que, bajo desigualdad estricta, la expresión (4.1) implica que Jes estrictamente convexo en U. Por otro lado, supongamos que Jes estrictamente convexo. En particular, Jes convexo, luego verifica (4.1). Tomamos ahora u, v ∈Utales que u=vyθ∈(0,1). Así, J(v)−J(u)>J(u+θ(v−u)) −J(u) θ≥⟨J′(u, θ(v−u))⟩ θ=⟨J′(u), v −u⟩. De esta forma se concluye el resultado. Observación 4.1.La condición de convexidad tiene una interpretación geométrica sencilla. Significa que el funcional convexo J(v)va a estar siempre por encima de su plano tangente a u. Lema 4.2.El funcional sobre el que se ha trabajado anteriormente J(u) = 1 2∥∇u∥2 L2(Ω) +ZΩ B(u)−ZΩ fu es estrictamente convexo en todo v∈H1 0(Ω). Demostración. Sean u, w ∈H1 0(Ω), ⟨J′(u), w⟩=ZΩ ∇u∇w+ZΩ β(u)w+ZΩ fw 4. minimización de funcionales convexos 33 ⇒J(v)−J(u)− ⟨J′(u), v −u⟩=1 2ZΩ |∇v|2−1 2ZΩ |∇u|2−ZΩ ∇u∇(v−u) + +ZΩ [B(v)−B(u)−B′(u)(v−u)] −ZΩ fv +ZΩ fu +ZΩ f(v−u) Considerando el desarollo en serie de B(v)hasta segundo orden B(v) = B(u)+(v−u)B′(u) + 1 2(v−u)2B′′(u+θx(v−u)), θx∈(0,1) y sustituyendo en la igualdad anterior, J(v)−J(u)− ⟨J′(u), v −u⟩=1 2ZΩ |∇(v−u)|2+ZΩ 1 2B′′(u+θx(v−u))(v−u)2 ≥1 2ZΩ |∇(v−u)|2>0, si u=v donde hemos usado que B′′ =β′≥0. Por tanto, la razón fundamental para que Jsea estrictamente convexo, además de la regularidad, es que βsea una función no decreciente, como se indica en las condiciones de la sección 1.1 . 4.1 Convexidad fuerte A continuación, vamos a introducir la idea de "convexidad fuerte", que será utilizada de aquí en adelante. Al igual que en la sección precedente, se han seguido principalmente [2], [5]y[8]. Definición 4.2. Sea Vun espacio de Hilbert y U⊂Vun subcojunto convexo no vacío. Se dirá que J:U⊂V−→ Res un funcional fuertemente convexo en Usi Jes continuo en Uy existe α > 0tal que para todo θ∈[0,1], se verifica J((1 −θ)u+θv)≤(1 −θ)J(u) + θJ(v)−αθ(1 −θ) 2∥u−v∥2,∀u, v ∈U(4.2) En tal caso, diremos que Jes un funcional α-convexo o elíptico con constante de elipticidad igual a α. 34 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Proposición 4.1.Sea Vun espacio de Hilbert, Ω⊂Vun abierto y U⊂Ωun convexo no vacío. Dado un funcional J: Ω −→ RG-diferenciable, las siguientes afirmaciones son equivalentes: Jes α-convexa.(4.3) J(v)≥J(u) + ⟨J′(u), v −u⟩+α 2∥v−u∥2,∀u, v ∈U. (4.4) ⟨J′(u)−J′(v), u −v⟩ ≥ α∥u−v∥2,∀u, v ∈U. (4.5) Demostración. Supongamos en primer lugar que tenemos (4.3). Como Jes α-convexo, si tomamos θ∈(0,1) y aplicamos (4.2), obtenemos J(u+θ(v−u)) −J(u) θ≤J(v)−J(u)−α 2(1 −θ)∥v−u∥2. Si tomamos límite en esta última expresión cuando θ→0+, deducimos (4.4). Veamos ahora (4.4)=⇒(4.5). Para ello, suponemos (4.4) y sean u, v ∈V. Si lo aplicamos a (u, v)y(v, u), J(v)≥J(u) + ⟨J′(u), v −u⟩+α 2∥v−u∥2, J(u)≥J(v) + ⟨J′(v), u −v⟩+α 2∥u−v∥2, Y sumando estas expresiones, obtenemos −⟨J′(u), v −u⟩−⟨J′(v), u −v⟩ ≥ α∥u−v∥2 de donde se deduce ⟨J′(u)−J′(v), u −v⟩ ≥ α∥u−v∥2. Por último, veamos que (4.5)=⇒(4.3). Sean u, v ∈Vy sea t∈R. Consideramos una nueva función φ(t) = J(u+t(v−u)). Esta función es derivable, con φ′(t) = ⟨J′(u+t(v−u)), v −u⟩, y por tanto continua. De (4.5), tenemos φ′(t)−φ′(s)≥α(t−s)∥v−u∥2,si t≥s. 4. minimización de funcionales convexos 35 Fijamos θ∈(0,1). Si integramos ahora esta última desigualdad desde t=θhasta t= 1, y desde s= 0 hasta s=θ, llegamos a θφ(1) + (1 −θ)φ(0) −φ(θ)≥αθ(1 −θ) 2∥v−u∥2, y por tanto Jes α-convexo en U. Observación 4.2.En la proposición anterior, en el caso de tener α= 0,Jes convexa sobre U. Lema 4.3.El funcional sobre el que se ha trabajado anteriormente J(u) = 1 2∥∇u∥2 L2(Ω) +ZΩ B(u)−ZΩ fu es α-convexo. Demostración. Es inmediato gracias a la la prueba del lema 4.2, donde se llega siguiente a la expresión: J(v)−J(u)− ⟨J′(u), v −u⟩ ≥ 1 2ZΩ |∇(v−u)|2=1 2∥v−u∥2 H1 0(Ω), y por tanto Jes α-convexo para α= 1. Observación 4.3.Es obvio que si un funcional Jes α-convexo en U, entonces, Jes estrictamente convexo en U. 4.2 Existencia y unicidad de óptimos En esta sección se dan una serie de resultados en relación con la existencia y unicidad de mínimos (óptimos) de una función. Principalmente, se han empleado [2] y [16]. También se pueden consultar los resultados en [8]. Definición 4.3. Un funcional J:U−→ R, con Uno acotado, se dice coercivo en el infinito si l´ım ∥v∥U→+∞ v∈U J(v) = +∞. Teorema 4.1. Sea U⊂Vun conjunto cerrado, convexo y no vacío. Sea J:U−→ R un funcional continuo, convexo y coercivo en el infinito. Entonces el problema de optimización (PO)Obtener u=arg m´ın v∈UJ(v) admite al menos una solución. 42 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Como Jes continuo, l´ım k→∞ J(ˆuk) = J(u). Es más, por ser uksolución del problema aproximado (Pk), vamos a tener J(uk)≤J(ˆuk),ˆuk∈Vk.(5.1) Por tanto, {J(uk)}k≥1está acotada. Luego, {uk}k≥1es una sucesión acotada. 2) Extracción de subsucesiones convergentes por compacidad. Por hipótesis del problema, Ves un espacio de Hilbert. Por tanto, existe una subsucesión {uk′}k′≥1débilmente convergente a u∗en V. Además, se tiene la inyección H1 0(Ω) ,→Lp(Ω), con 1≤p≤6y compacta para p < 6, por lo que podemos asumir que la sucesión converge fuertemente hacia u∗en cada uno de estos espacios. Teniendo en cuenta ahora el Teorema de Lebesgue inverso, podemos encontrar {uk′′ }k′′ ≥1 convergente c.p.d (a u∗) en Ωy una función g∈Lp(Ω), con 1≤p≤6tales que |uk′′ | ≤ gc.p.d en Ω. Al ser Buna función continua, se tiene B(uk′′ )−→ B(u∗)c.p.d en Ω. Suponiendo ahora |β′(t)| ≤ c1+c2|t|4−εyβ(0) = 0, entonces |B(t)| ≤ c3+c4|t|6−ε⇒ |B(uk′′ )| ≤ c3+c4|uk′′ |6−ε≤c3+c4|g|6−ε. Por el teorema de la convergencia dominada, existe l´ım k→∞ ZΩ B(uk′′ ) = ZΩ B(u∗).(5.2) Como uk′′ −→ u∗, entonces se tiene ZΩ f uk′′ −→ ZΩ f u∗.(5.3) Como además la norma es débilmente semicontinua inferiormente, ∥∇u∗∥L2(Ω) ≤l´ım inf k′′→∞ ∥∇uk′′ ∥L2(Ω) (5.4) Por tanto, J(u∗)≤l´ım inf k′′→∞ J(uk′′ )(5.5) De 5.1, como l´ım k→∞ J(ˆuk) = J(u), se va a tener l´ım inf k′′→∞ J(uk′′ )≤J(u). Esto nos lleva aJ(u∗)≤J(u). Como el problema de optimización admite unicidad de solución, al ser J(u) = ´ınf v∈VJ(v), se deduce inmediatamente que u∗=u. 5. condiciones de optimalidad del problema discreto 43 Se tiene que (PO)admite solución única, luego la sucesión {uk}k≥1converge débilmente a u. ∀v∈V , (∇uk,∇v)−→ (∇u, ∇v). En el caso de que esto no fuera cierto, existirían ciertos ˆv∈V,ε > 0y{uk′}k′≥1de forma que ε≤ |(∇uk,∇ˆv)−(∇u, ∇ˆv)|. Por el razonamiento anterior, existe {uk′′′ }k′′′ ≥1−→ u∗, lo que contradice el hecho de que aquí no se pueda pasar al límite. 3) Convergencia fuerte. Toda la sucesión converge débilmente a u.De5.1, se tiene que l´ım sup k→∞ J(uk)≤J(u). De 5.5,J(u)≤l´ım inf k→∞ J(uk). Así, va a existir l´ım k→∞ J(uk) = J(u). Por los razonamientos seguidos en la etapa anterior, l´ım k→∞ ZΩ B(uk) = ZΩ B(u)y l´ım k→∞ ZΩ f uk=ZΩ f uk. Se concluye así que l´ım k→∞ ∥∇uk∥L2(Ω) =∥∇u∥L2(Ω). Por lo tanto, ukconverge fuertemente a uen H1 0(Ω), como se quería probar. 5.2 Estimación del error a priori del problema elíptico Bajo las condiciones planteadas en el Capítulo 1para el problema que estamos estudiando en esta memoria, podemos obtener la siguiente cota de error. Teorema 5.5. Supongamos que β∈C2(R)y que |β′′(t)| ≤ M1+M2|t|3para todo t∈R. Entonces, existe c > 0que depende de Ωtal que ∥∇(u−uk)∥L2(Ω) ≤c1 + ∥∇u∥4 L2(Ω)d1(u, Vk), donde d1(u, Vk) = ´ınf v∈Vk ∥∇(u−v)∥L2(Ω) . Observación 5.1.La obtención a priori de estimaciones de error se transforma de esta forma en un nuevo problema de aproximación. 44 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Demostración. Dados u∈V, uk∈Vk, estos se caracterizan por (∇u, ∇v) + ZΩ β(u)v=ZΩ f v :⟨J′(u), v⟩= 0 , (∇uk,∇vk) + ZΩ β(uk)vk=ZΩ f vk:⟨J′ k(uk), vk⟩= 0 , siendo Jk=J|Vk. Se observa que J′es monótona: ⟨J′(v1)−J′(v2), v1−v2⟩ ≥ ZΩ |∇(v1−v2)|2,∀v1, v2∈V Sea ahora pkla proyección ortogonal de usobre Vk. Se tiene que ⟨J′(pk), v⟩=⟨J′(pk)−J′(u), v⟩=⟨δ, v⟩, siendo δ=J′(pk)−J′(u). Así, ∥∇(pk−uk)∥2 L2(Ω) ≤ ⟨J′(pk)−J′(uk), pk−uk⟩=⟨δ , pk−uk⟩. Entonces, ∥∇(u−uk)∥L2(Ω) ≤ ∥∇(u−pk)∥L2(Ω) +∥∇(pk−uk)∥L2(Ω) ≤ ≤d1(u, Vk) + ∥δ∥V′. Entonces, basta estimar ∥δ∥V′= sup w∈V w=0 ⟨δ, w⟩ ∥∇w∥L2(Ω) , ⟨δ, w⟩=fw(pk)−fw(u), donde para un cierto w∈Vfijado, se define el funcional fw:V−→ Rque a cada v∈Vle hace corresponder fw(v) = ⟨J′(v), w⟩. Si probamos que fwes F-derivable, tendremos que fw(pk)−fw(u) = ϕ(1) −ϕ(0) , con ϕ(t) = f(tpk+ (1 −t)u)) = ϕ′(θ), para θ∈(0,1) y ϕ′(t) = ⟨f′ w(tpk+ (1 −t)u), pk−u⟩. 5. condiciones de optimalidad del problema discreto 45 Así, se tiene ⟨δ, w⟩=⟨f′ w(σθ k), pk−u⟩, con σθ k=θpk+ (1 −θ)u. Por un lado, teníamos fw(v) = ⟨J′(v), w⟩=ZΩ ∇v∇w+ZΩ β(v)w−ZΩ fw . Por tanto, fw(v+h) = ZΩ ∇v∇w+ZΩ β(v)w−ZΩ fw +ZΩ ∇h∇w+ +ZΩ β′(v)h w +1 2ZΩ β′′(v+θλh)h2w−ZΩ fw . Luego, fw(v+h) = fw(v) + ZΩ ∇h∇w+ZΩ β′(v)h w | {z } ⟨f′ w(v),h⟩ + +1 2ZΩ β′′(v+θλh)h2w | {z } o(h) Si ahora usamos que |β′′(t)| ≤ M1+M2|t|3, obtenemos |σ(h)|=c∥w∥L3/2(Ω) +∥v∥3 L6(Ω)∥w∥L2(Ω) +∥h∥3 L6(Ω)∥w∥L6(Ω)∥∇h∥2 L2(Ω) . Luego fwes F-derivable, y ⟨δ, w⟩=ZΩ ∇(pk−u)∇w+ZΩ β′(σθ k)(pk−u)w ⇒ ⟨δ, w⟩ ≤ ∥∇(pk−u)∥L2(Ω)∥∇w∥L2(Ω) +ZΩc1+c3(|pk|4+|u|4)(pk−u)w. Si ahora usamos la desigualdad de Holder, podemos seguir acotando la expresión anterior: ⇒ ⟨δ, w⟩ ≤ ∥∇(pk−u)∥L2(Ω)∥∇w∥L2(Ω) +c3∥pk∥4 L6(Ω) +∥u∥4 L6(Ω)∥pk−u∥L6(Ω)∥w∥L6(Ω) ≤ ∥∇(pk−u)∥L2(Ω)∥(1 + c4(Ω)) ∥pk∥4 H1(Ω) +∥u∥4 H1(Ω)∥∇w∥L2(Ω) (5.6) 46 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa donde en la última desigualdad hemos hecho uso de H1 0(Ω) ,→L6(Ω). Por otro lado, en la segunda desigualdad hemos usado ZΩ |pk|4(pk−u)w≤ZΩ |pk|4α1 αZΩ |pk−u|β1 βZΩ |w|γ1 γ ≤ ≤ ∥pk∥4 L6(Ω) ∥pk−u∥L6(Ω) ∥w∥L6(Ω) . Volviendo a 5.6, ⟨δ, w⟩ ≤ h1+2c4∥∇u∥4 L2(Ω)∥∇w∥L2(Ω)i∥∇(pk−u)∥L2(Ω) . Así, ∥δ∥V′≤1+2c4∥∇u∥4 L2(Ω)∥∇(pk−u)∥L2(Ω) ⇒ ∥∇(uk−u)∥L2(Ω) ≤c51 + ∥∇u∥4 L2(Ω)d1(u, Vk). Observación 5.2.En el caso de Ω⊂R2,H1 0(Ω) ,→Lp(Ω), con 1≤p < +∞; todos los resultados vistos son ciertos si β, β′, β′′ tiene crecimiento polinómico cualquiera. 6 Métodos numéricos para la resolución del problema de optimización El objetivo de este capítulo es analizar y desarrollar algoritmos efectivos que nos permitan aproximar la solución del problema de optimización discreto que estamos estudiando: (Pk)Obtener uk=arg m´ın v∈Vk J(v). En particular, se van a desarrollar los algoritmos del gradiente (con paso fijo y paso óptimo), el algoritmo del gradiente conjugado (y sus variantes) y el método de Newton. Posteriormete, este último será utilizado de manera práctica en el Capítulo 7. También serán vistas las propiedades de convergencia de estos métodos con el fin de estimar la calidad de la solución. 6.1 Método del gradiente Los métodos que se van a estudiar son de carácter iterativo, es decir, parten de un punto inicial del espacio y van construyendo una sucesión {un}n∈Nque, bajo ciertas hipótesis, va a converger a la solución udel problema. En particular, suponiendo un funcional Jque sea F-diferenciable en un espacio de Hilbert V, vamos a plantear métodos del tipo descenso, caracterizados por tener la estructura del tipo un+1 =un−γdn, con γ∈Ry siendo dnla llamada dirección de descenso. Dentro de estos tipos de método, uno de los más conocidos es el método del gradiente, que es aquel en el que se toma como dirección de descenso dn=−J′(un). Este método y sus resultados han sido extraídos principalmente de [2], [16], [4]y[6]. 48 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa Algoritmo del método del gradiente con paso fijo 1. Inicialización: Se escoge u0∈Vk,γ > 0y fijamos ε > 0. 2. Iteraciones: Conocido previamente un∈Vk, a) Calcular J′(un)∈V. b) Calcular un+1 =un−γJ′(un). 3. Test de parada: Parar si ∥J′(un)∥V< ε. En general, se suele tomar como regla de parada cuando ∥un+1 −un∥< ε. Una opción para mejorar el método en cada iteración es ir actualizando nuestra longitude de paso γ(método del gradiente con paso óptimo). Para ello, podemos utilizar la regla de minimización, que consiste en elegir γktal que γk=arg m´ın γ>0J(un−γJ′(un)). La convergencia de dicho algoritmo puede estudiarse a partir del siguiente resultado. Teorema 6.1. Sea Jun funcional α-convexo y F-diferenciable y supongamos que J′ es globalmente Lipschitziana sobre los acotados de V. En tal caso, el método del gradiente con paso óptimo converge globalmente. Es decir, para cualquier inicialización del algoritmo u0, la sucesión (un)generada por el mismo converge a la solución udel problema de minimización. Demostración. Supongamos J′(un)= 0 ∀n≥0, pues en caso contrario se tendría J′(un)=0para cierto ny el método finaliza y converge con solución u∗=un. Introducimos la función f:γ∈R−→ f(γ) = J(un−γJ′(un)). Probemos que fes estrictamente convexa. Como Jes α-convexa, en particular es estrictamente convexa. Para cualquier valor de θ∈(0,1) yγ1=γ2, tenemos: f(θγ1+ (1 −θ)γ2) = J(un−(θγ1+ (1 −θ)2)J′(un)) =J(θ(un−γ1J′(un)) + (1 −θ)(un−γ2J′(un))) < θJ(un−γ1J′(un)) + (1 −θ)J(un−γ2J′(un)) =θf(γ1) + (1 −θ)f(γ2) y por tanto fes estrictamente convexa. 6. métodos numéricos para la resolución del problema de optimización 49 Por otro lado, si consideramos de nuevo que Jes α-convexa, por la Proposición 4.1, tenemos que J(u)≥J(v) + ⟨J′(v), u −v⟩+α 2∥u−v∥2,∀u, v ∈V, de donde tomando v=−γJ′(un), f(γ)≥J(un)−γ∥J′(un)∥2+α 2γ2∥J′(un)∥2=J(un) + α 2γ2−γ∥J′(un)∥2, donde tomando límite cuando |γ| → ∞, se llega a que fes coerciva en el infinito. Como fes estrictamente convexa y coercitiva, se deduce que posee un único mínimo γnen R f(γn) = ´ınf γ∈Rf(γ) por lo que: J(un−γnJ′(un)) = ´ınf γ∈RJ(un−γJ′(un)). y por tanto γnes el paso óptimo en nuestro algoritmo. En consecuencia, un+1 =un−γnJ′(un). Por ser un mínimo, se debe verificar f′(γn) = 0. Entonces: 0 = f′(γn) = −J′(un)′J′(un−γnJ′(un)) = −⟨J′(un+1), J′(un)⟩ y por tanto deducimos que dos direcciones de descenso consecutivas son ortogonales. Dado un+1 =un−γnJ′(un), se tiene que un+1 −un=γnJ′(un)y por tanto ⟨J′(un+1), un+1 −un⟩= 0. De nuevo, utilizando que Jes convexo: J(un)≥J(un+1) + ⟨J′(un+1), un+1 −un⟩+α 2∥un−un+1∥2 =J(un+1) + α 2∥un−un+1∥2, (6.1) de donde deducimos que la sucesión {J(un)}n≥0es monótona decreciente. Como {J(un)}n≥0es monótona decreciente y está acotada inferiormente por J(u), es convergente y de (6.1) se deduce que ∥un−un+1∥ → 0. De nuevo, por ser {J(un)}n≥0 50 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa monótona decreciente y Jcoerciva en el infinito, deducimos que la sucesión {un}n≥0 está acotada. Es decir, existe M > 0tal que ∥un∥ ≤ M, ∀n≥0.(6.2) Nos queda probar la convergencia de la sucesión. Debido a la ortogonalidad de dos direcciones consecutivas, tenemos: ∥J′(un)∥2=⟨J′(un), J′(un)⟩=⟨J′(un), J′(un)−J′(un+1)⟩≤∥J′(un)∥∥J′(un)−J′(un+1)∥, y si simplificamos, ∥J′(un)∥≤∥J′(un)−J′(un+1)∥. Si ahora utilizamos (6.2) junto con la hipótesis de que Jes globalmente lipschitziana, obtenemos lo siguiente: ∥J′(un)∥≤∥J′(un)−J′(un+1)∥ ≤ LM∥un−un+1∥, donde LMes la constante de Lipschitz de J′. Tomando límites en esta expresión y teniendo en cuenta que ∥un−un+1∥ → 0, se llega a que ∥J′(un)∥ → 0. Por último, utilizando de nuevo la Proposición 4.1 ya que Jes α-convexo, α∥un−u∥2≤ ⟨J′(un)−J′(u), un−u⟩=⟨J′(un), un−u⟩ ≤ ∥J′(un)∥∥un−u∥, lo que nos lleva a α∥un−u∥≤∥J′(un)∥=⇒l´ım n→∞ un=u. Lema 6.1.El método del gradiente converge para el funcional de energía sobre el que hemos estado trabajando. Demostración. Recordamos que nuestro funcional es J(v) = 1 2∥∇v∥2 L2(Ω) +ZΩ B(v)−ZΩ fv. Para ver la convergencia del algoritmo, tenemos que probar que se verifican las hipótesis del Teorema 6.1. En particular, tenemos que ver que J′es globalmente Lipschitziana en los acotados de V. 6. métodos numéricos para la resolución del problema de optimización 51 Basta aplicar el Teorema del Valor Medio en un espacio de Banach a J′y llegamos a que para cierto ψque se encuentra entre uyv, se tiene que: ∥J′(u)−J′(v)∥≤∥J′′(ψ)∥∥u−v∥ donde usando la desigualdad de Holder, así como teniendo en cuenta las inyecciones de Sobolev, se llega a que ∥J′′(ψ)∥ ≤ LM, con LM∈R, concluyéndose así la prueba. 6.2 Método del gradiente conjugado El método del gradiente, así como otras técnicas de descenso, utilizan en cada iteración del algoritmo una dirección localmente óptima en relación con el funcional Jy su derivada. Sin embargo, es posible encontrar direcciones de descenso mejoradas con respecto a este último algoritmo. Principalmente, hemos seguido [6], [8], [12], [4] y [14] para el desarrollo de la sección y sus subsecciones. 6.2.1 Método del gradiente conjugado para funcionales cuadráticos En el caso de tener un funcional J:RN−→ Rque sea cuadrático, el método del gradiente conjugado es uno de los métodos iterativos más eficientes para la resolución de sistemas lineales, que se caracteriza por converger en un número finito de etapas (menor o igual a k, la dimensión del espacio de dimensión finita en el que se busca la solución). Se trata de un método iterativo que es directo. Con este algoritmo, se minimiza Jen cada etapa en la expresión un+L{g1, ..., gn}, donde gi=J′(ui), con i= 1, ..., N. El nombre de conjugado se debe a la propiedad de conjugación que poseen las direcciones dk, que forman un sistema ortogonal y por tanto (di, dj)V= 0 para todo i=j. De manera general, los métodos de direcciones conjugadas minimizan un funcional Jpartiendo de un punto inicial, siguiendo de manera sucesiva una cantidad finita de direcciones linealmente independientes entre sí. Estas direcciones poseen la propiedad de conjugación con respecto a J. Más concretamente, si J:RN−→ Res un funcional cuadrático, J(v) = 1 2vtBv −ctv, con B∈ MN×N, c ∈RN 58 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa El siguiente Teorema, extraído de [4], proporciona una serie de condiciones suficientes que garantizan la existencia de un cero de Fen un entorno del punto inicial x0, así como la convergencia del algoritmo hacia este cero. Teorema 6.3. Sea Bun espacio de Banach y sea F: Ω ⊂B−→ Bun funcional F-derivable en Ω, un abierto no vacío. Supongamos que existen tres constantes δ, M, β tales que: δ > 0yB=B(x0, δ)⊂Ω, sup n≥0 sup x∈B ∥A−1 n(x)∥ ≤ M, (6.6) sup n≥0 sup x′∈B ∥F′(x)−An(x′)∥ ≤ β M, β < 1,(6.7) ∥F(x0)∥<δ M(1 −β).(6.8) Entonces, la sucesión (xn)n≥0que genera el método de Newton está contenida en la bola B y dicho método converge a un cero de F, que es el único contenido en dicha bola. Además, la convergencia es geométrica: ∥xn−x∗∥ ≤ βn 1−β∥x1−x0∥. Demostración. En primer lugar, vampos a demostrar que para todo n≥1, se tiene que: ∥xn−xn−1∥ ≤ M∥F(xn−1)∥, ∥xn−x0∥ ≤ δ , ∀xn∈B, ∥F(xn)∥ ≤ β M∥xn−xn−1∥. Vamos a aplicar para ello un razonamiento de inducción en n. Probemos, en primer lugar, que para n= 1 se verifican las desigualdades. De x1−x0=−A−1 0(x0)F(x0), deducimos de (6.6)y(6.8) tomando normas que: ∥x1−x0∥ ≤ M∥F(x0)∥ ≤ δ(1 −β)≤1. Además, podemos escribir F(x1) = F(x1)−F(x0)−A0(x0)(x1−x0). 6. métodos numéricos para la resolución del problema de optimización 59 Si ahora aplicamos el Teorema del Valor Medio a la función g(x) = F(x)−A0(x0)x, F(x1)≤sup x∈B ∥F′(x)−A0(x0)∥∥x1−x0∥ ≤ β M∥x1−x0∥,por (6.7). Aplicamos ahora inducción. Supongamos que las desigualdades vistas se cumplen para n−1, y veamos que son ciertas para n. Como xn−xn−1=−A−1 n−1(xn−1)F(xn−1), de (6.6) deducimos que ∥xn−xn−1∥ ≤ M∥F(xn−1)∥, y así queda probada la primera desigualdad para todo n≥1. De este modo y aplicando la hipótesis de inducción: ∥xn−xn−1∥ ≤ β∥xn−1−xn−2∥ ≤ ... ≤βk−1∥x1−x0∥, y así, usando (6.8) ∥xn−x0∥ ≤ n X j=1 ∥xj−xj−1∥ ≤ n X j=1 βj−1!∥x1−x0∥ ≤ ≤∥x1−x0∥ 1−β≤M 1−β∥F(x0)∥ ≤ r, y así se prueba que xn∈B. Por último, sabemos que F(xn) = F(xn)−F(xn−1)−An−1(xn−1)(xn−xn−1). Por lo tanto, teniendo en cuenta (6.7) y nuevamente el Teorema del Valor Medio para la función ˆg(x) = F(x)−An−1(xn−1)x, se llega a ∥F(xn)∥ ≤ sup x∈B ∥F′(x)−An−1(xn−1)∥∥xn−xn−1∥ ≤ β M∥xn−xn−1∥, y así queda probada la última desigualdad para todo n≥1. Demostremos ahora que existe un cero del funcional Fen la bola B. Como ∥xn+l−xn∥ ≤ l−1 X j=0 ∥xn+j+1 −xn+j∥ ≤ βn l−1 X j=0 βj∥x1−x0∥ ≤ βk 1−β∥x1−x0∥(6.9) 60 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa entonces para todo n≥0yl≥0, la sucesión (xn)n≥0es de Cauchy. Como se trata de una sucesión contenida en una bola cerrada de un espacio de Banach, entonces existe un punto x∗tal que l´ım n→∞ xn=x∗. La aplicación Fes continua en Ωpor ser F-derivable. En consecuencia, se tiene ∥F(x∗)∥= l´ım n→∞ ∥F(xn)∥ ≤ β Ml´ım n→∞ ∥xn−xn−1∥= 0, y por tanto F(x∗) = 0, como queríamos probar. Además, si ahora hacemos l→ ∞ en la expresión (6.9), obtenemos la desigualdad siguiente: ∥xn−x∗∥ ≤ βn 1−β∥x1−x0∥. Por último, queda probar que x∗es el único cero del funcional en la bola B. Por reducción al absurdo, supongamos que existe ˆxotro cero de F. Como F(x∗) = F(ˆx) = 0, podemos escribir x∗−ˆx=−A−1 0(x0)(F(ˆx)−F(x∗)−A0(x0)(ˆx−x∗)), de donde se deduce inmediatamente que ∥x∗−ˆx∥≤∥A−1 0(x0)∥sup x∈B ∥F′(x)−A0(x0)∥∥x∗−ˆx∥ ≤ β∥x∗−ˆx∥, y como β < 1, necesariamente x∗= ˆx, quedando así probada la unicidad. Observación 6.3.En general, el método de Newton plantea varios inconvenientes que dificultan su uso práctico. En particular, cuando tomamos como punto inicial del algoritmo un punto alejado de la solución, la convergencia no queda asegurada. Además, los elementos ∆nyAnpueden ser difíciles y/o costosos de calcular para valores de nmuy altos. En nuestro caso, la solución de (PO)va a estar caracterizada por J′(u) = 0, es decir, tendremos F=J′. De esta forma, F′(un) = J′′(un): J′′(un)(v, w) = ZΩ ∇v∇w+ZΩ β′(un)v w . Supongamos conocida un∈V. Queremos entonces calcular ∆n∈V, como indica el segundo paso del algoritmo, a partir de la solución de: ⟨J′(un), v⟩+J′′(un)(v, ∆n) = 0 6. métodos numéricos para la resolución del problema de optimización 61 ⇒ZΩ ∇un∇v+ZΩ β(un)v+ZΩ f v +ZΩ ∇v∇(∆n) + ZΩ β′(un)v∆n= 0 , para todo vde H1 0(Ω). Esta expresión admite una nueva interpretación, pues se puede entender la ecuación en un+1 como: ZΩ ∇un+1 ∇v+ZΩ (β(un) + β′(un)(un+1 −un)) v=ZΩ fv. (6.10) Nuestra idea desde un principio es resolver el problema ZΩ ∇u∇v+ZΩ β(u)v=ZΩ fv , por lo que (6.10) es similar a nuestro problema, solo que linealizando el factor no lineal β. En el próximo capítulo, veremos una implementación de este método para un problema particular. 7 Aplicación numérica: Resolución de una EDP elíptica no lineal En este último capítulo, vamos a aplicar los principales resultados que hemos ido demostrando a lo largo de este trabajo acerca del funcional de energía que queremos minimizar. Se recuerda en primer lugar que tenemos el funcional J:H1 0(Ω) →R dado por: J(v) = 1 2ZΩ |∇v|2+ZΩ B(v)−ZΩ fv. En el Capítulo 3hemos demostrado que Jes G-derivable y F-derivable. En el Capítulo 4hemos probado que es un funcional estrictamente convexo y en particular α-convexo. En ese mismo capítulo, también se prueba la coercividad en el infinito del funcional. En consecuencia, gracias a los Teoremas 4.1 y4.2, el siguiente problema de optimización (PO)Obtener u=arg m´ın v∈UJ(v) admite una única solución. En particular, (P O)es el problema de optimización asociada a la formulación variacional que se obtiene del siguiente problema de contorno, basado en una EDP no lineal con no linealidad suave: (P)−∆u+β(u) = f, en Ω u= 0 , sobre ∂Ω Ahora sí, para aproximar númericamente la solución del problema podemos utilizar algunos de los métodos que hemos visto en el Capítulo 6. 64 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa 7.1 Resolución con FreeFem++ Vamos a aproximar en particular el siguiente problema de contorno: Consideramos el dominio Ω = [0,1] ×[0,1] ⊂R2. Sea f= 40 y sea β(u) = u3. Queremos hallar u∈H1 0(Ω) tal que: (P)−∆u+u3= 40 , en Ω u= 0 , sobre ∂Ω Para ello, vamos a resolver el problema (PO)utilizando el método de Newton y el código que se encuentra en el Apéndice A. Para hayar el mínimo del funcional, se consideran los siguientes datos y notación: epsmn = 10−8como valor para el test de parada del método de Newton. i= 100, número máximo de iteraciones para el algoritmo. noruz, norma de la diferencia de las iteraciones, es decir, ∥un+1 −un∥ noru, norma de un+1 , es decir, ∥un+1∥ err equivale al error relativo. Como test de parada para el algoritmo hemos elegido noruz noru < epsmn, es decir, ∥un+1 −un∥ ∥un+1∥< epsmn = 10−8. Veamos en primer lugar el mallado de nuestro dominio: 7. aplicación numérica: resolución de una edp elíptica no lineal 65 Figura 7.1: Mallado del dominio En las figuras que se presentan a continuación, se pueden observar las sucesivas iteraciones del algoritmo. En particular, han sido necesarias 6 iteraciones del método para que se verifique nuestro test de parada. (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.2: Primera iteración del Método de Newton 66 aproximación numérica de ecuaciones en derivadas parciales elípticas formulables como problemas de optimización convexa (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.3: Segunda iteración del Método de Newton (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.4: Tercera iteración del Método de Newton (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.5: Cuarta iteración del Método de Newton 7. aplicación numérica: resolución de una edp elíptica no lineal 67 (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.6: Quinta iteración del Método de Newton (a) Solución de problema (b) Aplicación de beta en el resultado Figura 7.7: Sexta iteración del Método de Newton De esta forma, en la sexta iteración se llega a que ∥un+1 −un∥= 1.06358 ·10−16 y∥un+1∥= 40.0768, con lo cual se obtiene un error relativo de ε= 1.62907e·10−9. Entonces, se tiene que ε < epsmn , con lo cual se verifica nuestra regla de parada y podemos asegurar la convergencia del método que se ha planteado.