SOLUCIÓN NUMÉRICA DE UNA ECUACIÓN DIFERENCIAL

USANDO EL MÉTODO DE GALERKIN Y WAVELETS B-SPLINES

CARDINALES∗
RONALD LEÓN NAVARRO∗∗

Resumen. En el presente artı́culo, resolvemos numéricamente una ecuación diferencial ordinaria

elı́ptica con condiciones de frontera tipo Dirichlet. El tratamiento numérico se realiza usando el clásico
método de Galerkin y un tipo especial de bases wavelets; estas son wavelets B-splines cardinales. Los
resultados de la experimentación numérica realizada muestran que, aún considerando coeficientes, y
función en el segundo miembro, de la ecuación diferencial, como funciones discontinuas con salto finito
grande, el sistema que se genera es estable si se aplica precondicionamiento y la solución numérica
es suficientemente exacta para bajos niveles de aproximación.

Key words. Método Galerkin, bases, wavelets, B-splines, Análisis Multiresolución, ED elı́ptica.

1. Introducción. Una forma de definir una wavelet ([11]) es mediante una fun-
ción ψ de promedio cero:
Z +∞
ψ(t)dt = 0,
−∞

la cual es dilatada con un parámetro s, y trasladada por u:

1 t−u
ψu,s (t) = √ ψ( ).
s s

Existen muchas otras formas de definir tal función, en forma explı́cita como antes o
implı́citamente, o a través de otra función (“función escala”); lo cierto es que, funda-
mentalmente estas funciones deben ser bien localizadas en “tiempo” y “frecuencia”,
y dependiendo de la base de funciones que se desee construir para algún tipo de re-
presentación funcional que se quiera realizar, se adicionan otras propiedades. En la
historia, se conoce que Alfred Haar, en 1910, construyó la primera wavelet, después
aparecieron el análisis de Fourier como una alternativa a la transformada de Fourier
estándar, y en las últimas dos décadas del siglo pasado y el primero de este siglo, po-
demos encontrar a numerosos matemáticos, fı́sicos e ingenieros que han contribuido a
enriquecer la teorı́a de wavelets.
Por otro lado, el método de Galerkin es uno de los mejores métodos para resolver
numéricamente ecuaciones diferenciales. Este método consiste en hallar una base fun-
cional para el espacio solución de la ecuación, luego proyectar la solución sobre la base
funcional, y minimizar el “residual” respecto a la base funcional.
Las traslaciones de una wavelet para todas las dilataciones forman una base ortonor-
mal incondicional de L2 (R) ([6]) y las traslaciones de una función escala para todas
las dilataciones forman una base ortonormal incondicional para Vj ⊂ L2 (R), lo cual es
una gran mejora sobre bases de polinomios estándar o bases trigonométricas para el
método de Galerkin. Ası́ es explicable que, el uso de esquemas numéricos basados en
wavelets hayan tenido un crecimiento significativo en las últimas dćadas. Además, las
wavelets tienen varias propiedades que son especialmente útiles para representar so-
luciones de ecuaciones diferenciales (EDs), tal como ortogonalidad, soporte compacto

∗ Este trabajo es financiado con los recursos ordinarios asignados a la UNT.

∗∗ Departamento de Matemáticas, Universidad Nacional de Trujillo ([email protected]).
1
2

y representación exacta de polinomios de un cierto grado. Su capacidad de represen-

tar datos en diferentes niveles de resolución permiten al cálculo eficiente y estable de
funciones con gradientes grandes o singularidades, que necesitarı́an una malla densa
o elementos del mas alto orden en un análisis de Elemento Finito ([14]).
En el presente trabajo aplicamos una versión del clásico método de Galerkin usando
bases wavelets B-splines cardinales centralizadas para hallar una solución aproxima-
da de una ecuación diferencial ordinaria elı́ptica con condiciones de frontera tipo
Dirichlet. Las wavelets B-splines tienen una estructura muy simple y propiedades
especiales que nos permitieron obtener resultados con un alto orden de exactitud y
algoritmos asociados con una complejidad baja.
2. El método de Galerkin y eficiencia computacional. Los modelos ma-
temáticos que aparecen en la ciencia e ingenierı́a usualmente toman la forma de una
ecuación operador
(2.1) Tf = g
donde T : X −→ Y es un operador lineal entre ciertos espacios funcionales X y Y . Sin
embargo, no siempre es posible resolver (2.1) exactamente. Ası́, se tiene que depender
de métodos de aproximación que sean computacionalmente los mas simples de mane-
jar; el método de Galerkin tiene estas caracteristicas y es muy usado en ecuaciones
diferenciales.

2.1. El método de Galerkin. Como un primer paso para resolver (2.1) busca-
mos una aproximación f de f que se encuentre en cierto subespacio de dimensión finita
de X. Ası́, suponemos que X y Y son espacios producto interno con bases ortonor-
males enumerables {ϕ1 , ϕ2 , ...} y {ψ1 , ψ2 , ...}, respectivamente ([3, 8]). Considerando
que la ecuación (2.1) tiene una única solución, se deberı́a buscar una aproximación de
la solución f en
XM := gen{ϕ1 , ..., ϕM }.
De la ecuación (2.1) obtenemos:
(2.2) hT f, ψi i = hg, ψi i ∀i ∈ Z.
Sin embargo, necesitamos buscar un f ∈ XM tal que
(2.3) hT f , ψi i = hg, ψi i ∀i ∈ {1, ..., M }.
Escribiendo
M
X
(2.4) f := αj ϕj ,
j=1

obtenemos
M
X
(2.5) Tf = αj T ϕj
j=1

tal que (2.3) es igual a

M
X
(2.6) hT ϕj , ψi iαj = hg, ψi i, i ∈ {1, ..., M }
j=1
 3

Luego, llegamos al problema de resolver una ecuación matricial

(2.7) Ax = b

donde A := [aij ] es la matriz de M × M con ij-ésima entrada aij = hT ϕj , ψi i y

b = [bi ] es el M -vector con i-ésima entrada bi = hg, ψi i. Una vez que la ecuación (2.7)
es resuelta, f es tomada como en (2.4).
El procedimiento anterior de obtener una solución aproximada para (2.1) es llamado
el método de Galerkin.
Para problemas que se representan en ecuaciones diferenciales, usualmente se tiene
que resolver una ecuación de la forma

