Método POD y aplicaciones : caso discreto y continuo

universidad de sevilla Trabajo Fin de M´aster M´etodo POD y aplicaciones: caso discreto y continuo Presentado por: Dirigido por: Manuel Bernardino del Pino Macarena G´omez M´armol Samuele Rubino Departamento de Ecuaciones Diferenciales y An´ alisis Num´ erico MARZO DE 2021 A mi abuelo... Abstract The proper orthogonal decomposition (POD) is a technique that aims at obtaining a compact representation of the data reducing the computational cost. The key idea of the POD is to reduce a large number of interdependent variables to a much smaller number of uncorrelated variables while retaining as much as possible of the variation in the original variables. This method has been gaining popularity over the years and it is being used to detect the number of signals in a multichannel time-series, for forecasting in meteorology or for the design of gurney flaps, among others. In this work, a theoretical study is carried out on this method in both the continuous and the discrete cases, the effective calculation is exposed and two numerical simulations are made; one for each of the two cases mentioned above. ´ Indice general Introducci´on 3 1. Resultados preliminares 5 1.1. Optimizaci´on no lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. Espacios de Hilbert. Teor´ıa espectral . . . . . . . . . . . . . . . . . . . 8 2. El m´etodo POD en el espacio Rm12 2.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2. La descomposici´on en valores singulares (SVD) . . . . . . . . . . . . . 13 2.3. Relaci´on entre el m´etodo POD y la descomposici´on en valores singulares(SVD).................................. 14 2.4. Propiedades de la base POD . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5. El m´etodo POD para productos escalares con peso . . . . . . . . . . . 27 2.6. Algunas aplicaciones del m´etodo POD . . . . . . . . . . . . . . . . . . 34 2.6.1. Ecuaciones diferenciales ordinarias . . . . . . . . . . . . . . . . 34 2.6.2. Ecuaciones en derivadas parciales . . . . . . . . . . . . . . . . . 47 3. Aplicaciones pr´acticas del m´etodo POD 50 3.1. Casodiscreto ................................ 50 3.2. Casocontinuo................................ 53 Bibliograf´ıa 59 1 Introducci´on El m´etodo proper orthogonal descomposition (POD) es una t´ecnica que tiene como objetivo obtener una representaci´on compacta de una serie de datos reduciendo el coste computacional. La idea clave de este m´etodo es reducir un gran n´umero de variables interdependientes a un n´umero mucho menor de variables no correlacionadas mientras se retiene la “informaci´on esencial” de las variables originales. Se realiza una transformaci´on ortogonal a la base de los autovectores de una cierta matriz, y los datos se proyectan en el subespacio generado por los autovectores correspondientes a los autovalores m´as grandes. Esencialmente, proporciona una base ortogonal para representar un conjunto dado de datos en un cierto sentido ´optimo de m´ınimos cuadrados, es decir, ofrece formas de encontrar aproximaciones ´optimas de dimensiones inferiores para ese conjunto de datos. Este m´etodo, tambi´en conocido como la descomposici´on de Karhunen-Lo`eve (KLD), fue propuesto de forma independiente por varios cient´ıficos, entre los que destacan Karhunen, Kosambi, Lo`eve, Obukhov y Pougachev y fue concebido originalmente en el marco de procesos continuos de segundo orden. Cuando se restringe a un caso de dimensi´on finita y se trunca despu´es de algunos t´erminos, el POD es equivalente al an´alisis de componentes principales (ACP). Esta ´ultima metodolog´ıa se origin´o con el trabajo de Pearson como un medio de ajustar planos por m´ınimos cuadrados ortogonales, pero tambi´en fue propuesta por Hotelling. Con esto quiero enfatizar que en la literatura el m´etodo POD aparece con distintos nombres dependiendo del ´area de aplicaci´on, a saber, ACP en la literatura estad´ıstica, funci´on ortogonal emp´ırica en oceanograf´ıa y meteorolog´ıa y an´alisis factorial en psicolog´ıa y econom´ıa. En las referencias [8], [9] y [14] pueden encontrarse una discusi´on detallada sobre la equivalencia de POD, ACP y su conexi´on con la descomposici´on de valores singulares. 3 10 Manuel Bernardino del Pino Teorema 1.2.7. Un operador T∈ L(H1,H2) es compacto si y solo si su operador adjunto T∗∈ L(H2,H1) lo es. Teorema 1.2.8. Sea {Tn}una sucesi´on de operadores compactos en L(H1,H2) y T∈ L(H1,H2) tal que kTn−Tkn→∞ −→ 0.Entonces Tes un operador compacto. Teorema 1.2.9. Sean T:H → H un operador compacto y S∈ L(H). Entonces TS yST son compactos. Teorema 1.2.10. Sea T∈ L(H) un operador de rango finito, entonces Tes compacto. Se vio en Algebra Lineal que si Hes de dimensi´on finita, entonces todo operador sim´etrico es diagonalizable, es decir, tiene una base de autovectores: Teorema 1.2.11. Sea Hes un espacio eucl´ıdeo de dimensi´on nyT∈ L(H) un operador autoadjunto, entonces existe una base ortonormal {ψ1, . . . , ψn}de Hy unos n´umeros reales λ1, . . . , λntales que Tψi=λiψi,∀i= 1, . . . , n. El objetivo ahora ser´a generalizar este resultado al caso de un operador autoadjunto sobre un espacio de Hilbert que puede ser infinito dimensional. Definici´on 1.2.12. Sea T∈ L(H). Diremos que un escalar λ∈Res un autovalor si existe x6= 0 en Htal que Tx =λx. Al vector xse le llama autovector asociado al autovalor λ. Es claro que λ∈Res un autovalor de Tsi y solo si Ker (T−λId)6={0}, es decir, si y solo si T−λId no es inyectivo. Los autovalores de un operador autoadjunto verfican las siguientes propiedades: Proposici´on 1.2.13. Sea T∈ L(H) un operador autoadjunto, entonces sus autovalores (en caso de tenerlos) son reales. Proposici´on 1.2.14. Los autovalores asociados a distintos autovalores de un operador autoadjunto son ortogonales. Vamos a terminar con un resultado clave en la teor´ıa espectral: todo operador compacto y autoadjunto es diagonalizable. RESULTADOS PRELIMINARES 11 Teorema 1.2.15. Sea T∈ L(H) un operador compacto y autoadjunto. Entonces, existe un sistema ortonormal {xn}de autovectores de Ty su correspondiente sucesi´on de autovalores {λn}tal que para cada x∈ H se tiene que Tx =X n λnhx, xnixn. La sucesi´on {λn}es decreciente y si es infinita converge a 0. Cap´ıtulo 2 El m´etodo POD en el espacio Rm 2.1. Planteamiento del problema El m´etodo proper orthogonal descomposition, a partir de ahora, m´etodo POD, es una t´ecnica que permite calcular una aproximaci´on de unos vectores y1, . . . , yn∈Rm utilizando una serie de datos experimentales de los que se extrae su “informaci´on esencial” en t´erminos de unos pocos vectores, que llamaremos base POD. Es decir, se trata de sustituir las variables originales por un n´umero menor de variables que sean una combinaci´on lineal de las originales. Supongamos que y1, . . . , yn∈Rmson los vectores dados que queremos aproximar, entonces el problema de aproximar simult´aneamente estos vectores por `(< n) vectores normalizados y ortogonales dos a dos, puede ser planteado como el siguiente problema de maximizaci´on: m´ax ˜ ψ1,..., ˜ ψ`∈Rm ` X i=1 n X j=1 |hyj,˜ ψii|2sujeto a h˜ ψi,˜ ψji=δij para 1 ≤i, j ≤`, donde δij es la funci´on delta de Kronecker, es decir, toma el valor uno si i=jy cero en caso contrario. El objetivo, por tanto, ser´a resolver este problema de optimizaci´on no lineal, obtener algunas propiedades de la base POD y aplicar este m´etodo a alguna situaci´on de inter´es. Para ello, la descomposici´on en valores singulares ser´a una herramienta fundamental. Para el desarrollo de este cap´ıtulo hemos seguido la referencia de [12]. 12 EL M´ ETODO POD EN EL ESPACIO RM13 2.2. La descomposici´on en valores singulares (SVD) Sea Y= (y1|···|yn)∈Rm×nuna matriz de rango d≤m´ın{m, n}con las columnas yj∈Rmpara todo j= 1,...n. Recordemos que la descomposici´on en valores singulares, en adelante, SVD, de Y garantiza que existen unos n´umeros reales σ1≥σ2≥ ··· ≥ σd>0, una matriz ortogonal Ψ ∈Rm×mcon columnas {ψi}m i=1 y Φ ∈Rn×nortogonal con columnas {φi}n i=1 tales que: Ψ>YΦ =  DΘ Θ Θ  := Σ ∈Rm×n,(2.1) donde D=diag(σ1, . . . , σd)∈Rd×d, los ceros Θ denotan matrices nulas de dimensiones adecuadas y ‘>’ la traspuesta de una matriz (o vector) (ver [5], p´ag 76). Adem´as, se puede probar f´acilmente a partir de la definici´on (2.1) que los vectores {φi}d i=1 y{ψi}d i=1 verifican que Y φi=σiψiyY>ψi=σiφipara todo i= 1, . . . d. (2.2) Adem´as, utilizando estas igualdades obtenemos: Y>Y φi=σiY>ψi=σ2 iφipara todo i= 1, . . . d. Y Y >ψi=σiY φi=σ2 iψipara todo i= 1, . . . d. Por tanto, {φi}d i=1 y{ψi}d i=1 son el conjunto de autovectores de la matriz Y>Y eY Y >, respectivamente, asociado a los autovalores λi=σ2 i>0. Los vectores {φi}n i=d+1 ( si d < n) y {ψi}m i=d+1 ( si d < m) son autovectores de Y>YeY Y >asociados al autovalor nulo. Por otra parte, de (2.1) y de las propiedades de Ψ y Φ, se obtiene que Ypuede ser descompuesta como Y= ΨΣΦ>. Es m´as, Ypuede ser escrita como sigue: Y= ΨdD(Φd)>,(2.3) 14 Manuel Bernardino del Pino donde las matrices Ψd∈Rm×dy Φd∈Rn×dvienen dadas por Ψd ij = Ψij para 1 ≤i≤m; 1 ≤j≤d, Φd ij = Φij para 1 ≤i≤n; 1 ≤j≤d. Definiendo Bd:= D(Φd)>∈Rd×n, podemos reescribir (2.3) de la forma Y= ΨdBd. Por tanto, el espacio columna de Ypuede ser representado en t´erminos de las d columnas independientes de Ψd. Los coeficientes de la expansi´on para las columnas yj (j= 1, . . . , n) en la base {ψi}d i=1 vienen dados por la j-´esima columna de Bd. Adem´as, como Ψ es ortogonal, se tiene: yj= d X i=1 Bd ijΨd ·,i = d X i=1 (D(Φd)>)ijψi= d X i=1 ((Ψd)>Ψd | {z } =Id D(Φd)>)ijψi = d X i=1 ((Ψd)>Y)ijψi= d X i=1 m X k=1 Ψd kiYkj | {z } =ψ> iyj ψi= d X i=1hψi, yjiRmψi, donde Ides la matriz identidad de orden d, Ψd ·,i denota la i-´esima columna de la matriz Ψdyh·,·iRmrepresenta el producto escalar eucl´ıdeo en Rmque denotamos a partir de ahora simplemente por h·,·i para simplificar la notaci´on. Luego, hemos escrito los vectores columna de Ycomo combinaci´on lineal de {ψi}d i: yj= d X i=1hyj, ψiiψipara todo j= 1, . . . n. 2.3. Relaci´on entre el m´etodo POD y la descomposici´on en valores singulares (SVD) Supongamos que y1, . . . , yn∈Rmson los vectores dados que queremos aproximar. Para comenzar, vamos a considerar el problema de aproximar simult´aneamente todos estos vectores yjpor uno s´olo normalizado: (P1) m´ax ˜ ψ∈Rm n X j=1 |hyj,˜ ψi|2sujeto a k˜ ψk2= 1, EL M´ ETODO POD EN EL ESPACIO RM15 donde k˜ ψk=qh˜ ψ, ˜ ψipara ˜ ψ∈Rm. Notemos que (P1) es un problema de optimizaci´on no lineal con restricci´on de igualdad que puede ser resuelto considerando las condiciones necesarias de optimalidad de primer orden. Para ello, vamos a escribir nuestro problema en la forma est´andar definiendo la funci´on e:Rm→Rpor e(ψ) = 1 − kψk2para ψ∈Rm. Entonces, la restricci´on de (P1) puede ser expresada como e(ψ) = 0. Observemos que cualquier soluci´on ψde (P1) es no nula, ya que tiene que verificar kψk= 1. Por tanto, para cualquier soluci´on ψdel problema, el conjunto {∇e(ψ)}={−2ψ>}es linealmente independiente y por ende cualquier soluci´on de (P1) es un punto regular. Sea L:Rm×R→Rel funcional Lagrangiano asociado a (P1), es decir, L(ψ, λ) = n X j=1 |hyj, ψi|2+λ(1 −kψk2) para (ψ, λ)∈Rm×R. Supongamos que ψ∈Rmes una soluci´on de (P1). Como ψes un punto regular, se deduce del Teorema 1.1.4 que existe un ´unico multiplicador de Lagrange λ∈R verificando la condici´on necesaria de optimalidad de primer orden ∇L(ψ, λ) = 0 en Rm×R. Vamos a empezar calculando el gradiente de Lcon respecto a ψ: ∂L ∂ψi (ψ, λ) = ∂ ∂ψin X j=1  m X k=1 Ykjψk 2 +λ1− m X k=1 ψ2 k = 2 n X j=1 m X k=1 YkjψkYij −2λψi = 2 m X k=1 n X j=1 YijY> jk | {z } =(Y Y >)ik ψk−2λψi. Por tanto, ∇ψL(ψ, λ) = 2(Y Y >ψ−λψ) = 0 si y solo si Y Y >ψ=λψ. Ahora bien, como la matriz Y Y >es sim´etrica y semidefinida positiva, existen mautovalores λ1≥ λ2≥. . . ≥λm≥0 y autovectores asociados, que denotamos por {ψi}m i=1, que pueden ser elegidos de forma que el conjunto {ψ1, . . . , ψm}sea ortonormal. Es decir, se verifica 16 Manuel Bernardino del Pino que Y Y >ψi=λiψi(1 ≤i≤m) y adem´as, si ψ∈Rmes un vector arbitrario, entonces ψ= m X i=1hψ, ψiiψi. Por otra parte, la condici´on ∂L ∂λ (ψ, λ) = 0 nos conduce a kψk= 1. Es decir, la condici´on necesaria de optimalidad nos lleva a resolver el problema:    Y Y >ψ=λψ en Rm, kψk= 1 en R.(2.4) Debido a la SVD, en particular, el vector ψ1resuelve (2.4) y se tiene: n X j=1 |hyj, ψ1i|2= n X j=1hyj, ψ1ihyj, ψ1i= n X j=1hhyj, ψ1iyj, ψ1i=*n X j=1hyj, ψ1iyj, ψ1+ =*n X j=1 m X k=1 Ykj(ψ1)k!yj, ψ1+=*m X k=1   n X j=1 Y·jY> kj (ψ1)k , ψ1+ =DY Y >ψ1, ψ1E=λ1hψ1, ψ1i=λ1kψ1k2 |{z} =1 =λ1, donde (ψ1)kes la k-´esima componente del vector ψ1eY·jdenota la j-´esima columna de la matriz Y. Observemos que ∇ψψ L(ψ, λ) = 2(Y Y >−λIm)∈Rm×m. Si evaluamos en (ψ1, λ1) se concluye de λ1≥λ2≥. . . ≥λm≥0 que hψ, ∇ψψL(ψ1, λ1)ψi= 2 Dψ, (Y Y >−λ1Im)ψE = 2 m X i=1 m X j=1hψ, ψiihψ, ψjiDψ1,(Y Y >−λ1Im)ψjE = 2 m X i=1 m X j=1 (λj−λ1)hψ, ψiihψ, ψjihψi, ψji = 2 m X i=1 (λi−λ1)|hψ, ψii|2≤0. Por tanto, (ψ1, λ1) verifican las condiciones necesarias de optimalidad de segundo orden para un m´aximo, pero no las suficientes (Teoremas 1.1.6 y 1.1.7). Probemos a continuaci´on que en realidad ψ1resuelve el problema (P1). En efecto, supongamos que EL M´ ETODO POD EN EL ESPACIO RM17 ˜ ψ∈Rmes un vector arbitrario con k˜ ψk= 1. Como {ψi}m i=1 es una base ortonormal de Rm, tenemos: ˜ ψ= m X i=1h˜ ψ, ψiiψi. Por tanto: n X j=1 |hyj,˜ ψi|2= n X j=1 *yj, m X i=1h˜ ψ, ψiiψi+ 2 = n X j=1 m X i=1 m X k=1 Dyj,h˜ ψ, ψiiψiEDyj,h˜ ψ, ψkiψkE = n X j=1 m X i=1 m X k=1 hyj, ψiihyj, ψkih˜ ψ, ψiih˜ ψ, ψki = m X i=1 m X k=1  n X j=1hyj, ψiiyj | {z } =λiψi , ψkh˜ ψ, ψiih˜ ψ, ψki = m X i=1 m X k=1 hλiψi, ψki | {z } =λiδik h˜ ψ, ψiih˜ ψ, ψki = m X i=1 λi|h˜ ψ, ψii|2≤λ1 m X i=1 |h˜ ψ, ψii|2=λ1k˜ ψk2=λ1 = n X j=1 |hyj, ψ1i|2. Luego, ψ1resuelve (P1) y arg m´ax(P1) = σ2 1=λ1. Si ahora buscamos un segundo vector, ortogonal a ψ1que describa nuevamente el conjunto de datos {yi}n i=1 con el menor error posible, entonces tenemos que resolver el siguiente problema: (P2) m´ax ˜ ψ∈Rm n X j=1 |hyj,˜ ψi|2sujeto a k˜ ψk2= 1 y h˜ ψ, ψ1i= 0. La SVD implica que ψ2es una soluci´on de (P2) y arg m´ax(P2) = σ2 2=λ2. De hecho, ψ2resuelve las condiciones necesarias de optimalidad de primer orden y para ˜ ψ= m X i=2h˜ ψ, ψiiψi∈ hψ1i⊥ 18 Manuel Bernardino del Pino tenemos n X j=1 |hyj,˜ ψi|2≤λ2= n X j=1 |hyj, ψ2i|2. Claramente, este procedimiento se puede continuar por inducci´on finita. Vamos a resumir los resultados obtenidos con el siguiente teorema. Teorema 2.3.1. Sea Y= (y1|···|yn)∈Rm×nuna matriz de rango d≤m´ın{m, n} y denotemos por Y= ΨΣΦ>la descomposici´on en valores singulares de Y, donde Ψ = (ψ1|···|ψm)∈Rm×m,Φ = (φ1|···|φn)∈Rn×nson matrices ortogonales y la matriz Σ ∈Rm×nes: Σ =  DΘ Θ Θ  , donde D=diag(σ1, . . . , σd)∈Rd×d. Entonces, para cualquier `∈ {1, . . . , d}, la soluci´on de (P`) m´ax ˜ ψ1,..., ˜ ψ`∈Rm ` X i=1 n X j=1 |hyj,˜ ψii|2sujeto a h˜ ψi,˜ ψji=δij para 1 ≤i, j ≤` viene dada por los vectores {ψi}` i=1, es decir, por las primeras `columnas independientes de Ψ, donde δij es la delta de Kronecker. Adem´as, arg m´ax (P`) = ` X i=1 σ2 i= ` X i=1 λi. Antes de demostrar el anterior resultado, vamos a dar la siguiente definici´on: Definici´on 2.3.2. Para `∈ {1,...d}, los vectores {ψi}` i=1 se llaman base POD de rango `. Demostraci´on del Teorema. Como (P`) es un problema de optimizaci´on con restricci´on de igualdad vamos a definir la funci´on Lagrangiana L:Rm×. . . ×Rm | {z } `-veces ×R`×` por L(ψ1, . . . , ψ`,Λ) = ` X i=1 n X j=1 |hyj, ψii|2+ ` X i,j=1 λij (δij −hψi, ψji) EL M´ ETODO POD EN EL ESPACIO RM19 para ψ1, . . . , ψ`∈Rmy Λ = ((λij)) ∈R`×`.Las condiciones necesarias de optimalidad de primer orden para P`vienen dadas por ∂L ∂ψk (ψ1, . . . , ψ`,Λ) z= 0 para todo z∈Rmyk∈ {1, . . . , `}.(2.5) Ahora bien, ∂L ∂ψk (ψ1, . . . , ψ`,Λ) z= 2 ` X i=1 n X j=1 hyj, ψiihyj, ziδik − ` X i,j=1 λij hψi, ziδjk − ` X i,j=1 λij hz, ψjiδki = 2 n X j=1 hyj, ψkihyj, zi− ` X i=1 (λik +λki)hψi, zi =2 n X j=1 hyj, ψkiyj− ` X i=1 (λik +λki)ψi, z Luego (2.5), es equivalente a: n X j=1 hyj, ψkiyj=1 2 ` X i=1 (λik +λki)ψien Rmpara todo k∈ {1, . . . , `}.(2.6) Ahora bien, como Y Y >ψ= n X j=1 hyj, ψiyjpara ψ∈Rm, la condici´on (2.6) puede ser escrita de la forma: Y Y >ψk=1 2 ` X i=1 (λik +λki)ψien Rmpara todo k∈ {1, . . . , `}.(2.7) Ahora procedemos por inducci´on. Para `= 1 se escribe como Y Y >ψ1=λ1ψ1en Rm con λ1=λ11.Ahora supongamos que para `≥1 las condiciones de optimalidad de primer orden vienen dadas por Y Y >ψk=λkψken Rmpara todo k∈ {1, . . . , `}.(2.8) 26 Manuel Bernardino del Pino Demostraci´on. La prueba es una consecuencia inmediata de (2.14): n X j=1hyj, ψiihyj, ψki=n X j=1hyj, ψiiyj | {z } =Y Y >ψi , ψk=hσ2 iψi, ψki=σ2 iδik. C´alculo efectivo de la base POD Vamos a terminar con una observaci´on sobre la computaci´on pr´actica de una base POD de rango `. Si n<m, entonces se puede determinar la base POD de rango `de la siguiente forma: Se calculan los autovectores φ1, . . . , φ`∈Rnresolviendo el problema de autovalores n×n: Y>Y φi=λiφipara i= 1, . . . , ` (2.18) y se toma por (2.2), ψi=1 √λi Y φipara i= 1, . . . , `. (2.19) Si n>m, entonces se puede calcular la base POD resolviendo el problema de autovalores m×m: Y Y >ψi=λiψipara i= 1, . . . , `. (2.20) Por razones hist´oricas [11] este m´etodo de determinar la base POD es a menudo llamado m´etodo de los snapshots. Para la aplicaci´on de este m´etodo POD a problemas concretos, la elecci´on de `es de vital importancia. En general, no hay reglas definidas para ello. M´as bien, su elecci´on se basa en consideraciones heur´ısticas combinadas con la observaci´on de la relaci´on entre energ´ıa modelada y energ´ıa total contenida en el sistema Y, lo cual viene expresado como E(`) = ` X i=1 λi d X i=1 λi .(2.21) EL M´ ETODO POD EN EL ESPACIO RM27 Notemos que P` i=1 λi=Traza (Y Y >) =Traza (Y>Y). Mencionemos adem´as que el m´etodo POD se llama tambi´en An´alisis de Componentes Principales oDescomposici´on de Karhunen-Lo`eve. 2.5. El m´etodo POD para productos escalares con peso Dotemos a Rmdel siguiente producto escalar con peso: hψ, ˜ ψiW=ψ>W˜ ψ=hψ, W ˜ ψiRm=hWψ, ˜ ψiRmpara ψ, ˜ ψ∈Rm,(2.22) donde W∈Rm×mes una matriz sim´etrica y definida positiva. Adem´as, sea kψkW= phψ, ψiWpara ψ∈Rm, la norma inducida asociada. Observemos que para la elecci´on de W=Im,el producto escalar anteriormente definido coincide con el producto escalar eucl´ıdeo. Motivemos esta definici´on con el siguiente ejemplo. Ejemplo 2.5.1. Sea L2(0,1) el espacio de las funciones medibles de cuadrado integrable en (0,1), es decir, L2(0,1) = ϕ: (0,1) →Rmedibles : Z1 0|ϕ|2dx < +∞. Es bien conocido que L2(0,1) es un espacio de Hilbert con el producto escalar hϕ, ˜ϕiL2(0,1) =Z1 0 ϕ˜ϕ dx para ϕ, ˜ϕ∈L2(0,1) cuya norma inducida es kϕkL2(0,1) =qhϕ, ϕiL2(0,1) para ϕ∈L2(0,1). Nuestro objetivo es mostrar que el producto escalar en L2(0,1), espacio de dimensi´on infinita, se puede aproximar por un cierto producto escalar discreto en un espacio de dimensi´on finita. Para ello, consideramos una partici´on uniforme del intervalo (0,1): {0 = x1< x2<··· < xm−1< xm= 1}, donde xj= (j−1)hpara j= 1, . . . , m, h= 1/(m−1). 28 Manuel Bernardino del Pino Vamos a aproximar la integral Z1 0 ϕ(x) ˜ϕ(x)dx usando la f´ormula de cuadratura del trapecio: hϕ, ˜ϕiL2(0,1) =Z1 0 ϕ(x) ˜ϕ(x)dx = m−1 X j=1 Zxj+1 xj ϕ(x) ˜ϕ(x)dx ≈h 2ϕ(x1) ˜ϕ(x1) + h m−1 X j=2 ϕ(xj) ˜ϕ(xj) + h 2ϕ(xm) ˜ϕ(xm). Si tomamos ϕh=ϕh 1, . . . , ϕh m>∈Rmy ˜ϕh=˜ϕh 1,..., ˜ϕh m>∈Rm, donde ϕh i= ϕ(xi) y ˜ϕh i= ˜ϕ(xi) para todo j= 1, . . . , m, definimos el siguiente producto escalar discreto en L2(0,1): hϕ, ˜ϕiL2 h(0,1) := hϕh 1˜ϕh 1 2+ m−1 X j=2 ϕh j˜ϕh j+ϕh m˜ϕh m 2.(2.23) Si escribimos (2.23) en forma matricial, tenemos: hϕ, ˜ϕiL2 h(0,1) =ϕh 1. . . ϕh mW     ˜ϕh 1 . . . ˜ϕh m      =ϕhW˜ϕh, donde W:= diag(h/2, h, . . . , h, h/2) ∈Rm×m. Es decir, hemos escrito el producto escalar de la misma forma que aparec´ıa en (2.22): hϕ, ˜ϕiL2 h(0,1) =ϕh>W˜ϕh=ϕh,˜ϕhW. Ahora vamos a reemplazar (P1) por (P1 W) m´ax ˜ ψ∈Rm n X j=1 |hyj,˜ ψiW|2sujeto a k˜ ψk2 W= 1. An´alogamente a lo que hicimos en la Secci´on 2.3, vamos a tratar (P1 W) como un problema de optimizaci´on con restricci´on de igualdad. El Lagrangiano L:Rm×R→R asociado a (P1 W) viene dado por L(ψ, λ) = n X j=1 |hyj, ψiW|2+λ(1 −kψk2 W) para (ψ, λ)∈Rm×R. EL M´ ETODO POD EN EL ESPACIO RM29 Introducimos la funci´on e:Rm→Rdefinida por e(ψ) = 1 − kψk2 W= 1 −ψ>Wψ para ψ∈Rm. Entonces, la restricci´on del problema puede reescribirse como e(ψ) = 0. Observemos que cualquier soluci´on ψde (P1 W) es no nula, ya que tiene que verificar kψkW= 1. Por tanto, para cualquier soluci´on ψdel problema, el conjunto {∇e(ψ)}= {−2ψ>W}es linealmente independiente y por ende cualquier soluci´on de (P1 W) es un punto regular. En consecuencia, existe un ´unico multiplicador de Lagrange λ∈R asociado con la soluci´on ´optima ψverificando la condici´on necesaria de optimalidad de primer orden ∇L(ψ, λ) = 0 en Rm×R. Vamos a empezar calculando el gradiente de Lcon respecto a ψ. Como Wes una matriz sim´etrica, se tiene: ∂L ∂ψi (ψ, λ) = ∂ ∂ψi n X j=1  m X k=1 m X ν=1 Y> jν Wνkψk 2 +λ1− m X k=1 m X ν=1 ψνWνkψk! = 2 n X j=1 m X k=1 m X ν=1 Y> jν Wνkψk m X µ=1 Y> jµWµi −λm X ν=1 ψνWνi + m X k=1 Wikψk = 2 m X k=1 m X ν=1 m X µ=1 Wiµ n X j=1 YµjY> jν Wνkψk−2λm X k=1 Wikψk = 2 WY Y >Wψ −λWψi. Por tanto: ∇ψL(ψ, λ)=2WY Y >W ψ −λWψ. Si imponemos que ∇ψL(ψ, λ) = 0, entonces tenemos que resolver el siguiente problema de autovalores generalizado (WY )(WY )>ψ=λW ψ. (2.24) Como Wes sim´etrica y definida positiva, posee una descomposici´on de la forma W= QDQ>,donde D= diag (η1, . . . , ηm) que contiene los autovalores η1≥. . . ≥ηm>0 de WyQ∈Rm×mes una matriz ortogonal. Definamos: Wα:= Qdiag (ηα 1, . . . , ηα m)Q>para α∈R. 30 Manuel Bernardino del Pino Es f´acil comprobar las siguientes propiedades de la matriz Wα. (Wα)−1=W−αpara cualquier α∈R. Wα+β=WαWβpara cualesquiera α, β ∈R. Adem´as, tenemos hψ, ˜ ψiW=W1/2ψ, W1/2˜ ψRmykψkW=  W1/2ψ  Rmpara ψ, ˜ ψ∈Rm. Tomando ¯ ψ=W1/2ψ∈Rme¯ Y=W1/2Y∈Rm×ny multiplicando (2.24) por W−1/2 por la izquierda, llegamos al siguiente problema de autovalores ¯ Y¯ Y>¯ ψ=λ¯ ψen Rm. De ∂L ∂λ (ψ, λ) = 0 en Robtenemos la restricci´on kψkW= 1 que puede tambi´en ser expresada como k¯ ψkRm= 1. Por tanto, las condiciones de optimalidad de primer orden para P1 Wson an´alogas a las de P1(2.4), pero en este caso, la matriz Yy el vector ψtienen que estar pesados por la matriz W1/2. Notemos que ∇ψψL(ψ, λ)=2WY Y >Wψ −λWψ∈Rm×m.Sea ψ∈Rmun vector arbitrario. Como ¯ Y¯ Y>es sim´etrica, existen mautovectores ortonormales (con respec- to al producto escalar eucl´ıdeo) ¯ ψ1,..., ¯ ψm∈Rmde ¯ Y¯ Y>verificando ¯ Y¯ Y>¯ ψi=λi¯ ψi para 1 ≤i≤m. Tomemos ψi=W−1/2¯ ψi,1≤i≤m. Entonces, {ψi}m i=1 forman una base ortonormal (con respecto al producto escalar con peso) en Rmy WY Y >Wψi=λiWψi. Podemos escribir ψen la forma ψ= m X i=1 hψ, ψiiRmψi. EL M´ ETODO POD EN EL ESPACIO RM31 Veamos que ∇ψψL(ψ1, λ1) es semidefinida negativa: ψ>∇ψψL(ψ1, λ1)ψ= 2ψ>WY Y >W−λ1Wψ = 2 m X i=1 m X j=1 hψ, ψiiRmhψ, ψjiRmψ> iWY Y >W−λ1Wψj = 2 m X i=1 m X j=1 (λj−λ1)hψ, ψiiRmhψ, ψjiRmψ> iWψj = 2 m X i=1 (λj−λ1)|hψ, ψiiRm|2≤0. Por tanto, la matriz ∇ψψL(ψ1, λ1) es semidefinida negativa, lo cual implica que se verifica la condici´on necesaria de optimalidad de segundo orden. No obstante, haciendo un argumento an´alogo al de la Secci´on 2.3, se puede probar que ψ1=W−1/2¯ ψ1 resuelve P1 W,donde ¯ ψ1es un autovector de ¯ Y¯ Y>correspondiente al mayor autovalor λ1con  ¯ ψ1 Rm= 1. Debido a la descomposici´on en valores singulares (SVD), el vector ψ1puede tambi´en ser determinado resolviendo el problema de autovalores ¯ Y>¯ Y¯ φ1=λ1¯ φ1 donde ¯ Y>¯ Y=Y>WY , y tomando ψ1=W−1/2¯ ψ1=1 √λ1 W−1/2¯ Y¯ φ1=1 √λ1 Y¯ φ1. Como en la Secci´on 2.1, podemos continuar buscando un segundo vector ψ∈Rmcon hψ, ψ1iW= 0 que maximize n X j=1 |hyj, ψiW|2. Generalicemos a continuaci´on el Teorema 2.3.1: Teorema 2.5.2. Sea Y∈Rm×nuna matriz con rango d≤m´ın{m, n},Wuna matriz sim´etrica y definida positiva, ¯ Y=W1/2Yy`∈ {1, . . . , d}.Adem´as, sea ¯ Y=¯ ΨΣ¯ Φ>la descomposici´on en valores singulares de ¯ Y , donde ¯ Ψ = ¯ ψ1,..., ¯ ψm∈ Rm×m,¯ Φ = ¯ φ1,...,¯ φn∈Rn×nson matrices ortogonales y Σ una matriz de la forma ¯ Ψ>¯ Y¯ Φ =  DΘ Θ Θ  = Σ ∈Rm×n. 32 Manuel Bernardino del Pino Entonces, la soluci´on de (P` W) m´ax ˜ ψ1,..., ¯ ψ`∈Rm ` X i=1 n X j=1 |hyj,˜ ψiiW|2sujeto a h˜ ψi,˜ ψjiW=δij para 1 ≤i, j ≤` viene dada por los vectores ψi=W−1/2¯ ψi, i = 1, . . . , `. Adem´as, argmax (P` W) = ` X i=1 σ2 i= ` X i=1 λi. Demostraci´on. Usando argumentos similares de la prueba del Teorema 2.3.1, se puede demostrar que {ψi}` i=1 resuelve P` W. C´alculo efectivo de la base POD Debido a la SVD y ¯ Y>¯ Y=Y>WY , la base POD {ψi}` i=1 de rango `puede ser determinada por el m´etodo de los snapshots de la siguiente forma. Resolver el problema de autovalores: Y>WY ¯ φi=λi¯ φipara i= 1, . . . , ` y tomar ψi=W−1/2¯ ψi=1 √λi W−1/2¯ Y¯ φi=1 √λi W−1/2W1/2Y¯ φi=1 √λi Y¯ φi para i= 1, . . . , `. Notemos que hψi, ψjiW=ψ> iWψj=δijλj pλiλj para 1 ≤i, j ≤`. Para mnel m´etodo de los snapshots resulta ser m´as r´apido que calcular la base POD resolviendo el problema    ¯ Y¯ Y>¯ ψ=λ¯ ψen Rm, k¯ ψkRm= 1 en R. Nota 2.5.3. Siguiendo las mismas ideas que las ya desarrolladas en la Secci´on 2.4 podemos dar una reformulaci´on del problema (P` W). Nos planteamos, entonces, resolver el siguiente problema de optimizaci´on: (ˆ Pn,` W) m´ın ˜ ψ1,..., ˜ ψ`∈Rm n X j=1 αj  yj− ` X i=1hyj,˜ ψiiW˜ ψi   2 W s.t h˜ ψi,˜ ψjiW=δij para 1 ≤i, j ≤`, EL M´ ETODO POD EN EL ESPACIO RM33 donde los αjdenotan pesos no negativos que especificaremos m´as adelante. Observemos que si tomamos αj= 1 para todo j= 1, . . . , n yW=Im, el problema (ˆ Pn,` W) coincide con (ˆ P`). Para resolver ˆ Pn,` Wvamos a usar las t´ecnicas que hemos visto con anterioridad. Definimos la funci´on Lagrangiana L:Rm×. . . ×Rm |{z } `−veces ×R`×`→R por L(ψ1, . . . , ψ`,Λ) = n X j=1 αj   yj− ` X i=1 hyj, ψiiWψi    2 W + ` X i=1 ` X j=1 Λij 1−hψi, ψjiW para ψ1, . . . , ψ`∈Rmy Λ ∈R`×`con elementos Λij =δij,1≤i, j ≤`. Resulta que la soluci´on de ˆ Pn,` Wviene dada por las condiciones de optimalidad necesarias de primer orden: ∇ψiL(ψ1, . . . , ψ`,Λ) = 0 en Rm,1≤i≤`(2.25) y hψi, ψjiW=δij,1≤i, j ≤`. (2.26) De (2.25), haciendo c´alculos similares a los de la prueba del Teorema 2.3.1, se tiene: Y DY >Wψi=λiψipara i= 1, . . . , ` (2.27) donde D= diag (α1, . . . , αn)∈Rn×n.Escribiendo ψi=W−1/2¯ ψien (2.27) y multiplicando (2.27) por W1/2por la izquierda, llegamos a W1/2Y DY >W1/2¯ ψi=λi¯ ψi. Por otro lado, de (2.26), obtenemos: ¯ ψi,¯ ψjRm=¯ ψ> i¯ ψj=ψ> iWψj=hψi, ψjiW=δij,1≤i, j ≤`. Tomando ¯ Y=W1/2Y D1/2∈Rm×ny usando el hecho de que las matrices WyD son sim´etricas, inferimos que la soluci´on {ψi}` i=1 de ˆ Pn,` Wviene dada a trav´es del siguiente problema de autovalores sim´etrico m×m:    ¯ Y¯ Y>¯ ψi=λi¯ ψi,1≤i≤`, h¯ ψi,¯ ψjiRm=δij,1≤i, j ≤`. 34 Manuel Bernardino del Pino Notemos adem´as que ¯ Y>¯ Y=D1/2Y>WY D1/2∈Rn×n. Por tanto, la base POD de rango `puede tambi´en ser calculada por el m´etodo de los snapschots de la siguiente forma. C´alculo efectivo de la base POD Hay que resolver, en primer lugar, el problema de autovalores sim´etrico n×ndado por    ¯ Y>¯ Y¯ φi=λi¯ φi,1≤i≤` ¯ φi,¯ φjRn=δij,1≤i, j ≤` Despu´es, por la SVD, tomamos ψi=W−1/2¯ ψi=1 √λi W−1/2¯ Y¯ φi=1 √λi Y D1/2¯ φi,1≤i≤`. Observemos que hψi, ψjiW=ψ> iWψj=1 pλiλj ¯ ψ> iD1/2Y>WY D1/2 | {z } =¯ Y>¯ Y ¯ φj=λi pλiλj ¯ φ> i¯ φj=λiδij pλiλj para 1 ≤i, j ≤`, es decir, los vectores de la base POD ψ1, . . . , ψ`son ortonormales en Rmcon respecto al producto escalar h·,·iW. 2.6. Algunas aplicaciones del m´etodo POD 2.6.1. Ecuaciones diferenciales ordinarias Para T > 0 vamos a considerar el siguiente problema de valor inicial semilineal: (SL)   ˙y(t) = Ay(t) + f(t, y(t)) para t∈(0, T], y(0) = y0, donde y0∈Rmes una condici´on inicial elegida, A∈Rm×mes una matriz dada, f: [0, T]×Rm→Rmes una funci´on continua en ambos argumentos y Lipschitziana con respecto a la segunda variable. Es bien conocido que el problema (SL) tiene una ´unica soluci´on cl´asica y∈ C1(0, T;Rm)∩ EL M´ ETODO POD EN EL ESPACIO RM35 C([0, T]; Rm). Entonces, la soluci´on yde (SL) viene dada por la siguiente expresi´on integral en forma impl´ıcita: y(t) = etAy0+Zt 0 e(t−s)Af(s, y(s)) ds, con etA = ∞ X n=0 tnAn n!. Vamos a considerar ahora una partici´on del intervalo [0, T], es decir, {0 = t1< t2< ··· < tn=T}. Por simplificaci´on, vamos a escoger la partici´on de manera uniforme con paso de tiempo ∆t=T/(n−1), es decir, tj= (j−1)∆t. Supongamos que conocemos la soluci´on de (SL) en los instantes de tiempo tjpara j= 1, . . . , n. Nuestro objetivo es determinar una base POD de rango `≤m´ın{m, n}que describa el conjunto yj=y(tj) = etjAy0+Ztj 0 e(tj−s)Af(s, y(s)) ds, j = 1, . . . , n, con el menor error posible. Para ello, nos planteamos el problema (ˆ Pn,` W). Observemos que hemos pasado del caso continuo al discreto calculando yj=y(tj) para cada j= 1, . . . , n. El conjunto de snapshots {yj}n j=1 para ˆ Pn,` Wy por tanto el conjunto generado por {y1, . . . , yn}dependen de la elecci´on de los instantes de tiempo {tj}n j=1. Por consiguiente, los vectores de la base POD, {ψi}` i=1 y los correspondientes autovalores {λi}` i=1 dependen tambi´en de los instantes de tiempo, es decir, ψi=ψn iyλi=λn i,1≤i≤`. Adem´as, no hemos discutido hasta ahora cu´al es la motivaci´on para introducir los pesos no negativos {αj}n j=1 en ˆ Pn,` W. Por esta raz´on nos planteamos las siguientes cuestiones: ¿C´omo elegir unos instantes de tiempo que te permitan obtener una buena aproximaci´on en otros tiempos? ¿Cu´ales son los apropiados pesos no negativos {αj}n j=1 que debemos elegir para que el producto escalar con peso que genera estos valores {αj}n j=1 aproximen al producto escalar continuo? 42 Manuel Bernardino del Pino Si tomamos producto escalar, con respecto a Wcon ψi, obtenemos: hRψi, ψiiW=ZT 0hy(t), ψiiWy(t)dt> Wψi=ZT 0hy(t), ψiiWy(t)>Wψidt =ZT 0|hy(t), ψiiW|2dt. Ahora, sumamos desde i= 1, . . . , d: d X i=1 ZT 0|hy(t), ψiiW|2dt = d X i=1 hRψi, ψiiW= d X i=1hλiψi, ψiiW= d X i=1 λi= m X i=1 λi. Si desarrollamos y(t)∈Rmen t´erminos de {ψi}m i=1, tenemos y(t) = m X i=1 hy(t), ψiiWψi y por tanto, ZT 0 y(t) 2 Wdt = m X i=1 ZT 0|hy(t), ψiiW|2dt = m X i=1 λi Veamos, por ´ultimo, que se cumple b). En efecto, Rnψn i= n X j=1 αjhyj, ψn iiWyjpara todo i∈ {1, . . . , m}yn∈N. Multiplicamos escalarmente con respecto a W por ψn i: hRnψn i, ψn iiW=n X j=1 αjhyj, ψn iiWyj, ψn iW = n X j=1 αjhyj, ψn iiW2. Ahora, sumamos desde i= 1, . . . , d(n): d(n) X i=1 n X j=1 αjhyj, ψn iiW2= d(n) X i=1 hRnψn i, ψn iiW = d(n) X i=1hλn iψn i, ψn iiW = d(n) X i=1 λn i= m X i=1 λn i. EL M´ ETODO POD EN EL ESPACIO RM43 Si desarrollamos yj∈Rmen t´erminos de {ψn i}m i=1, tenemos yj= m X i=1 hyj, ψn iiWψn i, y por tanto, n X j=1 αjky(tj)k2 W= m X i=1 n X j=1 αjhyj, ψn iiW2= m X i=1 λn ipara todo n∈N. LLegados a este punto nos proponemos demostrar que el operador discreto Rnconverge al operador continuo Ry que las correspondientes sucesiones de autovalores y autovectores tambi´en lo hacen. Teorema 2.6.5. Supongamos que (SL) tiene una ´unica soluci´on y: [0, T ]→Rm,y∈ C1(0, T;Rm)∩C([0, T]; Rm) y sean {(ψn i, λn i)}m i=1 y{(ψi, λi)}m i=1 los pares autovectoresautovalores dados en (2.35) y (2.36). Supongamos que `∈ {1, . . . , m}es un valor fijado verificando λ`6=λ`+1 y m X i=`+1 λi6= 0, m X i=`+1 |hy0, ψiiW|26= 0. Entonces, se cumple que l´ım n→∞ kRn−RkL(Rm)= 0,(2.37) donde L(Rm) denota el espacio de Banach de todas las aplicaciones lineales y continuas de Rmen s´ı mismo y en consecuencia: l´ım n→∞ |λn i−λi|= l´ım n→∞ kψn i−ψikW= 0 para 1≤i≤`, l´ım n→∞ m X i=`+1 (λn i−λi)=0, l´ım n→∞ m X i=`+1 |hy0, ψn iiW|2= m X i=`+1 |hy0, ψiiW|2. Demostraci´on. Sea ψ∈Rmun vector arbitrario tal que kψkW= 1. Definamos la funci´on fψ: [0, T]→Rmpor: fψ(t) = hy(t), ψiWy(t) para t∈[0, T]. 44 Manuel Bernardino del Pino Entonces, tenemos que fψ∈C1([0, T ]; Rm) con ˙ fψ(t) = h˙y(t), ψiWy(t) + hy(t), ψiW˙y(t) para t∈[0, T]. Usando el desarrollo de Taylor para la funci´on fψ(t), existen τj1(t)∈(tj, t) y τj2(t)∈ (t, tj+1) dependientes de ttales que: fψ(t) = fψ(tj) + ˙ fψ(τj1(t)) (t−tj) y fψ(t) = fψ(tj+1) + ˙ fψ(τj2(t)) (t−tj+1). Entonces: Ztj+1 tj fψ(t)dt =1 2Ztj+1 tj fψ(tj) + ˙ fψ(τj1(t)) (t−tj)dt +1 2Ztj+1 tj fψ(tj+1) + ˙ fψ(τj2(t)) (t−tj+1)dt =∆t 2(fψ(tj) + fψ(tj+1)) + 1 2Ztj+1 tj ˙ fψ(τj1(t)) (t−tj)dt +1 2Ztj+1 tj ˙ fψ(τj2(t)) (t−tj+1)dt Por tanto,  Rnψ−Rψ W=    n X j=1 αjfψ(tj)−ZT 0 fψ(t)dt   W =    n−1 X j=1 ∆t 2(fψ(tj) + fψ(tj+1)) −Ztj+1 tj fψ(t)dt   W ≤1 2 n−1 X j=1 Ztj+1 tj ˙ fψ(τj1(t))  W|t−tj|+ ˙ fψ(τj2(t))  W|t−tj+1|dt ≤1 2m´ax t∈[0,T] ˙ fψ(t) W n−1 X j=1 (t−tj)2 2−(tj+1 −t)2 2! t=tj+1 t=tj =1 2m´ax t∈[0,T] ˙ fψ(t) W n−1 X j=1 (∆t)2=(∆t)2(n−1) 2m´ax t∈[0,T] ˙ fψ(t) W =∆tT 2m´ax t∈[0,T] ˙ fψ(t) W =∆tT 2m´ax t∈[0,T]kh˙y(t), ψiWy(t) + hy(t), ψiW˙y(t)kW = ∆tT m´ax t∈[0,T]k˙y(t)kWky(t)kW≤∆tTkyk2 C1([0,T];Rm). EL M´ ETODO POD EN EL ESPACIO RM45 Luego, kRn−RkL(Rm)= sup kψkW=1 kRnψ−RψkW≤∆tTkyk2 C1([0,T];Rm) ∆t→0 −→ 0. Vamos a probar a continuaci´on que m X i=1 λn i→ m X i=1 λicuando n→ ∞.(2.38) En efecto, por el Lema 2.6.4, la convergencia dada en (2.38) se puede reescribir como n X j=1 αjky(tj)k2 W→ZT 0ky(t)k2 Wdt cuando ∆t→0.(2.39) Luego, basta elegir los αjde manera que se garantice dicha convergencia. Para ello, si aproximamos esa integral usando la regla del trapecio, ZT 0ky(t)k2 Wdt = n−1 X j=1 Ztj+1 tjky(t)k2 Wdt ≈ n−1 X j=1 ∆t 2ky(tj)k2 W+ky(tj+1)k2 W, los pesos αj(j= 1, . . . , n) que tendr´ıamos que poner para tener la convergencia buscada ser´ıan los siguientes: α1=∆t 2, αj= ∆tpara 2 ≤j≤n−1, αn=∆t 2.(2.40) Ahora, fijemos `∈ {1, . . . , m}tal que λ`6=λ`+1. Entonces, por el an´alisis espectral de operadores compactos y (2.37), se sigue que λn i→λipara 1 ≤i≤`cuando n→ ∞.(2.41) Combinando (2.38) y (2.41), existe un n0∈Ntal que m X i=`+1 λn i≤2 m X i=`+1 λipara todo n≥n0(2.42) si m X i=`+1 λi6= 0. Adem´as, para `como antes, n0puede ser escogido de forma que se verifique: d(n) X i=`+1 |hy0, ψn iiW|2≤2 m X i=`+1 |hy0, ψiiW|2para todo n≥n0(2.43) 46 Manuel Bernardino del Pino siempre que Pm i=`+1 |hy0, ψiiW|26= 0 y se cumpla (2.37). En efecto, recordemos que el vector y0∈Rmes la condici´on inicial dada en el problema (SL). Por tanto, se verifica que ky0k2 W= m X i=1 |hy0, ψiiW|2.(2.44) Si t1= 0, tenemos que y0∈ hy1, . . . , ynipara cada ny ky0k2 W= d(n) X i=1 |hy0, ψn iiW|2.(2.45) Luego, para ` < d(n) y usando (2.44) y (2.45), es f´acil comprobar que d(n) X i=`+1 |hy0, ψn iiW|2= d(n) X i=1 |hy0, ψn iiW|2− ` X i=1 |hy0, ψn iiW|2+ ` X i=1 |hy0, ψiiW|2 + m X i=`+1 |hy0, ψiiW|2− m X i=1 |hy0, ψiiW|2 = ` X i=1 |hy0, ψiiW|2−|hy0, ψn iiW|2+ m X i=`+1 |hy0, ψiiW|2. Como consecuencia de (2.37) y de que `es tal que λ`6=λ`+1, tenemos que l´ımn→∞ kψn i−ψikW= 0 para i= 1, . . . , ` y por ende se verifica (2.43). Nota 2.6.6. Hemos realizado el estudio te´orico sobre el problema (SL), pero podr´ıamos haberlo hecho para cualquier otro problema realizando argumentos an´alogos a los ya descritos. C´alculo efectivo de la base POD Definamos el operador lineal y continuo Y:L2(0, T)→Rmde la siguiente forma. Yφ:= ZT 0 φ(t)y(t)dt para φ∈L2(0, T). Observemos que Yes un operador compacto ya que es de rango finito (Teorema 1.2.10). Definamos ahora su operador adjunto Y?:Rm→L2(0, T ) como hY?ψ, φiL2(0,T)=hψ, YφiWpara todo (ψ, φ)∈Rm×L2(0, T). EL M´ ETODO POD EN EL ESPACIO RM47 Usando esta definici´on, es f´acil comprobar que (Y?ψ) (t) = hψ, y(t)iWpara todo ψ∈Rmy para casi todo t∈[0, T]. Adem´as, Y∗es compacto al serlo Y(Teorema 1.2.7). Entonces, tenemos: YY?ψ=ZT 0hψ, y(t)iWy(t)dt =ZT 0hy(t), ψiWy(t)dt =Rψ para todo ψ∈Rm; es decir, R=YY?. M´as a´un, (Y?Yφ) (t) = ZT 0 φ(s)y(s)ds, y(t)W =ZT 0hy(s), y(t)iWφ(s)ds =: (Kφ)(t) para todo φ∈L2(0, T) y en casi todo t∈[0, T].Por tanto, K=Y?Y. Se puede probar que el operador K:L2(0, T)→L2(0, T) es lineal, acotado, autoadjunto y compacto (esto ´ultimo es trivial por el Teorema 1.2.9). Por tanto, la base POD puede tambi´en ser computada de la siguiente forma: Resolver Kφi=λiφipara 1 ≤i≤`, λ1≥. . . ≥λ`>0,ZT 0 φi(t)φj(t)dt =δij (2.46) y tomar ψi=1 √λiYφi=1 √λiZT 0 φi(t)y(t)dt para i= 1, . . . , `. Observemos que (2.46) es un problema de autovalores en el espacio de dimensi´on infinita L2(0, T). 2.6.2. Ecuaciones en derivadas parciales Vamos a considerar la ecuaci´on del calor uno-dimensional: (EC)       θt(t, x) = θxx(t, x) para todo (t, x)∈(0, T)×Ω, θx(t, 0) = θx(t, 1) = 0 para todo t∈(0, T), θ(0, x) = θ0(x) para todo x∈Ω, donde Ω = (0,1), T > 0 y θ0∈ C(Ω) es una condici´on inicial dada. Para resolver (EC) num´ericamente vamos a aplicar una aproximaci´on por diferencias finitas en la variable espacial x. Igual que hicimos en el Ejemplo 2.5.1, vamos a construir una partici´on uniforme {xi}m i=1 del intervalo [0,1]. Denotemos por yi: [0, T]→Ra la 48 Manuel Bernardino del Pino aproximaci´on num´erica de θ(·, xi) para i= 1, . . . , m. La derivada θxx y la condici´on de contorno son discretizadas usando diferencias finitas centradas de segundo orden; as´ı que obtenemos las siguientes ecuaciones diferenciales ordinarias para las funciones yi:             ˙y1(t) = −2y1(t)+2y2(t) h2, ˙yi(t) = yi−1(t)−2yi(t) + yi+1(t) h2, i = 2, . . . , m −1, ˙ym(t) = −2ym(t)+2ym−1(t) h2, para t∈(0, T]. Adem´as, de θ(0, x) = θ0(x), inferimos la condici´on inicial yi(0) = θ0(xi), i = 1, . . . , m. Vamos a expresar lo anterior como un sistema de ecuaciones diferenciales. Para ello, introducimos la matriz A=1 h2            −2 2 0 1−2 1 ......... 1−2 1 0 2 −2            ∈Rm×m y los vectores y(t) =      y1(t) . . . ym(t)      para t∈[0, T], y0=     θ0(x1) . . . θ0(xm)     ∈Rm, de modo que podemos escribir:    ˙y(t) = Ay(t) para t∈(0, T], y(0) = y0. Observemos que este problema lineal de valor inicial coincide con (SL) tomando f= 0. Adem´as ahora el vector y(t), t ∈[0, T],representa una funci´on en Ω evaluada en mpuntos. Por tanto, deber´ıamos dotar Rmpor un producto escalar con peso que represente un producto escalar discreto en un apropiado espacio de funciones. Aqu´ı, elegimos el producto escalar introducido en (2.23). Despu´es elegimos una partici´on EL M´ ETODO POD EN EL ESPACIO RM49 {tj}n j=1 del intervalo [0, T] y definimos yj=y(tj) para j= 1, . . . , n. Si estamos interesados en encontrar una base POD de rango `≤m´ın{m, n}que describa el conjunto de puntos {yj}n j=1, debemos plantear el problema ˆ Pn,` W. Cap´ıtulo 3 Aplicaciones pr´acticas del m´etodo POD En este cap´ıtulo vamos a aplicar los resultados te´oricos que hemos tratado con anterioridad en algunos casos pr´acticos, considerando un ejemplo para el caso discreto y otro para el caso continuo. Para ello, usaremos el software matem´atico MATLAB. 3.1. Caso discreto Vamos a considerar una serie de patios de los que tenemos mediciones experimentales de sus temperaturas internas cada hora en un determinado per´ıodo del a˜no. M´as concretamente, para que las temperaturas se mueven en un rango similar, hemos escogido patios en una misma zona clim´atica cuyas mediciones se corresponden a per´ıodos similares del a˜no (en concreto, verano - oto˜no). A continuaci´on, enumeramos los patios escogidos y las mediciones de las semanas consideradas. 1. Carlos Rubio en C´ordoba (26 de julio a 1 de agosto). 2. Carlos Rubio en C´ordoba (9-15 de agosto). 3. Hern´an Cort´es en Sevilla (5-11 de septiembre). 4. Hern´an Cort´es2 en Sevilla (5-11 de septiembre). 5. IMUS, Instituto de Matem´aticas de la Universidad de Sevilla (14-20 de agosto). 6. IMUS, Instituto de Matem´aticas de la Universidad de Sevilla (4-10 de septiembre). 50 APLICACIONES PR ´ ACTICAS DEL M´ ETODO POD 51 7. UNED, Universidad Nacional de Educaci´on a Distancia en Sevilla (14-20 de agosto). 8. UNED, Universidad Nacional de Educaci´on a Distancia en Sevilla (21-27 de agosto). Nuestro objetivo, por tanto, ser´a calcular una base POD para esta serie de datos experimentales que nos permita aproximar las temperaturas de un nuevo patio en un per´ıodo del a˜no similar. Para ello, vamos a guiarnos del c´alculo efectivo de la base POD descrito en la p´agina 26; concretamente usaremos (2.18) y (2.19). Es decir, vamos a construir la matriz Y= (y1|···|y8)∈R168×8, donde cada columna yjpara j= 1, . . . 8 representa las temperaturas de un patio medidas cada hora en la semana correspondiente, es decir, es un vector con 168 componentes. Posteriormente, calculamos los autovalores y autovectores de la matriz Y>Yy calculamos el valor de `observando la funci´on ε(`) definida en (2.21). En la siguiente tabla hemos recogido los valores de la funci´on ε(`) para cada `= 1,...,8. ε(1) ε(2) ε(3) ε(4) ε(5) ε(6) ε(7) ε(8) 99.689 % 99.863 % 99.930 % 99.972 % 99.990 % 99.995 % 99.998 % 100 % Tabla 1: Valores de la funci´on ε(`). Como podemos observar todos estos valores est´an muy pr´oximos al 100 %, lo que quiere decir que con muy pocas bases podemos capturar casi la totalidad de la energ´ıa del sistema. Por ello, vamos a elegir, por ejemplo, `= 2. Para terminar, vamos a considerar un nuevo patio (IMUS, 21-27 de agosto) y observemos qu´e aproximaci´on sobre las temperaturas internas obtenemos gracias a la base POD. Si yes el vector de las temperaturas internas de este nuevo patio, lo podemos aproximar como combinaci´on lineal de los vectores de la base POD de la forma: y≈ hy1, ψ1iψ1+hy2, ψ2iψ2, 58 Manuel Bernardino del Pino 15 d´ıas Valor aproximado N=1 N=2 N=3 y10.34912 0.35038 0.34704 0.34912 y2-0.95578 7.49848e-04 -0.02364 -0.95578 y30.96785 0.00119 0.99082 0.96785 Tabla 3: Aproximaci´on de la soluci´on a los 15 d´ıas seg´un el n´umero de funciones de base tomadas. 20 d´ıas Valor aproximado N=1 N=2 N=3 y1-0.54455 -0.53949 -0.54252 -0.54455 y20.88935 -0.00115 -0.02333 0.88935 y30.92059 -0.00184 0.89810 0.92059 Tabla 4: Aproximaci´on de la soluci´on a los 20 d´ıas seg´un el n´umero de funciones de base tomadas. 25 d´ıas Valor aproximado N=1 N=2 N=3 y10.47732 0.47856 0.47533 0.47732 y2-0.91582 0.00102 -0.02262 -0.91582 y30.93929 0.00163 0.96130 0.93929 Tabla 5: Aproximaci´on de la soluci´on a los 25 d´ıas seg´un el n´umero de funciones de base tomadas. Bibliograf´ıa [1] Barnston, A. G. & Ropelewski, C. F., Prediction of ENSO episodes using canonical correlation analysis, Journal of Climate 5, 1992, 1316–1345. [2] Bayly, P. V., Johnson, E. E., Wolf, P. D., Smith, W. M., & Ideker, R. E., Predicting patterns of epicardial potentials during ventricular fibrillation, IEEE Transactions on Biomedical Engineering 42, 1995, 898–907. [3] Chac´ on, T., Fern´ andez, U. & G´ omez M., Computational Modeling of Gurney Flaps and Microtabs by POD Method, Energies 11(8), 2018. [4] Epureanu, B. I., Hall, K. C., & Dowell, E. H., Reduced-order models of unsteady viscous flows in turbomachinery using viscous-inviscid coupling, Journal of Fluids and Structures 15, 2001, 255–273. [5] Golub, G.H. & Van Loan, C.F., Matrix Computations, The Johns Hopkins University Press, Baltimore, 1990. [6] Gohberg, I. & Goldberg, S., Basic operator theory, Birkh¨auser, 1980. [7] Leen, T. K., Rudnick, M., & Hammerstrom, R., Hebbian feature discovery improves classifier efficiency, in Proceedings of the IJCNN, IEEE, Piscataway, NJ, 1990, 51–56. [8] Liang, Y. C., Lee, H. P., Lim, S. P., Lin, W. Z., Lee, K. H., & Wu, C. G., Proper orthogonal decomposition and its applications, Part I: Theory’, Journal of Sound and Vibration 252, 2002, 527–544. [9] Mees, A. I., Rapp, P. E., & Jennings, L. S., Singular value decomposition and embedding dimension, Physical Review A 36, 1987, 340–346. 59 60 Manuel Bernardino del Pino [10] Nocedal, J. & Wright, S.J, Numerical Optimization, Springer Series in Operation Research, 2006. [11] Sirovich, L., Turbulence and the dynamics of coherent structures, Quarterly of Applied Mathematicss, XLV, 1987, 561-590. [12] Volkwein, S., Proper Orthogonal Descomposition: Theory and Reduced-Order Modelling, Lecture Notes, University of Konstanz, 2013. [13] Wax, M. & Kailath, T., Detection of signals by information theoretic criteria, IEEE Transactions on Acoustics, Speech and Signal Processing 33, 1985, 387–392. [14] Wu, C. G., Liang, Y. C., Lin, W. Z., Lee, H. P., & Lim, S. P., A note on equivalence of proper orthogonal decomposition methods, Journal of Sound and Vibration 265, 2003, 1103–1110.