Trabajo de Grado - Página...

REPÚBLICA BOLIVARIANA DE VENEZUELA

UNIVERSIDAD DEL ZULIA

FACULTAD DE INGENIERÍA

DIVISIÓN DE POSTGRADO

PROGRAMA DE POSTGRADO EN MATEMÁTICA APLICADA

SOLUCIÓN NUMÉRICA DE LA ECUACIÓN DE POISSON EN LA ESFERA S2

MEDIANTE LOS MÉTODOS DE GALERKIN Y DE FOURIER

Trabajo de Grado presentado ante la Ilustre Universidad del Zulia para optar al GradoAcadémico de:

MAGÍSTER SCIENTIARUM EN MATEMÁTICA APLICADA

Autor: Nithal El MejmissaniTutor: Jorge Guíñez Matamala

Maracaibo, Septiembre de 2006

APROBACIÓN

Este jurado aprueba el Trabajo de Grado titulado SOLUCIÓN NUMÉRICA DE LAECUACIÓN DE POISSON EN LA ESFERA S2 MEDIANTE LOS MÉTODOS DE GALER-KIN Y DE FOURIER que Nithal El Mejmissani, C.I.: 14.083.947 presenta ante el ConsejoTécnico de la División de Postgrado de la Facultad de Ingeniería en cumplimiento delArticulo 51, Parágrafo 51.6 de la Sección Segunda del Reglamento de Estudios para Gra-duados de la Universidad del Zulia, como requisito para optar al Grado Académico deMAGÍSTER SCIENTIARUM EN MATEMÁTICA APLICADA

Coordinador del JuradoJorge Guíñez

C. I.: 13.974.954

Atilio MorilloC. I.: 3.117.961

Robert QuinteroC.I.: 7.935.866

Directora de la División de PostgradoProfesora Cateryna Aiello

Maracaibo, Septiembre de 2006

El Mejmissani Sukkar, Nithal, SOLUCIÓN NUMÉRICA DE LA ECUACIÓN DE POISSONEN LA ESFERA S2 MEDIANTE LOS MÉTODOS DE GALERKIN Y DE FOURIER. (2006)Trabajo de Grado. Universidad del Zulia. Facultad de Ingeniería. División de Postgrado.Maracaibo, Tutor: Dr. Jorge Guíñez.

RESUMEN

En este trabajo de grado se aproxima la solución de la ecuación de Poisson en una es-fera de dimensión 2 usando el método de Galerkin de elementos finitos y el método deFourier fundado en la teoría de los valores propios del operador Laplace-Beltrami, dondecada método es implementado en una librería de software. Adicionalmente la tecnolo-gía de software desarrollada podrá ser usada con algunas adaptaciones al problema deaproximar numéricamente la solución de ecuaciones diferenciales que tienen dominio enuna geometría dada. Como reseña, existe un importante esfuerzo a nivel mundial parainvestigar la simulación numérica del clima y su pronóstico que gira en torno al desarrollode algoritmos y métodos numéricos para resolver Ecuaciones Diferenciales Parciales enuna geometría esférica.

Palabras Clave: Esfera, Triangulación, Método de Galerkin, Polinomios Armónicos, Méto-do de Fourier.

E-mail del Autor: [email protected]

El Mejmissani Sukkar, Nithal, NUMERICAL SOLUTION FOR THE POISON’S EQUA-TION ON S2 USING GALERKIN’S AND FOURIER’S METHOD. (2006) Trabajo de Gra-do. Universidad del Zulia. Facultad de Ingeniería. División de Postgrado. Maracaibo, Tu-tor: Dr. Jorge Guíñez.

ABSTRACT

In this grade work we approximate the solution of Poisson’s equation on the sphere of di-mension 2 using either the Galerkin’s finite element method as the Fourier method foundedon the theory of the Laplace-Beltrami operator, where both methods are implemented ina software library. Additionally the software technology developed will be used, with someadaptations, to the problem of approximate numerically the solution of partial differentialequations over any given geometry. As remark, there exist an important effort at worldlevel to investigate the climate numerical simulation and its forecast that evolve around ofthe development of algorithms and numerical methods to solve Partial Differential Equa-tions on spherical geometries.

Key Words: Sphere Triangulation, Galerkin’s Method, Harmonic Polynomials, FourierMethod.

Author’s e-mail: [email protected]

Índice general

RESUMEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

I. Introducción 11

II. Solución de la ecuación diferencial de Poisson en la esfera S2 por el Métodode elementos finitos de Galerkin 12II.1. Problema de Dirichlet y la ecuación de Poisson . . . . . . . . . . . . . . . . 12II.2. El método estándar de elementos finitos de Galerkin . . . . . . . . . . . . . 13

II.2.1. Método de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . 15II.3. Triangulación de la esfera . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

II.3.1. Estructura de los datos para un dominio discretizado . . . . . . . . . 17II.4. Programación de (∇Φi · ∇Φj) . . . . . . . . . . . . . . . . . . . . . . . . . . 19

II.4.1. Cálculo del gradiente ∇Φj de la función nodal para un triángulo τj. . 19II.4.2. Cálculo del producto escalar de los gradientes ∇F y ∇G . . . . . . 21II.4.3. Cálculo de (∇Φi · ∇Φj) . . . . . . . . . . . . . . . . . . . . . . . . . 22

II.5. Estimación del error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24II.6. Cálculo de (f, Φk) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25II.7. Ejemplo (comparando una solución real y una calculada) . . . . . . . . . . 25

II.7.1. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27II.8. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

II.8.1. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

III. Teoría de armónicos esféricos 33III.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33III.2. Ordenamiento de una base de Pk . . . . . . . . . . . . . . . . . . . . . . . 38III.3. Base para Hk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38III.4. Ejemplo (comparando una solución real y una calculada) . . . . . . . . . . 39

IV. Conclusiones 47

5

A. Implementación del método de elementos finitos 48

B. Implementación del método de Fourier para una base en el espacio de poli-nomios armónicos homogéneos restringidos a la esfera. 75B.1. Referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75B.2. Listados de la librería . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

C. Descripción y autirizaciones de uso de los programas informaticos emplea-dos en está tesis concedido por el(los) titular(es) del derecho de autor 93C.1. GNU GPL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93C.2. Sparse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93C.3. Gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94C.4. LYX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95C.5. LATEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Índice de figuras

II.1. Mallas esféricas de triángulos. Gen(k) para k = 0, 1, 2, 3, 4, 5 . . . . . . . 16II.2. Discretización del dominio Ω . . . . . . . . . . . . . . . . . . . . . . . . . . 17II.3. Discretización del dominio Ω: primera iteración. . . . . . . . . . . . . . . . . 18II.4. Transformación θ : Ω ⊂ R2 → R3 . . . . . . . . . . . . . . . . . . . . . . . . 20II.5. Caso particular de un elemento triangular en el plano. . . . . . . . . . . . . 22II.6. a) Solución real, b) Solución calculada para una malla de 1280 triángulos

generada con gen(3), c) Desviación entre las soluciones calculada y real. . 27II.7. a) Distribución sobre la esfera, b) Distribución con integral nula. . . . . . . . 29II.8. Solución calculada por el método de Galerkin para la ecuación diferencial

II.21. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

III.1. a) Solución real, b) Solución calculada para un polinomio de tercer gradosobre una malla de 5120 triángulos generada con gen(3), c) Desviaciónentre las soluciones calculada y real. . . . . . . . . . . . . . . . . . . . . . . 46

7

Índice de cuadros

II.1. Integración numérica cúbica. . . . . . . . . . . . . . . . . . . . . . . . . . . 25II.3. Resultados usando la función “solveGalerkinOnThe_tSurface” sobre mallas

de triángulos esférica, t≡triángulos y error por diferencia con la solución real. 28II.4. Solución de la ecuación II.21 por el método de elementos finitos de Galerkin. 30

8

Índice de listas

II.1. Programa en C para resolver la ecuación diferencial por el método de ele-mentos finitos de Galerkin. Este código está en el archivo fuente galerkin/-prueba01.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

II.2. Programa en C para resolver la ecuación diferencial por el método de ele-mentos finitos de Galerkin. Este código está en el archivo fuente galerkin/-prueba02.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

III.1. Programa en C para resolver una ecuación diferencial por el método Fourierpara una base de polinomios armónicos esféricos. Este código está en elarchivo fuente fourier/prueba01.c . . . . . . . . . . . . . . . . . . . . . . . . 39

III.2. Resultados producido por el program prueba01 . . . . . . . . . . . . . . . . 41A.1. Archivo de cabecera que contiene la declaración de las funciones para

el usuario pMedio() integrateOver() area(). Este código está en el archivofuente galerkin/tSurface.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

A.2. Módulo en C que define las funciones para el usuario pMedio() integrateO-ver() area(). Este código está en el archivo fuente galerkin/integrate.c . . . 49

A.3. Archivo de cabecera que contiene la declaración de las funciones para elusuario f() u() solveGalerkinOnThe tSurface(). Este código está en el archi-vo fuente galerkin/solveGalerkin.h . . . . . . . . . . . . . . . . . . . . . . . 51

A.4. Módulo en C para resolver una ecuación diferencial por el método de ele-mentos finitos de Galerkin. Este código está en el archivo fuente galerkin/-solveGalerkin.c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

A.5. Archivo de cabecera que contiene la declaración de las funciones parael usuario gen() tSurfacePlot() plotSolutionOver() y la estructura de datostSurface. Este código está en el archivo fuente galerkin/tSurface.h . . . . . 58

A.6. Módulo en C que contiene la definición de las funciones necesarias paragenerar una malla de triángulos esférica. Este código está en el archivofuente ../galerkin/tSurface.c . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

A.7. Archivo de proyecto galerkin/Makefile . . . . . . . . . . . . . . . . . . . . . 73

9

B.1. Archivo de cabecera que contiene la declaración de todos los datos y todaslas funciones usadas en la librería, así como, parámetros que sirven paraconfigurar el comportamiento de la librería durante el tiempo de compila-ción. El código está en el archivo fuente fourier/shDefs.h . . . . . . . . . . . 76

B.2. Módulo en C que define las funciones para construir la base. Este códigoestá en el archivo fuente fourier/shBase.c . . . . . . . . . . . . . . . . . . . 80

B.3. Módulo de ordenamimento contine las funciones indic() y elp(). Este códigoestá en el archivo fuente fourier/shIndic.c . . . . . . . . . . . . . . . . . . . 88

B.4. Módulo de aritmética, contiene las funciones coef() y factorial(). Este códigoestá en el archivo fuente fourier/shAritmetica.c . . . . . . . . . . . . . . . . 89

B.5. Módulo para el cálculo de los coeficientes de Fourier. Este código está enel archivo fuente fourier/shExp.c . . . . . . . . . . . . . . . . . . . . . . . . 90

B.6. Makefile. Este código está en el archivo fuente fourier/Makefile . . . . . . . 91

Capítulo I

Introducción

Sea Sn la esfera unitaria en Rn+1 de estructura interna Riemanniana. El problema deDirichlet sobre Sn consiste en resolver

−∆Snu = f sobre Sn (I.1)

donde ∆Sn es el operador de Laplace-Beltrami con respecto a la estructura Riemannianade la esfera Sn. Diremos así mismo que, I.1 es la ecuación de Poisson en Sn.

De modo explícito, el operador ∆Sn, para una función u diferenciable en una vecindadde Sn, se calcula así. Determinamos el gradiente

∇u(x) =

(∂u

∂x1

, · · · ,∂u

∂xn+1

)en Rn+1, que luego proyectamos sobre el espacio tangente al punto (x1, · · · , xn+1) ∈ Sn

(cf. Warner [6]). En seguida habiendo prolongado el campo así obtenido a un abierto quecontiene la esfera de modo que siga siendo tangente a las esferas interiores y exteriores,tomamos su divergencia, la que da en los puntos de la esfera Sn el valor del operador deLaplace-Beltrami en la función u.

En este trabajo de grado nos proponemos buscar soluciones aproximadas a la ecua-ción I.1 mediante dos métodos. El primero es el método de elementos finitos de Galerkin yel segundo el método de Fourier fundado en la teoría de los valores propios del operadorde Laplace-Beltrami. En cada uno de los casos desarrollamos sus fundamentos teóricosy construimos software apropiado para la solución del problema.

11

Capítulo II

Solución de la ecuación diferencial dePoisson en la esfera S2 por el Métodode elementos finitos de Galerkin

En esté capítulo se estudia la solución aproximada del Problema de Dirichlet para laecuación de Poisson por el métodos de elementos finitos de Galerkin, el método estándares presentado en la sección II.2. El procedimiento para la triangulación de la esfera y laestructura de datos relevante son presentados en la sección II.3. Resultados importantesque simplifican la programación del método de Galerkin son presentados en la secciónII.4. Finalmente, son presentados ejemplos en las secciones II.7 y II.8.

II.1. Problema de Dirichlet y la ecuación de Poisson

Sea dado un dominio n-dimensional Ω que es parte de Rn o parte de una variedad dedimensión n, sobre el cual esta definida la función f . Entonces el problema de Dirichletpara la ecuación de Poisson se establece como,

−∆u = f sobre Ω, con u = 0 sobre ∂Ω (II.1)

donde u ∈ C2 (Ω)∩C(Ω)

es llamada solución clásica del problema de Dirichlet. La ecua-ción diferencial ∆u = f es llamada ecuación de Poisson y ∆ es el operador de Laplace-Beltrami definido asi:

∆(f) = div (∇(f))

(cf. Warner [6]).Si ∂Ω = ∅, como será nuestro caso en la esfera S2 sustituimos la condición u = 0 en

∂Ω, por u (u0) = 0, para un punto de Ω.

12

II.2. El método estándar de elementos finitos de Galerkin

Antes de describir la solución aproximada por el método de elemento finito de (II.1),consideraremos brevemente la aproximación de funciones suaves en Ω que se anulan enla ∂Ω. Concretamente, se aproximara por funciones continuas lineales sobre Ω.

Tomamos un dominio Ω con frontera suave ∂Ω, y denotemos τh una partición de Ω entriángulos disjuntos τ de tal forma que ningún vértice de un triángulo este en el interiorde otro triángulo y de modo que la reunión de los triángulos determine un dominio poli-gonal convexo Ωh. Aqui h denota la máxima longitud de los lados de la triangulación τh.Las triangulaciones apropiadas para el método de Galerkin requiere que sus triángulostengan angulos acotados por abajo por una constante positiva, independiente de h y trian-gulaciones cuasiuniformes en el sentido de que los triángulos de τh son esencialmentedel mismo tamaño, lo cual expresamos demandando que el área de τ ∈ τh este acotadapor abajo por ch2, con c > 0, independientemente de h.

Denotemos Sh las funciones continuas en la cerradura Ω de Ω que son lineales encada triángulo de τh y que se anulan fuera de Ωh. Representemos con PjNh

j=1 al conjuntode vértice interiores de τh. Una función en Sh se determina únicamente por sus valores enel punto Pj y por lo tanto depende de Nh parámetros. Llamemos Φj a la función piramidalen Sh que toma el valor 1 en Pj pero se anula en los otros vértices. Por lo tanto ΦjNh

j=1

forma una base para Sh, y cada χ en Sh admite la representación

χ(x) =

Nh∑j=1

αjΦj(x), con αj = χ(Pj). (II.2)

Entonces dada una función suave ν sobre Ω que se anula en ∂Ω, ella puede seraproximada, por lo que llamaremos su interpolado Ihν en Sh, que definimos como elelemento de Sh que coincide con ν en los vértices interiores, es decir,

Ihν(x) =

Nh∑j=1

ν(Pj)Φj(x). (II.3)

En un contexto más general, para Ω ⊂ Rd la norma en L2 = L2(Ω) la denotaremos por‖·‖ y la norma en el espacio de Sobolev Hr = Hr(Ω) = W r

2 (Ω), la denotaremos por ‖·‖r,así para una función real ν,

‖ν‖ = ‖ν‖L2=

(∫Ω

ν2dx

) 12

(II.4)

13

y, para un entero positivo r,

‖ν‖r = ‖ν‖Hr =

∑|α|≤r

‖Dαν‖2

12

, (II.5)

donde, con α = (α1, . . . , αd), Dα =(

∂∂x1

)α1

· · ·(

∂∂xd

)αd

denota una derivada arbitraria con

respecto a x = (x1, . . . , xn) de orden |α| =∑d

j=1 αj, de tal forma que la suma en II.5contiene todas las derivadas con orden a lo sumo de r. Recordemos que para funcionesen H1

0 = H10 (Ω), es decir las funciones ν con ∇ν = grad ν en L2 y que se anulan sobre

∂Ω, ‖∇ν‖ y ‖ν‖1 son normas equivalentes (lema de Friedrichs, vea, p.ej. [8] o [7])

c ‖ν‖1 ≤ ‖∇ν‖ ≤ ‖ν‖1 ∀ν ∈ H10 , con c > 0. (II.6)

donde c y C denotan constantes positivas, las cuales son independientes de los paráme-tros y las funciones involucradas.

Asumamos que tenemos una familia Sh de subespacios H10 de dimensión finita de

tal forma, que para un entero dado r ≥ 2 y algún h pequeño,

ınfχ∈Sh

‖ν − χ‖+ h ‖ν − χ‖ ≤ Chs ‖ν‖s , for 1 ≤ s ≤ r, (II.7)

cuando ν ∈ Hs ∩ H10 . El número r se conocerá como orden de exactitud de la familia

Sh. Frecuentemente Sh consiste en polinomios seccionalmente continuos con ordende a lo sumo de r − 1 sobre una triangulación τh. Generalmente estimaciones comola (II.7) pueden obtenerse frecuentemente al considerar un operador de interpolaciónIh : Hr ∩H1