(2.8) Lf = g

donde L : X0 ⊆ X −→ X es un operador no acotado definido densamente. Usualmente

X0 es un espacio de Hilbert con respecto a otro producto interno h., .i0 el cual induce
una norma fuerte k . k0 ; es decir, existe c0 > 0 tal que

k f k≤ c0 k f k0 ∀f ∈ X0 .

También se tienen condiciones de la forma

|hLu, vi| ≤ c1 k u k0 k v k0 ∀u, v ∈ X0 ,

hLu, ui ≥ c2 k u k2 ∀u ∈ X0 ,

para algunas constantes positivas c1 , c2 . Observar que

a(u, v) := hLu, vi, u, v ∈ X0 ,

define una forma bilineal acotada sobre X0 × X0 , tal que las condiciones anteriores
asegurarán, usando el teorema de representación de Riesz, que para cada g ∈ X, existe
una única u ∈ X0 tal que

hLu, vi = hg, vi ∀v ∈ X0 .

La última ecuación es la formulación débil de (2.7).

ejemplo 1. Consideremos el problema de valor de frontera de Dirichlet:
 
d du
(Lu)(t) := − p(t) + q(t)u(t) = g(t), α ≤ t ≤ β,
dt dt

u(α) = 0 = u(β).

En este caso, X = L2 (α, β), y

X0 = H01 := {f ∈ L2 [α, β]/f es absolutamente continua, f 0 ∈ L2 , f (α) = f (β)}

con producto interno dado por

hu, vi0 := hu, vi + hu0 , v 0 i, u, v ∈ H01 ,

4
Z β
a(u, v) := hLu, vi = [p(t)u0 (t)v 0 (t) + q(t)u(t)v(t)]dt.
α

Suponiendo que p(.) y q(.) satisfacen

0 < a0 ≤ p(t) ≤ a1 , 0 < q(t) ≤ b0 ∀t ∈ [α, β],

puede ser probado ([10]) que

c20 a0 k u k2H 1 ≤ a(u, u), |a(u, v)| ≤ máx{a0 , b0 } k u kH01 k v kH01 ,

0

donde c0 > 0 es tal que k u k0 ≤ c0 k u k. Entonces, el método de Galerkin para el

problema (2.1) corresponde a la ecuación matricial (2.6) con A := [aij ] y b = [bi ],
donde

aij = a(ϕj , ψi ), bi = hg, ψi iL2

Además, también

(2.9) k f − f k≤ c ı́nf k f − ϕ k,
ϕ∈XM

donde c > 0 es una constante independiente de f y M . En consecuencia, la apro-

ximación Galerkin f es la mejor entre todas las aproximaciones de XM , hasta una
constante positiva.
2.2. Algunas propiedades deseables. Mientras investigamos el problema de
resolver la ecuación matricial (2.6), tomando aspectos computacionales en cuenta, no
es suficiente que pueda ser resuelto únicamente. Nos gustarı́a resolver la ecuación (2.6)
eficientemente lo que resultarı́a en pequeños errores computacionales, e incurrirı́a en
menor costo computacional. Tal deseable situación aparece si:
A tiene un número condición pequeño, y
A es dispersa.
El número condición de una matriz es medido en términos de la norma de la matriz.
Debemos recordar que, si k . k es una norma natural definida sobre Rn , entonces la
correspondiente norma natural de la matriz A, de n × n, es
 
k Ax k n
k A k:= sup : x ∈ R , x 6= 0 .
kxk

Una importante propiedad de esta norma es que

k Ax k≤k A kk x k ∀x ∈ Rn ,

con k I k= 1, donde I es la matriz identidad.

El número condición de una matriz invertible A, es dado por

KA :=k A kk A−1 k .

A diferencia del número condición, no hay una definición precisa para una matriz
dispersa. Sin embargo, significa que “la mayorı́a de sus entradas” son ceros. Es claro
que los cálculos realizados con una matriz “más dispersa” deberı́an ser mucho menores
que aquellos realizados con una matriz “menos dispersa”.
Veamos como un número condición alto influencia en errores de cálculo.
 5

Para vectores b, b̃ ∈ Rn , consideremos x y x̃ tales que Ax = b y Ax̃ = b̃. Entonces,

tenemos

x − x̃ = A−1 (b − b̃)

tal que

k x − x̃ k≤k A−1 kk b − b̃ k .

Además, tenemos que k b k=k Ax k≤k A kk x k, tal que

k x − x̃ k k b − b̃ k k b − b̃ k
(2.10) ≤k A kk A−1 k = κA .
kxk kbk kbk

En efecto, existen x, x̃, b, b̃ tales que la igualdad ocurre en (2.10). Para observar esto,
sean x y u vectores diferentes de cero tales que

k Ax k=k A kk x k, k A−1 u k=k A−1 kk u k

y sean

b = Ax, b = ˜b + u, x̃ = x + A−1 u.

Entonces, tenemos:

k x − x̃ k k A−1 u k k A−1 kk u k k A kk A−1 kk u k k b − b̃ k

= = = =k A kk A−1 k .
kxk kxk kxk k Ax k kbk
Por lo tanto, aún si el error relativo en b es pequeño, un número condición grande
puede conducir a un error relativo grande en la solución x.
3. Wavelets B-splines. Existen varias familias de wavelets que se pueden cla-
sificar según conformen bases ortonormales, semiortonormales o biortogonales. Ellas
pueden ser generadas desde un Análisis Multiresolución; a continuación estudiaremos
solo una clase de ellas.
3.1. Análisis Multiresolución.
Definición 1. Un Análisis Multiresolución (AMR) de L2 (R) es una sucesión de
subespacios cerrados enlazados (Vj )j∈Z de L2 (R) tales que
1. ∩j∈Z Vj = 0, y ∪j∈Z Vj es denso en L2 (R).
2. f (x) ∈ Vj ⇐⇒ f (2x) ∈ Vj+1 .
3. Existe una función φ ∈ V0 , tal que {φ(. − k)}k∈Z forma una base de Riesz de
V0 .
La función φ, en la definición anterior, es llamada generador o función escala del
AMR. Definamos ahora las funciones

φjk (x) = 2j/2 φ(2j x − k)

y φk = φ0k . Entonces, {φjk }k∈Z forma una base de Riesz de Vj . El subespacio Wj es

