Resolución numérica de problemas diferenciales por el Método de Soluciones Fundamentales

Universidad de Sevilla Facultad de Matem´aticas Trabajo Fin de Estudios Grado en Matem´aticas Resoluci´on num´erica de problemas diferenciales por el M´etodo de Soluciones Fundamentales David Pino Molero Tutor Anna Doubova Krasotchenko. Dpto. Ecuaciones Diferenciales y An´alisis Num´erico. Abstract This work examines the Method of Fundamental Solutions (MFS) together with other meshless numerical methods, such as the Method of Particular Solutions, for solving differential problems, both direct and inverse. The theoretical construction of the method is formalized in the context of partial differential equations (PDEs) using fundamental notions from Distribution Theory. The main objective has been to develop a practical approach, illustrated through classical examples that demonstrate the applicability and efficiency of the MFS to various types of problems (elliptic, parabolic and hyperbolic), as well as to provide MATLAB codes for effective resolution. A general description of the MFS is presented, supported by a collection of applied examples and their corresponding computational implementations, which can serve both as an introduction to meshless methods as a starting point for their extension to more complex problems. Keywords: Method of Fundamental Solutions, meshless method, partial differential equations, inverse problems, MATLAB. Resumen En este trabajo se estudia el M´etodo de Soluciones Fundamentales (MSF) junto con otros m´etodos num´ericos sin malla, como el M´etodo de Soluciones Particulares, para la resoluci´on de problemas diferenciales, tanto directos como inversos. Se formaliza la construcci´on te´orica del m´etodo en el contexto de ecuaciones en derivadas parciales (EDPs) usando nociones fundamentales de la Teor´ıa de Distribuciones. El objetivo principal ha sido desarrollar una aproximaci´on pr´actica, ilustrada mediante ejemplos cl´asicos que demuestran la aplicabilidad y eficiencia del MSF en distintos tipos de problemas (el´ıpticos, parab´olicos e hiperb´olicos), as´ı como ilustrar mediante c´odigos en MATLAB la resoluci´on efectiva. Se presenta una descripci´on general del MSF respaldada por una colecci´on de ejemplos y sus correspondientes desarrollos computacionales, que puede servir tanto de introducci´on a los m´etodos sin malla, como de punto de partida para su extensi´on a problemas de ´ındole m´as compleja. Palabras clave: M´etodo de Soluciones Fundamentales, m´etodo sin malla, ecuaciones en derivadas parciales, problemas inversos, MATLAB. ´ Indice Introducci´on 4 1. Conceptos previos 6 1.1. La distribuci´on de Dirac ................................ 7 1.2. Definici´on de soluci´on fundamental .......................... 14 1.3. Ejemplos en una, dos y tres dimensiones ....................... 15 1.3.1. Ejemplo de una EDO ............................. 15 1.3.2. La ecuaci´on de Laplace en dimensiones 2 y 3 ................ 16 2. Descripci´on general del m´etodo 19 2.1. Breve introducci´on hist´orica .............................. 19 2.2. Implementaci´on del m´etodo .............................. 20 2.2.1. Convergencia y estabilidad del m´etodo .................... 22 2.3. Combinaci´on del MSF con otros m´etodos sin malla ................. 22 3. Ejemplos de Aplicaciones 26 3.1. Problemas directos ................................... 26 3.1.1. Ecuaci´on de Laplace en una bola ....................... 26 3.1.2. Problema el´ıptico general ........................... 33 3.1.3. Ecuaci´on del Calor ............................... 40 3.1.4. Ecuaci´on de Ondas ............................... 42 3.2. Problemas Inversos ................................... 46 3.2.1. Problemas inversos para EDPs el´ıpticas ................... 46 3.2.2. Problemas inversos asociados a la ecuaci´on del calor ............ 50 4. Implementaci´on en MATLAB 53 4.1. C´odigo del problema el´ıptico con condiciones DirichletNeumann ........................................ 53 4.2. C´odigo de la ecuaci´on del calor con condiciones Dirichlet .............. 62 4.3. C´odigo del problema de ondas con condiciones Dirichlet .............. 65 4.4. C´odigo del problema el´ıptico inverso asociado a un dato en una subregi´on . . . . 68 4.5. C´odigo del problema el´ıptico inverso asociado a la ecuaci´on del calor ....... 71 Bibliograf´ıa 74 1 ´ Indice de figuras 2.1. Ejemplo de distribuci´on de las fuentes para el MSF en una elipse. ......... 21 3.1. Distribuci´on de puntos en la frontera y fuentes en una bola en R2......... 27 3.2. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF para la ecuaci´on de Laplace en una bola. ................................ 30 3.3. Error absoluto entre la soluci´on exacta y la obtenida con el MSF para la ecuaci´on de Laplace en una bola. ................................ 30 3.4. Distribuci´on de los puntos sobre la frontera y de las fuentes en el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. . . . . 31 3.5. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. ......... 31 3.6. Error absoluto entre la soluci´on exacta y la obtenida con el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. ......... 32 3.7. Evoluci´on del error para u(x) = ex1cos(x2) desde N= 50 hasta N= 100. . . . . 32 3.8. Evoluci´on del error para u(x) = cos(πx1) sinh(πx2) + sin(2πx1) cosh(2πx2) desde N= 50 hasta N= 100. ................................ 33 3.9. Distribuci´on de los puntos interiores, sobre la frontera y de las fuentes en el MSFMSP para el problema el´ıptico con condiciones Dirichlet-Neumann. ........ 38 3.10. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF-MSP para el problema el´ıptico con condiciones Dirichlet-Neumann. ............... 39 3.11. Error absoluto entre la soluci´on exacta y la obtenida con el MSF-MSP para el problema el´ıptico con condiciones Dirichlet-Neumann. ............... 39 3.12. Distribuci´on de los puntos interiores, sobre la frontera y de las fuentes en el MSFMSP para el problema el´ıptico con condiciones Dirichlet. .............. 40 3.13. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF-MSP para el problema el´ıptico con condiciones Dirichlet. ..................... 40 3.14. Error absoluto entre la soluci´on exacta y la obtenida con el MSF-MSP para el problema el´ıptico con condiciones Dirichlet. ..................... 41 3.15. Distribuci´on de los puntos interiores, sobre la frontera y de las fuentes en el MSFMSP para la ecuaci´on del calor con condiciones Dirichlet. ............. 42 3.16. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on del calor con condiciones Dirichlet. ..................... 43 2 3.17. Error absoluto entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on del calor en el tiempo final. ......................... 43 3.18. Comparaci´on entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on de ondas con condiciones Dirichlet. ..................... 45 3.19. Error absoluto entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on de ondas en el tiempo final. ......................... 45 3.20. Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de la subregi´on D................................ 47 3.21. Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de γ........................................ 49 3.22. Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de γ........................................ 51 3 Introducci´on El estudio y la resoluci´on de problemas asociados a ecuaciones en derivadas parciales (EDPs) es de vital importancia en el desarrollo de la Matem´atica Aplicada y la Ingenier´ıa, pues este tipo de problemas son los que permiten modelar multitud de fen´omenos f´ısicos y biol´ogicos, abarcando campos como Mec´anica, Termodin´amica, Electromagnetismo, Ciencias de la Salud, Econom´ıa, etc. Ante la patente imposibilidad de la resoluci´on exacta de la mayor parte de este tipo de problemas, los m´etodos num´ericos surgen como la alternativa para obtener una aproximaci´on fidedigna de la soluci´on de estos problemas, y de hecho son este tipo de m´etodos los que nos han permitido progresar y desarrollar las disciplinas antes mencionadas, y que hoy en d´ıa representan un pilar fundamental de la sociedad tecnol´ogica en la que estamos inmersos. Los procesos num´ericos que han ayudado a la resoluci´on de estos problemas hist´oricamente han necesitado de una discretizaci´on completa del dominio donde se plantea el problema. Ejemplo de este tipo de m´etodos pueden ser el M´etodo de Elementos Finitos (MEF) o de Vol´umenes Finitos (MVF). Pero, ¿qu´e ocurre si la geometr´ıa asociada al dominio es altamente compleja?, ¿y si tan s´olo conocemos parcialmente el interior de este dominio? Como respuesta a estas dificultades, surgen durante las ´ultimas d´ecadas los llamados m´etodos sin malla, como una alternativa m´as flexible e igualmente eficiente para resolver problemas asociados a EDPs. Dentro de esta familia de m´etodos, destaca como uno de sus principales exponentes el M´etodo de Soluciones Fundamentales (MSF), que brilla por su sencillez conceptual y por sus buenos resultados en los problemas que permiten su aplicaci´on. El presente trabajo se centra en la exploraci´on de las posibilidades de aplicaci´on del MSF a problemas de contorno el´ıpticos, as´ı como de estudiar formas de extender este m´etodo a problemas parab´olicos e hiperb´olicos mediante una discretizaci´on del intervalo temporal. Adem´as de centrarse en problemas directos, que pueden parecer a priori de mayor inter´es pr´actico, tambi´en se ejemplifica c´omo el MSF puede ser usado para resolver problemas inversos, que son de gran utilidad en ´areas tan diversas como Medicina, Geof´ısica, Ingenier´ıa de Materiales o Hidrolog´ıa, entre otros. Sin embargo, con un enfoque pr´actico centrado en las aplicaciones, de poco sirve dominar en profundidad el desarrollo te´orico si no se cuenta con capacidad de llevar a cabo, de facto, las aproximaciones que propone el MSF. Es por esto que se incluyen detallados desarrollos computacionales en MATLAB que permiten aplicar el m´etodo en todos los problemas que se 4 Introducci´on 5 tratan. As´ı, el objetivo principal de este trabajo es doble: por un lado, busca ofrecer una introducci´on rigurosa y accesible al MSF y a los fundamentos te´oricos necesarios para su aplicaci´on; por otro, ilustrar su implementaci´on computacional mediante ejemplos pr´acticos, empleando el entorno MATLAB para resolver problemas cl´asicos, tanto directos como inversos, en dominios bidimensionales. Asimismo, se deja la puerta abierta para los problemas en dominios tridimensionales, a los que el m´etodo es f´acilmente extensible, y que no se incluyen en el presente trabajo para no alargar en exceso el mismo. El desarrollo se divide en los siguientes cap´ıtulos: Cap´ıtulo 1: Conceptos previos de la Teor´ıa de Distribuciones para la formalizaci´on del m´etodo. Cap´ıtulo 2: Introducci´on hist´orica y descripci´on del m´etodo. Cap´ıtulo 3: Ejemplos de resoluci´on de problemas directos e inversos. Cap´ıtulo 4: Implementaci´on en MATLAB de ejemplos seleccionados y acceso al repositorio completo. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 12 A partir de esta funci´on, definimos ρε(x) = 1 εNρ(x/ε),∀x∈RNy∀ε > 0,que cumple          ρε≥0, sop(ρε)⊆B(0; ε), ZRN ρεdx = 1. Vamos a calcular el l´ımite de ρεen D′(Ω) cuando ε→0. Sea φ∈ D(RN), entonces, haciendo la identificaci´on dada por el Teorema 1.8, tenemos ⟨ρε, φ⟩=ZRN ρε(x)φ(x)dx = ±φ(0) ZRN ρε(x)φ(x)dx −ZRN ρε(x) | {z } =1 φ(0) dx +φ(0) = =ZRN ρε(x)(φ(x)−φ(0)) dx +φ(0). Veamos que el primer sumando tiende a cero. Como φes continua, para todo γ > 0, ∃ε0>0 tal que si |x|< ε0entonces |φ(x)−φ(0)|< γ. Entonces para todo ε≤ε0, se tiene ZRN ρε(x)(φ(x)−φ(0)) dx≤ZRN ρε(x)|φ(x)−φ(0)|dx = sop(ρε)⊆B(0;ε) ZB(0;ε) ρε(x)|φ(x)−φ(0)| | {z } |x|<ε0=⇒<γ dx < γ ZB(0;ε) ρε(x)dx = ρε=0 x/∈B(0;ε)γZRN ρε(x)dx | {z } =1 =γ. Luego ZRN ρε(x)(φ(x)−φ(0)) dx →0, y como ZRN ρε(x)(φ(x)−φ(0)) dx =ZRN ρε(x)φ(x)dx −φ(0) ZRN ρε(x)dx | {z } =1 , entonces, ⟨ρε, φ⟩ → φ(0) = ⟨δ0, φ⟩, por tanto, ρε⇀ δ0. Luego, podemos entender la distribuci´on δ0como el l´ımite de unas funciones cuyo soporte tiende a 0 y que cuanto menor es el soporte, m´as altos son los valores que toma. Entendemos que en el l´ımite, el soporte se reduce al punto 0, donde toma un valor infinito. Es decir, estar´ıamos aproximando una distribuci´on que no se puede identificar con una funci´on de L1 loc(Ω) por una sucesi´on de distribuciones que s´ı que pueden identificarse con una funci´on L1 loc(Ω) y que tienden a ella como distribuci´on. Esta ser´a la interpretaci´on que haremos de la delta de Dirac cuando trabajemos con ella. N´otese que, aunque hayamos tomado s= 0 por simplicidad en el desarrollo, el resultado es cierto para cualquier sbajo una traslaci´on. Observaci´on 1.14 La delta de Dirac tambi´en se suele denotar en la literatura como δs(x) = δ(x−s), que no hace sino resaltar la propiedad que cumple bajo traslaciones (ver [9]), es decir, δs(x) = δ0(x−s),∀x, s ∈RN. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 13 Teorema 1.15 (Propiedades fundamentales de δs) Dada f∈ D(Ω), se cumple lo siguiente 1. ZRN f(x)δs(x)dx =f(s)(Propiedad de muestreo). 2. ZRN δs(x)dx = 1 (Propiedad de normalizaci´on). Demostraci´on: Este resultado es inherente a la propia definici´on de δs, como se explica en la Secci´on 1.15 de [1]. Aunque se hable de δsunidimensional, el resultado es completamente generalizable a RN, ya que δsdefinida en RNse puede entender como un producto de Ndeltas unidimensionales: δs= N Y i=1 δsi. As´ı, utilizando el Teorema de Fubini (ver Teorema 8.8 de [18]), tenemos ZRN f(x)δs(x)dx =Z∞ −∞ . . . Z∞ −∞ f(x1, . . . , xN) N Y i=1 δsi(xi)dx1. . . dxN. Luego se tiene la generalizaci´on sin m´as que usar el resultado unidimensional. □ Observaci´on 1.16 En virtud de las propiedades vistas, entenderemos el producto de una funci´on φpor la delta de Dirac como: φ(x)δs(x) = φ(s)δs(x). Notemos una relaci´on interesante cuando trabajamos con una delta de Dirac definida sobre R. Veamos cu´al es su derivada distribucional, ya que, evidentemente, no es derivable en el sentido cl´asico. Definici´on 1.17 (La funci´on salto de Heaviside) Se define la funci´on salto de Heaviside, denotada por H, como H(x) = (1 si x≥0, 0 si x < 0. Se trata de una funci´on de L1 loc(R), y as´ı, H∈ D(Ω) se puede definir en el marco de la Teor´ıa de Distribuciones como ⟨H, φ⟩=Z+∞ −∞ H(x)φ(x)dx =Z+∞ 0 φ(x)dx TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 14 Teorema 1.18 (Relaci´on entre la delta de Dirac y la funci´on salto de Heaviside) Se cumple que ⟨H′, φ⟩=⟨δ0, φ⟩,∀φ∈ D(R). Demostraci´on: Basta calcular la derivada siguiendo la Definici´on 1.11: ⟨H′, φ⟩=−⟨H, φ′⟩=−Z+∞ −∞ H(x)φ′(x)dx =−Z+∞ 0 φ′(x)dx =−ZM 0 φ′(x)dx =φ(0), donde hemos usado que como φ∈ D(R), ∃M > 0 tal que sop(φ)⊂[−M, M]. As´ı, como se tiene que φ(0) = ⟨δ0, φ⟩, deducimos entonces que H′=δ0∈ D′(R), luego H′/∈L1 loc(R). □ 1.2. Definici´on de soluci´on fundamental Consideramos un problema del tipo (1.1), con f≡0, es decir, una EDP de la forma Lu= 0. Definici´on 1.19 (Soluci´on fundamental de una EDP lineal el´ıptica) Diremos que la distribuci´on en RN,G, es una soluci´on fundamental de (1.1) si verifica LG=δs. Se observa que la funci´on Ges singular en x=s. Se dice que ses la singularidad de la soluci´on fundamental. Las soluciones fundamentales adquieren especial relevancia, ya que a partir de ellas se puede construir cualquier otra soluci´on para un problema de la forma Lu=f, pues si Ges una soluci´on fundamental, entonces la convoluci´on de fcon Ges soluci´on de Lu=f (ver [9]), es decir, u(x)=(f∗G(x, s))(x) = ZΩ f(s)G(x, s)ds verifica Lu=f, donde denotamos por G(x, s) a la soluci´on de LG=δs. Notemos que la definici´on y el desarrollo que hemos hecho se encuadran en el contexto de un problema estacionario. Por tanto, la soluci´on fundamental para un problema el´ıptico representa la soluci´on de un problema con una fuente puntual en el punto s, si estamos considerando Lu=δs. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 15 Observaci´on 1.20 En el desarrollo de este trabajo, nos basta considerar las soluciones fundamentales de EDPs el´ıpticas, ya que en el caso de EDPs dependientes del tiempo, tomaremos una discretizaci´on del tiempo, de forma que trabajaremos con tantos problema el´ıpticos como elementos tenga la partici´on temporal. 1.3. Ejemplos en una, dos y tres dimensiones En lo que sigue, consideraremos por simplicidad δ0≡δla delta de Dirac centrada en s= 0. Los desarrollos siguientes ser´an v´alidos considerando tambi´en δspara cualquier fuente x=s, por traslaci´on. 1.3.1. Ejemplo de una EDO Para el siguiente desarrollo, nos basaremos en el inicio del Cap´ıtulo 4 de [23]. Consideremos la EDO lineal de primer orden, con coeficientes constantes, cuyo operador diferencial tiene la forma L=d dx −a, x ∈R, a ∈R. Estamos considerando una EDO, aunque esta idea se puede extender de forma natural a SDOs. Buscamos una funci´on G:R→Rque verifique la expresi´on dG dx −aG =δ. (1.5) Si a= 0, la respuesta es clara, G′(x) = δ=⇒G(x) = Zδ(x)dx =H(x) + C, donde Ces una constante. Si a= 0, para resolver (1.5) basta considerar como factor integrante e−ax, de modo que (1.5) equivale a e−axG′−ae−axG=e−axδ⇐⇒ e−axG′=e−axδ. En virtud de la Observaci´on 1.16, podemos considerar e−axδ=δ. Integrando, tenemos e−axG=H(x) + C=⇒G(x) = eaxH(x) + Ceax, donde C∈Res una constante arbitraria. Esta ser´ıa la soluci´on de nuestro problema (1.5). En particular, si consideramos C=−1/2, tenemos G(x) = eax H(x)−1 2=1 2sgn(x)eax. Con esto obtenemos una soluci´on sim´etrica de (1.5), siendo sgn(x) = (1 si x > 0, −1 si x < 0. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 16 Observaci´on 1.21 Si en lugar de considerar δ=δ0consideramos δs, con algunos simples cambios, tenemos que la expresi´on general de la soluci´on fundamental es: G(x, s) = ea(x−s)Hs(x) + Ceax, donde Hs(x) = (1 si x≥s, 0 si x < s. 1.3.2. La ecuaci´on de Laplace en dimensiones 2 y 3 Consideramos la EDP conocida como ecuaci´on de Laplace, que viene dada por ∆u= 0,(1.6) donde L=∆= ∂2 ∂x2 1 +... +∂2 ∂x2 N . Construcci´on de la soluci´on fundamental Nos basaremos en el desarrollo de la Secci´on 2.2 de [7]. Vamos a buscar la soluci´on de un problema de la forma: ∆G=δs. Se puede comprobar que la ecuaci´on de Laplace es invariante bajo rotaciones, parece por tanto razonable buscar soluciones con simetr´ıa radial. Por tanto buscaremos una soluci´on de (1.6) en Ω = RNque sea de la forma u(x) = v(r), donde r=|x|= (x2 1+...+x2 N)1/2yvuna funci´on que buscaremos que valga 0 si x= 0, es decir, si r= 0. Primero, notemos que para i= 1, ..., N, se tiene ∂r ∂xi =∂ ∂xi (x2 1+... +x2 N)1/2=1 2(x2 1+... +x2 N)−1/22xi=xi r. Por tanto, usando la Regla de la Cadena, tenemos para todo i= 1, ..., N uxi=v′(r)xi r=⇒uxixi=v′′(r)x2 i r2+v′(r)1 r−x2 i r3, en consecuencia ∆u=v′′(r) + N−1 rv′(r). TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 17 Si consideramos r= 0, e imponemos que la soluci´on cumpla ∆u= 0, obtenemos v′′ +N−1 rv′= 0. Esto no es m´as que una EDO de variables separables, luego asumiendo que v′= 0, tenemos v′′ v′=1−N r=⇒log v′= log r1−N+C=⇒v′=Cr1−N. Con lo que, finalmente, obtenemos v(r) =    blog(r) + csi N= 2, b rN−2+csi N≥3, donde bycson constantes. Aunque el proceso descrito es completamente general, nos restringiremos aqu´ı a los casos de N= 2 y N= 3. As´ı, tomando las constantes adecuadamente, se tiene la siguiente definici´on. Definici´on 1.22 (Soluci´on fundamental de la Ecuaci´on de Laplace para N= 2,3) Se define Soluci´on fundamental de la Ecuaci´on de Laplace bidimensional como G(x) = −1 2πlog(|x|). Soluci´on fundamental de la Ecuaci´on de Laplace para tridimensional como G(x) = 1 4π|x|. En ambos casos, |·|denota la norma eucl´ıdea en R2yR3respectivamente. Vamos a comprobar que, en efecto, estas soluciones verifican la ecuaci´on ∆G=δ. En primer lugar, como la funci´on Gque hemos obtenido s´olo depende de |x|, es interesante ver c´omo queda el Laplaciano en coordenadas polares o esf´ericas, para simplificar el proceso. Recordemos que tenemos la expresi´on Gxixi=d2G dr2 x2 i r2+dG dr 1 r−x2 i r3. Para el caso N= 2, el Laplaciano se escribe ∆G=d2G dr2 x2 1 r2+dG dr 1 r−x2 1 r3+d2G dr2 x2 2 r2+dG dr 1 r−x2 2 r3. Agrupando t´erminos, recordando que x2 1+x2 2=r2, se simplifica la expresi´on anterior, quedando ∆G=d2G dr2x2 1+x2 2 r2+dG dr 2 r−x2 1+x2 2 r3. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 1. Conceptos previos 18 Simplificando, obtenemos ∆G=d2G dr2+dG dr 1 r. Mientras que en el caso N= 3, tenemos ∆G=d2G dr2 x2 1 r2+dG dr 1 r−x2 1 r3+d2G dr2 x2 2 r2+dG dr 1 r−x2 2 r3+d2G dr2 x2 3 r2+dG dr 1 r−x2 3 r3. Agrupando t´erminos, recordando que x2 1+x2 2+x2 3=r2, se simplifica la expresi´on, teni´endose ∆G=d2G dr2x2 1+x2 2+x2 3 r2+dG dr 3 r−x2 1+x2 2+x2 3 r3, luego ∆G=d2G dr2+dG dr 2 r. Veamos c´omo queda entonces para N= 2, con G(r) = −1 2πlog(r). Por lo tanto, ∆G(r) = d2G dr2+dG dr 1 r=1 2πr2−1 2πr 1 r. Esta expresi´on es igual a 0 si r= 0 y tiende a infinito cuando r= 0, luego tenemos δ. Comprob´emoslo para N= 3, con G(r) = 1 4πr. Tenemos ∆G(r) = d2G dr2+dG dr 2 r=2 4πr3−1 4πr2 2 r. Esta expresi´on es igual a 0 si r= 0 y tiende a infinito cuando r= 0, luego de nuevo tenemos δ. Observaci´on 1.23 Los desarrollos anteriores son invariantes bajo traslaciones, esto es, como hemos considerado soluciones radiales, es indiferente considerarlas centradas en el 0 o en un punto s∈RN. Por tanto, sin m´as que considerar r=|x−s|, tenemos las soluciones fundamentales de la ecuaci´on de Laplace bidimensional y tridimensional centradas en una fuente s. Dichas soluciones son de la forma G(x) = −1 2πlog(|x−s|), G(x) = 1 4π|x−s|, para N= 2 y N= 3 respectivamente, y son soluciones de ∆G=δs. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 2 Descripci´on general del m´etodo 2.1. Breve introducci´on hist´orica El M´etodo de Soluciones Fundamentales (MSF), o tambi´en conocido por sus siglas en ingl´es como MFS (Method of Fundamental Solutions), es un m´etodo num´erico que se ide´o inicialmente para la aproximaci´on num´erica de problemas de contorno el´ıpticos, pero su uso se ha generalizado y extendido a otro tipo de problemas por su eficiencia y versatilidad, como veremos en este trabajo. Esta t´ecnica pertenece a la clase de m´etodos num´ericos conocida como m´etodos de frontera o m´etodos “meshless”(sin malla), cuyo nombre viene a decir que no necesitamos un mallado del dominio para la resoluci´on del problema en cuesti´on. Se puede consultar m´as sobre este tipo de t´ecnicas en [10]. Este hecho adquiere una gran utilidad e importancia cuando trabajamos con dominios geom´etricamente complejos, o cuyo interior nos es parcial o totalmente desconocido. En cuanto al recorrido hist´orico que nos lleva a este m´etodo tiene su precedente, como se detalla en [5], en los m´etodos de Trefftz, presentados por Erich Trefftz en 1926 (ver [22]). Esta t´ecnica parte de un problema de contorno el´ıptico, y aproxima su soluci´on mediante una combinaci´on lineal de soluciones particulares de la ecuaci´on. Todo esto se hace bajo la hip´otesis de que dichas combinaciones lineales son densas en el espacio de soluciones. M´as adelante, en 1952, Mergelyan demostr´o (ver [15]) que las funciones holomorfas en dominios simplemente conexos de Cpueden ser aproximadas por polinomios, mientras que si el dominio es conexo pero no simplemente conexo (por ejemplo, m´ultiplemente conexo), se pueden aproximar por funciones racionales. Este trabajo aparece como culminaci´on de toda la labor realizada por Runge, Walsh, Lavrent’ev y Keldysh sobre aproximaciones polin´omicas y racionales desde finales del siglo XIX. As´ı, el M´etodo de Soluciones Fundamentales, presentado en primera instancia por Kupradze y Aleksidze en 1963 (ver [14]), y conocido al principio como M´etodo de Series Generalizadas de Fourier, aparece como una continuaci´on de todo lo antes mencionado, donde las soluciones particulares de la ecuaci´on consideradas para la combinaci´on lineal vienen dadas por soluciones 19 Cap´ıtulo 2. Descripci´on general del m´etodo 20 fundamentales de la misma. Algunos ejemplos de la aplicaci´on del m´etodo m´as all´a de los problemas de contorno el´ıpticos pueden ser los problemas de contorno parab´olicos o hiperb´olicos. Por ejemplo, se puede usar para la ecuaci´on de ondas, como podemos ver en [11]o[12]. Otro ejemplo son los problemas inversos, como se puede ver en [3]. En el presente trabajo, adem´as de trabajar los problemas de contorno el´ıpticos, veremos c´omo se puede generalizar el m´etodo para estas aplicaciones. 2.2. Implementaci´on del m´etodo En esta secci´on realizaremos una descripci´on general del m´etodo en el marco del siguiente problema que se considera de base; es un problema de contorno el´ıptico: (Lu= 0 en Ω, Bu=gsobre ∂Ω,(2.1) donde Les un operador diferencial el´ıptico, ues la funci´on inc´ognita y Ω ⊂RNes un abierto acotado no vac´ıo, cuya frontera es ∂Ω. Por otro lado, la condici´on de contorno viene dada por el operador diferencial lineal By la funci´on gdada. El operador Bexpresa las condiciones de contorno vistas en (1.3), pudiendo ser una combinaci´on de estas, estableci´endolas en zonas disjuntas de la frontera. El objetivo de este m´etodo ser´a obtener una aproximaci´on de la soluci´on de (2.1) mediante una expresi´on de la forma uN(x) = N X j=1 cjG(x, sj),con cj∈R,para j= 1, ..., N. (2.2) Es decir, una combinaci´on lineal de soluciones fundamentales G(x, sj) de la ecuaci´on Lu= 0 centradas en sj, esto es, la funci´on tal que LG(x, sj) = δsj(x), j = 1, ..., N. Seguiremos la referencia [5] para describir la implementaci´on de la primera etapa del m´etodo por su simplicidad. Este texto se basa eminentemente a su vez en [24]. Para implementar el MSF, lo primero es ubicar las llamadas fuentes puntuales sj, con j= 1, ..., N, que son las singularidades de las soluciones fundamentales. Las dos opciones m´as usadas para la ubicaci´on de estas fuentes son las siguientes: Uniformemente distribuidas sobre la frontera de un c´ırculo que contenga a Ω. Sobre una frontera virtual, que denotaremos por ∂Ω′, situada de manera equidistante a ∂Ω. Esta segunda forma generalmente proporciona mejores resultados en dominios con formas b´asicas (ver [24]). TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 2. Descripci´on general del m´etodo 21 Consideraremos la segunda propuesta, de forma que, si Ω es el dominio tratado y ∂Ω su frontera, seguiremos los siguientes pasos: 1. Sean xj,j= 1, ..., N,Npuntos distribuidos uniformemente sobre la frontera ∂Ω. 2. Localizamos el centro geom´etrico o centroide (ver [20]), xc, del dominio Ω. En los casos de dominios regulares, no es m´as que la noci´on cl´asica de centro de una figura geom´etrica. 3. Emplazamos las fuentes puntuales sj,j= 1, ..., N, siguiendo la siguiente expresi´on sj=xj+λ(xj−xc), j = 1, ..., N, donde xjson los puntos tomados sobre la frontera, y λ > 0 un par´ametro escalar que determina la distancia entre la frontera real del dominio, ∂Ω, donde est´an los xj, y la frontera virtual ∂Ω′, donde ubicamos los sj. Un ejemplo de la colocaci´on de las fuentes en un dominio Ω bidimensional dado por una elipse para N= 10 y λ= 1 se puede observar en la Figura 2.1. Figura 2.1: Ejemplo de distribuci´on de las fuentes para el MSF en una elipse. De esta forma, la soluci´on que proporciona el m´etodo, que viene dada por (2.2), queda determinada una vez calculados los coeficientes escalares cj,j= 1, ..., N. Para determinar estos coeficientes, impondremos la condici´on de contorno sobre los puntos xi, con i= 1, ..., N, que hemos escogido sobre la frontera: BuN(xi) = g(xi),∀i= 1, ..., N. N´otese que hemos cambiado ligeramente la notaci´on de los puntos en la frontera, cambiando el sub´ındice jpor ipara poder reescribir lo anterior, en virtud de la linealidad del operador B, TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 28 Hacemos las sustituciones teniendo en cuenta (3.2) y obtenemos lo siguiente: uN(xb1) = −1 2πhc2log pR2+ (R+ 1)2+c3log (2R+ 1) + c4log pR2+ (R+ 1)2i= 1, uN(xb2) = −1 2πhc1log pR2+ (R+ 1)2+c3log pR2+ (R+ 1)2+c4log (2R+ 1)i= 1, uN(xb3) = −1 2πhc1log (2R+ 1) + c2log pR2+ (R+ 1)2+c4log pR2+ (R+ 1)2i= 1, uN(xb4) = −1 2πhc1log pR2+ (R+ 1)2+c2log (2R+ 1) + c3log pR2+ (R+ 1)2i= 1. Si denotamos α:= log pR2+ (R+ 1)2yβ:= log (2R+ 1), el sistema queda como sigue: Ac =b, con A=       0α β α α0α β β α 0α α β α 0        , c ∈R4, b =−2π       1 1 1 1        . La matriz Atiene inversa, luego el sistema es compatible determinado, cuya soluci´on es cj=−2π 2α+β, para todo j= 1,2,3,4. Luego, la soluci´on num´erica que nos da el m´etodo viene dada por: uN(x) = N X j=1 −2π 2α+βG(x, sj),con N= 4. Este ejemplo es extremadamente simple en aras de la f´acil comprensi´on y la posibilidad de realizarlo a mano. Sin embargo, este procedimiento es totalmente generalizable para una funci´on g(x) arbitraria. Podemos tomar tambi´en como par´ametros el n´umero de puntos en la frontera, el centro y radio de la bola y el valor de λ. No obstante, es importante notar que la soluci´on fundamental G(x, s) tiene una singularidad cuando x=s, por tanto hay que tomar λde forma que en todo momento |x−s|no tome valores cercanos a 0. El proceso anterior, que hemos simplificado, se puede generalizar f´acilmente para un problema de la forma (3.1), donde consideraremos Ω = B(xc;R). M´as concretamente, tomamos Npuntos xb1, ..., xbN,distribuidos en la frontera y las correspondientes Nfuentes dadas por sj=xbj+λ(xbj−xc), donde λ > 0 es un par´ametro previamente fijado. Considerando la soluci´on fundamental del Laplaciano para cada fuente puntual, G(x, sj) = −1 2πlog(|x−sj|), j = 1, ..., N, tenemos que la soluci´on dada por el MSF viene dada por uN= N X j=1 cjG(x, sj), TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 29 donde c1, ..., cNson los coeficientes que hacen que uN(xbj) = g(xbj)∀j= 1, ..., N, (3.3) pues el hecho de que G(x, sj) sea soluci´on fundamental ya nos asegura que ∆uN= 0 en Ω, pues ∆G(x, sj) = δsj, y como estamos evaluando en puntos de Ω, con sj/∈Ω, para todo j= 1, ..., N, entonces ∆uN= N X j=1 cj∆G(x, sj) = 0. Por tanto, si escribimos las condiciones (3.3) en forma de sistema lineal, tenemos Ac =b, (3.4) siendo A= (G(xbi, sj))i,j=1,...N , b = (g(xbi))i=1,...N . As´ı, obtenemos los coeficientes dados por el vector c= (c1, ...cN), y con ello obtenemos la soluci´on del problema dada por el MSF. Un problema que aparece en este m´etodo es que la matriz Aque define el sistema (3.4) es una matriz mal condicionada, luego la precisi´on va disminuyendo al resolver el sistema de forma num´erica cuando el n´umero de puntos crece. Esto se puede solucionar con alg´un m´etodo de precondicionamiento de la matriz. Podemos usar, por ejemplo, la regularizaci´on de Tikhonov (ver [16]), es decir, consideramos una peque˜na variaci´on ε > 0, y resolvemos el sistema (AtA+εI)c=Atb, en vez de Ac =b, de forma que podemos considerar Nun n´umero elevado. No obstante, hacer crecer Nprovoca que aumente la complejidad computacional del m´etodo. En este sentido, podemos ver sendos ejemplos en (3.7) y (3.8) de un caso en el que aumentar Nrepercute significativamente en la bondad de la aproximaci´on, hasta un punto en el que se estabiliza, y otro donde no se mejora demasiado aumentando N. Por tanto, dependiendo del caso, es conveniente encontrar un equilibrio. Ejemplo 3.1 Para ver c´omo se comporta el m´etodo para (3.1), consideremos la funci´on arm´onica u(x) = ex1cos(x2). Tomemos N= 100, R= 1, xc= (0,0) y λ= 1.75 como par´ametro para definir las fuentes. Las soluciones aproximada y exacta pueden verse en la Figura 3.2. Asimismo, podemos ver en la Figura 3.3 que el error entre ambas, dado por ue−uN, siendo uela soluci´on exacta y uNla proporcionada por el m´etodo, es considerablemente peque˜no. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 30 (a) Soluci´on con el MSF. (b) Soluci´on exacta. Figura 3.2: Comparaci´on entre la soluci´on exacta y la obtenida con el MSF para la ecuaci´on de Laplace en una bola. Figura 3.3: Error absoluto entre la soluci´on exacta y la obtenida con el MSF para la ecuaci´on de Laplace en una bola. Observaci´on 3.2 Por simplicidad, hemos considerado el problema (3.1) definido en una bola. Sin embargo, el m´etodo y su implementaci´on es totalmente generalizable a cualquier dominio. Para ilustrar este hecho, vamos a considerar un dominio cuyas ecuaciones param´etricas vienen dadas por (x1(θ) = (1 + 0.3 cos(5θ)) cos(θ), x2(θ) = (1 + 0.3 cos(5θ))sen(θ), donde θtoma valores en [0,2π). Consideramos de nuevo g(x) = ex1cos(x2), de forma que la soluci´on del problema es u=g, y podemos implementar el MSF en la regi´on antes descrita. Tomaremos N= 100 y λ= 1.75. Tenemos as´ı una representaci´on de la distribuci´on de las fuentes en la Figura 3.4. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 31 Figura 3.4: Distribuci´on de los puntos sobre la frontera y de las fuentes en el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. Las soluciones aproximada y exacta pueden verse en la Figura 3.5, mientras que el error entre ambas aparece ilustrado en la Figura 3.6. (a) Soluci´on con el MSF. (b) Soluci´on exacta. Figura 3.5: Comparaci´on entre la soluci´on exacta y la obtenida con el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. En todos los ejemplos posteriores que aparecen en el presente trabajo nos restringiremos a tomar como dominio bolas, aunque con sencillas modificaciones del c´odigo MATLAB usado se puede generalizar a cualquier dominio deseado. Estos c´odigos pueden encontrarse en el repositorio enlazado al comienzo del Cap´ıtulo 4. El motivo de esta decisi´on es ilustrar el m´etodo en los casos m´as sencillos, as´ı como de simplificar los procesos computacionales. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 32 Figura 3.6: Error absoluto entre la soluci´on exacta y la obtenida con el MSF para el problema de Laplace en un dominio dado por sus ecuaciones param´etricas. Observaci´on 3.3 Como el tama˜no de Ninfluye considerablemente en la complejidad computacional del m´etodo, es interesante realizar un peque˜no estudio de c´omo var´ıa el error en funci´on del Nelegido. En general, el error decrece conforme Naumenta. Podemos ilustrarlo tomando la funci´on del ejemplo anterior, u(x) = ex1cos(x2), que sabemos que es soluci´on del siguiente problema, construido ad hoc a partir de u: (∆u= 0 en B((0,0); 1), u=ex1cos(x2) sobre ∂B((0,0); 1). Podemos ver la evoluci´on del error entre N= 50 y N= 100 en la Figura 3.7. Figura 3.7: Evoluci´on del error para u(x) = ex1cos(x2) desde N= 50 hasta N= 100. Como se puede observar, el error sigue decreciendo conforme Nva aumentando. Sin embargo, puede haber funciones para las cuales haya un decrecimiento considerable a partir de un N TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 33 cr´ıtico, y luego se produzca una estabilizaci´on del error. En estos casos, es interesante saber d´onde se produce esto, ya que hay un salto considerable en la bondad de la soluci´on con esfuerzo peque˜no a nivel computacional. Para presentar un ejemplo de este hecho, podemos proceder como antes para la funci´on arm´onica u(x) = cos(πx1)senh(πx2) + sen(2πx1) cosh(2πx2). En este caso, vemos que entre N= 50 y N= 55 se produce una reducci´on dr´astica del error, mientras que este apenas decrece entre N= 55 y N= 100. El resultado se muestra en la Figura 3.8. Figura 3.8: Evoluci´on del error para u(x) = cos(πx1) sinh(πx2) + sin(2πx1) cosh(2πx2) desde N= 50 hasta N= 100. N´otese que hemos elegido funciones arm´onicas para poder contar con la soluci´on exacta del problema. 3.1.2. Problema el´ıptico general Consideremos de nuevo una regi´on Ω ⊂R2, y el siguiente problema:          −∆u+a(x)u=h(x) en Ω, u=g(x) sobre Γ1⊂∂Ω, ∂u ∂n(x) = ξ(x) sobre Γ2⊂∂Ω, (3.5) donde nes el vector normal unitario externo a ∂Ω, y Γ1∪Γ2=∂Ω, con Γ1∩Γ2=∅. Tomaremos como Ω una bola, de forma que el vector normal en la frontera est´a perfectamente definido. Buscaremos una soluci´on de la forma u=uP+uH, siendo uPuna soluci´on particular de −∆u+au =hyuHes una soluci´on de la EDP homog´enea ∆u+au = 0. Espec´ıficamente, buscamos combinaciones lineales de funciones radiales en uPy de soluciones TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 34 fundamentales en uH, es decir, u(x) = uP(x) + uH(x) = Nf X j=1 βjF(|x−rj|) + Nb X k=1 αkG(x, sk),(3.6) donde Nfes el n´umero de puntos en el interior del dominio (asociados a la base radial) y Nb es el n´umero de fuentes (asociadas a las soluciones fundamentales) y situadas en el exterior del dominio. Tengamos en cuenta que Nb=Nb1+Nb2, siendo Nb1el n´umero puntos que tomaremos sobre Γ1, para la condici´on de tipo Dirichlet y Nb2el n´umero de puntos que tomaremos sobre Γ2, para la condici´on de tipo Neumann. En (3.6), G(x, sk) es la soluci´on fundamental del Laplaciano, dada por G(x, sk) = −1 2πlog (|x−sk|). Por otro lado, asumiremos que Fes una funci´on de base radial integrada, es decir, una funci´on obtenida por integraci´on a partir de ∆F=f, donde f=f(r) es en este caso una funci´on radial con soporte compacto dada por f(x) = f(|x|) = f(r) =   1−r λ2si r≤λ, 0 si r > λ, (3.7) con λun factor escalar. En nuestro caso, tomemos (ver [13]) F(r) =        r4 16λ2−2r3 9λ+r4 4si r≤λ, 13λ2 144 +λ2 12 log r λsi r > λ. (3.8) Veamos c´omo hemos llegado a esta expresi´on. En R2, para una funci´on radial F(r), el operador Laplaciano tiene la forma ∆F(r) = ∂2 ∂x1 F(r) + ∂2 ∂x2 F(r),donde r=|x| y ∂F ∂xi =F′(r)∂r ∂xi =xi rF′(r), ∂2F ∂x2 i =F′′(r)xi r xi r+∂ ∂xixi r=F′′(r)x2 i r2+F′(r)r−x2 i r r2. Sumando y simplificando, tenemos entonces la siguiente expresi´on del laplaciano de F: ∆F(r) = F′′(r) + 1 rF′(r) = 1 r d dr rdF dr .(3.9) TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 35 Si r≤λ, tenemos f(r) = 1−r λ2, y como ∆F=f, usando la expresi´on de (3.9), tenemos que d dr rdF dr =rf(r) = r1−r λ2=r−2r2 λ+r3 λ2 =⇒rdF dr =Zr−2r2 λ+r3 λ2dr =r2 2−2r3 3λ+r4 4λ2+C1 =⇒dF dr =1 rr2 2−2r3 3λ+r4 4λ2+C1, de donde F(r) = Zr 2−2r2 3λ+r3 4λ2+C1 rdr =r2 4−2r3 9λ+r4 16λ2+C1log(r) + C2. Para que Fno sea singular en r= 0, tomamos C1= 0. La constante C2la tomaremos igualmente nula por simplicidad. Con esto, tenemos F(r) = r2 4−2r3 9λ+r4 16λ2si r≤λ. Si r > λ, entonces f(r) = 0 y d dr rdF dr =rf(r) = 0 =⇒rdF dr =C3=⇒dF dr =C3 r=⇒F(r) = C3log(r) + C4. Por tanto, F(r) = C3log(r) + C4si r > λ. Imponemos continuidad de Fen r=λ: l´ım r→λ−F(r) = λ2 4−2λ3 9λ+λ4 16λ2=λ2 4−2λ2 9+λ2 16 =13λ2 144 , l´ım r→λ+F(r) = C3log(λ) + C4. Por tanto, imponemos C3log(λ) + C4=13λ2 144 . Buscamos que la derivada sea igualmente continua, luego l´ım r→λ−F′(r) = 4λ3 16λ2−6λ2 9λ+λ 2=λ 12, l´ım r→λ+F′(r) = C3 λ. Imponemos entonces C3 λ=λ 12, de donde C3=λ2 12. Sustituyendo, obtenemos λ2 12 log(λ) + C4=13λ2 144 =⇒C4=13λ2 144 −λ2 12 log(λ). Si unimos todo lo que hemos ido deduciendo, llegamos a (3.8). Vamos a ir imponiendo las condiciones que tiene que cumplir la soluci´on que estamos construyendo. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 36 •La EDP en los puntos interiores del dominio ri, con i= 1, ..., Nf: Nf X j=1 βj[−f(|ri−rj|) + a(ri)F(|ri−rj|)] + Nb X k=1 αka(ri)G(ri, sk) = h(ri). Esta expresi´on se tiene, ya que ∆G=δsk, y como todos los puntos riest´an dentro de Ω y las fuentes skfuera, entonces ∆G= 0. Adem´as, ∆F=fpor construcci´on. Luego, en efecto, tomando la aproximaci´on de u(x) que estamos manejando, la expresi´on anterior es equivalente a−∆u+au =hen los puntos interiores. •La condici´on Dirichlet en Γ1⊂∂Ω para los puntos mi∈Γ1⊂∂Ω, donde i= 1, ..., Nb1: Nf X j=1 βjF(|mi−rj|) + Nb X k=1 αkG(mi, sk) = g(mi). •La condici´on Neumann en los puntos zi∈Γ2⊂∂Ω para i= 1, ..., Nb2: Nf X j=1 βj ∂ ∂nF(|zi−rj|) + Nb X k=1 αk ∂ ∂nG(zi, sk) = ξ(zi), donde estamos entendiendo G(x, sk) como funci´on de x, luego ∂ ∂nG(x, sk) = ∇xG·n. Observaci´on 3.4 Para calcular ∂F ∂n =∇F·no∂G ∂n =∇G·n, basta tener en cuenta que la normal exterior a una bola B(xc;R) viene dada por n=x−xc |x−xc|. Por tanto, como las funciones FyGson radiales, tenemos que ∇F=F′(r)x−xc ry∇G=G′(r)x−xc r, luego ∂F ∂n =∇F·n=F′(r)·(x−xc)2 |x−xc|2=F′(r), ∂G ∂n =∇G·n=G′(r)·(x−xc)2 |x−xc|2=G′(r). Obtenemos entonces el siguiente sistema lineal: TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 37 Nf X j=1 βj cij z }| { [−f(|ri−rj|) + a(ri)F(|ri−rj|)] + Nb X k=1 αk dik z }| { a(ri)G(|ri−sk|) = h(ri), i = 1, ..., Nf, Nf X j=1 βj eij z }| { F(|mi−rj|) + Nb X k=1 αk lik z }| { G(|mi−sk|) = g(mi), i = 1, ..., Nb1, Nf X j=1 βj pij z }| { ∂ ∂nF(|zi−rj|) + Nb X k=1 αk qik z }| { ∂ ∂nG(|zi−sk|) = ξ(zi), i = 1, ..., Nb2. (3.10) En consecuencia, el sistema lineal en forma matricial es                      c11 · · · c1Nfd11 · · · d1Nb . . .. . .. . .. . . cNf1· · · cNfNfdNf1· · · dNfNb e11 · · · e1Nfl11 · · · l1Nb . . .. . .. . .. . . eNb11· · · eNb1NflNb1· · · lNb1Nb p11 · · · p1Nfq11 · · · q1Nb . . .. . .. . .. . . pNb21· · · pNb2NfqNb21· · · qNb2Nb                                   β1 . . . βNf α1 . . . αNb              =                      h(r1) . . . h(rNf) g(m1) . . . g(mNb1) ξ(z1) . . . ξ(zNb2)                      ⇐⇒     C D E L P Q    "β α#=    h(r) g(m) ξ(z)    , donde hemos usado la notaci´on de (3.10) para los elementos de la matriz. Observaci´on 3.5 Aunque el procedimiento anterior est´a expresado para un problema el´ıptico con condiciones Dirichlet-Neumann, con simples modificaciones podemos adaptar el desarrollo a un problema con condiciones s´olo de tipo Dirichlet, es decir (−∆u+a(x)u=h(x) en Ω, u=g(x) sobre ∂Ω.(3.11) Todo el desarrollo es v´alido, el ´unico cambio ser´ıa no incluir en el sistema lineal resultante la ´ultima condici´on, es decir, en la matriz por bloques que hemos definido, habr´ıa que quitar los bloques PyQy los puntos donde se eval´ua gestar´ıan distribuidos sobre toda la frontera, quedando "C D E L#"β α#="h(r) g(m)#, siendo la soluci´on de la misma forma u(x) = Nf X j=1 βjF(|x−rj|) + Nb X j=1 αkG(x, sk), TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 44 entonces aproximaremos la segunda derivada, utt =∂ut ∂t , como utt(x, tn+1)≃ut(x, tn+1)−ut(x, tn) k, donde, usando la aproximaci´on que ya hab´ıamos hecho para la derivada primera, tenemos que esta ´ultima expresi´on es equivalente a un+1 −un k−un−un−1 k k=un+1 −2un+un−1 k2, n = 1, ..., M −1. Este esquema se conoce como un esquema de diferencias finitas centradas. Para inicializar el proceso tenemos que tener en cuenta que la expresi´on buscada para utt depende de la iteraci´on actual y de la anterior. Por tanto, necesitamos partir de t2, y conocer u0yu1. En cuanto a u0, tenemos por hip´otesis en (3.16) que u0=u0. Para aproximar u1, nos valdremos de la derivada temporal en el tiempo inicial y de un esquema de Euler expl´ıcito para la derivada temporal, teni´endose as´ı v0=ut(x, t0)≃u1−u0 k=⇒u1≃u0+kv0=u0+kv0. Consideramos entonces la aproximaci´on de utt −∆u=h(x, t) dada por el esquema anterior, es decir, un+1 −2un+un−1 k2−∆un+1 =h(x, tn+1), n = 1, ..., M −1. Despejando, obtenemos la expresi´on 1 k2un+1 −∆un+1 =h(x, tn+1) + 2un−un−1 k2, n = 1, ..., M −1. Si ahora imponemos la condici´on de contorno, conocida la soluci´on del problema unen el instante de tiempo tn, la soluci´on en el instante tn+1 viene dada por el problema    1 k2un+1 −∆un+1 =h(x, tn+1) + 2un−un−1 k2en Ω, un+1 =g(x, tn+1) sobre ∂Ω. (3.17) Este problema no es m´as que un caso espec´ıfico de (3.11), luego ya sabemos resolverlo. Como la soluci´on en el instante inicial es conocida, u(x, 0) = u0(x), y tambi´en es conocida (por aproximaci´on) la soluci´on en t=t1,u1≃u0+kv0, podemos iterar este proceso desde n= 1 en adelante para conseguir la soluci´on del problema en cada instante de tiempo. Ejemplo 3.10 Vamos a probar el m´etodo antes descrito para un problema de tipo (3.16). Consideremos como las funciones que describen el problema h(x, t) = 2π2 100(t+1) cos πx1 10cos πx2 10, g(x, t) = (t+ 1) cos πx1 10cos πx2 10,u0(x) = cos πx1 10cos πx2 10, y v0(x) = u0(x), donde estamos considerando Ω = B((0,0); 1). La soluci´on exacta para este problema es u(x, t) = cos πx1 10cos πx2 10 (t+ 1). Tomamos la distancia entre los puntos del mallado interior como TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 45 dist = 0.2, y Nb= 50. Por otro lado, se ha tomado λ0= 1.75 como par´ametro para definir las fuentes y λ1= 10 como par´ametro para definir la funci´on fen (3.7). El mallado y la distribuci´on de fuentes son los mismos que los antes vistos para el ejemplo de la ecuaci´on del calor, de forma que su representaci´on puede encontrarse en la Figura 3.15. Supondremos de nuevo un horizonte temporal de T= 1, y una partici´on en M= 50 subintervalos. Las soluciones en el tiempo final t= 1 pueden verse en la Figura 3.18. Por otro lado, el error en t= 1, dado por la diferencia entre soluci´on exacta y aproximaci´on, aparece en la Figura 3.19. (a) Soluci´on con el MSF-MSP en t= 1. (b) Soluci´on exacta en t= 1. Figura 3.18: Comparaci´on entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on de ondas con condiciones Dirichlet. Figura 3.19: Error absoluto entre la soluci´on exacta y la obtenida con el MSF-MSP para la ecuaci´on de ondas en el tiempo final. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 46 3.2. Problemas Inversos En esta secci´on vamos a tratar brevemente algunos problemas inversos asociados a los problemas directos vistos en la secci´on anterior. Mientras que los problemas directos consisten en obtener la soluci´on de un problema dado, los problemas inversos tratan de, a partir de un dato adicional de la soluci´on, obtener alg´un par´ametro o dato que desconocemos. Utilizaremos el MSF combinado con t´ecnicas de minimizaci´on y optimizaci´on num´erica como las que podemos encontrar en [21]. 3.2.1. Problemas inversos para EDPs el´ıpticas Consideraremos el problema de contorno el´ıptico con condiciones Dirichlet visto en (3.11), esto es: (Pa)(−∆u+a(x)u=h(x) en Ω, u=g(x) sobre ∂Ω. Llamaremos uaa la soluci´on del problema anterior asociada a la funci´on a=a(x), que hemos denotado como (Pa). La funci´on a(x) es el par´ametro del problema que buscaremos a partir de un dato que suponemos que podemos medir. Como el espacio de las funciones es demasiado grande para buscar, asumiremos que conocemos un conjunto donde se encuentra a, que denotaremos como Aad, es decir, el conjunto de las funciones admisibles. Seg´un el dato de la soluci´on que conozcamos, podemos plantear distintos problemas inversos. Soluci´on conocida en una subregi´on Vamos a suponer que el abierto Ω donde est´a planteado el problema (3.11) es una bola de centro xcy radio R, esto es, Ω = B(xc;R). Si uAes la soluci´on de (PA), donde Aes el par´ametro que queremos calcular, consideraremos conocido el valor de esta funci´on en la bola de centro xcy radio r, con 0 < r < R. Es decir, si D=B(xc;r), consideraremos que η(x) = uA|D(x) es una funci´on conocida. Consideraremos el funcional J=J(a), con J:Aad →R, que trataremos de minimizar para hallar la soluci´on. En este caso, tomaremos como funcional J(a) = 1 2ZD |η(x)−ua(x)|2dx, donde uaes la soluci´on de (Pa), que calcularemos usando el MSF-MSP visto en la Subsecci´on 3.1.2. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 47 As´ı, podremos buscar Acomo J(A) = m´ın a∈Aad J(a). Este problema se puede resolver con t´ecnicas de minimizaci´on num´ericas como la funci´on fmincon de Optimization Toolbox de MATLAB. Para estudiar la robustez del m´etodo, as´ı como ajustarlo m´as a la realidad, podemos plantear el caso en el que el dato η(x) conocido se encuentra levemente perturbado por, por ejemplo, ruido en la medici´on del mismo. Ejemplo 3.11 Tomemos como espacio de las posibles funciones los polinomios de grado menor o igual que uno, esto es, Aad ={a(x) = a1+a2x1+a3x2:aj∈R, j = 1,2,3}. Por simplicidad, identificaremos el polinomio a(x) con el vector a= (a1, a2, a3). Tomaremos como centro xc= (0,0), R= 2 como radio de Ω, y r= 1 como radio de D. Por otro lado, para implementar el MSF-MSP para un problema el´ıptico, tomaremos los par´ametros dist = 0.4 para el mallado interior, Nb= 50 para los puntos sobre la frontera, λ0= 1.75 para la distancia de las fuentes a la frontera y λ1= 10 como par´ametro para definir la funci´on radial f en (3.7). Podemos encontrar una representaci´on del dominio y de las fuentes en la Figura 3.20. Figura 3.20: Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de la subregi´on D. Las funciones que definen el problema (3.11) ser´an h(x) = p|x1+x2|ex2sen(x1) y g(x) = ex2sen(x1). El polinomio que buscamos ser´a A(x) = x1+x2. Tomando a0= (0,0,0) como el vector para inicializar fmincon, tenemos los resultados de la Tabla 3.1 para distintas perturbaciones aleatorias. Observamos que se sigue obteniendo convergencia hacia el valor buscado hasta un ruido aleatorio del 0.1 %. No obstante, si llegamos hasta un 1 %, el polinomio ´optimo que devuelve el TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 48 m´etodo difiere ligeramente del esperado, aunque no demasiado, lo que demuestra que el m´etodo es bastante robusto. No obstante, hemos de tener en cuenta que el a0elegido para inicializar fmincon est´a ya cerca del buscado, lo que mejora la convergencia. Tabla 3.1: Resultados de optimizaci´on con fmincon para diferentes niveles de ruido aleatorio en el problema inverso el´ıptico asociado a un dato en una subregi´on D. Ruido aleatorio ( %) Iteraciones de fmincon Valor objetivo Polinomio ´optimo 0.01 16 1.755547e−10 0.00 + 1.00x1+ 1.00x2 0.1 16 1.471382e−08 0.00 + 1.00x1+ 1.00x2 1 16 1.229420e−06 0.02 + 1.01x1+ 1.01x2 Derivada normal de la soluci´on conocida en un subconjunto de la frontera Consideramos de nuevo el problema (3.11) planteado en Ω = B(xc;R). Tomaremos como γ⊂∂Ω la semicircunferencia superior de ∂Ω. Si uAes la soluci´on de (PA), siendo este el valor que queremos calcular, en este caso consideraremos como funci´on conocida η(x) = ∂uA ∂n (x), x ∈γ, con γ⊂∂Ω, y nel vector normal unitario exterior a ∂Ω. El funcional considerado en este caso para resolver el problema ser´a J(a) = 1 2Zγη(x)−∂ua ∂n (x) 2 dσ, siendo uala soluci´on de (Pa), que podemos calcular con lo visto en la Subsecci´on 3.1.2. Como en el caso anterior, podemos buscar Acomo J(A) = m´ın a∈Aad J(a). Este problema es resoluble fmincon de MATLAB. En este caso, consideraremos de nuevo un cierto porcentaje de ruido aleatorio para comprobar la fiabilidad del m´etodo. Ejemplo 3.12 Tomemos de nuevo como espacio de las posibles funciones los polinomios de grado menor o igual que uno, esto es, Aad ={a(x) = a1+a2x1+a3x2:aj∈R, j = 1,2,3}. Identificamos nuevamente el polinomio a(x) con el vector a= (a1, a2, a3). TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 49 Tomaremos como centro xc= (1,2), R= 1 como radio de Ω. Por otro lado, para implementar el MSF-MSP para un problema el´ıptico, tomaremos los par´ametros dist = 0.2 para el mallado interior, Nb= 50 para los puntos sobre la frontera, λ0= 1.75 para la distancia de las fuentes a la frontera y λ1= 10 como par´ametro para definir fen (3.7). La representaci´on del dominio junto con las fuentes se encuentra en la Figura 3.21. Figura 3.21: Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de γ. Las funciones que definen el problema (3.11) ser´an de nuevo h(x) = p|x1+x2|ex2sen(x1) y g(x) = ex2sen(x1). Consideremos de nuevo A(x) = x1+x2. Tomando a0= (0,0,0) como el vector para inicializar fmincon, tenemos los resultados de la Tabla 3.2 para distintas perturbaciones aleatorias. Tabla 3.2: Resultados de optimizaci´on con fmincon para A(x) = x1+x2y diferentes niveles de ruido aleatorio en el problema inverso el´ıptico asociado a un dato de la derivada normal sobre γ. Ruido aleatorio ( %) Iteraciones de fmincon Valor objetivo Polinomio ´optimo 0.01 16 2.611114e−07 0.00 + 1.00x1+ 1.00x2 0.1 17 1.203192e−04 −0.01 + 1.00x1+ 1.00x2 1 15 9.467925e−03 0.11 + 1.02x1+ 0.93x2 Por otro lado, considerando A(x) = 1+2x1+4x2, de nuevo tomando a0= (0,0,0), obtenemos los resultados de la Tabla 3.3 tomando distintos porcentajes de ruido. Observamos como en ambos casos se obtiene convergencia hacia el valor esperado para un ruido aleatorio del 0.01 %, si consideramos un ruido del 0.1 %, para A(x) = x1+x2se obtiene convergencia a una funci´on que difiere m´ınimamente de la original, mientras que para A(x) = 1+2x1+ 4x2s´ı que hay una diferencia no despreciable. Esto se debe eminentemente a que el a0inicial que hemos elegido est´a mucho m´as cerca del primer polinomio que del segundo. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 50 Finalmente, para una perturbaci´on del 1 %, la diferencia entre el polinomio devuelto por el m´etodo y el esperado es considerable en ambos, aunque con una diferencia m´as marcada en el caso de A(x) = 1 + 2x2+ 4x4por el fen´omeno antes mencionado. Tabla 3.3: Resultados de optimizaci´on con fmincon para A(x) = 1 + 2x1+ 4x2y diferentes grados de ruido aleatorio en el problema inverso el´ıptico asociado a un dato de la derivada normal sobre γ. Ruido aleatorio ( %) Iteraciones de fmincon Valor objetivo Polinomio ´optimo 0.01 18 3.644543e−06 1.00 + 2.00x1+ 4.00x2 0.1 17 2.965464e−04 1.11 + 2.00x1+ 3.95x2 1 18 1.644015e−01 0.86 + 1.97x1+ 4.05x2. 3.2.2. Problemas inversos asociados a la ecuaci´on del calor Vamos a considerar como ejemplo de problema inverso parab´olico el asociado al problema de contorno de la ecuaci´on del calor con condiciones de tipo Dirichlet visto en la Secci´on 3.1.3. Por tanto, consideraremos el problema (3.12), es decir: (Ph)       ut−∆u=h(x, t) en Ω ×(0, T), u=g(x, t) sobre ∂Ω×(0, T), u(x, 0) = u0en Ω. Denotaremos por uha la soluci´on de este problema asociada a h=h(x, t), que ser´a el par´ametro buscado. Al problema asociado a un cierto valor de la funci´on hlo llamaremos (Ph). An´alogamente al caso anterior, denotaremos como Aad al conjunto de los posibles valores que puede tomar h. Es decir, buscaremos h∈Aad. Para este problema volveremos a tomar Ω = B(xc;R), y consideraremos que el dato conocido es la derivada normal de la soluci´on en un subconjunto de la frontera en cada instante de tiempo. Consideraremos γ⊂∂Ω como la semicircunferencia superior. As´ı, dado Hel par´ametro buscado, consideramos uHla soluci´on de (PH), y el dato η(x, t) = ∂uH ∂n (x, t),(x, t)∈γ×(0, T). El funcional considerado para minimizar ser´ıa J(h) = 1 2ZT 0Zγη(x, t)−∂uh ∂n (x, t) 2 dσ dt, donde htoma valores en Aad. Sin embargo, como el c´alculo de uhlo realizaremos mediante el m´etodo de MSF-MSP para la ecuaci´on del calor visto en la Secci´on 3.1.3, lo que conocemos es, para una partici´on temporal TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 51 de (0, T) en M+ 1 instantes de tiempo, con M∈N, la aproximaci´on para cada instante tm, es decir, um h: Ω →R,para m= 0, ..., M, donde um h(x)≃uh(x, tm). Por tanto, el funcional que consideraremos en la pr´actica ser´a J(h) = 1 2 M X m=0 Zγη(x, tm)−∂um h ∂n (x) 2 dσ. De esta forma, como en los casos anteriores, podemos buscar Hcomo J(H) = m´ın h∈Aad J(h), problema que resolveremos con fmincon de MATLAB. Al igual que en las situaciones planteadas anteriormente vamos a suponer un cierto porcentaje de ruido aleatorio para comprobar la robustez del m´etodo. Ejemplo 3.13 Tomemos como espacio de las posibles funciones como las constantes, es decir, Aad ={h∈R}. Tomaremos como centro xc= (0,0), R= 1 como radio de Ω. Por otro lado, para implementar el MSF-MSP para un problema el´ıptico, tomaremos los par´ametros dist = 0.4 para el mallado interior, Nb= 50 para los puntos sobre la frontera, λ0= 1.75 para la distancia de las fuentes a la frontera y λ1= 10 como par´ametro para definir las funciones radiales. Podemos encontrar una representaci´on del dominio espacial y las fuentes en la Figura 3.20. Figura 3.22: Distribuci´on de los puntos interiores, sobre la frontera, de las fuentes y representaci´on de γ. En cuanto al tiempo, consideramos un horizonte temporal de T= 1, y como par´ametro de la partici´on tomamos M= 20. TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 3. Ejemplos de Aplicaciones 52 Las funciones que definen el problema (3.12) ser´an g(x, t) = e−tsen(πx1)sen(πx2) y u0(x) = sen(πx1)sen(πx2). El valor de la constante hbuscada ser´a H= 10. Tomando h0= 0 como el valor para inicializar fmincon, obtenemos los resultados de la Tabla 3.4 para distintas perturbaciones aleatorias. Observamos como se sigue obteniendo convergencia hacia el valor buscado hasta un ruido aleatorio del 1 %, en este caso tenemos que llegar hasta una perturbaci´on del 10 % para encontrar que el valor difiere del esperado, pero hemos de tener en cuenta que el espacio Aad es bastante reducido en este caso, y que la distancia entre h0yHno es demasiado elevada. Tabla 3.4: Resultados de optimizaci´on con fmincon para diferentes ruidos aleatorios en el problema inverso parab´olico de la ecuaci´on del calor. Ruido aleatorio ( %) Iteraciones de fmincon Valor objetivo Valor ´optimo 0.01 3 3.109938e−06 10.00 0.1 4 2.582204e−04 10.00 1 4 2.809573e−02 10.00 10 4 2.772585e+ 00 9.96 TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4 Implementaci´on en MATLAB En este cap´ıtulo se exponen algunos de los c´odigos que se han usado para resolver los problemas con los que se ha ido ilustrando este trabajo. La explicaci´on de los mismos se encuentra dentro del c´odigo como comentarios. Todos los c´odigos de MATLAB de los m´etodos utilizados en este trabajo, as´ı como algunos ejemplos de aplicaci´on de los mismos, se pueden encontrar en el siguiente repositorio en Github: https://github.com/dpino03/Metodo-de-Soluciones-Fundamentales-con-MATLAB.git. A continuaci´on se exponen algunos de mayor importancia. 4.1. C´odigo del problema el´ıptico con condiciones DirichletNeumann Implementaci´on del MSF-MSP para el problema el´ıptico bidimensional en una bola con condiciones Dirichlet-Neumann function [u] = MSF_MSPDirNeu(R,xc,dist,Nb1,Nb2,lambda0,lambda1,a,h,g,xi,u_exacta) % ------------------------------------------------------------------------- % Consideramos un problema de la forma: % % -Lap(u) + a*u = h en B(0;R) % u = g sobre Gamma1 contenida en la frontera de B(0;R) % grad(u)*n = xi sobre Gamma2 contenida en la frontera de B(0;R) % % Tomaremos como Gamma1 la semicircunferencia superior de la frontera de la % bola y como Gamma2 la semicircunferencia inferior. % % Entradas: % % R = Radio de la bola. % xc = Centro de la bola (vector 2x1). % dist = Distancia entre los puntos del mallado regular interior para el 53 Cap´ıtulo 4. Implementaci´on en MATLAB 60 for j = 1:Ntheta x_eval = [X1(i,j);X2(i,j)]; % Punto en el que evaluamos u U(i,j) = u(x_eval); end end % ------------------------------------------------------------------------- % Evaluamos la soluci´on exacta en el mallado. % ------------------------------------------------------------------------- ue = zeros(size(X1)); for i = 1:Nr for j = 1:Ntheta x_eval = [X1(i,j);X2(i,j)]; % Punto en el que evaluamos u_exacta ue(i,j) = u_exacta(x_eval); end end % ------------------------------------------------------------------------- % Graficamos la soluci´on del m´etodo y exacta, y tomamos el error absoluto % y relativo en cada punto. El error que tomaremos como referencia ser´a la % norma infinito del error relativo. % ------------------------------------------------------------------------- % ----- Gr´aficas ----- figure(1) surf(X1,X2,U); xlabel(’x_1’); ylabel(’x_2’); zlabel(’u(x_1,x_2)’); title(’Soluci´on con MSF-MSP’) shading interp; % Para suavizar la visualizaci´on de la superficie colorbar; % A~nadir barra de colores figure(2) surf(X1,X2,ue); xlabel(’x_1’); ylabel(’x_2’); zlabel(’u(x_1,x_2)’); title(’Soluci´on exacta’) shading interp; colorbar; figure(3) surf(X1,X2,ue-U); xlabel(’x_1’); ylabel(’x_2’); zlabel(’Error absoluto’); title(’Error absoluto del MSF-MSP respecto de la soluci´on exacta’) TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 61 shading interp; colorbar; err_abs_max = max(max(abs(ue-U))); % Como es una matriz, hay que aplicar dos veces m´aximo. disp([’El error absoluto m´aximo es ’,num2str(err_abs_max)]); % ------------------------------------------------------------------------- % Graficar puntos del mallado, Gamma1, Gamma2 y fuentes. % ------------------------------------------------------------------------- figure(4); hold on; % 1. Puntos del mallado. scatter(puntos_campo(1,:), puntos_campo(2,:), 10, ’b’, ’filled’); % 2. Puntos de Gamma1. scatter(puntos_Gamma1(1,:), puntos_Gamma1(2,:), 30, ’r’, ’filled’); % 3. Puntos de Gamma2. scatter(puntos_Gamma2(1,:), puntos_Gamma2(2,:), 30, ’g’, ’filled’); % 4. Fuentes fuera de la bola. scatter(fuentes_Gamma1(1,:), fuentes_Gamma1(2,:), 50, ’m’, ’x’, ’LineWidth’, 1.5); scatter(fuentes_Gamma2(1,:), fuentes_Gamma2(2,:), 50, ’m’, ’x’, ’LineWidth’, 1.5); % Dibujar la frontera de la bola theta_circ = linspace(0, 2*pi, 100); plot(R*cos(theta_circ) + xc(1), R*sin(theta_circ) + xc(2), ’k-’, ’LineWidth’, 1.5); % Configuraci´on de la gr´afica xlabel(’x_1’, ’Color’, ’k’, ’FontWeight’, ’bold’); ylabel(’x_2’, ’Color’, ’k’, ’FontWeight’, ’bold’); title(’Distribuci´on de puntos y fuentes’, ’Color’, ’k’, ’FontWeight’, ’bold’); legend(’Puntos interiores r_i ’, ’Puntos m_i sobre \Gamma_1 (Dirichlet)’, ... ’Puntos z_i sobre \Gamma_2 (Neumann)’, ’Fuentes’, ’Location’, ’best’); grid on; axis equal; hold off; end TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 62 4.2. C´odigo de la ecuaci´on del calor con condiciones Dirichlet Implementaci´on del MSF para el problema del calor bidimensional en una bola con condiciones Dirichlet function [u] = MSFCalor(R,xc,T,M,u0,dist,Nb,lambda0,lambda1,h,g,u_exacta) % ------------------------------------------------------------------------- % Consideramos un problema de la forma: % % u_t - Lap(u) = h(x,t) en B(xc;R)x(0,T] % u = g(x,t) sobre fr(B(xc;R))x(0,T] % u(x,0) = u0 en B(xc;R) % % Entradas: % % R = Radio de la bola. % xc = Centro de la bola (vector 2x1). % T = Tiempo final. % M = N´umero de subintervalos en los que dividiremos [0,T]. % u0 = Valor de la funci´on en el instante inicial t = 0. % dist = Distancia entre los puntos del mallado regular interior para el % c´alculo de la soluci´on particular u_P en cada problema el´ıptico. Se usar´a % al invocar la funci´on MSF_MSPEliptDir. % Nb = N´umero de puntos de la frontera para aplicar el MSF en cada problema % el´ıptico. Se usar´a al invocar la funci´on MSF_MSPEliptDir. % lambda0 = Factor para ubicar las fuentes fuera de la bola. Se usar´a al % invocar la funci´on MSF_MSPEliptDir. % lambda1 = Constante que aparece en la definici´on de las funciones % radiales f y F, relacionadas por Lap(F) = f. Se usar´a al invocar la funci´on % MSF_MSPEliptDir. % h = Funci´on que describe la EDP. % g = Funci´on que describe las condiciones de Dirichlet en la frontera para % cada tiempo t. % u_exacta = Soluci´on exacta para realizar la comparaci´on. % % Salida: % % Aproximaci´on de la soluci´on u(x) por el MSF-MSP y representaci´on gr´afica. % Comparaci´on de la soluci´on del m´etodo con la soluci´on exacta en cada % tiempo. % % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- % En primer lugar, vamos a realizar la partici´on uniforme del intervalo % [0,T] en los N subintervalos, es decir, tomar N+1 puntos. % ------------------------------------------------------------------------- t = linspace(0,T,M+1); % Dividimos [0,T] en N+1 puntos k = T/M; % Este es el paso, pues la partici´on es uniforme TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 63 % ------------------------------------------------------------------------- % Definimos ahora un bucle que nos va a devolver la soluci´on en cada % instante de tiempo, resolviendo el problema el´ıptico visto en la parte % te´orica en cada instante de tiempo. % ------------------------------------------------------------------------- % Definimos u como una celda que ir´a almacenando en cada entrada la % soluci´on en el instante de tiempo correspondiente como funci´on an´onima. u = cell(1,M+1); % Inicializamos con la soluci´on en el instante inicial t=0. u{1} = u0; % Definimos el siguiente valor como escalar y como funci´on an´onima, ya que % en diferentes contextos en el bucle lo necesitaremos de una forma u otra. k1 = 1/k; a = @(x) k1; % --- Desarrollo del bucle --- for n = 1:M un = u{n}; tn1 = t(n+1); h1 = @(x) h(x,tn1); h_elipt = @(x) h1(x) + k1*un(x); g1 = @(x) g(x,tn1); [uu] = MSF_MSPEliptDir(R,xc,dist,Nb,lambda0,lambda1,a,h_elipt,g1); u{n+1} = uu; end % ------------------------------------------------------------------------- % Seguiremos con la representaci´on gr´afica de la soluci´on. Lo haremos de % forma que veamos c´omo var´ıa la soluci´on a lo largo del tiempo. % ------------------------------------------------------------------------- % --- Construcci´on del mallado --- % N´umero de puntos de evaluaci´on en el radio y en el ´angulo. Nr = 70; Ntheta = 70; % Creamos una malla en coordenadas polares dentro de la bola. r = linspace(0,R,Nr); % Radios desde 0 hasta R theta_mesh = linspace(0,2*pi,Ntheta); % ´ Angulos de 0 a 2pi [Rad,Theta] = meshgrid(r,theta_mesh); % Malla en coordenadas polares TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 64 % Convertimos la malla a coordenadas cartesianas. X = Rad.*cos(Theta) + xc(1); Y = Rad.*sin(Theta) + xc(2); % ------------------------------------------------------------------------- % Vamos a ir construyendo a la vez la representaci´on de la soluci´on, de la % soluci´on exacta y del error, ya que el proceso es an´alogo. Es por esto, % para no doblar los bucles, haremos los desarrollos simult´aneamente. % ------------------------------------------------------------------------- % Inicializamos la soluci´on, el error y la soluci´on exacta. U = nan(size(X)); E = nan(size(X)); Ue = nan(size(X)); % Graficamos la soluci´on, soluci´on exacta y el error. for n = 1:M+1 for i = 1:size(X,1) for j = 1:size(X,2) U(i,j) = u{n}([X(i,j);Y(i,j)]); % Soluci´on en cada punto de la malla E(i,j) = abs(U(i,j)-u_exacta([X(i,j);Y(i,j)],t(n))); % Error en cada punto de la malla Ue(i,j) = u_exacta([X(i,j);Y(i,j)],t(n)); % Soluci´on exacta en cada punto de la malla end end % Gr´afica de la la soluci´on y el error para el tiempo t(n). figure(1) surf(X,Y,U,’EdgeColor’,’none’); % Gr´afica 3D sin bordes colorbar; % Barra de color title([’Soluci´on del m´etodo en t = ’,num2str(t(n))]); % T´ıtulo xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); % Etiquetas de los ejes axis(’equal’,’tight’); % Configuraci´on de los ejes view(3); % Vista en 3D shading(’interp’); % Sombreado por interpolaci´on pause(0.1) % Pausa de 0.1 segundos figure(2) surf(X,Y,Ue,’EdgeColor’,’none’); colorbar; title([’Soluci´on exacta en t = ’,num2str(t(n))]); xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); axis(’equal’,’tight’); view(3); shading(’interp’); pause(0.1); figure(3) TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 65 surf(X,Y,E,’EdgeColor’,’none’); colorbar; title([’Error absoluto en t = ’,num2str(t(n))]); xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); axis(’equal’,’tight’); view(3); shading(’interp’); pause(0.1) end end 4.3. C´odigo del problema de ondas con condiciones Dirichlet Implementaci´on del MSF para el problema de ondas bidimensional en una bola con condiciones Dirichlet function [u] = MSFOndas(R,xc,T,M,u0,v0,dist,Nb,lambda0,lambda1,h,g,u_exacta) % ------------------------------------------------------------------------- % Consideramos un problema de la forma: % % u_tt - Lap(u) = h(x,t) en B(xc;R)x(0,T) % u = g(x,t) sobre fr(B(xc;R))x(0,T) % u(x,0) = u0 en B(xc;R) % u_t(x,0) = v0 en B(xc;R) % % Entradas: % % R = Radio de la bola. % xc = Centro de la bola (vector 2x1). % T = Tiempo final. % M = N´umero de subintervalos en los que dividiremos [0,T]. % u0 = Valor de la funci´on en el instante inicial t = 0. % v0 = Valor de la derivada temporal en el instante inicial t = 0. % dist = Distancia entre los puntos del mallado regular interior para el % c´alculo de la soluci´on particular u_P en cada problema el´ıptico. Se usar´a % al invocar la funci´on MSF_MSPEliptDir. % Nb = N´umero de puntos de la frontera para aplicar el MSF en cada problema % el´ıptico. Se usar´a al invocar la funci´on MSF_MSPEliptDir. % lambda0 = Factor para ubicar las fuentes fuera de la bola. Se usar´a al % invocar la funci´on MSF_MSPEliptDir. % lambda1 = Constante que aparece en la definici´on de las funciones % radiales f y F, relacionadas por Lap(F) = f. Se usar´a al invocar la funci´on % MSF_MSPEliptDir. % h = Funci´on que describe la EDP. % g = Funci´on que describe las condiciones de Dirichlet en la frontera para % cada tiempo t. % u_exacta = Soluci´on exacta para realizar la comparaci´on. % TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 66 % Salida: % % Aproximaci´on de la soluci´on u(x) por el MSF-MSP y representaci´on gr´afica. % Comparaci´on de la soluci´on del m´etodo con la soluci´on exacta. % % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- % En primer lugar, vamos a realizar la partici´on uniforme del intervalo % [0,T] en los N subintervalos, es decir, tomar N+1 puntos. % ------------------------------------------------------------------------- t = linspace(0,T,M+1); % Dividimos [0,T] en N+1 puntos k = T/M; % Este es el paso, pues la partici´on es uniforme % ------------------------------------------------------------------------- % Definimos ahora un bucle que nos va a devolver la soluci´on en cada % instante de tiempo, calculado mediante MSF_MSPEliptico. % ------------------------------------------------------------------------- % Definimos u como una celda que ir´a almacenando en cada entrada la % soluci´on en el instante de tiempo correspondiente como funci´on an´onima. u = cell(1,M+1); % Inicializamos con la soluci´on en el instante inicial t=0. u{1} = u0; % La soluci´on en t_2 la aproximamos usando la derivada en el instante % inicial seg´un se ha visto en el desarrollo te´orico. u2 = @(x) u0(x) + k*v0(x); u{2} = u2; % Definimos el siguiente valor como escalar y como funci´on an´onima, ya que % en diferentes contextos en el bucle lo necesitaremos de una forma u otra. k2 = 1/k^2; a = @(x) k2; % --- Desarrollo del bucle --- for n = 2:M un_1 = u{n-1}; un = u{n}; tn1 = t(n+1); h1 = @(x) h(x,tn1); h_elipt = @(x) h1(x) + k2*(2*un(x) - un_1(x)); g1 = @(x) g(x,tn1); TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 67 [uu] = MSF_MSPEliptDir(R,xc,dist,Nb,lambda0,lambda1,a,h_elipt,g1); u{n+1} = uu; end % ------------------------------------------------------------------------- % Seguiremos con la representaci´on gr´afica de la soluci´on. Lo haremos de % forma que veamos c´omo var´ıa la soluci´on a lo largo del tiempo. % ------------------------------------------------------------------------- % --- Construcci´on del mallado --- % N´umero de puntos de evaluaci´on en el radio y en el ´angulo Nr = 70; Ntheta = 70; % Creamos una malla en coordenadas polares dentro de la bola r = linspace(0,R,Nr); % Radios desde 0 hasta R theta_mesh = linspace(0,2*pi,Ntheta); % ´ Angulos de 0 a 2pi [Rad,Theta] = meshgrid(r,theta_mesh); % Malla en coordenadas polares % Convertimos la malla a coordenadas cartesianas. X = Rad.*cos(Theta) + xc(1); Y = Rad.*sin(Theta) + xc(2); % ------------------------------------------------------------------------- % Vamos a ir construyendo a la vez la representaci´on de la soluci´on y del % error, ya que el proceso es an´alogo. Es por esto, para no doblar los % bucles, haremos los dos desarrollos simult´aneamente. % ------------------------------------------------------------------------- % Inicializamos la soluci´on, el error y la soluci´on exacta. U = nan(size(X)); E = nan(size(X)); Ue = nan(size(X)); % Graficamos la soluci´on, soluci´on exacta y el error. for n = 1:M+1 for i = 1:size(X,1) for j = 1:size(X,2) U(i,j) = u{n}([X(i,j);Y(i,j)]); % Soluci´on en cada punto de la malla E(i,j) = abs(U(i,j)-u_exacta([X(i,j);Y(i,j)],t(n))); % Error en cada punto de la malla Ue(i,j) = u_exacta([X(i,j);Y(i,j)],t(n)); % Soluci´on exacta en cada punto de la malla end end % Gr´afica de la la soluci´on y el error para el tiempo t(n). TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 68 figure(1) surf(X,Y,U,’EdgeColor’,’none’); % Gr´afica 3D sin bordes colorbar; % Barra de color title([’Soluci´on del m´etodo en t = ’,num2str(t(n))]); % T´ıtulo xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); % Etiquetas de los ejes axis(’equal’,’tight’); % Configuraci´on de los ejes view(3); % Vista en 3D shading(’interp’); % Sombreado por interpolaci´on pause(0.1); % Pausa de 0.1 segundos figure(2) surf(X,Y,Ue,’EdgeColor’,’none’); colorbar; title([’Soluci´on exacta en t = ’,num2str(t(n))]); xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); axis(’equal’,’tight’); view(3); shading(’interp’); pause(0.1); figure(3) surf(X,Y,E,’EdgeColor’,’none’); colorbar; title([’Error absoluto en t = ’,num2str(t(n))]); xlabel(’x_1’);ylabel(’x_2’);zlabel(’u(x)’); axis(’equal’,’tight’); view(3); shading(’interp’); pause(0.1) end end 4.4. C´odigo del problema el´ıptico inverso asociado a un dato en una subregi´on Implementaci´on del MSF para el problema del calor bidimensional en una bola con condiciones Dirichlet function [a_opt] = PbInversoD(R,r,xc,dist,Nb,lambda0,lambda1,A,h,g,pert,a0) % ------------------------------------------------------------------------- % Consideramos un problema (P), de la forma: % % -Lap(u) + a*u = h en B(0;R) % u = g sobre fr(B(0;R)) % % Sea u_a el valor de la soluci´on del problema asociado al valor a = a(x). % TFG David Pino Molero. Grado en Matem´aticas Cap´ıtulo 4. Implementaci´on en MATLAB 69 % Consideraremos a = A(x) como el valor de a buscado. Dada la soluci´on u_A % asociada a la soluci´on de (P_a) para este valor de a, realizaremos una % medici´on en D, la bola de centro xc y radio r, del valor de u_A. % % El espacio donde consideraremos a(x) ser´an los polinomios de grado 1. % % Entradas: % % R = Radio de la bola. % r = Radio del dominio donde se realiza la medici´on. % xc = Centro de la bola (vector 2x1). % dist = Distancia entre los puntos del mallado regular interior para el % c´alculo de la soluci´on particular u_P. % Nb = N´umero de puntos sobre la frontera para la condici´on Dirichlet % lambda0 = Factor para ubicar las fuentes fuera de la bola. % lambda1 = Constante que aparece en la definici´on de las funciones % radiales f y F, relacionadas por Lap(F) = f. % A = Funci´on que describe la EDP para la que realizaremos la medici´on. % h = Funci´on que describe la EDP. % g = Funci´on que describe las condiciones Dirichlet en la frontera. % pert = Porcentaje de ruido en la medici´on. Por ejemplo, pert = 1.e-2 es % una perturbaci´on aleatoria del 1%. % a0 = Aproximaci´on de los coeficientes del polinomio buscado, para poder % aplicar de forma eficiente fmincon. % % Salida: % % Valor del polinomio a(x) ´optimo. % % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- % Vamos a resolver en primer lugar el problema (P) para el valor a = A(x) % sobre el que realizaremos la medici´on. % ------------------------------------------------------------------------- [u_A] = MSF_MSPEliptDir(R,xc,dist,Nb,lambda0,lambda1,A,h,g); % Definimos ahora la funci´on eta(x) como el valor de u_A restringido a la % bola B(xc;r). Como u_A est´a definido en todo B(xc;R), basta considerar % eta(x) como u_A, y evaluarla en puntos de B(xc;r): eta = @(x) u_A(x); % Esta ser´a nuestra medici´on en D. % Vamos a considerar una perturbaci´on de la medici´on, es decir, % consideraremos la funci´on eta modificada aleatoriamente, lo que % representa el ruido en la medici´on. Para considerar una medici´on exacta, % basta considerar pert = 0. El par´ametro pert determinar´a el porcentaje de % ruido aleatorio. TFG David Pino Molero. Grado en Matem´aticas