0 → Sh tal como

‖Ihν − ν‖+ h ‖∇ (Ihν − ν)‖ ≤ Chs ‖ν‖s para 1 ≤ s ≤ r. (II.8)

El orden óptimo con el que pueden ser aproximadas las funciones y sus gradientesbajo nuestra suposición (II.7), es O (hr) y O (hr−1) , respectivamente, y a continuaciónse construirán aproximaciones de ese orden para la solución del problema de Dirichlet(II.1). Para este propósito escribiremos este problema en su forma variacional o débil:Multiplicaremos la ecuación elíptica por una función suave ϕ que se anula en la frontera∂Ω (es suficiente requerir ϕ ∈ H1

0 ), luego integramos sobre Ω, y aplicamos la fórmula deGreen en el lado izquierdo, para obtener

(−∆u, ϕ) = (∇u,∇ϕ) = (f, ϕ) , ∀ϕ ∈ H10 , (II.9)

14

en donde hemos usado el producto interno en L2,

(v, w) =

∫Ω

vwdx, (∇v,∇w) =

∫Ω

d∑j=1

∂v

∂xj

∂w

∂xj

dx. (II.10)

II.2.1. Método de Galerkin

En el método de elemento finito nos planteamos resolver el problema buscando unaaproximación uh ∈ Sh de tal manera que

(−∆uh, χ) = (∇uh,∇χ) = (f, χ), ∀χ ∈ Sh, (II.11)

definir así la solución aproximada en términos de la formulación variacional del proble-ma es conocida como método de Galerkin, en honor al matemático aplicado Ruso BorisGrigorievich Galerkin (1871-1945).

Note que de las ecuaciones (II.9) y (II.11) resulta

(∇ (uh − u) ,∇χ) = 0, ∀χ ∈ Sh (II.12)

es decir, que el error en la solución discreta es ortogonal a Sh con respecto al productointerno de Dirichlet (∇v,∇w).

En términos de una base ΦjNh

1 para el espacio de elementos finitos Sh, planteamosnuestro problema discreto: Encontrar los coeficientes αj en uh(x) =

∑Nh

j=1 αjΦj(x) de talforma que

Nn∑j=1

αj (∇Φj,∇Φk) = (f, Φk) , para k = 1, . . . , Nh. (II.13)

lo que podemos expresar en notación matricial como

Bα = f ,

donde B = (bjk) es la matriz de rigidez con elementos bjk = (∇Φj,∇Φk), f = (fk) esun vector donde fk = (f, Φk), y α el vector de las incógnitas αj. La dimensión de estosarreglos es Nh, que es también la dimensión de Sh (el cual es igual al número de vérticesinteriores)

El siguiente teorema (probado en Thomée [1]) proporciona un medio para estimar elerror entre la solución continua y discreta del problemas.

Teorema II.1 Asuma que se cumple II.7, y sean uh y u las soluciones de II.11 y II.1,

15

respectivamente. Entonces para 1 ≤ s ≤ r,

‖uh − u‖ ≤ C hs ‖u‖s y ‖∇uh −∇u‖ ≤ C hs−1 ‖u‖s .

II.3. Triangulación de la esfera

Consideremos la esfera S2 nuestro dominio, τh denotará una partición de la esfera entriángulos disjuntos τ de tal forma que ningún vértice de cada triángulo este en el ladointerior de otro triángulo y tal que la unión de los triángulos determine un dominio poligo-nal S2h

⊂ S2. Estas triangulaciones son generadas recursivamente, empezando con unicosaedro inscrito en la esfera unitaria, luego subdividiendo recursivamente cada triángu-lo en cuatro al dividir cada lado en dos por el medio y luego proyectando este punto ala esfera unitaria en la dirección normal. Este proceso genera una secuencia de triangu-laciones cuasiuniformes que tienen 20, 80, 320, 1280, 5120, 20480, 81820,..., triángulosrespectivamente. En efecto, los triángulos tienen ángulos acotados inferiormente por 3

10π.

Denotaremos estas triangulaciones como Gen(0), Gen(1),...,. La figura II.1 muestra lasprimeras seis mallas esféricas de triángulos Gen(k) para k = 0, 1, 2, 3, 4, 5.

Figura II.1: Mallas esféricas de triángulos. Gen(k) para k = 0, 1, 2, 3, 4, 5

16

II.3.1. Estructura de los datos para un dominio discretizado

Para propósitos de programación establecemos la siguiente estructura de dato. Seaun dominio Ω y una partición τh de Ω que genera un conjunto de nodosPjNh

j=1. El do-minio plano mostrado en la figura II.2 servirá de ejemplo. Podemos reconocer en τh lossiguientes elementos de información:

nodos que identificamos con enteros positivo y guardan información sobre las coor-denadas. (Ej. nodo 8 coordenadas (x, y, z) = (1

4; 1

2;√

114

)). También tienen datobandera que indica si el nodo esta en el borde, esta condición no se usaráen esta tesis, pero es necesaria para un futuro, cuando consideremos unacondición mixta en el borde;

triángulos que identificamos con enteros positivos y guardan información de los nodos.(Ej. triángulo 6 nodos 3, 4, 2.);

lados que identificamos con enteros positivo y guardan información de los nodosque forman un lado de un triángulo dado;

En la lista A.5 se muestra la implementación de la estructura de datos tSurface.

Ωh

Ω

τ

NODO

Lados

Figura II.2: Discretización del dominio Ω

El procedimiento recursivo necesita que las dimensiones de la estructura de datosaumente en cada iteración. Para ilustrar el proceso usemos el dominio mostrado en lafigura II.2 subdividamos cada triángulo en cuatro al dividir cada lado en dos por el mediocomo en la figura II.3. Para este procedimiento tenemos las siguientes relaciones:

17

Ωh

Ω

τ

NODO

Figura II.3: Discretización del dominio Ω: primera iteración.

li = 2× li−1

ni = 2× li−1 + 3× ti−1

ti = 4× ti−1

donde para una iteración i dada:

li es el número de lados,

ni es el número de nodos,

ti es el número de triángulos.

En la lista A.6 se muestra la implementación de una función gen(i) que genera una esferaconstruida con una malla de triángulos con el método recursivo antes mencionado. Encaso de un dominio no esférico aún existe el problema de redistribuir de forma adecuadalos nuevos nodos generados PjNh

j=1.

18

II.4. Programación de (∇Φi · ∇Φj)

II.4.1. Cálculo del gradiente ∇Φj de la función nodal para un triángu-lo τj.

Para realizar este cálculo se uso la fórmula expuesta en un trabajo de Zeeman [2] parael ∇ de una función u : X → R,

(∇u)i =n∑

j=1

(g−1)

ij

∂u

∂xj

(II.14)

donde X es una variedad Riemanniana n-dimensional cuya estructura está determinadapor una forma simétrica positiva definida g sobre los planos tangentes de X. En coorde-nadas locales g esta dada por una matriz m × n positiva definida (gij) que depende dex ∈ X, y que determina la siguiente métrica Riemanniana

dx2 =n∑

i,j=1

gijdxidxj.

Queremos tener una función Φ(x, y, z) definida en una vecindad del triángulo4 (P1, P2, P3), cuyos vértices se encuentran en el espacio tri-dimensional y determinarla expresión de su gradiente en coordenadas locales, se tiene,

(∇Φ)i = g−1 (∇ (Φ θτ )) donde i = 1, 2, 3

donde,

g =

(g11 g12

g12 g22

)donde g11 =

∂θτ

∂u

∂θτ

∂u; g12 =

∂θτ

∂u

∂θτ

∂v; g22 =

∂θτ

∂u

∂θτ

∂v

para una parametrización θτ del triángulo P1, P2 , P3.Si tomamos la parametrización lineal que lleva el triángulo A1, A2, A3 del plano R2 en

el triángulo P1, P2 , P3, se tiene θτ (u, v) = α(P2−P1)+β(P3−P1) siendo α y β las funcionescaracterizadas por (u, v) = α(A2 − A1) + β(A3 − A1) (vea la figura II.4).

de donde resultan

α =

∣∣∣∣∣ (u, v)(A2 − A1) (A3 − A1)(A2 − A1)

(u, v)(A3 − A1) (A3 − A1)2

∣∣∣∣∣area(A1, A2, A3)

19

A2

P3

P1

z

y

x

vA3

P3

A1v

Figura II.4: Transformación θ : Ω ⊂ R2 → R3

β =

∣∣∣∣∣ (A2 − A1)2 (u, v)(A2 − A1)

(A2 − A1)(A3 − A1) (u, v)(A3 − A1)

∣∣∣∣∣area(A1, A2, A3)

area(A1, A2, A3) =

∣∣∣∣∣ (A2 − A1)2 (A2 − A1) (A3 − A1)

(A2 − A1) (A3 − A1) (A3 − A1)2

∣∣∣∣∣∂θτ

∂u=

∂α

∂u(P2 − P1) +

∂β

∂u(P3 − P1) y

∂θτ

∂u=

∂α

∂v(P2 − P1) +

∂β

∂v(P3 − P1)

g12 =

(∂α

∂u

∂α

∂v

)(−−→P1P2

)2

+∂β

∂u

∂β

∂v

(−−→P1P3

)2

+

(∂α

∂u

∂β

∂v+

∂α

∂v

∂β

∂u

)−−→P1P2

−−→P1P3 (II.15)

g11 =

(∂α

∂u

)2 (−−→P1P2

)2

+

(∂β

∂u

)2 (−−→P1P3

)2

+ 2

(∂α

∂u

∂β

∂u

)−−→P1P2

−−→P1P3 (II.16)

g22 =

(∂α

∂u

)2 (−−→P1P2

)2

+

(∂β

∂u

)2 (−−→P1P3

)2

+ 2

(∂α

∂v

∂β

∂v

)−−→P1P2

−−→P1P3 (II.17)

∂α

∂u=

1

4 (area (A1, A2, A3))2

[(u2 − u1)

(−−−→A1A3

)2

− (u3 − u1)(−−−→A1A3 −

−−−→A1A2

)]∂α

∂v=

1

4 (area (A1, A2, A3))2

[(v2 − v1)

(−−−→A1A3

)2

− (v3 − v1)(−−−→A1A3 −

−−−→A1A2

)]

20

∂β

∂u=

1

4 (area (A1, A2, A3))2

[(u3 − u1)

(−−−→A1A2

)2

− (u2 − u1)(−−−→A1A3 −

−−−→A1A2

)]∂β

∂v=

1

4 (area (A1, A2, A3))2

[(v3 − v1)

(−−−→A1A2

)2

− (u2 − u1)(−−−→A1A3 −

−−−→A1A2

)]Para una función F definida en el triángulo (o la superficie) se tiene

∇F = (∇F )u

∂θτ

∂u+ (∇F )v

∂θτ

∂v

y las componentes (∇F )u, (∇F )v se obtienen((∇F )u

(∇F )v

)= g−1

(∂(Fθτ )

∂u∂(Fθτ )

∂v

)

siendo g−1 la matriz inversa de la matriz

g =

(g11 g12

g12 g22

)

de donde sigue que

g−1 =1

det(g)

(g22 −g12

−g12 g11

)

II.4.2. Cálculo del producto escalar de los gradientes ∇F y ∇G

Para dos funciones F y G definidas sobre la imagen de θτ se tendrá[(∇F )u

(∇F )v

]= [∇F ] = g−1

(∂(Fθτ )

∂u∂(Fθτ )

∂v

)= g−1

(∂F∂u∂F∂v

)[

(∇G)u

(∇G)v

]= [∇G] = g−1

(∂(Gθτ )

∂u∂(Gθτ )

∂v

)= g−1

(∂G∂u∂G∂v

)

∇F = (∇F )u

∂θτ

∂u+ (∇F )v

∂θτ

∂v

∇G = (∇G)u

∂θτ

∂u+ (∇G)v

∂θτ

∂v

por lo tanto, el producto escalar de estos gradientes sobre la superficie da como resultado

21

∇F · ∇G = [∇F ]T g [∇G]

=(

∂F∂u

, ∂F∂v

)g−1g g−1

(∂G∂u∂G∂v

)

=(

∂F∂u

, ∂F∂v

)g−1

(∂G∂u∂G∂v

)

= 1det(g)

(∂F∂u

, ∂F∂v

)[ g22 −g12

−g12 g11

](∂G∂u∂G∂v

)= 1

det(g)

[∂F∂u

∂G∂u

g22 −(

∂F∂u

∂G∂v

+ ∂F∂v

∂G∂u

)g12 + ∂F

∂v∂G∂v

g11

]II.4.3. Cálculo de (∇Φi · ∇Φj)

El cálculo del producto interno de los gradientes de funciones nodales

(∇Φi · ∇Φj) =

∫4(P1,P2,P3)

(∇Φi · ∇Φj) (II.18)

se reduce a un cálculo algebraico que considera solo la geometría del triángulo involu-crado. Sean los índices i, j, k una permutación del conjunto 1, 2, 3 que se correspondecon los vértices del triángulo 4 (P1, P2, P3). Tomemos por ejemplo el caso mostrado en lafigura II.5, para el cual se cumple lo siguiente

u

v

A2

Aj=3 Ai=1

Figura II.5: Caso particular de un elemento triangular en el plano.

−−−→A1A3 = (−h, 0) ,

−−−→A1A2 = (0, h)

u2 − u1 = 0, v2 − v1 = 0, u3 − u1 = −h, v3 − v1 = 0.

22

De alli se obtiene

∂α

∂u=

1

h4(0 + h · 0) = 0,

∂α

∂v=

1

h4

(h · h2

)=

1

h,

∂β

∂u= −h · h2

h4= −1

h,

∂β

∂v= 0

∇ΦP1 = ∇ΦA1 =

(1

h,−1

h

), ∇ΦP2 = ∇ΦA2 =

(0,

1

h

), ∇ΦP3 = ∇ΦA3 =

(−1

h, 0

)de las ecuaciones II.15, II.16, II.17, seguidamente podemos escribir

g11 =1

h2

(−−→P1P3

)2

, g12 = − 1

h2

−−→P1P2 ·

−−→P1P3, g22 =

1

h2

(−−→P1P2

)2

det(g) =1

h4

∣∣∣∣∣∣(−−→P1P3

)2 −−→P1P2 ·

−−→P1P3

−−→P1P2 ·

−−→P1P3

(−−→P1P2

)2

∣∣∣∣∣∣=

1

h4

((−−→P1P3

)2

·(−−→P1P2

)2

−(−−→P1P2 ·

−−→P1P3

)2)

Ahora podemos resolver

(∇Φi · ∇Φj) =

∫4(A1,A2,A3)

J

det(g)

[∂Φi

∂u

∂Φj

∂ug22 −

(∂Φi

∂u

∂Φj

∂v+

∂Φi

∂v

∂Φj

∂u

)g12 +

∂Φi

∂v

∂Φj

∂vg11

]

En primer lugar consideremos i 6= j, en particular i = 1 y j = 2

(∇Φi · ∇Φj) =

(−−→P1P2 ·

−−→P1P3 −

(−−→P1P3

)2)

1h4

1h4

((−−→P1P3

)2

·(−−→P1P2

)2

−(−−→P1P2 ·

−−→P1P3

)2) ∫

4(A1,A2,A3)

(∇Φi · ∇Φj) =

(−−→P1P2 ·

−−→P1P3 −

(−−→P1P3

)2)

12√(−−→

P1P3

)2

·(−−→P1P2

)2

−(−−→P1P2 ·

−−→P1P3

)2

sacando el factor común−−→P1P3, y tomando en cuenta que

−−→P1P2 −

−−→P1P3 =

−−→P1P3 ·

−−→P3P2,

conduce a la siguiente forma general

23

(∇Φi,∇Φj)4(P1,P2,P3) = −

(−−→PkPj ·

−−→PiPk

)4 · area (P1, P2, P3)

note que las consideraciones particulares hechas en el plano desaparecen. Finalmentepara i = j obtenemos

(∇Φi · ∇Φj) =

((−−→P1P2

)2

− 2(−−→P1P2 ·

−−→P1P3

)+(−−→P1P3

)2)

1h4

1h4

((−−→P1P3

)2

·(−−→P1P2

)2

−(−−→P1P2 ·

−−→P1P3

)2) ∫

4(A1,A2,A3)

(∇Φi · ∇Φj) =

((−−→P1P2

)2

− 2(−−→P1P2 ·

−−→P1P3

)+(−−→P1P3

)2)

12√(−−→

P1P3

)2

·(−−→P1P2

)2

−(−−→P1P2 ·

−−→P1P3

)2

la expresión((−−→

P1P2

)2

− 2(−−→P1P2 ·

−−→P1P3

)+(−−→P1P3

)2)

se transforma en(−−→P1P2 −

−−→P1P3

)2

que es(−−→P3P2

). Lo que lleva a escribir la siguiente expresión general

(∇Φi,∇Φi)4(P1,P2,P3) =

(−−→PkPj

)2

4 · area (P1, P2, P3)

en donde también desaparecen las consderaciones generales que se hicieron en el planopara realizar este cálculo. Ahora es posible independizarse de una parametrización de lasuperficie, es decir, se puede calcular de forma precisa los elementos de la matriz queresulta en el problema de Galerkin. Así reducimos el esfuerzo de computo numérico alevitar una integración numérica.

II.5. Estimación del error

El error en la norma L2 (S2) entre la solución exacta y la calculada se estima del modosiguiente

‖uh − u‖L2≤√

2π

3× n

n∑i=1

|sci − sri| . (II.19)

donde n es el número de triángulos de la triangulación considerada, sci y sri son lassolución calculada y la real para el punto i respectivamente.