definido como el complemento ortogonal de Vj con respecto a Vj+1 ,

(3.1) Vj+1 = Vj ⊕ Wj , Wj ⊥ Vj .

El espacio L2 (R), por tanto, tiene una descomposición L2 (R) = ⊕j∈Z Wj . Ha sido
probado ([5]) que existe una función ψ ∈ W0 tal que {ψjk }k∈Z forma una base de Wj
6

y por lo tanto {ψjk }(j,k)∈Z2 es una base de L2 (R). La función ψ es llamada wavelet
semi-ortogonal dado que

hψjk , ψj 0 k0 i = 0, j 6= j 0 .

Llamaremos {ψjk }(j,k)∈Z2 una Base Wavelet de L2 (R).

Análogamente, ψ es llamada una wavelet ortonormal si

hψjk , ψj 0 k0 i = δjj 0 δkk0 ,

y la base {ψjk }(j,k)∈Z 2 generada por una wavelet ortonormal ψ es llamada una Base
Wavelet Ortonormal de L2 (R).
Observación 1. Algunas veces, el subespacio Wj definido por (3.1) no necesita
satisfacer la condición Wj ⊥ Vj . Entonces, la correspondiente wavelet ψ no es semi-
ortogonal. En este caso, una función ψ ? ∈ L2 (R) es llamada wavelet dual de ψ con
respecto a L2 (R) si

ψjk , ψj∗0 k0 = δjj 0 δkk0 .



El par de (ψ, ψ ∗ ) es llamado biortogonal.

Observación 2. Un análisis multiresolución y bases wavelet para otros espacios
diferentes a L2 (R) pueden ser definidos en forma similar a la definición anterior. Es
claro que la función escala φ satisface la función de refinamiento
X
φ (x) = pk φ (2x − k)
k∈Z

donde
P la sucesión (pk ) es llamada la máscara de φ y la serie de Laurent p (z) =
1 k
2 k∈Z pk z es llamada el sı́mbolo de φ. En la mayorı́a de aplicaciones, se necesita
que φ sea compactamente soportada o tenga un decaimiento exponencial.
3.2. Wavelets splines. Los splines son buenos ejemplos de funciones escala,
especialmente los B-splines cardinales que definimos a continuación.
Sea χ la función caracterı́stica del intervalo unitario [0, 1).
Definición 2. Sea m ∈ Z+ . La función Nm : R −→ R recursivamente definida
por
Z 1
Nm (x) := (Nm−1 ∗ N1 ) (x) = Nm−1 (x − t) dt
0

donde N1 (x) = χ[0,1) (x), es llamada B-spline cardinal de orden m.

En aplicaciones, es común usar los splines con orden par. En este trabajo, suponemos
también que el spline siempre es de orden par. Algunas veces omitiremos el subı́ndice
m de las notaciones como Nm si eso no induce a alguna confusión.
Varias de las propiedades de las B-splines se resumen en la siguiente proposición, cuya
demostración se puede ver en [5, 16].
Proposición 1. Los B-splines Cardinales tienen las siguientes propiedades:
a) Ellos son compactamente soportados con sopNm ⊂ [0, m] .
b) Los B-splines son no-negativos, i.e. Nm (x) ≥ 0 y Nm (x) > 0 para x ∈ (0, m) .
c) Ellos forman una partición de la unidad; es decir,
Z X
Nm (x) dx = 1, Nm (x − k) = 1.
R k∈Z
 7

0
d) Nm ∈ C m−2 (R) y Nm (x) = Nm−1 (x) − Nm−1 (x − 1) .  
1−m
Pm m
e) Los B-splines Cardinales son refinables con Nm (x) = 2 k=0 Nm (2x − k) .
k
m m



f) Nm 2 + x = Nm 2 − x .
x
g) Nm (x) = m−1 Nm−1 (x) + m−x
m−1 Nm−1 (x − 1) .
h) Se cumple:
m  
1 X
k m
Nm (x) = (−1) (x − k)m−1
+ ,
(m − 1)! k
k=0

donde xn+ := (x+ )n = χ[0,∞) (x)xn .

Figura 3.1. B-splines cardinales de orden 1, 2 y 4

La transformada de Fourier de N (x) es

m
1 − e−iω

N̂ (ω) = .
ω

Se verifica que N (x) satisface la ecuación de refinamiento

m  
X m
(3.2) N (x) = 21−m N (2x − k).
k
k=0

El B-spline N (x) genera un AMR de L2 (R), es decir, los espacios splines Vk :=

closL2 (gen{Nm;k,j : j ∈ Z}), k ∈ Z, constituyen un análisis multiresolución de L2
en el sentido que
8

(i) ... ⊂ V−1 ⊂ V0 ⊂ V1 ⊂ ...;

(ii) closL2 (∪k∈Z Vk ) = L2 ;
(iii) ∩k∈Z Vk = {0};
(iv) para cada k, {Nm;k,j : j ∈ Z} es una base incondicional de Vk ;
(v) Vk+1 = Vk ⊕ Wk , para toda k ∈ Z;
(vi) Wk ⊥ Wj , para toda k 6= j;
(vii) L2 = ⊕k∈Z Wk .
Observación 3. La importancia de las propiedades de aproximación (ii)-(iii) y
las propiedades (v)-(vi) es que cada función f en L2 puede ser aproximada tan cerca
como se desee por algún fk ∈ Vk , para un valor suficientemente grande de k, y que fk
tiene una (única) descomposición ortogonal

(3.3) fk = gk−1 + ... + gk−l + fk−l

donde gj ∈ Wj , j = k − l, ..., k − 1, y fk−l ∈ Vk−l , y l es un entero positivo arbitra-

riamente grande, tal que k fk−l k2 es tan pequeño como se desee. Llamaremos a (3.3)
una descomposición wavelet de fk .
3.3. Wavelets splines sobre intervalos finitos. Ahora construimos una base
semi-ortogonal del espacio L2 (I), donde I es un intervalo finito.
Sin pérdida de generalidad, suponemos que I = [0, M ], M ∈ N y que {Vj }j∈Z es un
AMR generado por el B-spline cardinal N de orden m. Además, sea VjI ⊂ L2 (R) el
subespacio que contiene las restricciones sobre I de las funciones en Vj . Luego {VjI }j≥0
es un AMR de L2 (I) y dimVj = m + 2j M − 1. Para construir la base de Vj , hacemos
I
Njk = Njk χI , y Sj = {2−j k : sop{Njk } ∩ I 6= ∅}. Ası́ se tiene que |Sj | = dimVj y
algunos resultados que mostramos a continuación.
I
Proposición 2. {Njk }2−j k∈Sj forma una base de VjI , y existen dos constantes
positivas C1 y C2 independientes de j tales que
X X X
C1 |cjk |2 ≤ k I
cjk Njk k2 ≤ C2 |cjk |2 .
2−j k∈Sj 2−j k∈Sj 2−j k∈Sj