24

Orden Figura Error Puntos CoordenadasTriangulares Pesos

Cúbica R = O (h3)abc

12, 1

2, 0

0, 12, 1

212, 0, 1

2

131313

Cuadro II.1: Integración numérica cúbica.

En efecto, considerando el teorema II.1 y (uh − u)Pi=∑n

j=1 (sci − sri) tenemos

‖uh − u‖L2≤

n∑i=1

|sci − sri| ‖Φi‖L2

y de la cuasiuniformidad de la triangulación

‖Φi‖L2≤√

2π

3n

II.6. Cálculo de (f, Φk)

Para este cálculo usaremos la fórmula cúbica de integración numérica para triángulospresentada en Zienkiewics y Taylor [5] y que se reproduce en la tabla II.1 que sigue,obviamente por su simplicidad. Ciertamente podemos esperar mejores resultados parainterpolaciones de orden superior.

II.7. Ejemplo (comparando una solución real y una cal-

culada)

Sea la Ecuación Diferencial

∆S2u = 8z − 48x2z + 20y − 120z2y (II.20)

donde el lado derecho se obtuvo aplicando la ecuación III.3 al polinomio u = 4x2z +10z2y.En la lista II.1 implementamos un programa en lenguaje C que calcula la solución apro-ximada por el método de Galerkin para una triangulación dada generada con la funcióngen(#).

25

Lista II.1: Programa en C para resolver la ecuación diferencial por el método de elementosfinitos de Galerkin. Este código está en el archivo fuente galerkin/prueba01.c

#include " tSur face . h "#include " so lveGa le rk in . h "#include <math . h>#include < s t d i o . h>

doublef ( vec to r p )

double x = p . x , y = p . y , z = p . z ;return (−(8 − 24 ∗ pow ( x , 2) + 20 ∗ y − 120 ∗ pow ( z , 2) ∗ y ) )

;

doubleu ( node n )

double x = n . x , y = n . y , z = n . z ;return (4 ∗ pow ( x , 2) + 10 ∗ pow ( z , 2) ∗ y ) ;

i n tmain ( )

i n t i ;tSur face esfera ;for ( i = 0 ; i <= 3 ; i ++)

es fera = gen ( i ) ;tSu r faceP lo t ( esfera , " good " ) ;es fera = solveGalerkinOnThe_tSurface ( esfera , 12) ;i f ( i == 4)

p lo tSo lu t i onOver ( esfera , " solRealPrueba01 . dat " , ’ r ’ ) ;p lo tSo lu t i onOver ( esfera , " solCalcPrueba01 . dat " , ’ c ’ ) ;p lo tSo lu t i onOver ( esfera , " solErrorPrueba01 . dat " , ’ e ’ ) ;

26

a) b)

c)

‘

Figura II.6: a) Solución real, b) Solución calculada para una malla de 1280 triángulosgenerada con gen(3), c) Desviación entre las soluciones calculada y real.

tSur faceDest roy ( es fera ) ;p r i n t f ( " \ n \ n \ n " ) ;

II.7.1. Resultados

La tabla II.3 muestra la solución real y calculada en los doce nodos iniciales del icosae-dro inscrito en la esfera, es decir, los nodos de la esfera generada con la función gen(0).Notamos que medida que el número de triángulos aumenta y la malla se hace más fina,entonces, la solución calculada tiende a igualarse a la solución real. La escala de coloresen la figura II.6 indica que las soluciones real II.6.a) y la calculada II.6.b) son similares.Además la figura II.6.c) muestra que el error es pequeño y varía en función de los cambiosde u.

27

NODO 1 2 3 4 5 6SOL.REAL 2.894 2.894 2.894 2.894 1.106 1.106

20 t. 3.378 2.942 2.942 3.378 0.8897 0.8897err. 0.4838 0.04731 0.04731 0.4838 0.2159 0.215980 t. 1.396 1.396 1.396 1.396 -0.7457 -0.7457err. 1.498 1.498 1.498 1.498 1.851 1.851

320 t. 2.463 2.461 2.461 2.463 0.5985 0.5985err. 0.4315 0.4334 0.4334 0.4315 0.5071 0.5071

1280 t. 2.778 2.778 2.778 2.778 0.9766 0.9766err. 0.1164 0.1166 0.1166 0.1164 0.1289 0.1289

5120 t. 2.863 2.863 2.863 2.863 1.073 1.073err. 0.03109 0.0311 0.0311 0.03109 0.03253 0.03253

20480 t. 2.886 2.886 2.886 2.886 1.097 1.097err. 0.008385 0.008386 0.008386 0.008385 0.008319 0.008319

81920 t. 2.892 2.892 2.892 2.892 1.103 1.103err. 0.002271 0.002271 0.002271 0.002271 0.002148 0.002148

NODO 7 8 9 10 11 12SOL.REAL 1.106 1.106 3.804 -3.804 -3.804 1.04e-312

20 t. 3.336 3.336 1.767 0.5819 1.164 1.889e-312err. 2.23 2.23 2.037 4.386 4.968 1.616e-31080 t. 0.5561 0.5561 1.573 -3.692 -3.692 1.573err. 0.5495 0.5495 2.231 0.1123 0.1123 2.231

320 t. 0.9152 0.9152 3.153 -3.753 -3.745 3.155err. 0.1904 0.1904 0.6513 0.05125 0.05876 0.6492

1280 t. 1.049 1.049 3.627 -3.783 -3.783 3.627err. 0.05637 0.05637 0.1771 0.02086 0.0217 0.1769

5120 t. 1.089 1.089 3.757 -3.797 -3.797 3.757err. 0.01637 0.01637 0.04747 0.006755 0.006828 0.04745

20480 t. 1.101 1.101 3.791 -3.802 -3.802 3.791err. 0.004795 0.004795 0.01276 0.001904 0.00191 0.01276

81920 t. 1.104 1.104 3.801 -3.804 -3.804 3.801err. 0.001397 0.001397 0.003433 0.0005041 0.0005045 0.003433

Cuadro II.3: Resultados usando la función “solveGalerkinOnThe_tSurface” sobre mallasde triángulos esférica, t≡triángulos y error por diferencia con la solución real.

28

a) b)

Figura II.7: a) Distribución sobre la esfera, b) Distribución con integral nula.

II.8. Ejemplo

Cosideremos una ecuación diferencial con la siguiente forma

∆S2u = exp

(20×

(z − 1

z2

))+ k (II.21)

donde la forma exponencial en el lado derecho representa una incursión violenta distri-buida sobre la esfera como se muestra en la figura II.7.a. k es una constante que anula laintegral sobre la esfera de todo el lado derecho de la expresión mostrada, esto equivale a,una fina película de espesor k dispersa en el interior esfera capaz de absorber la interven-ción representada por la distribución mencionada. La figura II.7.b muestra la distribucióncon integral nula.

II.8.1. Resultados

En la tabla II.4 se muestra la solución de para los 42 nodos que se obtienen despuésde afinar una vez el icosaedro inicial para una malla esférica. La figura II.8 muestra enescala de colores la solución sobre la esfera. El código del programa usado para resolvereste problema es mostrado en la lista II.2.

Lista II.2: Programa en C para resolver la ecuación diferencial por el método de elementosfinitos de Galerkin. Este código está en el archivo fuente galerkin/prueba02.c

29

Solución calculada según el número de triángulos (t.) de la malla.N 20 t. 80 t. 320 t. 1280 t. 5120 t. 20480 t. 81920 t.1 5.49506e-08 -2.37411e-06 0.0002841 0.0007424 0.0009815 0.001091 0.0011422 1.2039e-07 4.17927e-05 0.003026 0.006467 0.007785 0.008196 0.0083243 1.2039e-07 4.17927e-05 0.003026 0.006467 0.007785 0.008196 0.0083244 5.49506e-08 -2.37411e-06 0.0002841 0.0007424 0.0009815 0.001091 0.0011425 8.13785e-08 1.71766e-05 0.001312 0.00283 0.003447 0.003661 0.0037396 8.13785e-08 1.71766e-05 0.001312 0.00283 0.003447 0.003661 0.0037397 9.39625e-08 9.82469e-06 0.001231 0.002798 0.003439 0.003659 0.0037398 9.39625e-08 9.82469e-06 0.001231 0.002798 0.003439 0.003659 0.0037399 2.17292e-07 0.000105235 0.005802 0.01214 0.01452 0.01522 0.0154310 2.21487e-07 0.000101926 0.005774 0.01213 0.01451 0.01522 0.0154311 4.61413e-08 -1.18194e-05 -0.0002492 -0.0002748 -0.0001732 -9.347e-05 -4.805e-0512 -8.50992e-06 -0.0001514 -0.000196 -0.000141 -8.269e-05 -4.471e-0513 1.74627e-05 0.001307 0.002825 0.003445 0.00366 0.00373914 4.6938e-05 0.002944 0.006225 0.007482 0.007876 0.00815 4.6938e-05 0.002944 0.006225 0.007482 0.007876 0.00816 7.7122e-05 0.005157 0.01091 0.01308 0.01372 0.0139117 2.95867e-05 0.002158 0.004628 0.005595 0.005908 0.00601118 1.37852e-05 0.00128 0.002815 0.003442 0.00366 0.00373919 4.73443e-06 0.0006595 0.0015 0.001875 0.002023 0.00208420 -8.44278e-06 -0.0001128 -9.367e-05 -1.162e-05 5.478e-05 9.491e-0521 3.65749e-07 0.0003501 0.0008377 0.001084 0.001195 0.00124722 3.65749e-07 0.0003501 0.0008377 0.001084 0.001195 0.00124723 -1.02728e-05 -0.0001565 -0.0001201 -1.974e-05 5.266e-05 9.438e-0524 -1.31913e-05 -0.0003722 -0.0006027 -0.0006022 -0.0005573 -0.000522525 -8.44278e-06 -0.0001128 -9.367e-05 -1.162e-05 5.478e-05 9.491e-0526 -1.02728e-05 -0.0001565 -0.0001201 -1.974e-05 5.266e-05 9.438e-0527 4.73443e-06 0.0006595 0.0015 0.001875 0.002023 0.00208428 2.95867e-05 0.002158 0.004628 0.005595 0.005908 0.00601129 1.37852e-05 0.00128 0.002815 0.003442 0.00366 0.00373930 7.7122e-05 0.005157 0.01091 0.01308 0.01372 0.0139131 0.000158326 0.0122 0.02707 0.03238 0.03389 0.0343132 7.5292e-05 0.005141 0.0109 0.01307 0.01372 0.0139133 7.5292e-05 0.005141 0.0109 0.01307 0.01372 0.0139134 2.63819e-05 0.002124 0.004614 0.005591 0.005907 0.00601135 4.04942e-05 0.00289 0.006204 0.007477 0.007875 0.00836 1.52968e-06 0.0006104 0.001477 0.001869 0.002021 0.00208337 -6.07804e-06 0.0002019 0.000773 0.001067 0.001191 0.00124638 -6.07804e-06 0.0002019 0.000773 0.001067 0.001191 0.00124639 1.52968e-06 0.0006104 0.001477 0.001869 0.002021 0.00208340 2.63819e-05 0.002124 0.004614 0.005591 0.005907 0.00601141 4.04942e-05 0.00289 0.006204 0.007477 0.007875 0.00842 0.001197 0.002784 0.003435 0.003658 0.003738

Cuadro II.4: Solución de la ecuación II.21 por el método de elementos finitos de Galerkin.

30

Figura II.8: Solución calculada por el método de Galerkin para la ecuación diferencialII.21.

#include " tSur face . h "#include " so lveGa le rk in . h "#include <math . h>#include < s t d i o . h>


double x = p . x , y = p . y , z = p . z ;i f ( z >= 0)

return ( exp (20 ∗ (pow ( z , 2) − 1 / pow ( z , 2) ) ) ) ;else

return 0;

doubleu ( node n )

double x = n . x , y = n . y , z = n . z ;return ( 0 ) ;

i n t

31

main ( )

i n t i ;tSur face esfera ;for ( i = 4 ; i <= 4 ; i ++)

es fera = gen ( i ) ;/ / tSu r faceP lo t ( esfera , " good " ) ;esfera = solveGalerkinOnThe_tSurface ( esfera , 42) ;i f ( i == 4)

p lo tSo lu t i onOver ( esfera , " solCalcPrueba02 . dat " , ’ c ’ ) ;tSur faceDest roy ( es fera ) ;p r i n t f ( " \ n \ n \ n " ) ;

32

Capítulo III

Teoría de armónicos esféricos

III.1. Introducción

Sea Pk el espacio de los polinomios homogéneos de grado k sobre Rn+1, Sn la esferaunitaria en Rn+1 de estructura interna Riemanniana, y sea

Hk = P ∈ Pk : ∆P = 0 , (III.1)

Hk = P |Sn : P ∈ Hk .

Es decir, Hk es el espacio de los polinomios armónicos homogéneos de grado k y Hk

es el espacio de su restricción a la esfera unitaria. Los elementos de Hk son llamadosarmónicos esféricos de grado k. El mapeado de Hk a Hk es un isomorfismo; cuyo inversoes el mapa Y → P donde P (x) = |x|kY (|x|−1x).

Usemos r2 para denotar a la función x →∑n

1 x2j sobre Rn, para un elemento de P2.

Proposición III.1 Pk = Hk ⊕ r2Pk−2, donde r2Pk−2 = r2P : P ∈ Pk−2. Es decir, Hk esel complemento ortogonal de r2Pk−2.

Prueba: cf. Folland [4].

de III.1 se sigue el siguiente corolario que se puede probar por inducción de k, (cf. Folland[4]).

Corolario III.2 Pk = Hk ⊕ r2Pk−2 ⊕ r4Pk−4 ⊕ · · · , (cf. Folland [4]).

Corolario III.3 La restricción a la esfera unitaria de cualquier elemento de Pk es unasuma de armónicos esféricos con grado k a lo sumo, (cf. Folland [4]).

33

Teorema III.4 L2 (Sn) =⊕∞

0 Hk, donde la expresión en el lado derecho es una sumaortogonal directa con respecto al producto escalar en el espacio de Hilbert L2(Sn). Esdecir se tiene,

(f, g) = 0 , f ∈ Hi , g ∈ Hj para i 6= j

(cf. Folland [4]).

Proposición III.5

dimPk =(k + n− 1)!

k! (n− 1)!.

(cf. Folland [4]).

Teorema III.6 Sea f una función polinomial en Sn, se tiene

∆Sn (f) =n+1∑i=1

(1− x2

i

) ∂2f

∂x2i

− 2∑i<j

∂2f

∂xi∂xj

xixj − nn+1∑i=1

∂f

∂xi

xi. (III.2)

(cf. Guíñez y Rueda [3]).

En el caso n = 2 tenemos

∆S2 (f) =(1− x2

) ∂2f

∂x2+(1− y2

) ∂2f

∂y2+(1− z2

) ∂2f

∂z2

−2

[∂2f

∂x∂yxy +

∂2f

∂x∂zxz +

∂2f

∂y∂zyz

](III.3)

−2

[∂f

∂xx +

∂f

∂yy +

∂f

∂zz

]Aplicando (III.2) a un monomio f tenemos:

Lema III.7 Sea f un monomio con grado k, entonces

∆Sn (f) = (−k (k − 1)− nk) · f +n+1∑i=1

∂2f

∂x2i

(III.4)

con H(k, n) = (−k(k − 1)− nk) (cf. Guíñez y Rueda [3]).

Prueba: apliquemos el la ecuación III.2 a un momio f con grado k,

∆Sn (f) =n+1∑i=1

(1− x2

i

) ∂2f

∂x2i

− 2∑i<j

∂2f

∂xi∂xj

xixj − nn+1∑i=1

∂f

∂xi

xi

34

para f = xα11 xα2

2 · · ·xαn+1

n+1 , resulta

∆Sn (f) =n+1∑i=1

f

x2i

(αi) (αi − 1)−

n+1∑i=1

(αi) (αi − 1) f − 2∑i<j

(αi) (αj) f − nn+1∑i=1

αif

=n+1∑i=1

f

x2i

(αi) (αi − 1)−

f

[n+1∑i=1

(αi) (αi − 1)− 2∑i<j

(αi) (αj)− n

n+1∑i=1

αi

]

=n+1∑i=1

f

x2i

(αi) (αi − 1)− f

×

[n+1∑i=1

(αi) (αi) + 2∑i<j

(αi) (αj)−n+1∑i=1

(αi) + nn+1∑i=1

αi

]

=n+1∑i=1

∂2f

∂x2i

− f[(q)2 − q + nq

]= (−k (k − 1)− nk) f +

n+1∑i=1

∂2f

∂x2i

si el monomio f ∈ Hk, entonces (III.4) se reduce a

∆Sn (f) = (−k (k − 1)− nk) · f

Puesto que la esfera Sn no tiene borde, el problema de Dirichlet para la ecuación dePoisson sobre Sn se puede escribir,

−∆Snu = f sobre Sn

donde ∆Sn es el operador de Laplace-Beltrami definido por Guíñez y Rueda [3] con res-pecto a la estructura Riemanniana de la esfera Sn, f es una función definida en la esferaSn. La función u ∈ C2 (Ω) ∩ C

(Ω)

es la solución clásica del problema de Dirichlet.Por el teorema (III.4) para f ∈ L2 (Sn) escribimos

f =∞∑0

fk

con fk ∈ Hk.

35

Por su parte, si Y1, ..., Ydkes una base de Hk donde dk = dimPk − dimPk−2es la

dimensión de esta base,

f =

dk∑j=1

Ak,jYj

en donde Ak,j verificaM (A) = F

donde M es una matriz de dk × dk

Mij = (Yi, Yj)