En la resolución de ecuaciones diferenciales parciales (EDP), las soluciones algunas

veces deben satisfacer ciertas condiciones de frontera homogéneas. Ası́, necesitamos
que las bases wavelets estén en el espacio

H0s (I) = {f ∈ C0s−1 (I); f (s−1) ∈ L2 (I)}, 0 ≤ s ≤ m − 1

Convenimos que H00 (I) = L2 (I). Dado que los B-splines truncados no están en el es-
pacio H0s (I), s ≥ 1, tenemos que adaptar aquellos construidos sobre todo R a aquellos
definidos sobre I tal que pertenezcan al espacio H0s (I). Los B-splines con múltiples
nodos son las herramientas apropiadas para tal adaptación.
Sea [x1 , ..., xn ]f la diferencia dividida de f en x1 , ..., xn . Denotamos la sucesión de
nodos
m m
z }| { z }| {
0, 0, ..., 0, 1, 2, ..., M, M, ..., M
−1
por {ti }m+M
i=m+1 , y definimos

m m
(3.4) φIk (x) = [tk− m2 , tk+1 , ..., tk+ m2 ](. − x)m−1
+ , − +1≤k ≤M −1+
2 2
 9

El soporte de φIk es [tk− m2 , tk+ m2 ], y φIk (x) = Nk− m2 (x), m2 ≤ k ≤ M −

m
2. Por
conveniencia, el AMR de H0s (I) aún es denotado por {VjI }∞ j=0 .
M −1−s+ m
Proposición 3. Φ0 := {φIk }k=s+1− m2 es una base de V0I , 0 ≤ s ≤
2
mı́n{m − 1, [ m+M
2 ] − 1}. Cuando m
2 ≤ k ≤ M − m
2, φIk
es la traslación del B-
spline central; es decir, φIk = N (x − k + m2 ).
Para construir las bases de VjI , dilatamos las funciones escala en Φ0 . Definimos
φIjk (x) = 2j/2 φIk (2j x) para s+1− m m j m j
2 ≤ k ≤ 2 −1 y 2 M − 2 +1 ≤ k ≤ 2 M −1−s+ 2 ;
m

I j/2 j m m j m
definimos φjk (x) = 2 N (2 x − k − 2 ) para 2 ≤ k ≤ 2 M − 2 . Entonces,

2j M −1−s+ m
(3.5) Φj := {φIjk }k=s+1− m 2
2

es una base de VjI . Escribimos SjL = 2−j {s+1− m m C −j m j

2 , ..., 2 −1}, Sj = 2 { 2 , ..., 2 M −
m R −j j m j m L I −j L
2 }, Sj = 2 {2 M − 2 + 1, ..., 2 M − 1 −s + 2 }. Entonces, Φj := {φjk : 2 k ∈ Sj }
R I
contiene las funciones escala con frontera izquierda en Φj , mientras que Φj := {φjk :
2−j k ∈ SjR } contiene las funciones escala de frontera derecha en Φj . Las funciones
−j
escala centrales están en ΦC I C L
j := {φjk : 2 k ∈ Sj }. Cuando s = m − 1, Φj y Φj son
R
L C R
vacı́os. Haciendo Sj = Sj ∪ Sj ∪ Sj , obtenemos una aplicación inyectiva Sφ de Φj a
Sj tal que Sφ (φIjk ) = 2−j k ∈ Sj .

Bases wavelets semi-ortogonales del espacio H0s (I). Sea WjI el subespacio
wavelet de H0s (I) tal que WjI ⊕ VjI = Vj+1I
y WjI ⊥ VjI . Haciendo S j = (Sj+1 \ Sj ) ∩
[0, M ], obtenemos dimWjI = |S j |. El conjunto S j es simétrico con respecto a M 2 , es
decir, si d ∈ S j , entonces M − d ∈ S j . Sin pérdida de generalidad, solo discutiremos
la construcción de la base wavelet de W I . (Por conveniencia, cualquier subindice 0
será omitido). Se verifica que dimW I = min(M, (2M + m − 2s − 1)+ ). La base de
W I no es única. En aplicaciones, dos propiedades principales son usualmente necesa-
rias. Primera, los soportes de los elementos en bases wavelets de W I son requeridos
tan pequeños como sea posible. Segunda, las bases son esperadas que tengan cierta
simetrı́a.
Ahora descomponemos el conjunto S j en la forma

L C R
(3.6) Sj = Sj ∪ Sj ∪ Sj ,

C L R
donde S j = {m − 21 , m + 12 , ..., M − m + 12 }, S j = S j ∩ (0, m − 1], y S j = S j ∩ [M −
C
m + 1, M ). Observar que si M < 2m − 1, entonces S j = ∅. Puede ser verificado que
la matriz

A = (< φIl , φI1k >)l∈S, k ∈S1

2

es una matriz de rango completo. La ecuación lineal Ad = 0 tiene M soluciones

linealmente independientes d1 , ..., dM . Ahora escribimos S = {k1 , ..., kM }, donde
k1 < ... < kM , y hacemos ψkI l = (dl )T Φ1 . Entonces {ψdI }d∈S es una base de W I . Selec-
cionando cuidadosamente vectores d1 , ..., dM , podemos hacer ψdI (x) = ψM I
−d (M − x)
L I
√ 1 C
para d ∈ S y ψd (x) = 2w(2(x − d) + m − 2 ) para d ∈ S donde w es definido por
 −iω/2
m Q
−iω/2 b
w(ω)
b = 1−e 2 m (e )N (ω/2)
10
P2m−2
con m (z) := k=0 N2m (k + 1)z k y N es un B-spline de orden m; w ası́ definido
Q
es una wavelet semi-ortogonal con respecto a N .

ejemplo 2. Considerar bases wavelet del espacio H0m−1 (I). En este caso, las
bases de W I pueden ser construidas de la siguiente forma.
1. Si M < m − 1, entonces V I = {0} y V1I = W I . Una base de V1I es también
una base de W I .
2. Si m ≤ M ≤ 2m − 2, entonces, para cualquier d ∈ S con d ≤ M 2 , podemos
−m+1+k
hallar (qlk )M
l=k tal que la función

M −m+k−1
X
(3.7) ψdI = qlk φ1l , k = 2d
l=k

es ortogonal al espacio V1I y k ψdI k2 = 1. Definimos ψdI (x) = ψM

I
−d (M − x)
M I I
para d ∈ S con d > 2 . Entonces {ψd }d∈S es una base de W .
L
3. Si M ≥ 2m − 1, tenemos S = { m m 1 m 1 m+1 m+3
4 , 4 + 2 , ..., 2 − 2 } ∪ { 2 , 2 , ...,
2m−3
2 }.
Hagamos K = min{M, 3m − 3}. Es posible verificar que la matriz
K− m ,K+ m
2 −1
A = (< φk , φ1l >)k,l= 2m
2

es de rango completo. Por lo tanto, la ecuación Ac = 0 tiene m−1 soluciones

independientes cl , 1 ≤ l ≤ m − 1. Podemos elegir cl tal que sus últimas
(m − 1 − l) componentes sean ceros. Sea ψdI = Φ1 cl , donde
 m l−1 m
d= 4 + 2 , 1≤l ≤ 2
1 m
l − 2, 2 +1≤l ≤m−1

Entonces el soporte de ψdI es

K +l K −1+m
sopψdI = [0, ] ⊂ [0, ].
2 2
R √
Consideremos ψdI (x) = ψM
I
−d (M − x) para d ∈ S y sea ψdI (x) = 2w(2(x −
1 C
d) + m − 2 ), para d ∈ S . Ası́, el conjunto {ψdI }d∈S forma una base de W I .

Observación 4. Si I es un intervalo diádico, por ejemplo I = [2−J k, 2−J (k +

l)], para J fijo, la construcción de la base de funciones escala de VjI y la base de
wavelets de WjI es análogo. Definimos los conjuntos Sj y S j de la misma forma como
anteriormente lo hicimos. Consecuentemente, tenemos aún que dimWjI = |S j |. Una
evidente modificación se necesita realizar para la construcción de la base de wavelets
de WjI .
4. Solución Numérica de una ecuación diferencial mediante el método
Galerkin y wavelets B-splines cardinales. Es de conocido que en los métodos
de Galerkin el sistema discreto para la resolución numérica de problemas elı́pticos en
un dominio acotado es mal condicionado, si se usan métodos de elementos finitos o
diferencias finitas. Usualmente, el número condición del sistema discreto para un pro-
blema elı́ptico de segundo orden en dos dimensiones es del orden O(1/h2 ); utilizando
métodos de pre-condicionamiento, podemos reducir el número condición a O(1/h).
 11

Sin embargo, si las bases wavelets son utilizadas, el precondicionamiento alcanza un

número condición de orden O(1), por lo cual, los métodos wavelet conducen a estabi-
lizaciones numéricas en la resolución de tales EDs. Además, los algoritmos iterativos
son muy populares en la resolución numérica de las ecuaciones diferenciales, especial-
mente en las PDEs, y el método wavelet puede acelerar la convergencia de algoritmos
iterativos.

A continuación realizamos una aplicación del método de Galerkin en la resolu-

ción de una ecuación diferencial elı́ptica unidimensional con condiciones de frontera
tipo Dirichlet utilizando bases wavelet B-splines cardinales. Consideramos el siguiente
Problema de Dirichlet.

(4.1) −u00 (x) + cu(x) = f (x), x ∈ (0, 1), c ≥ 0, u(0) = u(1) = 0,

donde f ∈ L2 (0, 1).

En sección “Método de Galerkin” establecimos, como un primer paso del procedi-
miento, elegir el espacio que contenga la solución del problema planteado, al que
llamaremos espacio muestral (al espacio de llegada del operador que define la ecua-
ción diferencial se le llama espacio test), luego formulamos el problema variacional
asociado. Realizada esta formulación, debemos garantizar que existe una solución del
problema. Finalmente, utilizamos bases wavelet B-splines cardinales para construir
soluciones aproximadas del problema variacional y consecuentemente del problema
original.
4.1. Formulación Variacional. Sea el espacio de funciones Test:

C0∞ ((0, 1)) := {v ∈ C ∞ ((0, 1)) : v(0) = v(1) = 0}.

Multiplicando ambos lados de la ecuación diferencial con una función test φ ∈ C0∞ ((0, 1))
y aplicando integración por partes, considerando que φ(0) = φ(1) = 0, obtenemos:
Z 1 Z 1 Z 1
(4.2) f (x)φ(x)dx = u0 (x)φ0 (x)dx + c u(x)φ(x)dx
0 0 0

Observamos que esta última ecuación está bien definida si u y φ están en el espacio
de Sobolev H01 ((0, 1)), definido por:

H01 ((0, 1)) = {v ∈ H 1 ((0, 1)), v(0) = v(1) = 0}

donde H 1 ((0, 1)) = {v ∈ L2 ((0, 1)) : v 0 ∈ L2 ((0, 1))}. Utilizando como espacio mues-
tral y Test a X := H01 ((0, 1)), la formulación variacional (o débil) se enuncia como
sigue:
Hallar u ∈ X tal que

(4.3) a(u, v) = (f, v)0;(0,1) , v ∈ X,

donde a : X × X −→ R es la forma bilineal definida por

(4.4) a(u, v) := (u0 , v 0 )0;(0,1) + c(u, v)0;(0,1)

Z 1
(4.5) = [u0 (x)v 0 (x) + cu(x)v(x)]dx.
0
12

Definición 3. Sea u ∈ L2 ((0, 1)). Una función v ∈ L2 ((0, 1)) es llamada derivada
débil de u, es decir, v = u0 si
Z 1 Z 1
(4.6) v(x)φ(x)dx = − u(x)φ0 (x)dx
0 0

para todas las funciones test φ ∈ C0∞ ((0, 1)).

Definición 4. Sea m ∈ N.
a) El espacio de Sobolev de orden m es definido por
H m ((0, 1)) := {v ∈ L2 ((0, 1)) : v (k) ∈ L2 ((0, 1)), 1 ≤ k ≤ m},
donde las derivadas son entendidas en sentido débil. Una norma sobre H m ((0, 1))
es definida por
m
X
kukm;(0,1) := ( ku(k) k20;(0,1) )1/2 ,
k=0