y F es la columnaFj = (f, Yi) , j = 1, . . . , dk

Por lo tanto, con el fin de tener un algoritmo explícito, que nos permita un tratamientode la expresión de fk, es necesario obtener una base de Hk. Para ello determinamosuna base Pk y una expresión para la matricial del operador ∆ =

∑n+1i=1

∂2

∂x2i

y seguido unalgoritmo preciso para determinar una base adecuada en el espacio nulo de tal matriz. Elcálculo de la matriz M , obliga entonces al cálculo de integrales

∫Sn

f ds, para polinomios,lo que nos lleva a establecer una fórmula para

∫Sn

x2α11 x2α2

2 · · ·x2αn+1

3 ds

Teorema III.8 Sea f una función polinomial de la forma x2α11 x2α2

2 · · ·x2αn+1

3 con 2q = 2α1 +

2α2 + · · ·+ 2αn+1en Sn, se tiene∫S2

f =(2α1)! (2α2)! · · · (2αn+1)!

H(2q, n)×H(2q − 1, n) · · ·H(2, n)

(q)!

(α1)! (α2)! · · · (αn+1)!× vol (Sn) (III.5)

donde H(k, n) = −k(k − 1)− n k. (cf. Guíñez y Rueda [3]).

Prueba La fórmula es válida para q = 1. En efecto por razones de simetría de Sn∫Sn

x2i =

1

n + 1

∫Sn

(n+1∑i=1

x2i

)=

1

n + 1vol (Sn) .

El multi-índice que corresponde a x2i es (2α1, · · · , 2αn+1) con αj = 0 si i 6=

i , αi = 1. Para este multi-índice el segundo miembro de la ecuación III.5da

(2α1)! (2α2)! · · · (2αn+1)!

H(2q, n)×H(2q − 1, n) · · ·H(2, n)

(q)!

(α1)! (α2)! · · · (αn+1)!× vol (Sn) =

(−1)2!

H(2, n)

1!

1!× vol (Sn) =

(−1) · 2−2 · (2− 1)− 2n

vol (Sn) =1

n + 1vol (Sn)

36

Suponemos ahora que la fórmula III.5 es válida para x2β1

1 x2β2

2 · · ·x2βn+1

3 , dondeq − 1 = β1 + β2 + · · · + βn+1, nuestra hipótesis de inducción, y tomemos f =

x2α11 x2α2

2 · · ·x2αn+1

3 con q = α1 + α2 + · · · + αn+1. Aplicando a este f la fórmula III.4conseguimos

∆Sn (f) = H(2q, n) · f +n+1∑i=1

(2αi) (2αi − 1)f

x2i

.

Integrando sobre la esfera Sn ambos miembros de esta ecuación sale

0 = H(2q, n)

∫Sn

f +n+1∑i=1

(2αi) (2αi − 1)

∫Sn

f

x2i

de donde,

∫Sn

f =−1

H(2q, n)

(n+1∑i=1

(2αi) (2αi − 1)

∫Sn

f

x2i

). (III.6)

Por la hipótesis de inducción∫Sn

f

x2i

= (−1)q−1 (2α1)! · (2α2)! · · · (2αn+1)!

2αi · (2αi − 1) ·H(2(q − 1), n) · · ·H(2, n)

× (q − 1)!αi

α1!α2! · · ·αn+1!.

Reemplazando en la fórmula III.6 resulta ’

∫Sn

f = (−1)q

n+1∑i=1

(2α1)! · · · (2αn+1)! (q − 1)! αi

H(2q, n) · · ·H(2, n) (α1)! (α2)! · · · (αn+1)!

= (−1)q (2α1)! · · · (2αn+1)! (q − 1)! q!

H(2q, n) · · ·H(2, n) (α1)! (α2)! · · · (αn+1)!

En particular cuando n = 2, es decir f = x2α11 x2α2

2 x2α33 resulta∫

Sn

f =(2α1)! (2α2)! · · · (2αn+1)!

(q + 1)!

(q)!

(α1)! (α2)! · · · (αn+1)!× vol (Sn)

lo que nos permite escribir la siguiente relación de recurrencia para coef (2α1, 2α2, 2α3) =∫Sn

f

coef (2α1, 2α2, 2α3) =2α1 − 1

q + 1× coef (2α1 − 2, 2α2, 2α3)

37

III.2. Ordenamiento de una base de Pk

Sabemos que dimP,k = (k+n−1)!k!(n−1)!

sea Jk,n el intervalo de enteros [1, dimPk] y definimosinductivamente

indic(k, n) : Jk,n → Ukn

siendo Ukn el conjunto de los multi-índices (α1, . . . , αn) con αi ≥ 0 y α1 + · · · + αn = k,

estipulando para 1 ≤ p ≤ dimPk

αn = sup j|p− dimPk + dimPk−j ≥ k − c + 1

(α1, . . . , αn−1) = indic (k− αn, n− 1) (p− dimPk + dimPk−αn − 1)

y poniendoindic(1, n)(p) = p

para1 ≤ p ≤ n = dimPk

esta función resulta entonces una biyección, por razonamiento inductivo. Su inversa

elp : Ukn → Jk,n

eselp (α1, . . . αn) = dimPk − dimPk−αn + elp (α1, . . . αn−1)

lo que se comprueba fácilmente.

III.3. Base para Hk

Corolario III.9 La dimensión de Hk es

dk =(2k + n− 1) (k + n− 2)!

k! (n− 1)!.

Prueba: por la proposición III.1, dk = dimPk − dimPk−2.

La matriz de ∆ : Pk → Pk−2 respecto de la base

(ui)dim(Pk)i=1 , (vj)

dim(Pk−2)j=1

38

donde

ui = xα11 · · ·xαn+1

n+1 , α = (α1, . . . , αn+1) = indic(k, n + 1)(i)

vi = xβ1

1 · · ·xβn+1

n+1 , β = (β1, . . . , βn+1) = indic(k− 2, n + 1)(i)

es la matriz M , que verifica

Mi,j =n+1∑t=1

αt (αt − 1) δ (α1, β1) · · · δ (αt−1, βt−1) δ (αt−2, αt)

× δ (αt+1, βt+1) · · · δ (αn+1, βn+1)

siendo δ el delta de Dirac. En efecto

∆(xα1

1 · · ·xαn+1

n+1

)=

n+1∑t=1

(αt − 1) αtxα11 · · ·xαt−2

t xαt+1

t+1 · · ·xαn+1

n+1 .

Habiendo determinado M , mediante el algoritmo de Gauss, obtenemos la correspondien-te matriz escalonada reducida M equivalente por filas. Así M esta estrictamente deter-minada. Designemos ip, p = 1, . . . , dimPk−2, el índice de columnas donde aparece elprimer término no nulo en M en la fila p. Sean jh, h = 1, . . . , dk los índices de columnasrestantes. para cada h = 1, . . . , dk, definimos

Yh =

dimPk∑i=1

[Yh]i xindic(k,i)

siendo[Yh]ip = −Mip,h, p = 1, . . . , dimPk−2

[Yh]jt= δ(t, h), t = 1, . . . , dk

III.4. Ejemplo (comparando una solución real y una cal-

culada)

Resolvemos el problema planteado en la sección II.7 por el método de Fourier con unabase de polinomios armónicos, y así comparar con la solución real. Para este propósito,se siguen las instrucciones enumeradas en el apéndice B para escribir un programa enlenguaje C que es mostrado en la lista III.1. Este programa genera las imágenes mostra-das en la figura III.1 y una salida con los resultados mostrada en la lista III.2.

39

Lista III.1: Programa en C para resolver una ecuación diferencial por el método Fourierpara una base de polinomios armónicos esféricos. Este código está en el archivo fuentefourier/prueba01.c

#include " shDefs . h "#include <math . h>#include < s t d i o . h>


double x = p . x , y = p . y , z = p . z ;return (−(8 − 24 ∗ pow ( x , 2) + 20 ∗ y − 120 ∗ pow ( z , 2) ∗ y ) )

;/ / . . . mas pol inomios para probarreturn −(12 ∗ pow ( x , 2) − 20 ∗ pow ( x , 4) ) ;return −(6 ∗ x − 12 ∗ pow ( x , 3) ) ;return −pow ( x , 2) ;

doubleu ( node n )

double x = n . x , y = n . y , z = n . z ;return (4 ∗ pow ( x , 2) + 10 ∗ pow ( z , 2) ∗ y ) −1.35;/ / . . . mas pol inomios para probarreturn (pow ( x , 4) ) ;return pow ( x , 3) ;return ( ( pow ( x , 2) − 1 / 3) / 6) ;

i n tmain ( )

tSur face esfera ;i n t i = 4 ;es fera = gen ( i ) ;/ / tSu r faceP lo t ( esfera , " good " ) ;

40

esfera = shSolve (3 , esfera , 0 .1 ) ;s h P r i n t S o l u t i o n ( esfera , 42) ;i f ( i == 4)

p lo tSo lu t i onOver ( esfera , " solRealPrueba01 . dat " , ’ r ’ ) ;p lo tSo lu t i onOver ( esfera , " solCalcPrueba01 . dat " , ’ c ’ ) ;p lo tSo lu t i onOver ( esfera , " solErrorPrueba01 . dat " , ’ e ’ ) ;

tSur faceDest roy ( es fera ) ;

Lista III.2: Resultados producido por el program prueba01

∗∗∗ BASE_1 ∗∗∗

MULTI−INDDICESELP ( A, B, C)

1 ( 1 , 0 , 0)2 ( 0 , 1 , 0)3 ( 0 , 0 , 1)

VECTORES BASE ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003

1 0 00 1 00 0 1

VECTORES BASE ORTHO−NORMALIZADOS ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003

1.73 0 00 1.73 00 0 1.73

COEFICIENTES DE FOURIERA_1,1= 6.62063e−16A_1,2= 2.28341A_1,3= −1.64108e−15

∗∗∗ BASE_2 ∗∗∗

MULTI−INDDICES

41

ELP ( A, B, C)1 ( 2 , 0 , 0)2 ( 1 , 1 , 0)3 ( 0 , 2 , 0)4 ( 1 , 0 , 1)5 ( 0 , 1 , 1)6 ( 0 , 0 , 2)

VECTORES BASE ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003 VEC # 0004 VEC # 0005

0 −1 0 0 −11 0 0 0 00 1 0 0 00 0 1 0 00 0 0 1 00 0 0 0 1

VECTORES BASE ORTHO−NORMALIZADOS ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003 VEC # 0004 VEC # 0005

0 −1.94 0 0 −1.123.87 0 0 0 0

0 1.94 0 0 −1.120 0 3.87 0 00 0 0 3.87 00 0 0 0 2.24

COEFICIENTES DE FOURIERA_2,1= −7.16453e−16A_2,2= −6.17155A_2,3= −1.18525e−15A_2,4= 1.48205e−15A_2,5= −3.56314

∗∗∗ BASE_3 ∗∗∗

MULTI−INDDICESELP ( A, B, C)

1 ( 3 , 0 , 0)2 ( 2 , 1 , 0)3 ( 1 , 2 , 0)4 ( 0 , 3 , 0)5 ( 2 , 0 , 1)6 ( 1 , 1 , 1)

42

7 ( 0 , 2 , 1)8 ( 1 , 0 , 2)9 ( 0 , 1 , 2)

10 ( 0 , 0 , 3)

VECTORES BASE ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003 VEC # 0004 VEC # 0005 VEC # 0006

−0.333 0 0 0 −0.333 00 −3 0 0 0 −11 0 0 −1 0 00 1 0 0 0 00 0 0 0 0 00 0 1 0 0 00 0 0 1 0 00 0 0 0 1 00 0 0 0 0 10 0 0 0 0 0

VEC # 000700

−30000001

VECTORES BASE ORTHO−NORMALIZADOS ∗∗∗VEC # 0001 VEC # 0002 VEC # 0003 VEC # 0004 VEC # 0005 VEC # 0006

−2.09 0 0 −1.21 −1.46 00 −6.27 0 0 0 −1.62

6.27 0 0 −1.21 −1.46 00 2.09 0 0 0 −1.620 0 0 0 0 00 0 10.2 0 0 00 0 0 4.83 −0.728 00 0 0 0 6.55 00 0 0 0 0 6.480 0 0 0 0 0

43

VEC # 0007−0.703

0−0.703

000

−4.22−0.703

02.58

COEFICIENTES DE FOURIERA_3,1= 1.14674e−15A_3,2= 0.000551746A_3,3= 1.71479e−15A_3,4= −1.91606e−15A_3,5= −3.00607e−15A_3,6= 14.7338A_3,7= 3.54941e−15

NUMERO DE NODOS=2562NUMERO DE TRIANGULOS=5120NUMERO DE Lados=7680

NODO SOL.REAL SOL.CAL ERROR1 1.544 1.555 0.010312 1.544 1.555 0.010313 1.544 1.555 0.010314 1.544 1.555 0.010315 −0.2444 −0.2369 0.007546 −0.2444 −0.2369 0.007547 −0.2444 −0.2168 0.027658 −0.2444 −0.2168 0.027659 2.454 2.45 0.00446

10 −5.154 −5.106 0.0486511 −5.154 −5.106 0.0486512 2.454 2.45 0.0044613 −1.35 −1.34 0.0103814 1.055 1.055 0.000190215 1.055 1.055 0.000190216 1.673 1.676 0.0035317 1.745 1.748 0.00289418 2.65 2.656 0.0058119 1.745 1.748 0.00289420 1.673 1.676 0.00353

44

21 1.055 1.055 0.000190222 1.055 1.055 0.000190223 −2.373 −2.34 0.0325224 −1.35 −1.328 0.0220925 1.673 1.676 0.0035326 −2.373 −2.34 0.0325227 1.745 1.748 0.00289428 1.745 1.748 0.00289429 2.65 2.656 0.0058130 1.673 1.676 0.0035331 −1.35 −1.328 0.0220932 −2.373 −2.34 0.0325233 −2.373 −2.34 0.0325234 0.7906 0.8106 0.0199835 −2.991 −2.95 0.0408936 0.7906 0.8106 0.0199837 −2.991 −2.95 0.0408938 −2.991 −2.95 0.0408939 0.7906 0.8106 0.0199840 0.7906 0.8106 0.0199841 −2.991 −2.95 0.0408942 −1.35 −1.316 0.03381

45

a) b)

c)

Figura III.1: a) Solución real, b) Solución calculada para un polinomio de tercer grado so-bre una malla de 5120 triángulos generada con gen(3), c) Desviación entre las solucionescalculada y real.

46

Capítulo IV

Conclusiones

Demostramos que es posible aproximar la solución de una ecuación diferencial dePoisson en una esfera tanto por el método de Galerkin como el método de Fourier. Loslogros en teoría más notables son una fórmula para calcular el producto interno de los gra-dientes de funciones nodales lo que hace posible independizarse de una parametrizaciónde la superficie y nos permite calcular de forma precisa los elementos de la matriz queresultan en el problema de Galerkin cuya utilidad se extendería al caso de una ecuaciónde Poisson no lineal. El método de Fourier nos llevo a demostrar una fórmula para calcu-lar la integral de un polinomio sobre la esfera que no hemos encontrado en la literatura yque es válida para una esfera de dimensión arbitraria. La librería de software desarrolladapermite resolver con poco esfuerzo una ecuación diferencial sobre la esfera en el caso li-neal y tan solo con modificaciones obvias se puede extender al caso no lineal. Además lalibrería es un punto de partida para el desarrollo de un software que pueda funcionar so-bre otras geometrías y resolver problemas de ecuaciones diferenciales de orden elípticoy parabólico con otras condiciones de contorno. El trabajo realizado generó necesidadesque marcan la pauta a seguir en el futuro inmediato. Es importante mejorar la precisiónde la integracción numérica sobre el elemento considerado, así como, diseñar un interfazde programación de aplicación adecuada para el manejo de mallas de triángulos másamistosa al usuario.

47

Apéndice A

Implementación del método deelementos finitos

Se implementa una librería que pone a posposición del usuario las siguientes funcio-nes:

SOLUCIÓN POR EL MÉTODO DE GALERKIN SOBRE LA SUPERFICIE TRIAN-GULADA

solveGalerkinOnThe_tSurface () Solución del problema de Direchlet ∆S2u = f

por el método de elementos finitos de Galerkin. Debe definir las funciones f yu antes.

Definición: tSurface solveGalerkinOnThe_tSurface (tSurface esfera, int lastNode-ToBePrinted);

Regresa: La solución por el método de elementos finitos de Galerkin en forma deun arreglo contenido dentro de la estructura de datos de una malla de triángulos(tSurface).

Argumentos:

esfera <input> (tSurface) Una estructura de datos de una malla de triángulos(tSurface).

lastNodeToBePrinted <input> (int) Imprime la solución hasta alcanzar el nodolastNodeToBePrinted.

LAS FUNCIONES F Y U

Definición: double f (vector p) double u (node n)

Regresa: el valor correspondiente en una posición (x,y,z)

48

p <input> (vector) estructura de datos que contiene las coordenadas (p.x,p.y,p.z)

n <input> (node) estructura de datos que contiene además de las coordenadas(n.x,n.y,n.z) el índice del nodo (n.index).

GENERANDO LA ESTRUCTURA DE LA MALLA DE TRIÁNGULOS PARA UNASUPERFICIE

gen() genera una primera superficie triangulada. Luego toma cada triángulo en estasuperficie y lo divide en cuatro nuevos triángulos el numero de veces especi-ficado por el usuario. Así que si empezamos con 20 triángulos y queremosrefinar esta estructura 3 veces obtendremos 20*4*4*4.

tSurface gen(n)

Regresa

output (tSurface) una estructura de datos de la superficie triangulada

Argumentos

n input (int) el número de veces que quiere refinar la superficie triangulada. Cadavez que se refine un triángulo se obtienen cuatro nuevos triángulos.

Lista A.1: Archivo de cabecera que contiene la declaración de las funciones para el usua-rio pMedio() integrateOver() area(). Este código está en el archivo fuente galerkin/tSurfa-ce.h

vec to r pMedio ( node , node ) ;double area ( node , node , node ) ;double i n tegra teOver ( tSur face ) ;double c u b i c I n t e g r a t i o n ( node , node , node , int , double ) ;

Lista A.2: Módulo en C que define las funciones para el usuario pMedio() integrateOver()area(). Este código está en el archivo fuente galerkin/integrate.c

#include " tSur face . h "#include " so lveGa le rk in . h "#include < s t d i o . h>#include <math . h>

doublearea ( node n1 , node n2 , node n3 )

49

vec to r v23 = getVector ( n2 , n3 ) , v21 = getVector ( n2 , n1 ) ;double uv = ( v23 . x ∗ v21 . x + v23 . y ∗ v21 . y + v23 . z ∗ v21 . z ) ,

uu = ( v21 . x ∗ v21 . x + v21 . y ∗ v21 . y + v21 . z ∗ v21 . z ) ,vv = ( v23 . x ∗ v23 . x + v23 . y ∗ v23 . y + v23 . z ∗ v23 . z ) ;

return ( s q r t ( uu ∗ vv − pow ( uv , 2) ) / 2) ;

vec to rpMedio ( node n1 , node n2 )

vec to r pMedio ;pMedio . x = ( n1 . x + n2 . x ) / 2 ;pMedio . y = ( n1 . y + n2 . y ) / 2 ;pMedio . z = ( n1 . z + n2 . z ) / 2 ;return ( pMedio ) ;

doublei n tegra teOver ( tSur face esfera )

double i n t e g r a l , A ;i n t i ;node n1 , n2 , n3 ;i n t e g r a l = 0 ;for ( i = 0 ; i < es fera . nTr iang les ; i ++)

n1 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 0 ] , es fera ) ;n2 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 1 ] , es fera ) ;n3 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 2 ] , es fera ) ;A = area ( n1 , n2 , n3 ) ;

i n t e g r a l += ( f ( pMedio ( n1 , n3 ) ) +f ( pMedio ( n1 , n2 ) ) + f ( pMedio ( n2 , n3 ) ) ) ∗ A

/ 3 . 0 ;

p r i n t f ( " \ nINTEGRAL DE f SOBRE LA ESFERA= %g \ n " , i n t e g r a l ) ;return i n t e g r a l ;

50

tSur facec u b i c I n t e g r a t i o n ( node n1 , node n2 , node n3 , i n t maxit , double es )

tSur face t r i a n g u l o = T_SURFACE_EMPTY;

/ / t r i a n g u l o = t S u r f a c e I n i t i a t e ( t r i a n g u l o , 3 , 1 , 3) ;

nodeAdd (1 , t r i a n g u l o . nodes , n1 . x , n1 . y , n1 . z ) ; / / ∗∗∗∗∗∗nodeAdd (2 , t r i a n g u l o . nodes , n2 . x , n2 . y , n2 . z ) ;nodeAdd (3 , t r i a n g u l o . nodes , n3 . x , n3 . y , n3 . z ) ;

t r i ang leAdd (1 , t r i a n g u l o . t r i a n g l e s , 1 , 2 , 3) ;

return ( t r i a n g u l o ) ;

doublei n t eg ra teOverT r i ang le ( node n1 , node n2 , node n3 )

double i n t e g r a l , A ;

Lista A.3: Archivo de cabecera que contiene la declaración de las funciones para el usua-rio f() u() solveGalerkinOnThe tSurface(). Este código está en el archivo fuente galerkin/-solveGalerkin.h

tSur face nlso lveGalerk inOnThe_tSurface ( tSur face , int , i n t ) ;tSur face solveGalerkinOnThe_tSurface ( tSur face , i n t ) ;double f ( vec to r ) ;double u ( node ) ;double k ( double ) ;

Lista A.4: Módulo en C para resolver una ecuación diferencial por el método de elementosfinitos de Galerkin. Este código está en el archivo fuente galerkin/solveGalerkin.c

#include < s t d i o . h>#include " tSur face . h "

51

#include " i n t e g r a t e . h "#include " so lveGa le rk in . h "#include <spConfig . h>#include <spmatr ix . h>#include <math . h>

doubleiP roduc t ( vec to r v1 , vec to r v2 )

return ( v1 . x ∗ v2 . x + v1 . y ∗ v2 . y + v1 . z ∗ v2 . z ) ;

/∗SOLUCIÓN POR EL MÉTODO DE GALERKIN SOBRE LA SUPERFICIE

TRIANGULADA

solveGalerkinOnThe_tSurface ( )Soluc ión de l problema de D i r i c h l e t \ Delta_S_2 U=F por e l

método deelementos f i n i t o s de Galerk in . Debe d e f i n i r l as func iones F y

U antes .

>>> D e f i n i c i ó n :tSur face solveGalerkinOnThe_tSurface ( tSur face esfera , i n t

lastNodeToBePrinted ) ;

>>> Regresa :La so luc ión por e l método de elementos f i n i t o s de Galerk in en

formade un a r reg lo contenido dentro de l a e s t r u c t u r a de datos de

unamal la de t r i á n g u l o s ( tSur face ) .

>>> Argumentos :

es fera < input > ( tSur face )Una e s t r u c t u r a de datos de una mal la de t r i á n g u l o s ( tSur face )

.

52

lastNodeToBePrinted < input > ( i n t )Imprime l a so luc ión hasta a lcanzar e l nodo

lastNodeToBePrinted .

LAS FUNCIONES F Y U

>>> D e f i n i c i ó n :double f ( vec to r p )double u ( node n )

>>> Regresa :e l va l o r cor respondiente en una pos ic ión ( x , y , z )

p < input > ( vec to r )e s t r u c t u r a de datos que cont iene las coordenadas ( p . x , p . y , p . z

)

n < input > ( node )e s t r u c t u r a de datos que cont iene además de las coordenadas( n . x , n . y , n . z ) e l í nd i ce de l nodo ( n . index ) .

∗ /tSur facesolveGalerkinOnThe_tSurface ( tSur face esfera , i n t

lastNodeToBePrinted )

char ∗M;unsigned long nodeOmitted ;i n t e r r o r ;unsigned long i , j ;node n1 , n2 , n3 ;vec to r v23 , v32 , v31 , v13 , v12 , v21 ;

53

double A, k ;double normaL2 ;spREAL rhs [ es fera . nNodes − 1 ] ;spREAL rhs_ rea l [ es fera . nNodes − 1 ] ;spREAL so l_ca l c [ es fera . nNodes − 1 ] ;spREAL s o l _ r e a l [ es fera . nNodes − 1 ] ;

for ( i = 0 ; i < es fera . nNodes − 1; i ++) / / . . . I n i c i a n d ovectores

rhs [ i ] = 0 ;rhs_ rea l [ i ] = 0 ;so l_ca l c [ i ] = 0 ;s o l _ r e a l [ i ] = 0 ;

/ / p r i n t f ( " Número de nodos= %d \ n " , es fera . nNodes−1) ;M = spCreate ( es fera . nNodes − 1 , 0 , &e r r o r ) ;spClear (M) ;

nodeOmitted = esfera . nNodes ;k = in tegra teOver ( es fera ) / (4 ∗ M_PI ) ;

for ( i = 0 ; i < es fera . nTr iang les ; i ++)

n1 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 0 ] , es fera ) ;n2 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 1 ] , es fera ) ;n3 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 2 ] , es fera ) ;A = area ( n1 , n2 , n3 ) ;

v13 = getVector ( n1 , n3 ) ;v31 = getVector ( n3 , n1 ) ;v32 = getVector ( n3 , n2 ) ;v23 = getVector ( n2 , n3 ) ;v21 = getVector ( n2 , n1 ) ;v12 = getVector ( n1 , n2 ) ;

i f ( ! ( n1 . index == nodeOmitted ) )

54

spADD_REAL_ELEMENT ( spGetElement (M, n1 . index , n1 . index ) ,iP roduc t ( v23 , v23 ) / (4 ∗ A) ) ;

i f ( ! ( n1 . index == nodeOmitted | | n2 . index == nodeOmitted ) )spADD_REAL_ELEMENT ( spGetElement (M, n1 . index , n2 . index ) ,

−iP roduc t ( v31 , v32 ) / (4 ∗ A) ) ;i f ( ! ( n1 . index == nodeOmitted | | n3 . index == nodeOmitted ) )

spADD_REAL_ELEMENT ( spGetElement (M, n1 . index , n3 . index ) ,−iP roduc t ( v21 , v23 ) / (4 ∗ A) ) ;


−iP roduc t ( v32 , v31 ) / (4 ∗ A) ) ;i f ( ! ( n2 . index == nodeOmitted ) )



−iP roduc t ( v12 , v13 ) / (4 ∗ A) ) ;i f ( ! ( n3 . index == nodeOmitted | | n1 . index == nodeOmitted ) )

spADD_REAL_ELEMENT ( spGetElement (M, n3 . index , n1 . index ) ,−iP roduc t ( v23 , v21 ) / (4 ∗ A) ) ;


−iP roduc t ( v13 , v12 ) / (4 ∗ A) ) ;i f ( ! ( n3 . index == nodeOmitted ) )


/∗ . . . RHS element∗ numer ica l i n t e g r a t i o n ∗ /

i f ( ! ( n1 . index == nodeOmitted ) )rhs [ n1 . index − 1] += ( f ( pMedio ( n1 , n3 ) ) +

f ( pMedio ( n1 , n2 ) ) − 2 ∗ k ) ∗ A /6 . 0 ;

i f ( ! ( n2 . index == nodeOmitted ) )rhs [ n2 . index − 1] += ( f ( pMedio ( n2 , n3 ) ) +

f ( pMedio ( n2 , n1 ) ) − 2 ∗ k ) ∗ A /6 . 0 ;

i f ( ! ( n3 . index == nodeOmitted ) )

55

rhs [ n3 . index − 1] += ( f ( pMedio ( n3 , n1 ) ) +f ( pMedio ( n3 , n2 ) ) − 2 ∗ k ) ∗ A /

6 . 0 ;

/ / . . . Store the r e a l s o l u t i o nesfera . so lReal = vCreate_dbl ( es fera . nNodes − 1) ;for ( i = 1 ; i <= esfera . nNodes − 1; i ++)

s o l _ r e a l [ i − 1] = esfera . so lReal [ i − 1] = u ( nodeGet ( i ,es fera ) ) ;

#define MULTIPLY_TO_SEE_WHAT_HAPEN# i f d e f MULTIPLY_TO_SEE_WHAT_HAPEN

Crea te In te rna lVec to rs (M) ;s p M u l t i p l y (M, rhs_rea l , s o l _ r e a l ) ;

#endif

/ / spP r i n t (M,0 ,1 ,1 ) ;

#define FACTOR_AND_SOLVE# i f d e f FACTOR_AND_SOLVE

e r r o r = spFactor (M) ;spSolve (M, rhs , so l_ca l c ) ;

/ / The next f i v e l i n e s w r i t e s the l a s t nodes coord ina tesp r i n t f ( "NUMERO DE TRIÁNGULOS= %d \ n " , es fera . nTr iang les ) ;p r i n t f ( "NUMERO DE LADOS= %d \ n " , es fera . nSides ) ;double valorUlt imoNodo ;n1 = nodeGet ( es fera . nNodes , es fera ) ;p r i n t f ( " Ul t imo Nodo= %d \ tU ( x= %g , y= %g , z= %g )= %g \ n " ,

n1 . index , n1 . x , n1 . y , n1 . z , valorUl t imoNodo = u ( n1 ) ) ;

for ( i = 1 ; i <= esfera . nNodes − 1; i ++)so l_ca l c [ i − 1] += valorUl t imoNodo ;

es fera . so lCalc = vCreate_dbl ( es fera . nNodes − 1) ;for ( i = 1 ; i <= esfera . nNodes − 1; i ++)

56

esfera . so lCalc [ i − 1] = so l_ca l c [ i − 1 ] ;#endif

esfera . e r r o r = vCreate_dbl ( es fera . nNodes − 1) ;for ( i = 1 ; i <= esfera . nNodes − 1; i ++)

es fera . e r r o r [ i − 1] =fabs ( fabs ( es fera . so lReal [ i − 1 ] ) − fabs ( es fera . so lCa lc [ i

− 1 ] ) ) ;

/ / . . . E r ro r en Norma L2normaL2 = 0;for ( i = 1 ; i <= esfera . nNodes − 1; i ++)

normaL2 += esfera . e r r o r [ i − 1 ] ;normaL2 ∗= s q r t (2 ∗ M_PI / 3 / es fera . nNodes ) ;

p r i n t f ( "ERROR EN NORMA L2= %g \ n " , normaL2 ) ;

#define PRINT_SOLUTION# i f d e f PRINT_SOLUTION

i f ( lastNodeToBePrinted >= esfera . nNodes − 1)

p r i n t f ( " %5s %12s %12s %12s %12s %12s \ n " , " NODO" ,"SOL REAL" , "SOL CAL" , "RHS_INT" , "RHS_REAL" , "

ERROR" ) ;for ( i = 1 ; i <= esfera . nNodes − 1; i ++)

p r i n t f ( " %5d %12g %12g %12g %12g %12g \ n " , i ,

es fera . so lReal [ i − 1 ] , es fera . so lCalc [ i − 1 ] ,rhs [ i − 1 ] ,

rhs_ rea l [ i − 1 ] , es fera . e r r o r [ i − 1 ] ) ;

else

p r i n t f ( " %5s %12s %12s %12s %12s %12s \ n " , " NODO" ,

"SOL.REAL" , "SOL.CAL" , "RHS_INT" , "RHS_REAL" , "ERROR" ) ;

for ( i = 1 ; i <= lastNodeToBePrinted ; i ++)

57

p r i n t f ( " %5d %12.4g %12.4g %12.4g %12.4g %12.4g \ n "

, i ,es fera . so lReal [ i − 1 ] , es fera . so lCalc [ i − 1 ] ,

rhs [ i − 1 ] ,rhs_ rea l [ i − 1 ] , es fera . e r r o r [ i − 1 ] ) ;

#endifreturn ( es fera ) ;

Lista A.5: Archivo de cabecera que contiene la declaración de las funciones para el usua-rio gen() tSurfacePlot() plotSolutionOver() y la estructura de datos tSurface. Este códigoestá en el archivo fuente galerkin/tSurface.h

typedef double ∗∗mMatr ix_dbl ;typedef unsigned long ∗∗mMatr ix_ in t ;typedef double ∗ vec to rAr ray_o f_db l ;

typedef struct

mMatr ix_dbl nodes ;mMat r ix_ in t t r i a n g l e s ;mMat r ix_ in t s ides ;vec to rAr ray_o f_db l so lReal ;vec to rAr ray_o f_db l so lCalc ;vec to rAr ray_o f_db l e r r o r ;unsigned long nSides ;unsigned long nTr iang les ;unsigned long nNodes ;

tSur face ;

#define T_SURFACE_EMPTY NULL,NULL,NULL,NULL,NULL,NULL,0 ,0 ,0

typedef struct

unsigned long index ;double x , y , z ;

58

node ;

typedef struct

node n1 ;node n2 ;node n3 ;

t r i a n g l e ;

typedef struct

double x , y , z ; vec to r ;

tSur face gen ( i n t ) ;tSur face f s t ( void ) ;tSur face gntSurface ( tSur face ) ;void tSuf raceDest roy ( tSur face ) ;void tSu r faceP lo t ( tSur face , char ∗ ) ;node nodeGet ( unsigned long , tSur face ) ;vec to r getVector ( node , node ) ;void p lo tSo lu t i onOver ( tSur face , char ∗ , char ) ;vec to rAr ray_o f_db l vCreate_dbl ( unsigned long ) ;

#define HELP_TO_PLOT " #GNUPLOT f i l e \ n \#∗∗∗ EJEMPLO ∗∗∗ \ n \# . . . l a es fera es guardada en e l a rch ivo \ " hsur face \ " \ n \#set hidden noo f f se t \ n \#set isosamples 100\n \#set s ize square \ n \#set t i c s l e v e l 0 \ n \#set s t y l e data l i ne \ n \#set cbrange [ . 5 : 2 ] \ n \#set xrange [ −1 :1 ] \ n \#set yrange [ −1 :1 ] \ n \#set view 60 ,30 ,1 ,1 .5 \ n \#set t e rm ina l X11 \ n \

59

#unset pm3d \ n \# s p l o t \ " hsur face \ " n o t i t l e w l i n e lw 2 l t 3 , exp (20∗ ( s q r t (1−x

∗∗2−y ∗∗2)−1/(1−x∗∗2−y ∗∗2) ) ) + s q r t (.95−x∗∗2−y ∗∗2)−0.05 n o t i t l e wl i n e l t 2 \ n \

#pause −1 \ " H i t r e t u r n to cont inue \ " \ n \#set pm3d \ n \# s p l o t \ " hsur face \ " n o t i t l e w l i n e lw 2 l t 3 , exp (20∗ ( s q r t (1−x

∗∗2−y ∗∗2)−1/(1−x∗∗2−y ∗∗2) ) ) + s q r t (1−x∗∗2−y ∗∗2) n o t i t l e w l i n elw 2 l t 0 \ n \