la cual es inducida por el producto interno

m
X
hu, vim;(0,1) := hu(k) , v (k) i0;(0,1) .
k=0

b) El espacio de Sobolev con condiciones de frontera de Dirichlet homogéneas

generalizadas es definido por
H0m ((0, 1)) := closk.km;(0,1) (C0∞ ((0, 1))).
4.2. Existencia y unicidad. Para averiguar si existe una única solución del
problema planteado, debemos analizar las propiedades de la forma bilineal a(., .), an-
tes definida, como veremos a continuación.

Primero, usando la desigualdad de Hölder se demuestra que a(., .) es acotada; es

decir,
(4.7) a(u, v) ≤ Ckuk1;Ω kuk1;Ω .
Segundo, debemos demostrar que la forma bilineal a es coerciva; es decir,
(4.8) a(u, v) ≥ αkuk21;Ω
Las constantes C y α son también conocidas como constantes de continuidad y coer-
cividad, respectivamente. Una forma bilineal que es coerciva y acotada es llamada
elı́ptica. Para demostrar que se cumple la última desigualdad necesitamos el siguiente
resultado, conocido como la desigualdad de Poincaré - Friedrichs.
Proposición 4. Para cualquier intervalo acotado Ω := [a, b], existe una constante
CΩ = (b − a) > 0 tal que
(4.9) kvk0;Ω ≤ CΩ kv 0 k0;Ω
para toda v ∈ H01 (Ω). Demostración: Consideremos una función φ ∈ C0∞ (Ω). Por el
teorema fundamental del cálculo, tenemos
Z x Z x
φ(x) = φ(0) + φ0 (ξ)dξ = φ0 (ξ)dξ.
a a
 13

Utilizando la desigualdad de Cauchy-Schwarz se obtiene

Z b Z b Z x Z x
2
|φ(x)| dx ≤ [ 2
1 dξ][ |φ0 (ξ)|2 dξ]dx
a a a a
Z b
(4.10) ≤ (b − a)2 |φ0 (ξ)|2 dξ = (b − a)2 kφ0 k20;Ω .
a

Desde que, por definición, el espacio de Sobolev H01 (Ω) es la clausura de C0∞ (Ω), el
espacio C0∞ (Ω) es denso en H01 tal que la desigualdad anterior se extiende también
para φ ∈ H01 (Ω) lo cual prueba la proposición.
Con esta preparación, (3,8) es fácilmente probado. En efecto,

1 1 1 1
a(u, u) = ku0 k20;Ω ≥ 2 kuk20;Ω + ku0 k20;Ω ≥ mı́n{ , }kuk21;Ω ,
2CΩ 2 2 2CΩ2

i.e., se obtiene (4.8) con α := mı́n{ 21 , 2C1 2 }.

Ω
Con estas propiedades y aplicando el Teorema de Lax-Milgram aseguramos la exis-
tencia de una única solución del problema variacional formulado.
4.3. Discretización. Consideremos subespacios finito-dimensionales XΛ ⊂ X,
dimXΛ = |Λ| < ∞, donde Λ es el conjunto de grados de libertad, el cual es supuesto
a ser finito para poder ser capaces de tratar el correspondiente problema numérico.
Luego, la versión discreta de (4.3) se enuncia:
Hallar uΛ ∈ XΛ tal que

(4.11) a(uΛ , vΛ ) = hf, vΛ i0;(0,1) , vΛ ∈ XΛ .

Aunque este problema no es completamente accesible en una computadora, cuenta el

hecho de tener una base de XΛ para manipular, digamos

XΛ = genΘΛ , ΘΛ = {ϑλ : λ ∈ Λ},

y ΘΛ es una familia de elementos linealmente independientes de X. Entonces la solu-

ción de (4.11) puede ser escrita como
X
(4.12) uΛ = uλ ϑλ =: uTΛ ΘΛ , uΛ := (uλ )λ∈Λ ,
λ∈Λ

donde ahora el vector de coeficientes uΛ necesita ser determinado por

X
(4.13) a(uΛ , ϑµ ) = uλ a(ϑλ , ϑµ ) = (f, ϑµ )0;(0,1) , µ ∈ Λ
λ∈Λ

Esto puede ser formulado equivalentemente en la forma matriz-vector

(4.14) AΛ uΛ = fΛ

donde AΛ := (a(ϑλ , ϑµ ))λ,µ∈Λ es la matriz de rigidez y el vector fΛ := ((f, ϑµ )0;(0,1) )µ∈Λ

es parte derecha de la igualdad.
Observar en (4.11), que el espacio test y el espacio muestral coinciden, es decir, ambos
son XΛ ; esto es común en el método clásico de Galerkin.
14

4.4. Estabilidad. Antes de tratar con la solución numérica de (4.11), investi-

gamos sobre la estabilidad y propiedades de aproximación del método de Galerkin. A
continuación observamos un estimado.
Proposición 5. La solución de (3.11) es estable independiente de la elección del
subespacio XΛ , es decir,
1
kuΛ kX ≤ kf kX 0 .
α
Demostración: Elegir vΛ = uΛ en (3.3) y usar la elipticidad para obtener

αkuΛ k2X ≤ a(uΛ , uΛ ) = hf, uΛ i ≤ kf kX 0 kuΛ kX

lo que prueba la aserción.

4.5. Estimaciones de error. Ahora, sea u ∈ V la solución de (4.3) y uΛ ∈ VΛ
la solución de (4.11), es decir,

a(u, v) = hf, vi ∀v ∈ X,
a(uΛ , vΛ ) = hf, vΛ i ∀vΛ ∈ XΛ .

Dado que XΛ ⊂ X, podemos también satisfacer la primera ecuación para vΛ ∈ XΛ .

Ası́, restando miembro a miembro ambas ecuaciones obtenemos:

(4.15) a(u − uΛ , vΛ ) = a(u, vΛ ) − a(uΛ , vΛ ) = hf, uΛ i − hf, uΛ i = 0.