# \ n \# # . . . GUARDANDO COMO PNG\ n \#set t e rm ina l png s ize 1000,1000 crop enhanced \ n \#set output \ " d i s t r i b u t i o n −med. png \ " \ n \# s p l o t \ " hsur face \ " n o t i t l e w l i n e lw 2 l t 3 , exp (20∗ ( s q r t (1−x

∗∗2−y ∗∗2)−1/(1−x∗∗2−y ∗∗2) ) ) + s q r t (1−x∗∗2−y ∗∗2) n o t i t l e w l i n elw 2 l t 0 \ n \

# \ n \#set t e rm ina l png s ize 1000,1000 crop enhanced \ n \#set output \ " d i s t r i b u t i o n −b ig . png \ " \ n \# s p l o t \ " hsur face \ " n o t i t l e w l i n e lw 2 l t 3 , exp (20∗ ( s q r t (1−x

∗∗2−y ∗∗2)−1/(1−x∗∗2−y ∗∗2) ) ) + s q r t (1−x∗∗2−y ∗∗2) n o t i t l e w l i n elw 2 l t 0 \ n \ "

Lista A.6: Módulo en C que contiene la definición de las funciones necesarias para gene-rar una malla de triángulos esférica. Este código está en el archivo fuente ../galerkin/tSur-face.c

#include <stdarg . h>#include < s t d i o . h>#include < s t d l i b . h>#include <math . h>#include " tSur face . h "

vec to rAr ray_o f_db lvCreate_dbl ( unsigned long m)

vec to rAr ray_o f_db l vec to r ;

60

unsigned long i ;/ / Reserva memoria para e l vec to rvec to r = c a l l o c (m, sizeof ( double ) ) ;return ( vec to r ) ;

mMatr ix_dblmCreate_dbl ( unsigned long m, unsigned long n )

mMatr ix_dbl mat r i x ;unsigned long i ;/ / Reserva memoria para e l vec to r de punteros : mat r i x [ ]mat r i x = c a l l o c (m, sizeof ( double ∗ ) ) ;/ / Reserva memoria para toda l a mat r i z : mat r i x [ ] [ ]mat r i x [ 0 ] = c a l l o c (m ∗ n , sizeof ( double ) ) ;/ / Se asigna e l va l o r para los elementos de l vec to r de punteros

mat r i x [ ]for ( i = 1 ; i < m; i ++)

mat r i x [ i ] = mat r i x [ i − 1] + n ;return ( mat r i x ) ;

mMat r ix_ in tmCreate_int ( unsigned long m, unsigned long n )

mMat r ix_ in t mat r i x ;unsigned long i ;/ / Reserva memoria para e l vec to r de punteros : mat r i x [ ]mat r i x = c a l l o c (m, sizeof ( unsigned long ∗ ) ) ;/ / Reserva memoria para toda l a mat r i z : mat r i x [ ] [ ]mat r i x [ 0 ] = c a l l o c (m ∗ n , sizeof ( unsigned long ) ) ;/ / Se asigna e l va l o r para los elementos de l vec to r de punteros

mat r i x [ ]for ( i = 1 ; i < m; i ++)

mat r i x [ i ] = mat r i x [ i − 1] + n ;return ( mat r i x ) ;

61

voidnodeAdd ( unsigned long nFi la , mMatr ix_dbl node , double x , double

y , double z )

double norm = s q r t (pow ( x , 2) + pow ( y , 2) + pow ( z , 2) ) ;−−nF i l a ;node [ nF i l a ] [ 0 ] = x / norm ;node [ nF i l a ] [ 1 ] = y / norm ;node [ nF i l a ] [ 2 ] = z / norm ;

/∗nodeCopy ( )copy a vec to r row node from mOldNode mat r i x to mNewNode mat r i x

∗ /voidnodeCopy ( unsigned long nFi la , mMatr ix_dbl mOldNode , mMatr ix_dbl

mNewNode)−−nF i l a ;unsigned long i = 0 ;for ( i = 0 ; i < 3 ; i ++)

mNewNode[ nF i l a ] [ i ] = mOldNode [ nF i l a ] [ i ] ;

voidt r i ang leAdd ( unsigned long nFi la ,

mMat r ix_ in t t r i a n g l e ,unsigned long i , unsigned long j , unsigned long k )

−−nF i l a ;t r i a n g l e [ nF i l a ] [ 0 ] = i ;t r i a n g l e [ nF i l a ] [ 1 ] = j ;t r i a n g l e [ nF i l a ] [ 2 ] = k ;

/∗t S u r f a c e I n i t i a t e ( surface , nNodes , nTr iangles , nSides )

62

i n i t i a t e s an empty t r i a n g u l a t e d sur face data s t r u c t u r e ( sur face) wh i t

the number o f nodes ( nNodes ) the number o f t r i a n g l e s (nTr iang les ) and

the number o f s ides ( nSides ) . The s ize arguments nNodes ,nTr iang les and

nSides are ignored i f the t r i a n g u l a t e d sur face data s t r u c t u r ei s not

empty and only the ( surface , argument i s used to re− i n i t i a t e toa new

s ize the t r i a n g u l a t e d sur face data s t r u c t u r e .

t S u r f a c e I n i t i a t e ( sur face )re− i n i t i a t e s an non empty t r i a n g u l a t e d sur face data s t r u c t u r e (

sur face )using the o ld s ize to accomplish i t .

tSur face t S u r f a c e I n i t i a t e ( surface , nNodes , nTr iangles , nSides )tSur face t S u r f a c e I n i t i a t e ( sur face )

RETURNED ( tSur face )a t r i a n g u l a t e d sur face data s t r u c t u r e .

ARGUMENTS

sur face inpu t ( tSur face )a t r i a n g u l a t e d sur face data s t r u c t u r e

nNodes inpu t ( unsigned long )the number o f nodes used to i n i t i a t e the sur face data

s t r u c t u r e .

nTr iang les i npu t ( unsigned long )the number o f t r i a n g l e s used to i n i t i a t e the sur face data

s t r u c t u r e .

nSides i npu t ( unsigned long )

63

the number o f s ides used to i n i t i a t e the sur face datas t r u c t u r e .

∗ /tSur facet S u r f a c e I n i t i a t e ( tSur face old , . . . )

unsigned long i ;tSur face new ;

i f ( o ld . nNodes == 0 && old . nTr iang les == 0 && old . nSides == 0)

v a _ l i s t args ;va_s ta r t ( args , o ld ) ;new . nNodes = va_arg ( args , i n t ) ;new . nTr iang les = va_arg ( args , i n t ) ;new . nSides = va_arg ( args , i n t ) ;va_end ( args ) ;new . nodes = mCreate_dbl (new . nNodes , 3) ;new . t r i a n g l e s = mCreate_int (new . nTr iangles , 3) ;new . s ides = mCreate_int (new . nSides , 2) ;

else

new . nSides = 3 ∗ o ld . nTr iang les + 2 ∗ o ld . nSides ;new . nTr iang les = 4 ∗ o ld . nTr iang les ;new . nNodes = o ld . nNodes + o ld . nSides ;

new . nodes = mCreate_dbl (new . nNodes , 3) ;for ( i = 1 ; i <= o ld . nNodes ; i ++)

nodeCopy ( i , o ld . nodes , new . nodes ) ;

new . t r i a n g l e s = mCreate_int (new . nTr iangles , 3) ;new . s ides = mCreate_int (new . nSides , 2) ;

new . solReal=NULL ;new . so lCalc=NULL ;new . e r r o r =NULL ;

64

return (new) ;

voidtSur faceDest roy ( tSur face toBeDestroyed )

f r ee ( toBeDestroyed . nodes [ 0 ] ) ;f r ee ( toBeDestroyed . nodes ) ;f r ee ( toBeDestroyed . t r i a n g l e s [ 0 ] ) ;f r ee ( toBeDestroyed . t r i a n g l e s ) ;f r ee ( toBeDestroyed . s ides [ 0 ] ) ;f r ee ( toBeDestroyed . s ides ) ;i f ( toBeDestroyed . solReal != NULL)

f r ee ( toBeDestroyed . solReal ) ;i f ( toBeDestroyed . so lCalc != NULL)

f r ee ( toBeDestroyed . so lCalc ) ;i f ( toBeDestroyed . e r r o r != NULL)

f r ee ( toBeDestroyed . e r r o r ) ;

/∗locateNodeBeween ( ) checks whether a s ide was inc luded i n the

sur faceand inc ludes i t i f i t i sn ’ t . The n1 and n2 nodes both belong a

group ofth ree nodes stored i n a t r i a n g l e ar ray . This f u n c t i o n i s used

bygenerateSides ( ) .

tSur face locateNodeBetween ( n1 , n2 , old , new)

RETURNEDA tSur face s t r u c t u r e data .

ARGUMENTS

n1 , n2 inpu t ( unsigned long )

65

two of the th ree nodes inc luded i n the row of a t r i a n g l emat r i x t h a t

w i l l be s tored independent ly as a s ide i n a row of anothermat r i x

o f s ides . Both the t r i a n g l e mat r i x [ ] [ 3 ] and de s ide mat r i x[ ] [ 2 ] are

pa r t o f a tSur face s t r u c t u r e passed as argument .

o ld i npu t ( tSur face )a s t r u c t u r e t h a t conta ins the o ld t r i a n g u l a t e d sur face .

∗ /unsigned longlocateNodeBetween ( unsigned long n1 ,

unsigned long n2 , tSur face old , tSur face new)

double x , y , z ;unsigned long j = 0 ;while ( 1 )

i f ( ( n1 == old . s ides [ j ] [ 0 ] && n2 == old . s ides [ j ] [ 1 ] ) | |

( n2 == old . s ides [ j ] [ 0 ] && n1 == old . s ides [ j ] [ 1 ] ) )return ( j + 1 + o ld . nNodes ) ;

else i f ( o ld . s ides [ j ] [ 0 ] == 0 && old . s ides [ j ] [ 1 ] == 0)

o ld . s ides [ j ] [ 0 ] = n1 ;o ld . s ides [ j ] [ 1 ] = n2 ;x = ( o ld . nodes [ n1 − 1 ] [ 0 ] + o ld . nodes [ n2 − 1 ] [ 0 ] ) / 2 ;y = ( o ld . nodes [ n1 − 1 ] [ 1 ] + o ld . nodes [ n2 − 1 ] [ 1 ] ) / 2 ;z = ( o ld . nodes [ n1 − 1 ] [ 2 ] + o ld . nodes [ n2 − 1 ] [ 2 ] ) / 2 ;nodeAdd ( j + 1 + o ld . nNodes , new . nodes , x , y , z ) ;return ( j + 1 + o ld . nNodes ) ;

j ++;

tSur facef s t ( void )

66

/∗ Bu i ld an icosahedron from three mutua l l y perpend icu la r

golden rec tang les∗ /

tSur face icosahedron=T_SURFACE_EMPTY;double phi = (1 + pow (5 , 0 .5 ) ) / 2 . 0 ;

icosahedron = t S u r f a c e I n i t i a t e ( icosahedron , 12 , 20 , 30) ;

nodeAdd (1 , icosahedron . nodes , −phi , 0 , −1) ; / / ∗∗∗∗∗∗nodeAdd (2 , icosahedron . nodes , −phi , 0 , 1) ;nodeAdd (3 , icosahedron . nodes , phi , 0 , 1) ;nodeAdd (4 , icosahedron . nodes , phi , 0 , −1) ;nodeAdd (5 , icosahedron . nodes , −1, phi , 0) ; / / ∗∗∗∗∗∗nodeAdd (6 , icosahedron . nodes , 1 , phi , 0) ;nodeAdd (7 , icosahedron . nodes , 1 , −phi , 0) ;nodeAdd (8 , icosahedron . nodes , −1, −phi , 0) ;nodeAdd (9 , icosahedron . nodes , 0 , 1 , ph i ) ; / / ∗∗∗∗∗∗nodeAdd (10 , icosahedron . nodes , 0 , −1, ph i ) ;nodeAdd (11 , icosahedron . nodes , 0 , −1, −phi ) ;nodeAdd (12 , icosahedron . nodes , 0 , 1 , −phi ) ;

t r i ang leAdd (1 , icosahedron . t r i a n g l e s , 5 , 6 , 9) ;t r i ang leAdd (2 , icosahedron . t r i a n g l e s , 5 , 9 , 2) ;t r i ang leAdd (3 , icosahedron . t r i a n g l e s , 5 , 2 , 1) ;t r i ang leAdd (4 , icosahedron . t r i a n g l e s , 5 , 1 , 12) ;t r i ang leAdd (5 , icosahedron . t r i a n g l e s , 5 , 12 , 6) ;t r i ang leAdd (6 , icosahedron . t r i a n g l e s , 1 , 11 , 12) ;t r i ang leAdd (7 , icosahedron . t r i a n g l e s , 11 , 12 , 4) ;t r i ang leAdd (8 , icosahedron . t r i a n g l e s , 12 , 4 , 6) ;t r i ang leAdd (9 , icosahedron . t r i a n g l e s , 4 , 6 , 3) ;t r i ang leAdd (10 , icosahedron . t r i a n g l e s , 6 , 3 , 9) ;t r i ang leAdd (11 , icosahedron . t r i a n g l e s , 3 , 9 , 10) ;t r i ang leAdd (12 , icosahedron . t r i a n g l e s , 9 , 10 , 2) ;t r i ang leAdd (13 , icosahedron . t r i a n g l e s , 10 , 2 , 8) ;t r i ang leAdd (14 , icosahedron . t r i a n g l e s , 2 , 8 , 1) ;t r i ang leAdd (15 , icosahedron . t r i a n g l e s , 8 , 1 , 11) ;t r i ang leAdd (16 , icosahedron . t r i a n g l e s , 7 , 11 , 4) ;

67

t r i ang leAdd (17 , icosahedron . t r i a n g l e s , 7 , 4 , 3) ;t r i ang leAdd (18 , icosahedron . t r i a n g l e s , 7 , 3 , 10) ;t r i ang leAdd (19 , icosahedron . t r i a n g l e s , 7 , 10 , 8) ;t r i ang leAdd (20 , icosahedron . t r i a n g l e s , 7 , 8 , 11) ;

return ( icosahedron ) ;

tSur facegntSurface ( tSur face o ld )

unsigned long i = 0 , j = 1 ;unsigned long n1 , n2 , n3 , n4 , n5 , n6 ;tSur face new ;

new = t S u r f a c e I n i t i a t e ( o ld ) ;

for ( i = 1 ; i <= o ld . nTr iang les ; i ++)

n1 = o ld . t r i a n g l e s [ i − 1 ] [ 0 ] ;n3 = o ld . t r i a n g l e s [ i − 1 ] [ 1 ] ;n5 = o ld . t r i a n g l e s [ i − 1 ] [ 2 ] ;n2 = locateNodeBetween ( n1 , n3 , old , new) ;n4 = locateNodeBetween ( n3 , n5 , old , new) ;n6 = locateNodeBetween ( n1 , n5 , old , new) ;

t r i ang leAdd ( j , new . t r i a n g l e s , n1 , n2 , n6 ) ;t r i ang leAdd (++ j , new . t r i a n g l e s , n2 , n3 , n4 ) ;t r i ang leAdd (++ j , new . t r i a n g l e s , n6 , n2 , n4 ) ;t r i ang leAdd (++ j , new . t r i a n g l e s , n6 , n4 , n5 ) ;++ j ;

tSur faceDest roy ( o ld ) ;

return (new) ;

68

/∗GENERANDO LA ESTRUCTURA DE LA MALLA DE TRIÁNGULOS PARA UNA

SUPERFICIE

gen ( )genera una pr imera s u p e r f i c i e t r i angu lada . Luego toma cada

t r i á n g u l oen esta s u p e r f i c i e y l o d i v i d e en cuat ro nuevos t r i á n g u l o s e l

númerode veces espec i f i cado por e l usuar io . Así que s i empezamos con

20t r i á n g u l o s y queremos r e f i n a r esta e s t r u c t u r a 3 veces

obtendremos20∗4∗4∗4.

tSur face gen ( n )

RETURNS output ( tSur face )una e s t r u c t u r a de datos de l a s u p e r f i c i e t r i angu lada

ARGUMENTSn inpu t ( i n t )e l número de veces que qu iere r e f i n a r l a s u p e r f i c i e

t r i angu lada . Cadavez que se r e f i n e un t r i á n g u l o se obt ienen cuat ro nuevos

t r i á n g u l o s .∗ /

tSur facegen ( i n t n )

tSur face new ;i n t i ;

new = f s t ( ) ;for ( i = 1 ; i <= n ; i ++)

new = gntSurface (new) ;return (new) ;

69

vec to rgetVector ( node n1 , node n2 )

vec to r d i f f ;d i f f . x = n2 . x − n1 . x ;d i f f . y = n2 . y − n1 . y ;d i f f . z = n2 . z − n1 . z ;return ( d i f f ) ;

nodenodeGet ( unsigned long nodeNumber , tSur face sur face )

node n ;n . index = nodeNumber ;n . x = sur face . nodes [ nodeNumber − 1 ] [ 0 ] ;n . y = sur face . nodes [ nodeNumber − 1 ] [ 1 ] ;n . z = sur face . nodes [ nodeNumber − 1 ] [ 2 ] ;return n ;

voidtSu r faceP lo t ( tSur face toBePloted , char ∗ f i leName )

unsigned long i , j ;node n1 , n2 , n3 ;FILE ∗ f i l e ;f i l e = fopen ( fi leName , "w" ) ;f p r i n t f ( f i l e , " #GNUPLOT f i l e \ n " ) ;f p r i n t f ( f i l e , HELP_TO_PLOT) ;

for ( i = 0 ; i < toBePloted . nTr iang les ; i ++)

n1 = nodeGet ( toBePloted . t r i a n g l e s [ i ] [ 0 ] , toBePloted ) ;n2 = nodeGet ( toBePloted . t r i a n g l e s [ i ] [ 1 ] , toBePloted ) ;n3 = nodeGet ( toBePloted . t r i a n g l e s [ i ] [ 2 ] , toBePloted ) ;

70

f p r i n t f ( f i l e , " %12g %12g %12g \ n \%12g %12g %12g \ n \ n \%12g %12g %12g \ n \%12g %12g %12g \ n \ n \ n " ,n1 . x , n1 . y , n1 . z ,n2 . x , n2 . y , n2 . z ,n1 . x , n1 . y , n1 . z ,n3 . x , n3 . y , n3 . z ) ;

f c l o s e ( f i l e ) ;

voidp lo tSo lu t i onOver ( tSur face esfera , char ∗ f i leName , char

inputRealOrCalcu la ted )

unsigned long i , j ;double n1Relocate , n2Relocate , n3Relocate ;double max, min ;vec to rAr ray_o f_db l so lu t ionToBePr in ted ;node n1 , n2 , n3 ;FILE ∗ f i l e ;f i l e = fopen ( fi leName , "w" ) ;f p r i n t f ( f i l e , HELP_TO_PLOT) ;

i f ( inputRealOrCalcu la ted == ’ r ’ )so lu t ionToBePr in ted = esfera . so lReal ;

else i f ( inputRealOrCalcu la ted == ’ c ’ )so lu t ionToBePr in ted = esfera . so lCalc ;

else i f ( inputRealOrCalcu la ted == ’ e ’ )so lu t ionToBePr in ted = esfera . e r r o r ;

for (max=0 , i =0; i <es fera . nNodes−1; i ++) max= (max>so lu t ionToBePr in ted [ i ] ) ? max : so lu t ionToBePr in ted [

i ] ;

71

for ( min=0 , i =0; i <es fera . nNodes−1; i ++) min= ( min<so lu t ionToBePr in ted [ i ] ) ? min : so lu t ionToBePr in ted [

i ] ;

for ( i = 0 ; i < es fera . nTr iang les ; i ++)

n1 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 0 ] , es fera ) ;n2 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 1 ] , es fera ) ;n3 = nodeGet ( es fera . t r i a n g l e s [ i ] [ 2 ] , es fera ) ;