Esta ecuación es llamada Ortogonalidad de Galerkin.

A continuación enunciamos, sin demostración, el Lema de Céa, el cual muestra que el
orden de aproximación de un método de Galerkin depende de la potencia de aproxi-
mación de los espacios de prueba XΛ . Ası́, necesitamos construir apropiados espacios
XΛ de tal forma que se pueda alcanzar un orden óptimo de aproximación con un
mı́nimo número de grados de libertad N.
Teorema 1. Sea a(., .) elı́ptica. Entonces, se cumple
C
(4.16) ku − uΛ kX ≤ ı́nf ku − vΛ kX .
α vΛ ∈XΛ

El término

(4.17) EΛ (u) := ı́nf ku − vΛ kX

vΛ ∈XΛ

es conocido como el error de la mejor aproximación en XΛ .

A partir del Lema de Céa se derivan una serie de estimaciones del error de aproxima-
ción para el método de Galerkin que no los usamos en el presente trabajo.
5. Resultados. En esta parte, utilizamos un Análisis Multiresolución (AMR)
para la discretización Galerkin de problemas con valor de frontera. En la literatura,
esto es algunas veces llamado Método Wavelet-Galerkin.
5.1. Discretización multiescala y precondicionamiento. Es conocido hace
algún tiempo que las soluciones de muchos problemas de valor en la frontera exhiben
comportamientos que se detectan en escalas de longitud ligeramente diferentes. Se ha
demostrado que esquemas numéricos que usan esta información pueden ser altamente
eficientes [urban]. En consecuencia, nos interesamos en discretizaciones multiescala.
 15

Estas son usualmente realizadas considerando una familia completa de espacios mues-
trales y Test
X0 ⊂ X1 ⊂ . . . ⊂ Xj ⊂ Xj+1 ⊂ . . . ⊂ X.
El intercambio entre estos espacios para diferentes valores de j es la clave para la
eficiencia del método resultante.
En métodos clásicos tales como discretizaciones con diferencias finitas o elemento
finito, la jerarquia multiescala puede fácilmente ser formada mediante un parámetro de
tamaño de malla hj → 0 para j → ∞ (por ejemplo, hj = 2−j ) y Xj como un espacio de
elemento finito correspondiente a una malla de tamaño hj . Luego, métodos conjugado
gradiente (cg) tipo precondicionado BPX o métodos multigrid tienen un rendimiento
óptimo, lo que significa que el problema de Galerkin sobre Xj puede ser resuelto en
O(dimXj ) operaciones; es decir, complejidad lineal. En este trabajo usamos un AMR
generado por un B-spline cardinal centralizado S = {Vj }j∈N0 para la discretización,
esto es, Xj := Vj . El esquema resultante es llamado Método multiresolución - Galerkin,
también abreviado como MGM o método MG. Describimos el método con un poco
de mas detalle. Suponemos que los espacios Test y muestral son dados en términos de
bases, es decir,
Xj = Vj = gen{Φj }, Φj := {ϕj,k : k ∈ Ij }.
Ası́, la sucesión Ij representa conjuntos de grados de libertad, |Ij |, tal que |Ij | −→ ∞
para j −→ ∞, tı́picamente en una forma exponencial, es decir,
|Ij | ∼ 2j .
Por lo tanto, el número de operaciones para la multiplicación matriz-vector con
Aj = AVj crece (exponencialmente) cuando crece j. Usando un esquema iterati-
vo, esto tı́picamente significa que la cantidad de trabajo por paso es O(|Ij |). En orden
de tener un mejor esquema numérico posible en mano, es deseable que el número de
iteraciones del esquema sea acotado independiente de j, tal que aumente el mismo
número de iteraciones para reducir el error mediante un factor dado independiente
de cuan fina sea la discretización. Teniendo el estimado del error, por ejemplo para
el método CG, esto significa que el número condición deberı́a también ser acotado
independiente de j, esto es, cond2 (Aj ) = O(1), j → ∞.
Para una base de funciones escalas, esto no es alcanzable. Se necesitará un precondi-
cionador Dj a Aj tal que
Dj sea simétrico y definido positivo;
−1/2
la aplicación de Dj a un vector dado vj puede ser realizado en O(|Ij |)
operaciones;
−1/2 −1/2
cond2 (Dj Aj Dj ) = O(1) para j → ∞.
5.2. Una aplicación Galerkin - Wavelets B-splines. En esta parte presen-
tamos los resultados de experimentación numérica realizada. Hemos aplicado una dis-
cretización multiescala multiresolución a casos especiales; utilizamos funciones escala
y bases wavelets B-splines cardinales centralizaods para realizar las correspondientes
discretizaciones. Si bien la teorı́a desarrollada nos garantiza únicas soluciones para el
problema de Dirichlet homogéneo cuando c es una constante positiva, consideramos
el caso cuando c es una función constante por tramos con soporte, el intervalo [0, 1].
Esta experimetación numérica nos permite verificar algunos de los resultados obteni-
dos y discutidos en las secciones anteriores. El problema que se eligió para realizar la
experimentación numérica es el siguiente:
16


−u”(x) + c(x)u(x) = f (x), x ∈ Ω = (0, 1)
(5.1)
u(0) = u(1) = 0
donde:

 1, si x ∈ [0; 0,5)
c(x) = −R, si x ∈ [0,5; 1)
0 en caso contrario


y además


 −4x2 + 4x − 8, si x ∈ [0; 0,5)
f (x) = 4Rx2 − 4Rx − 8, si x ∈ [0,5; 1)
0 en caso contrario


Se tiene entonces un Problema de Dirichlet donde el coeficiente c(x) y la función

f (x) son funciones discontinuas. Se trata de hacer variar el parámetro R y ver como
es el comportamiento de la solución. Como es de nuestro interés, las bases wavelets
elegidas para aproximar la solución son las wavelets B-splines cardinales centralizadas
de orden 2 y 4, respectivamente. Se utilizó el método de Galerkin para discretizar el
problema y se determinó si los sistemas resultantes son bien condicionados o no, a
través de los números condición de la matrices de coeficientes. Mostramos resultados
para los sistemas sin precondicionamiento y con precondicionamiento. El método que
se utilizó para resolver los sistemas de ecuaciones lineales generados fue el de Con-
jugado Gradiente y Conjugado Gradiente Precondicionado, dadas las caracterı́sticas
de los sistemas a resolver; es decir, sistemas con matrices de coeficientes simétricas y
definidas positivas.
A continuación presentamos un listado de los resultados obtenidos (ver Cuadro 5.1,
Cuadro 5.2 y Cuadro 5.3); solo presentamos para los casos de R = 10 y R = 500,
los cuales son casos extremos. Dado que se cuenta con la solución explı́cita ha sido
posible obtener información que permite comparar el comportamiento de la solución
aproximada obtenida con la exacta y los errores correspondientes.

Nivel j Error res. Error l2 Nro. iter

4 6.48982e-14 9.72701e-02 15
6 6.55582e-13 2.55338e-02 63
8 9.74582e-11 6.46227e-03 255
10 2.83888e-09 1.61521e-03 1023
Cuadro 5.1
Errores de aproximación residual y l2 usando B-splines de orden 2, R=10, sin precondiciona-
miento

6. Conclusiones. De lo investigado, podemos establecer las siguientes conclu-

siones:
El método de Galerkin involucra la elección de una base del espacio que con-
tiene la solución, la formulación del Problema Variacional, la determinación
de la existencia de una solución del problema y la formulación del Problema
Discreto.
 17

Nivel j Error res. Error l2 Nro. iter

4 5.80149e-11 7.57413e-01 17
6 4.60217e-10 3.51286e+00 68
8 1.50453e-09 3.28961e-01 262
10 7.99947e-09 7.65207e-02 1048
Cuadro 5.2
Errores de aproximación residual y l2 usando B-splines de orden 2, R=500, sin precondiciona-
miento

Nivel j Error res. Error l2 Nro. iter

4 5.5351e-14 7.57413e-01 1
6 4.52399e-10 3.51286e+00 65
8 1.38431e-09 3.28961e-01 109
10 3.11206e-09 7.65207e-02 131
Cuadro 5.3
Errores de aproximación residual y l2 usando B-splines de orden 2, R=500, con precondicio-
namiento

Un B-spline cardinal genera un Análisis Multiresolución.

Es posible construir una base con B-splines cardinales de L2 (I) adaptando la
base con B-splines cardinales de L2 (R).
Aplicando el método de Galerkin al problema (4.1), se determinó que el Pro-
blema Variacional obtenido tiene única solución y es estable.
El orden de aproximación correspondiente a un método de Galerkin depende
de la potencia de aproximación de los espacios Test; en nuestro caso, los
espacios wavelets.
Usando bases wavelet con B-splines cardinales centralizados de orden 2, y
por la experimentación numérica realizada, se determinó que las aproxima-
ciones tienen bajo error y tienden a decrecer cuando se incrementa el nivel
de resolución j; es decir, cuando se trabaja en espacios de mayor resolución.
Dado que la matriz de coeficientes del sistema, correspondiente al problema
discreto, es mal condicionada, los resultados con alto orden de exactitud se
obtienen con un nivel j relativamente alto; sin embargo, cuando se usa pre-
condicionamiento, el error decrece exponencialmente, logrando altos niveles
de exactitud con un nivel de resolución bajo.

Referencias

[1] K. Amaratunga, J. Williams, S. Qian y J. Weiss, Wavelet-Galerkin solutions for one di-
mensional partial differential equations, IESL Technical Report No. 92-05, Intelligent En-
gineering Systems Laboratory, M.I.T., 1992.
[2] G. Beylkin, R. Coifman y V. Rokhlin, Fast wavelet transforms and numerical algorithms I,
Comm. Pure Appl. Math., 44(1991) pp. 141–183.
[3] H. Brézis, Análisis Funcional, Alianza Editorial S.A., Madrid, 1984.
[4] R. Burgos, R. Silva y M. Santos, Direct solution of differential equations using the wavelet-
Galerkin method, Mecánica Computacional, 29(2010) pp. 4573–4584.
[5] C. Chui y J-Z Wang, On Compactly Spline Wavelets and a Duality Principle, Transactions
of the American Mathematical Society, 330 (1992) pp. 903–915.
[6] I. Daubechies, Ten Lectures on Wavelets, SIAM Publications, Philadelphia, 1992.
[7] W. R. Glowinski, W. Lawton, M. Ravachol y E. Tenenbaum, Wavelet solution of linear
and nonlinear elliptic, parabolic and hyperbolic problems in one space dimension, Proc.
18

9th International Conference on Numerical Methods in Applied Sciences and Engineering,

SIAM, Philadelphia, 1990.
[8] S. Kesavan, Topics in Functional Analysis and Applications, Wiley Eastern Limited, New
Delhi, 1989.
[9] W. Lawton, W. Morrell, E. Tenebaum, y J. Weiss, The wavelet galerkin method for partial
differential equations, Technical Report AD901220, Aware, Inc., 1990.
[10] A. Louis, D. Maass y A. Rieder, Wavelets: Theory and Applications, Johb-Wiley & Sons,
USA, 1997.
[11] S. Mallat, A wavelet tour of signal processing, Academic Press, San Diego, 1999.
[12] M. Nair, Wavelet-Galerkin Method, Talk at QIP-short Term Course, Department of Mathe-
matics, IIT Madras, Diciembre, 2004.
[13] H. Nguyen y G. Evangelista, A continuous wavelet-Galerkin method for the linear wave
equation [Elektronisk resurs], SIAM Journal on Scientific Computing., 2007
http://liu.diva-portal.org/smash/get/diva2:24195/FULLTEXT01
[14] S. Qian y J. Weiss, Wavelets and the numerical solution of partial differential equations,
Journal of Computational Physics, 106 (1992) pp. 155–175.
[15] K. Urban, Wavelets in Numerical Simulation Problem Adapted Construction and Applications,
Ed. Springer, Heilderberg, 2002.
[16] K. Urban, Wavelet Methods for Elliptic Partial Differential Equations, Oxford Science Publi-
cations, Oxford, 2009.
[17] J. Wang, Cubic spline wavelet Bases of Sobolev spaces and multilevel interpolation, Appl.
Compu. Harm. Anal., 3 no.2 (1996) pp. 154–163.
[18] R. Wells y X. Zhou, Wavelet Solutions for the Dirichlet Problem, AFOSR Technical Report
No. 90-0334, Aware, Inc., 1993.