i f ( so lu t ionToBePr in ted [ n1 . index ] >=0)n1Relocate =1; / / +so lu t ionToBePr in ted [ n1 . index ] / max

/ 2 0 ;else

n1Relocate =1;i f ( so lu t ionToBePr in ted [ n2 . index ] >=0)

n2Relocate =1; / / +so lu t ionToBePr in ted [ n2 . index ] / max/ 2 0 ;

elsen2Relocate =1;

i f ( so lu t ionToBePr in ted [ n3 . index ] >=0)n3Relocate =1; / / +so lu t ionToBePr in ted [ n3 . index ] / max

/ 2 0 ;else

n3Relocate =1;

f p r i n t f ( f i l e , " %12g %12g %12g %12g \ n \%12g %12g %12g %12g \ n \ n \%12g %12g %12g %12g \ n \%12g %12g %12g %12g \ n \ n \ n " ,n1 . x∗n1Relocate , n1 . y∗n1Relocate , n1 . z∗n1Relocate ,so lu t ionToBePr in ted [ n1 . index ] ,n2 . x∗n2Relocate , n2 . y∗n2Relocate , n2 . z∗n2Relocate ,so lu t ionToBePr in ted [ n2 . index ] ,

n1 . x∗n1Relocate , n1 . y∗n1Relocate , n1 . z∗n1Relocate ,so lu t ionToBePr in ted [ n1 . index ] ,

72

n3 . x∗n3Relocate , n3 . y∗n3Relocate , n3 . z∗n3Relocate ,so lu t ionToBePr in ted [ n3 . index ]) ;

f c l o s e ( f i l e ) ;

Lista A.7: Archivo de proyecto galerkin/Makefile

CC=gcc# − l g s l − l g s l c b l a sLIB=−lm −l sparse#−mtune=pentium2LDFLAGS=CFLAGS=−gHEADERS=tSur face . hMODULES=tSur face . c poisson . cOFILES=tSur face . o so lveGa le rk in . o i n t e g r a t e . o n l so l veGa le rk i n . oLIBRARY= l i b s g . aLOADLIBES=$ (LIBRARY) −lm −l sparse

a l l : $ (LIBRARY) prueba01 prueba02 nlprueba01

prueba01 : prueba01 . o $ (LIBRARY)

prueba02 : prueba02 . o $ (LIBRARY)

nlprueba01 : nlprueba01 . o $ (LIBRARY)

$ (LIBRARY) : $ (OFILES )ar r $@ $?r a n l i b $@

so lveGa le rk in . o : n l so l veGa le rk i n . c so lveGa le rk in . c so lveGa le rk in. h tSur face . h

tSur face . o : tSur face . c tSur face . h

73

i n t e g r a t e . o : i n t e g r a t e . c tSur face . h so lveGa le rk in . h

clean :rm − f prueba0? ∗ . bak $ (LIBRARY) $ (OUTPUT) ∗ . o core ∗ . swp

∗~

c leanP lo ts :rm − f ∗ . dat

74

Apéndice B

Implementación del método de Fourierpara una base en el espacio depolinomios armónicos homogéneosrestringidos a la esfera.

Escribimos una librería de funciones para resolver una ecuación diferencial de Poissonsobre la esfera por el método de Fourier de modo que el usuario tenga que tomar encuenta unas pocas consideraciones que en seguida se enumeran:

1. escribir la definición de la solución real f y calculada u. u debe retornar cero si nose conoce.

2. usar la función gen() para discretizar la esfera usando una malla de triángulos.

3. usar la función shSolve() para calcular la solución del problema. Aquí podemosindicar hasta que grado queremos la base.

4. escribimos la solución con la función shPrintSolution().

5. podemos dibujar en escala de colores la solución calculada, real y el error sobre laesfera con la función plotSolutionOver()

6. finalmente liberamos la memoria que ocupa la estructura tSurfaceDestroy().

B.1. Referencia

SOLUCIÓN POR EL MÉTODO DE FOURIER SOBRE LA SUPERFICIE TRIANGU-LADA

75

shSolve() Solución del problema de Direchlet \Delta_S_2U=F por el método Fou-rier usando polinomios armónicos homogéneos. Debe definir las funciones F yU antes.

Definición: tSurface shSolve (int grado, tSurface,double);

Regresa: La solución por el método de Fourier para una base de polinomios ar-mónicos esféricos en forma de un arreglo contenido dentro de la estructura dedatos de una malla de triángulos (tSurface).

Argumentos:

grado <input> (int) grado máximo a usar para aproximar la solución.

error <input> (double) Criterio de parada usando la igualdad de Parseval.

ESCRIBIENDO LA SOLUCIÓN

Definición: void shPrintSolution (tSurface esfera, int lastNodeToBePrinted);

Regresa: void.

Argumentos:

esfera <input> (tSurface) Una estructura de datos de una malla de triángulos(tSurface).

lastNodeToBePrinted <input> (int) Imprime la solución hasta alcanzar el nodolastNodeToBePrinted.

B.2. Listados de la librería

Lista B.1: Archivo de cabecera que contiene la declaración de todos los datos y todaslas funciones usadas en la librería, así como, parámetros que sirven para configurar elcomportamiento de la librería durante el tiempo de compilación. El código está en elarchivo fuente fourier/shDefs.h

/ / Imprime in formac ión re levan te para l a func ióncreateBaseVectors ( )

/ / d e f i n i d a en shBase . c#define PRINT_MULTI_INDICE 1#define SHOW_LEFT_MATRIX 0#define SHOW_USED_INDICES 0#define WRITE_BASE_MATRIX 1/ / Imprime l a mat r i z base . Cada columna es un vec to r base de l

espacio

76

/ / de pol inomios armónicos homogéneos r e s t r i n g i d o s a l a/ / es fera .#define WRITE_BASE_MATRIX_ORTHO_NORMALIZED 1/ / Ac t i va código usado para probar las func iones d e f i n i d a s en e l

a rch ivo/ / shAr i tme t i ca . c#define COEF_TEST 0/ / CONSTANTES DE LA MAQUINA#define DBL_EPSILON 8.9e−15/∗

SPHERE HARMONIC DATA SHEET STRUCTURE

>>> St ruc tu re f i e l d s :

k ( i n t )The polynomia ls degree .

dimPk ( i n t )I s the dimension o f the polynomia ls o f degree k .

dimPkmenos2 ( i n t )I s the dimenssion o f the po l i nom ia l wh i t order k−2.

matr izBase ( double ∗ )A mat r i x ar ray whose columns vec to r are a base f o r the

po l i nom ia lexpansion .

∗ /typedef struct

i n t k ;i n t dimPk ; /∗ i s the dimension o f the

polynomia ls o fdegree k .

∗ /i n t dimPkmenos2 ; /∗ i s the dimension o f the

polynomia ls o f

77

degree k−2.∗ /

double ∗matrizBase ;struct

i n t m, n ; matr izBaseSize ;

char ∗shFurierSystem ; /∗ i t ’ s a p o i n t e r to a sparsemat r i x

o f the second member o f thenext equat ion

f ∗ Y_ i = (sum_ i =0 ( A_ k ,j ∗ Y_ j ) )

∗ Y_ i where : 1<= i <=dimPk and 1<= j

<=dimPk∗ /

char ∗M_Pk; /∗ i t ’ s untended t h a t t h i s shouldbe

a p o i n t e r a sparse mat r i x t h a tconta ins the t rans form mat r i x

f o rthe Laplace opera tor∗∗∗ NOT USED YET ∗ /

shDataSheet ;

/∗ the next s t r u c t u r e i s used to handle a mu l t i−index vec to r ∗ /typedef struct

i n t a , b , c ; m u l t i _ i n d i c e ;

m u l t i _ i n d i c e i n d i c ( int , i n t ) ;i n t elp ( int , int , i n t ) ;

typedef shDataSheet ∗shDataSheetPtr ;

78

i n t j ( int , i n t ) ;void shBui ldBaseVectors ( shDataSheet ∗ ) ; / / . . . d e f i n i d a en

shBase . cvoid shOrthoNormalizeBaseVectors ( shDataSheet ∗ ) ; / / . . .

d e f i n i d a en shBase . c

double ∗baseMatrixGetElement ( shDataSheet ∗ , int , i n t ) ;

double f a c t o r i a l ( int , i n t ) ;double coef ( int , int , i n t ) ;

double i n t e r n a l P r o d u c t ( shDataSheet ∗ , int , i n t ) ;void shBui ldFour ierSystem ( shDataSheet ∗ ) ; / / . . . d e f i n i d a en

shExp . cshDataSheet ∗shCreateDataSheet ( i n t ) ;

/∗∗ MEMORY ALLOCATION∗ /

#include < s t d l i b . h>

#define ALLOC( type , number ) ( ( type ∗ ) mal loc ( sizeof ( type ) ∗ ( number )) )

#define REALLOC( p t r , type , number ) \p t r = ( type ∗ ) r e a l l o c ( ( char ∗ ) p t r , ( unsigned ) ( sizeof ( type )

∗ ( number ) ) )

#define FREE( p t r ) i f ( ( p t r ) != NULL) f ree ( ( char ∗ ) ( p t r ) ) ; ( p t r )= NULL ;

/∗ Cal loc t h a t p rope r l y handles a l l o c a t i n g a c leared vec to r . ∗ /#define CALLOC( p t r , type , number ) \ i n t i ; p t r = ALLOC( type , number ) ; \

i f ( p t r != ( type ∗ )NULL) \for ( i =(number )−1; i >=0; i −−) p t r [ i ] = ( type ) 0 ; \

79

typedef double mat r i x ;

#define CREATE_MATRIX( type , s ize ) ( ( type ∗ ) c a l l o c ( ( s i ze ) , ( sizeof( type ) ) ) )

/∗ADD_REAL_ELEMENTAdds a r e a l number i n the d i r e c t i o n po in ted by p t r .p t r i s the d i r e c t i o n o f an element o f a mat r i x Row and Col .

∗ /#define ADD_REAL_ELEMENT( p t r , r e a l ) ∗ ( p t r ) +=(double ) ( r e a l )#define EQUAL_REAL_ELEMENT( p t r , r e a l ) ∗ ( p t r ) =( double ) ( r e a l )#include " . . / g a l e r k i n / tSur face . h "#include " . . / g a l e r k i n / i n t e g r a t e . h "double f ( vec to r ) ;double u ( node ) ;tSur face shBu i l dSo lu t i on ( shDataSheet ∗ , tSur face , double ∗ ) ;tSur face shSolve ( int , tSur face , double ) ;void s h P r i n t S o l u t i o n ( tSur face , i n t ) ;void shDataSheetDestroy ( shDataSheet ∗ ) ;

Lista B.2: Módulo en C que define las funciones para construir la base. Este código estáen el archivo fuente fourier/shBase.c

/∗∗ BUILD AND ORTHONORMALICE BASE VECTORS MODULE∗∗ >>> Funct ions access ib le out o f t h i s f i l e∗ shBui ldBaseVectors∗ shOrthoNormalizeBaseVectors∗∗ >>> Funct ions access ib le i n t h i s f i l e∗ baseMatrixGetElement∗ /

#include < s t d i o . h>#include < s t d l i b . h>#include <spConfig . h>#include <spmatr ix . h>#include <math . h>

80

#include " shDefs . h "

/∗ELEMENT ADDITION TO BASE MATRIX BY INDEX

baseMatrixGetElement ( )Finds element [Row, Col ] and re tu rns a p o i n t e r to i t . I f the

elementi s not found then an e r r o r message i s d isp layed and a n u l l

p o i n t e ri s re turned . The mat r i x as i t s s ize i s s tored i n t o a data

sheets t r u c t u r e which i s passed as an argument .

>>> D e f i n i t i o n :double ∗ baseMatrixGetElement ( shDataSheet ∗ ds , i n t i , i n t j )

>>> Returns :Returns a p o i n t e r to the element . This p o i n t e r i s then used

to d i r e c t l yaccess the element dur ing successive b u i l d s .

>>> Arguments :

ds < input > ( shDataSheet ∗ )p o i n t e r to the sphere harmonic data sheet where the base

mat r i x i ss tored .

i < input > ( i n t )Row index f o r element . Must be i n the range of [ 1 . .m] m=dimPk

. I f i >mthen an e r r o r message i s d isp layed and n u l l p o i n t e r i s

re turned .

j < input > ( i n t )

81

Col index f o r element . Must be i n the range of [ 1 . . n ] n=2∗k+1. I f j >n

then an e r r o r message i s d isp layed and n u l l p o i n t e r i sre turned .

∗ /double ∗baseMatrixGetElement ( shDataSheet ∗ ds , i n t i , i n t j )

i f ( i > ds−>matr izBaseSize .m)

p r i n t f( "ERROR: i n f i l e \ " shBase . c \ " c a l l to \ "

baseMatrixGetElement ( ds , i , j ) \ " wh i t i >m ( %d> %d ) \ n " ,i , ds−>matr izBaseSize .m) ;

return ( double ∗ ) NULL ;

i f ( j > ds−>matr izBaseSize . n )

p r i n t f( "ERROR: i n f i l e \ " shBase . c \ " c a l l to \ "

baseMatrixGetElement ( ds , i , j ) \ " wh i t j >n ( %d> %d ) \ n " ,j , ds−>matr izBaseSize . n ) ;

return ( double ∗ ) NULL ;

return (&ds−>matrizBase [ 0 ] + ( j − 1) ∗ ds−>matr izBaseSize .m + (i − 1) ) ;

i n tj ( i n t k , i n t i )

m u l t i _ i n d i c e v ;v = i n d i c ( k , i ) ;return e lp ( v . a + 2 , v . b , v . c ) ;

82

voidshPr intBaseVectors ( shDataSheet ∗ ds )

i n t step = 5; / / . . . nÃomero de columnas amostrar por l Ã n e a

i n t j = 1 , i , nColumna ;i f ( step >= ds−>matr izBaseSize . n )

step = ds−>matr izBaseSize . n − 1;do

p r i n t f ( " \ n " ) ;for ( nColumna = j ; nColumna <= j + step ; nColumna++)

p r i n t f ( "VEC # %04d " , nColumna ) ;p r i n t f ( " \ n " ) ;for ( i = 1 ; i <= ds−>matr izBaseSize .m; i ++)

for ( nColumna = j ; nColumna <= j + step ; nColumna++)

p r i n t f ( " %10.3g " , ds−>matrizBase [ ( nColumna − 1)

∗ds−>

matr izBaseSize.m +

( i − 1) ] ) ;

p r i n t f ( " \ n " ) ;

j += ( step + 1) ;i f ( j + step > ds−>matr izBaseSize . n )

step = ds−>matr izBaseSize . n − j ;

while ( j <= ds−>matr izBaseSize . n ) ;

voidshBui ldBaseVectors ( shDataSheet ∗ ds )

char ∗M;

83

i n t er ro r , i , nColumna ;i n t h ;m u l t i _ i n d i c e v ;spREAL rhs [ ds−>dimPkmenos2 ] ;spREAL s o l u t i o n [ ds−>dimPkmenos2 ] ;i n t J [ ds−>dimPk + 1 ] ;ds−>matr izBaseSize .m = ds−>dimPk ;ds−>matr izBaseSize . n = 2 ∗ ds−>k + 1;ds−>matrizBase =

CREATE_MATRIX ( double , ds−>matr izBaseSize .m ∗ ds−>matr izBaseSize . n ) ;

# i f PRINT_MULTI_INDICEi n t p ;p r i n t f ( " \ nMULTI−INDDICES \ n " ) ;p r i n t f ( " \ t %3s ( %3s, %3s, %3s ) \ n " , "ELP" , "A" , "B" , "C" ) ;for ( p = 1 ; p <= ds−>dimPk ; p++)

v = i n d i c ( ds−>k , p ) ;p r i n t f ( " \ t %3d ( %3d, %3d, %3d ) \ n " , e lp ( v . a , v . b , v . c ) , v . a ,

v . b , v . c ) ;

#endifi f ( ds−>k == 1)

for ( nColumna = 1; nColumna <= 3; nColumna++)

for ( i = 1 ; i <= 3 ; i ++)

i f ( nColumna == i )ADD_REAL_ELEMENT ( baseMatrixGetElement ( ds , i ,

nColumna ) , 1) ;else

ADD_REAL_ELEMENT ( baseMatrixGetElement ( ds , i ,nColumna ) , 0) ;

else

84

M = spCreate ( ds−>dimPkmenos2 , 0 , &e r r o r ) ;spClear (M) ;for ( i = 1 ; i <= ds−>dimPkmenos2 ; i ++)

v = i n d i c ( ds−>k − 2 , i ) ;spADD_REAL_ELEMENT ( spGetElement (M, i , i ) , ( v . a + 2) ∗

( v . a + 1) ) ;i f ( v . b >= 2)

spADD_REAL_ELEMENT ( spGetElement(M, e lp ( v . a + 2 , v . b − 2 , v . c ) ,

i ) ,( v . b ) ∗ ( v . b − 1) ) ;

i f ( v . c >= 2)spADD_REAL_ELEMENT ( spGetElement

(M, e lp ( v . a + 2 , v . b , v . c − 2) ,i ) ,

( v . c ) ∗ ( v . c − 1) ) ;

# i f SHOW_LEFT_MATRIXspPr i n t (M, 0 , 1 , 1) ;

#endife r r o r = spFactor (M) ;/ / . . . Creando un Ã n d i c e para descar ta r l as pos ic iones ya

usadas .for ( i = 1 ; i <= ds−>dimPk ; i ++)

J [ i ] = i ;for ( i = 1 ; i <= ds−>dimPkmenos2 ; i ++)

J [ j ( ds−>k − 2 , i ) ] = 0 ;# i f SHOW_USED_INDICES

p r i n t f ( " \nUSED INDEXES %s: %d \ n " , __FILE__ , __LINE__ ) ;for ( i = 1 ; i <= ds−>dimPk ; i ++)

p r i n t f ( " \ t J[ %d]= %d \ n " , i , J [ i ] ) ;#endif

nColumna = 1;for ( h = 1 ; h <= ds−>dimPk ; h++)

i f ( J [ h ] == h )

85

for ( i = 0 ; i < ds−>dimPkmenos2 ; i ++)

rhs [ i ] = 0 ; / / . . . i n i c i a n d o a cerov = i n d i c ( ds−>k , h ) ;i f ( v . a >= 2)

rhs [ e lp ( v . a − 2 , v . b , v . c ) − 1] = −v . a ∗ ( v . a −1 . ) ;

i f ( v . b >= 2)rhs [ e lp ( v . a , v . b − 2 , v . c ) − 1] = −v . b ∗ ( v . b −

1 . ) ;i f ( v . c >= 2)

rhs [ e lp ( v . a , v . b , v . c − 2) − 1] = −v . c ∗ ( v . c −1 . ) ;

spSolve (M, rhs , s o l u t i o n ) ;

for ( i = 1 ; i <= ds−>dimPkmenos2 ; i ++)ADD_REAL_ELEMENT ( baseMatrixGetElement ( ds , i ,

nColumna ) ,s o l u t i o n [ i − 1 ] ) ;

ADD_REAL_ELEMENT ( baseMatrixGetElement ( ds , h ,nColumna ) , 1) ;

nColumna++;

spDestroy (M) ;

# i f WRITE_BASE_MATRIX

p r i n t f ( " \ n \nVECTORES BASE ∗∗∗ " ) ;shPr intBaseVectors ( ds ) ;

#endif

doublei n t e r n a l P r o d u c t ( shDataSheet ∗ ds , i n t vect1 , i n t vect2 )

double i n t e r n a l P r o d u c t = 0 ;i n t k , l ;

86

m u l t i _ i n d i c e mi1 , mi2 ;for ( k = 1 ; k <= ds−>matr izBaseSize .m; k++)

mi1 = i n d i c ( ds−>k , k ) ;for ( l = 1 ; l <= ds−>matr izBaseSize .m; l ++)

mi2 = i n d i c ( ds−>k , l ) ;i n t e r n a l P r o d u c t +=∗ ( baseMatrixGetElement ( ds , k , vect1 ) ) ∗( ∗ ( baseMatrixGetElement ( ds , l , vect2 ) ) ) ∗coef ( mi1 . a + mi2 . a , mi1 . b + mi2 . b , mi1 . c + mi2 . c ) ;

return i n t e r n a l P r o d u c t ;

voidshOrthoNormalizeBaseVectors ( shDataSheet ∗ ds )

i n t nColumna , i , j ;double f ac t , norm ;/∗ . . . La columna 1 permanece i n t a c t a ∗ /for ( nColumna = 2; nColumna <= ds−>matr izBaseSize . n ; nColumna

++)for ( j = 1 ; j < nColumna ; j ++)

f a c t =− i n t e r n a l P r o d u c t ( ds , nColumna , j ) / i n t e r n a l P r o d u c t (

ds , j , j ) ;for ( i = 1 ; i <= ds−>matr izBaseSize .m; i ++)

ADD_REAL_ELEMENT ( baseMatrixGetElement ( ds , i , nColumna) ,

∗ ( baseMatrixGetElement ( ds , i , j ) ) ∗f a c t ) ;

#define NORMALIZANDO# i f d e f NORMALIZANDO

/ / . . . NORMALIZANDO

87

for ( j = 1 ; j <= ds−>matr izBaseSize . n ; j ++)

norm = s q r t ( i n t e r n a l P r o d u c t ( ds , j , j ) ) ;for ( i = 1 ; i <= ds−>matr izBaseSize .m; i ++)

EQUAL_REAL_ELEMENT ( baseMatrixGetElement ( ds , i , j ) ,∗ ( baseMatrixGetElement ( ds , i , j ) ) /

norm ) ;

#endif

/ / . . . IMPRIMIENDO LA MATRIZ BASE# i f WRITE_BASE_MATRIX_ORTHO_NORMALIZED

p r i n t f ( " \ n \nVECTORES BASE ORTHO−NORMALIZADOS ∗∗∗ " ) ;shPr intBaseVectors ( ds ) ;

#endif

Lista B.3: Módulo de ordenamimento contine las funciones indic() y elp(). Este código estáen el archivo fuente fourier/shIndic.c

#include " shDefs . h "

m u l t i _ i n d i c ei n d i c ( i n t k , i n t p )

i n t a , b , c ;m u l t i _ i n d i c e v ;

c = 0 ;for ( c = 0 ;

( p − ( k + 2) ∗ ( k + 1) / 2 + ( k − c + 2) ∗ ( k − c + 1) / 2>

k − c + 1) ; c++) ;

b = p − ( k + 2) ∗ ( k + 1) / 2 + ( k − c + 2) ∗ ( k − c + 1) / 2 −1;

88

a = k − b − c ;

v . a = a ;v . b = b ;v . c = c ;

return v ;

i n telp ( i n t a , i n t b , i n t c )

i n t k = a + b + c ;return ( b + ( k + 2) ∗ ( k + 1) / 2 − ( k − c + 2) ∗ ( k − c + 1) /

2 + 1) ;

Lista B.4: Módulo de aritmética, contiene las funciones coef() y factorial(). Este códigoestá en el archivo fuente fourier/shAritmetica.c

#include <math . h>#include " shDefs . h "# i f COEF_TEST#include < s t d i o . h>#endif

doublef a c t o r i a l ( i n t n , i n t m)

double f = ( double ) n ;i f ( n <= ( double ) m + DBL_EPSILON)

return 1 . 0 ;f = f ∗ f a c t o r i a l ( n − 1 , m) ;return f ;

doublecoef ( i n t a , i n t b , i n t c )

89

double coef ;i n t q ;i f ( a % 2 != 0)

return 0;i f ( b % 2 != 0)

return 0;i f ( c % 2 != 0)

return 0;q = a + b + c ;coef = pow (1 , q / 2) ∗ f a c t o r i a l ( a , a / 2) ∗

f a c t o r i a l ( b , b / 2) ∗ f a c t o r i a l ( c , c / 2) /f a c t o r i a l ( q , q / 2) / ( q + 1 .0 ) ;

# i f COEF_TESTp r i n t f ( " \nCOEF( %d, %d, %d )= %g \ n " , a , b , c , coef ) ;

#endifreturn ( coef ) ;

Lista B.5: Módulo para el cálculo de los coeficientes de Fourier. Este código está en elarchivo fuente fourier/shExp.c

/∗∗ FOURIER EXPANSION MODULE∗∗ This f i l e con ta in rou t i nes to perform the Four ie r expansion∗∗ >>> Funct ions conta ined i n t h i s f i l e∗ i n t e r n a l P r o d u c t∗ shBui ldFour ierSystem∗ /

#include < s t d i o . h>#include <spConfig . h>#include <spmatr ix . h>#include " shDefs . h "

voidshBui ldFour ierSystem ( shDataSheet ∗ ds )

90

i n t i , j , e r r o r ;ds−>shFurierSystem = spCreate ( ds−>matr izBaseSize . n , 0 , &e r r o r )

;spClear ( ds−>shFurierSystem ) ;for ( i = 1 ; i <= ds−>matr izBaseSize . n ; i ++)

for ( j = 1 ; j <= ds−>matr izBaseSize . n ; j ++)

spADD_REAL_ELEMENT ( spGetElement ( ds−>shFurierSystem , i ,j ) ,

i n t e r n a l P r o d u c t ( ds , i , j ) ) ;

p r i n t f ( " ∗∗∗ SUM Aki ∗ Yi ∗ Yj ∗∗∗ \ n " ) ;spPr i n t ( ds−>shFurierSystem , 0 , 1 , 1) ;e r r o r = spFactor ( ds−>shFurierSystem ) ;

Lista B.6: Makefile. Este código está en el archivo fuente fourier/Makefile

CC=gcc# − l g s l − l g s l c b l a sCFLAGS=−g −O2LDFLAGS=HFILES=shDefs . hCFILES= sh Ind ic . c shBui ld . c shBase . c shExp . c shAr i tme t i ca . cOFILES= sh Ind ic . o shBui ld . o shBase . o shExp . o shAr i tme t i ca . oLIBRARY= l i b s h . aTSURFACELIB = . . / g a l e r k i n / l i b s g . aLOADLIBES= $ (TSURFACELIB) −lm −l sparse#$ (LIBRARY)

a l l : $ (LIBRARY) $ (TSURFACELIB) prueba01 prueba02

prueba01 : prueba01 . o $ (LIBRARY)prueba02 : prueba02 . o $ (LIBRARY)

$ (TSURFACELIB) :cd . . / g a l e r k i n && make ; cd . . / f o u r i e r

91

$ (LIBRARY) : $ (OFILES )ar r $@ $?r a n l i b $@

$ (OFILES ) : $ ( HFILES )

clean :rm − f . ∗ . swp ∗ . bak prueba0? ∗~ ∗ . o ∗ . a core t e s t

c leanP lo ts :rm − f ∗ . dat

92

Apéndice C

Descripción y autirizaciones de uso delos programas informaticos empleadosen está tesis concedido por el(los)titular(es) del derecho de autor

C.1. GNU GPL

La GNU GPL (General Public License o licencia pública general) es una licencia crea-da por la Free Software Foundation a mediados de los 80, y está orientada principalmentea proteger la libre distribución, modificación y uso de software. Su propósito es declararque el software cubierto por esta licencia es software libre y protegerlo de intentos deapropiación que restrinjan esas libertades a los usuarios.

Existen varias licencias "hermanas" de la GPL, como la licencia de documentación li-bre GNU (GFDL), la Open Audio License, para trabajos musicales, etcétera, y otras menosrestrictivas, como la MGPL, o la LGPL (Lesser General Public License o Library GeneralPublic License), que permiten el enlace dinámico de aplicaciones libres a aplicaciones nolibres.

C.2. Sparse

Sparse 1.3 es un paquete de subrutinas flexible escrito en C usado para resolver rá-pida y exactamente sistemas de ecuaciones lineales dispersas. El paquete es capaz detrabajar con matrices cuadradas reales o complejas. Además es capaz de resolver sis-temas lineales, también puede resolver rápidamente sistemas transpuestos, encontrardeterminantes, y estimar errores debido a mal condicionamiento en el sistema de ecua-

93

ciones e inestabilidad en los cálculos. Sparse proporciona un programa de prueba que escapaz de leer ecuaciones matriciales de un archivo, resolverlas, e imprimir informaciónútil sobre la ecuación y su solución.

A continuación reproducimos la nota de derechos de autor del paqueteCopyright (c) 1985,86,87,88by Kenneth S. Kundert and the University of California.

Permission to use, copy, modify, and distribute this software and its documentationfor any purpose and without fee is hereby granted, provided that the copyright noticesappear in all copies and supporting documentation and that the authors and the Universityof California are properly credited. The authors and the University of California make norepresentations as to the suitability of this software for any purpose. It is provided ‘as is’,without express or implied warranty.

C.3. Gnuplot

‘gnuplot’ es un programa de comandos interactivo para dibujar funciones y datos. ‘gnu-plot’ se empleo en este trabajo para realizar dibujos en tres dimensiones de esferas y defunciones definidas sobre la esfera. A continuación se reproduce la nota de derechos deautor.

Copyright (C) 1986 - 1993, 1998, 2004 Thomas Williams, Colin KelleyPermission to use, copy, and distribute this software and its documentation for any

purpose with or without fee is hereby granted, provided that the above copyright noticeappear in all copies and that both that copyright notice and this permission notice appearin supporting documentation.

Permission to modify the software is granted, but not the right to distribute the com-plete modified source code. Modifications are to be distributed as patches to the releasedversion. Permission to distribute binaries produced by compiling modified sources is gran-ted, provided you 1. distribute the corresponding source modifications from the releasedversion in the form of a patch file along with the binaries, 2. add special version identifica-tion to distinguish your version in addition to the base release version number, 3. provideyour name and address as the primary contact for the support of your modified version,and 4. retain our contact information in regard to use of the base software. Permission todistribute the released version of the source code along with corresponding source modi-fications in the form of a patch file is granted with same provisions 2 through 4 for binarydistributions.

This software is provided "as is" without express or implied warranty to the extentpermitted by applicable law.

94

AUTHORSOriginal Software: Thomas Williams, Colin Kelley.Gnuplot 2.0 additions: Russell Lang, Dave Kotz, John Campbell.Gnuplot 3.0 additions: Gershon Elber and many others.Gnuplot 4.0 additions: See list of contributors at head of this document.

C.4. LYX

LYX es un sistema para preparar documentos. Que permite crear artículos de matemá-ticas y científicos complejos manejando referencias cruzadas, bibliografías, índices, etc.A continuación se reproduce la nota de derechos de autor.

LYX is Copyright (C) 1995 by Matthias Ettrich, 1995-2001 LYX TeamThis program is free software; you can redistribute it and/or modify it under the terms

of the GNU General Public License as published by the Free Software Foundation; eitherversion 2 of the License, or (at your option) any later version.

LYX is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICU-LAR PURPOSE. See the GNU General Public License for more details. You should havereceived a copy of the GNU General Public License along with this program; if not, writeto the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

C.5. LATEX

LATEXes un procesador de textos formado por un gran conjunto de macros de TEX,escritas inicialmente por Leslie Lamport (LamportTEX) en 1984, con la intención de facilitarel uso del lenguaje de composición tipográfica creado por Donald Knuth. Es muy utilizadopara la composición de artículos académicos, tesis y libros técnicos, dado que la calidadtipográfica de los documentos realizados con LATEX es comparable a la de una editorialcientífica de primera línea. LATEX es software libre bajo licencia LPPL. La Licencia Públicadel Proyecto LATEX (LPPL) no contiene todos los términos de la distribución de LATEX. Segúnlo que se puede leer, se trata de una licencia de software libre, pero incompatible con laGPL porque tiene muchos requisitos que no están en ésta.

95

Bibliografía

[1] Thomée, V. (1997) Galerkin finite element method for parabolic problems. Springer-Verlag. New York, Inc.

[2] E. C. Zeeman. Stability of dynamical systems. Nonlinearity. Vol 1 (1998). 1,115-155

[3] Guíñez J. y Rueda A. D. (2002) Steady states for a Flokker-Planck equation on Sn.Acta Math. Hungar, 94 (3), 211-221.

[4] Folland, Gerald B. (1995) Introduction to Partial Differential Equations. Princeton Uni-versity Press, United States of America.

[5] Zienkiewicz ,O. C. y Taylor, R. L. (2000) The Finite Element Method, Volume 1: TheBasis. Butterworth Heinemann, Woburn.

[6] Warner, Frank W. (1987) Foundations of Differentiable Manifolds and Lie Groups.Springer-Verlag. United States of America.

[7] Berner S. C. and Scott L. R. (1994) The Mathematical teory of Finite Element Met-hods, Springer-Verlag, New York.

[8] Ciarlet P. G. (1978) The Finite Element Method for Elliptic Problems, North Holland,Amsterdam.

96

Trabajo de Grado - Página...

Documents

Transcript of Trabajo de Grado - Página...