3 ALGORITMOS MATEMÁTICOS

3.1 INTRODUCCIÓN
Con el planteamiento físico del problema vimos que la mayor complicación surge de la
necesidad de calcular y representar la función b(V,s) de cada modo, así como sus derivadas
primera y segunda respecto de V. Para ello debemos calcular la constante de propagación
axial correspondiente a la frecuencia adimensional V, lo que nos lleva a encontrar los ceros
de la ecuación de dispersión. Esta es una ecuación algebraica que implica el cálculo de las
funciones de Bessel de primera especie. Con esta constante obtenemos el valor de b(V,s) y
sólo nos queda calcular, numéricamente, sus derivadas. Matemáticamente, todo esto nos
genera, seis problemas fundamentales:

1) Cálculo de las funciones de Bessel de primera especie, Jn(u) y Kn(w) ∀n∈{0}∪N 8,

∀u, w ∈ ℜ+. Serán los algoritmos más internos y más utilizados, por lo que
incidirán de forma notable en el rendimiento y precisión de los resultados.

2) Cuando nos acercamos al corte, w → 0+ y las funciones de Bessel modificadas de

primera especie tienden a infinito, por lo que es difícil evaluar la ecuación de
dispersión. Hay que dar mecanismos para resolver este problema.

3) Encontrar los ceros de distintas ecuaciones algebraicas. Es necesario utilizar el

algoritmo adecuado para cada caso, sopesando la robustez, rapidez y fiabilidad del
mismo9.

4) Al no existir una representación analítica, las derivadas de b(V,s) deben obtenerse

numéricamente. En general los sistemas de derivación numérica necesitan un gran
esfuerzo computacional10, por lo que sólo la usaremos para cálculos puntuales como
las condiciones de contorno en los extremos.

5) El cálculo de la función b(V,s) es muy laborioso de por sí, por lo que es necesario
un sistema de interpolación que nos permita obtener los valores de la función de
forma rápida, con precisión adecuada. Si se requiere el uso intensivo de la derivada,
o es preciso obtener derivadas de orden superior a uno, se hace necesario utilizar las
derivadas analíticas de la función de interpolación en lugar de la original. En este
caso, hemos de calcular con frecuencia tanto la primera, como la segunda derivada
de b(V,s). Por ello, el sistema de interpolación debe facilitar su cálculo.

8
En fibras monomodo, el valor de Vmax está limitado, teóricamente, a 2,4048. Para tratar las fibras multimodo
a salto de índice hemos ampliado a 20 el rango de V. Por tanto, grosso modo, n ∈ [0, Vmax], con Vmax ≤ 20.
9
Por robustez entendemos que garantice la convergencia y que avise adecuadamente si esta no se produce. La
rapidez dependerá del orden de convergencia, entendido como la relación entre los dígitos de precisión
obtenidos respecto del número de evaluaciones realizadas en cada iteración. La fiabilidad implica que, en caso
de convergencia, la solución obtenida es la que estábamos buscando y no otra.
10
La cantidad de cálculo varía exponencialmente con el orden de derivada.

31
 6) Definir el criterio al elegir los puntos que se utilizarán en la representación de cada
función, para obtener, en el menor tiempo posible, imágenes fieles de las mismas.

3.2 FUNCIONES DE BESSEL DE PRIMERA ESPECIE,

ORDEN ENTERO Y ARGUMENTO REAL
−1
∀n ∈ {0} U N, ∀u ∈ ℜ + ; Para n = 0, H 0 (u) =
J n (u)
Definimos Hn(u) = .
J n −1 (u) H 1 (u)

∀n ∈ {0} U N, ∀w ∈ ℜ + ; Para n = 0, G 0 (w) =

K (w) 1
Igualmente: Gn(w) = n
K n −1 (w) G 1 (w)

Para resolver la ecuación de dispersión necesitamos conocer el valor de Hn(u) y Gn(w),

mientras que el valor de los campos eléctricos y magnético viene dado por Jn(u) y Kn(w).
Vamos a concretar los rangos de cálculo, pues ello nos dará los criterios para decidirnos por
el mejor método. También hay que decidir la precisión de cálculo para estas funciones,
teniendo en cuenta la precisión de los resultados finales.

En los apartados 2.2.4, “Familias de modos reales y su denominación” y 2.2.6,

“Frecuencias adimensionales de corte”, vimos que u se encuentra acotada por:
A) Modos TE0,k, TM0,k, EHn,k : u∈ [Cn,k, Cn+1,k]
B) Modo HEn,k :
n = 1, k = 1: u ∈ [0, C0,1]
n = 1, ∀k ≥ 2: u ∈ [C1,k-1, C0,k]
∀n ≥ 2, ∀k ≥1: u ∈ [C*n-2,k, Cn-1,k]
C*n-2,k ≡ k-ésimo cero positivo de [Link]-2(x) + (s2-1).(n-1).Jn-1(x) = 0
C) Modo LPm,k :
m = 0, k = 1: u ∈ [0, C0,1]
m = 0, ∀k ≥ 2: u ∈ [C1,k-1, C0,k]
∀m ≥ 1, ∀k ≥1: u ∈ [Cm-1,k, Cm-1,k]

La ecuación de dispersión de los modos TE0,k y EHn,k depende de Hn+1(u), la de los

modos TM0,k, la de LP0,k y la de HE1,k de H1(u), la de los HEn,k, con n ≥ 2, de 1/Hn(u) y la
de los modos LPm,k , con m ≥ 1, de 1/Hm(u). Por otro lado, usando el desarrollo asintótico
para n grande del primer cero:
C0,1 = 2,4048 > 2 = 0 + 2 n=0
Cn,1 ∼ n +1.8557571.n1/3 + 1.033150.n-1/3 + O(n-1) > n +2 ∀n ≥ 1

Esto nos permite acotar más estrechamente el rango de n y u

Hn(u): n∈N con 1≤ n ≤ Vmax y n ≤ u ≤ Vmax.
Jn(u): n∈{0}∪N con 0 ≤ n ≤ Vmax y 0 ≤ u ≤ Vmax.
Gn(w): n∈N con 1≤ n ≤ Vmax y 0 ≤ w ≤ Vmax.
Kn(w): n∈{0}∪N con 0 ≤ n ≤ Vmax y 0 ≤ w ≤ Vmax

32
 Por último, hemos de decidir la precisión requerida, simple o doble11. En el primer caso
trabajaríamos con variables tipo Single, en el segundo usaríamos el tipo Double. El primer
tipo de variables ocupa la mitad de memoria y las operaciones son muchísimo más rápidas.
Además, los algoritmos serán más sencillos y rápidos de ejecutar. Necesitamos unos tres
dígitos de precisión en los resultados.

Ya vimos que Hn(u) puede estar en la base de dos etapas de cálculo. Primero, solución
numérica de la ecuación de dispersión para obtener las curvas de propagación axial
adimensional. Segundo, y más desfavorable, cálculo de la dispersión intramodal derivando
dos veces las curvas anteriores. Suponiendo que los algoritmos de cálculo de Hn(u) nos
ofrecen resultados exactos hasta el último bit. Si en la solución de la ecuación “gastamos”
un dígito de precisión y en las dos derivaciones “consumimos” otros tres, podríamos usar
variables tipo Single. De todas formas, la cosa va muy justa, por lo que es prudente estudiar
también los algoritmos en doble precisión12.

El cálculo de Jn(u) sólo es necesario, en principio, para tener los valores de los campos
en el núcleo, por lo que su uso es directo, y la precisión simple es suficiente.

Vamos a considerar los diversos métodos para el cálculo de Hn(u) y Jn(u):

k
 u2 
− 
u  ∞  4 
n
1) 
Utilizando la definición en serie de Jn(u) =   .∑ ∀n ∀u.
 2  k = 0 k! (n + k)!
Converge rápidamente cuando u << n. Si u es del orden de n o mayor, la
convergencia es lenta y aun peor, tenemos cancelaciones y pérdidas de precisión.

2 .( n − 1)
2) Recurrencia. J n ( u ) = J n −1 ( u ) − J n − 2 ( u ) ∀ n
u
Estudiaremos el sentido en que la relación es estable13. Hay dos tramos:

0 ≤ u ≤ n: la relación es estable en sentido descendente. Usamos el algoritmo de

Miller. Consiste en empezar por un orden m: m>>n. Empezamos con un valor
arbitrario14, J*m(u) = 0 y J*m-1(u) = 1. Utilizamos la relación de recurrencia hasta
llegar a J*0(u). Al ser lineal la relación de recurrencia, J*k(u) = [Link](u) ∀k. Basta
con obtener el valor de la constante. Para ello usamos la propiedad:
∞ ∞
1 = J 0 (u) + 2.∑ J 2k (u) , o en nuestro caso: C = J *0 (u) + 2.∑ J *2k (u)
k =1 k =1

11
Precisión simple significa unos 7 dígitos de precisión. En doble precisión tenemos unos 15 dígitos.
12
Empezamos con algoritmos de precisión simple (variables Single). Todo fue bien hasta que llegamos a las
curvas de dispersión intramodal. Tuvimos que cambiarlo todo a doble precisión (variables Double).
13
En este caso, estabilidad significa que los errores de truncamiento no se amplifican al recorrer la relación.
Cuando una relación de recurrencia es inestable en un sentido significa que estable en el otro.
14
Tomamos estos valores iniciales pues lim J k (u) = 0 , ∀u ≤ n.
k→∞

33
 El único problema es elegir adecuadamente m. Para ello tenemos una regla
heurística, por la que m = n + P.n1/2, donde P es una constante cuyo valor es igual
al número de dígitos de precisión que deseamos alcanzar. En este caso tenemos
que realizar m iteraciones. Si queremos calcular Hn(u) no hace falta normalizar
pues al dividir la función de orden n por la de orden n-1 desaparece la constante
de normalización. En este caso, realizaremos P.n1/2 iteraciones.

n ≤ u: aquí la relación es estable en sentido ascendente. Empezamos calculando

J0(u) y J1(u) y usando la relación de recurrencia llegamos a Jn-1(u) y Jn(u).

El cálculo de J0(u) y J1(u) depende de la precisión que se requiera:

Para precisión simple tomaremos un valor de referencia uo >>1. Si 0≤ u ≤ uo, las
aproximaremos mediante funciones racionales. Para uo < u, nos basaremos en el
desarrollo en series de Hankel. Las funciones Pn(u) y Qn(u) las aproximaremos
mediante polinomios que den la precisión deseada con ordenes bajos. Aun así, es
necesario evaluarlos internamente en doble precisión para evitar errores.
Para precisión doble este método no es adecuado, ya que, al no ser uniforme el
desarrollo de Hankel, necesitaríamos que uo fuese muy grande, Vmax < uo, por lo
que tendríamos que trabajar sólo con funciones racionales. Pedir precisión doble
en un intervalo grande nos llevaría a usar polinomios de grado elevado, sensibles
al redondeo y a las cancelaciones.

3) Fracciones continuas: este método converge en nuestro rango de trabajo. Sea el

punto de inflexión up = n.(n + 1) ≈ n, a partir del cual Jn(u) se hace oscilatoria.
Si 0 ≤ u < up, entonces el tiempo de convergencia crece exponencialmente con u
y casi no depende del orden n . Para up < u, cada iteración es equivalente a hacer
que el orden aumente en uno, hasta que u < up, y volvemos al caso anterior. Para
evaluar las fracciones continuas usaremos dos métodos: el algoritmo de Steed y
el modificado de Lentz. El primero es más rápido, sin embargo, en el cálculo de
las funciones de Bessel de 1ª especie utilizaremos el segundo para evitar pérdidas
de precisión.

Sea una función f(x) dada mediante una expresión en fracciones continuas y sea
fn el resultado de evaluar, de izquierda a derecha, la expresión hasta los
coeficientes an y bn. Por ejemplo, para f3, evaluaríamos hasta b3, esperando que el
resto de la expresión, los puntos suspensivos, sea despreciable frente a b3.
a1
f(x) = b 0 +
a2
b1 +
a3
b2 +
b3 + ⋅⋅ ⋅
El método de Steed consiste en ir calculando los incrementos de fj, ∆fj = fj – fj-1.
Por inducción se pueden demostrar las siguientes relaciones:
∆fj = ([Link] –1).∆fj-1 y Dj = 1/(bj + [Link]-1), ∀j = 2,...,n. Con D1 = 1/b1, ∆f1 = a1/b1 y
n
f0 = b0. Evidentemente fn = f0 + ∑ ? f j
j=1

34
La condición de finalización será: |(∆fj)/fj-1| < ε, con ε la precisión requerida en
los cálculos. Esta condición no es de validez general, pero sirve perfectamente
para los casos que vamos a estudiar.

El algoritmo presenta el inconveniente de que si el denominador de la expresión

de Dj se hace pequeño entonces Dj y ∆fj se hacen muy grandes. En la siguiente
iteración ∆fj+1 se hará pequeño, pero se pierde precisión al evaluar los fj como el
sumatorio de los incrementos. Esto es muy difícil de evitar, por lo que el método
de Steed sólo se utiliza cuando se sabe por adelantado que no se dará este caso.

El algoritmo modificado de Lentz es una evolución sobre el de Steed para evitar

el problema anterior. Es algo más complicado y, por tanto, más lento. Por ello,
sólo será utilizado en aquellos casos en que el uso del de Steed no sea seguro. La
implementación genérica de este algoritmo sería:

• f0 = b0
• Si f0 = 0, entonces f0 = δ
• C0 = f0
• D0 = 0
• Bucle para j =1, 2,...
o Dj = bj + [Link]-1
o Si Dj = 0, entonces Dj = δ
o Cj = bj + aj/Cj-1
o Si Cj = 0, entonces Cj = δ
o Dj = 1/Dj
o ∆j = [Link]
o fj = fj-1.∆j
o Condición de salida del bucle
Si |∆j –1| < ε, entonces:
Fin
• Repetición del bucle

Donde: ε la precisión requerida, con 0 < ε << 1

δ un número real, con 0 < δ < ε.|bj| ∀j

Vamos a particularizarlo para el cálculo de la función Hn(u):

1
Partiendo de la relación de recurrencia llegamos a: H n (u) =
2.n
− H n +1 (u)
u
Por tanto: b0 = 0, bi = 2.(n-1+i)/u ∀i ≥ 1 a1 = 1, ai = -1 ∀i ≥ 2

Finalmente, hemos de destacar que, aunque no hay estudios teóricos sobre la

propagación de errores en ambos algoritmos, las pruebas muestran un buen
comportamiento.

35
 4) Expansiones en series.
Expansión asintótica de Hankel para argumentos grandes ≡ u >> n:
Jn(u) ≈
2
{Pn (u).cos( ? ) − Q n (u).sen( ? )}, con ? = u − p  n + 1  .
pu 2 2
∞
(n,2k) ∞
(n,2k + 1)
donde Pn(u) ≈ ∑ ( − 1) k y Qn(u) ≈ ∑ ( − 1) k .
k =0 (2.u) 2k
k=0 (2.u) 2k +1
(u) los términos k-ésimos de Pn(u) y Qn(u) y µ = 4.n :
2
Sean TPnk (u) y TQn
k

TPnk (u) = -
(µ - (4.k - 1) ).(µ - (4.k - 3) ) T
2 2
k -1
(u) ∀k ≥ 2
(2k − 2 )(. 2k − 3 )(. 8.u ) 2 Pn

k
TQn (u) = -
(µ - (4.k + 1)2 ).(µ - (4.k - 1)2 ) T k -1 (u) ∀k ≥ 2
(2k − 1).2 k.(8.u )2 Qn

Cuanto mayor es u con respecto al orden, más rápido converge. Si evaluamos

Pn(u) y Qn(u) hasta el término k-ésimo: k > n/2, verificaran que el resto es menor,
en valor absoluto, que el valor absoluto del término k+1-ésimo y del mismo
signo. Esto es interesante pues tengo que evaluar sólo n/2 términos y tengo
garantizada la precisión, excepto errores de truncamiento. El problema es que la
precisión alcanzable depende de u. A mayor u, mayor precisión.

Expansión asintótica de Debye para ordenes grandes. Es complicada y no vale

para ordenes pequeños.

Expansiones asintóticas uniformes. Son complicadas y lentas.

5) Métodos integrales que son muy lentos y complicados de implementar.

Por rapidez y sencillez, nos decantamos por el método de fracciones continuas y el de

recurrencia. Entre ellos, la elección depende de la función a calcular, del rango de trabajo y
de la precisión requerida.

Precisión simple.
En primer lugar, el método de las fracciones continuas calcula Hn(u) directamente, pero
necesitamos otra relación para calcular Jn(u). El método de recurrencia va calculando una
sucesión de ordenes hasta llegar a n, por lo que obtenemos Jn-1(u) y Jn(u) simultáneamente y
sólo con dividir obtenemos Hn(u). Además, en este caso no hace falta normalizar, pues al
dividir eliminamos la constante de normalización., y por tanto, no es necesario seguir hasta
n = 0, ahorrándonos n iteraciones.

Vamos a comparar los tiempos de ejecución de ambos métodos para el cálculo de

Hn(u). Vimos que para las fracciones continuas, el tiempo de convergencia varía de forma
aproximadamente exponencial con u y casi no depende de n. Para la recurrencia el tiempo
no depende de u y varía linealmente con n. También observamos que la recurrencia hacia
abajo es algo más lenta que hacia arriba, ello es debido, a que en el primer caso realizamos
P.n1/2 iteraciones y en el segundo n. Pero en precisión simple P = 6 > Vmax1/2 ≈ 4.

36
 1

FC ( u)

Rec( 13 , u)
Tiempos

0.5
Rec( 10 , u)

Rec( 5 , u) Recurrencia Recurrencia

hacia abajo hacia arriba

0 5 10 15
u
Evidentemente, dado un n fijo existirá un valor de u, al que denominaremos Un(n), para
el que el tiempo de ejecución de las fracciones continuas es igual al de recurrencia.
Variando n dentro del intervalo máximo de trabajo, [0, Vmax = 20], e imponiendo un ajuste
exponencial llegamos a:
na
Un(n) = C a −1 , con a = 2.2, b = 20 y C = 1.1
n +b
Para sacar el máximo rendimiento usaríamos fracciones continuas para 0 < u < Un(n) y
recurrencia para Un(n) < u. Dado que Un(n) < n, ∀n: 0 < n ≤ Vmax ≤ 20 y que trabajaremos
con u > n, el método más rápido es el de recurrencia ascendente. Finalmente, puesto que la
diferencia entre n y Un(n) es grande, asumir que el resultado experimental obtenido es un
comportamiento intrínseco a los algoritmos utilizados, tiene validez en general y que no
depende del hardware utilizado.
15

Recurrencia
hacia arriba

10

Un ( n)

n Recurrencia
hacia abajo
5

Fracciones
continuas

0 2 4 6 8 10 12 14
n

37
 Resumiendo:
Hn(u):
- Recurrencia ascendente.
Jn(u):
- ∀n: si 0 ≤ u < n, recurrencia descendente. Si n ≤ u, recurrencia ascendente.

Precisión doble.
Hn(u): fracciones continuas en todo el rango utilizando el método de Lentz modificado.
Jn(u): ya vimos que no será necesario su cálculo con tanta precisión.

Por último destacar que obtenemos precisión relativa lejos de los ceros de Jn(u) = 0.
Cuando nos acercamos a esos puntos la precisión obtenida es absoluta. Lo mismo sucede
con Hn(u) pero con los ceros de Jn(u) = 0 y los de Jn-1(u) = 0.

Pasamos ahora a considerar los diversos métodos para el cálculo de Gn(w) y Kn(w):

1) Expansión en series ascendentes.

1 2 n −1 (n − k − 1)! w 
n 2k

Kn(w) =   ∑  
w
+ ( − 1) n +1 ln   I n (w) +
2  w  k = 0 (-1) k!  2 
k
 2
(−1) n  w 
n ∞
[? (k + 1) + ? (n + k + 1)] w  2k , ∀n∈{0}∪ N, ∀w, con:
+
2 2
  ∑
k =0 k!(n + k)!
 
2
n −1
1
ψ(n) = -γ + ∑ ∀n∈N, siendo γ = 0.57721..., la constante de Euler.
k =1 k
n 2k
w ∞ 1 w
In(w) =   ∑   ,∀n∈{0}∪ N, ∀w
 2  k =0 k!(n + k)! 2 
La expansión es complicada, sobre todo para n medio o grande. Sólo vale para w
pequeños. La usaremos para evaluar K0(w) y K1(w), con 0< w < wo y wo ≈ 1.

2 .( n − 1)
2) Recurrencia. K n ( w ) = K n −1 ( w ) + K n − 2 ( w ) ∀ n
w
Hay que estudiar el sentido de la estabilidad:

La relación es estable en sentido descendente ∀w ∈ ℜ+ y ∀n∈{0}∪ N. De todas

formas, Kn(w) creciente con n para un w cualquiera fijo, por lo que, también
puede usarse en sentido ascendente, ya que se mantiene estable el error relativo.

En resumen, podemos usar la relación en cualquier sentido. Sin embargo, en este

caso el algoritmo de Miller no es de aplicación, por tanto usaremos la relación de
recurrencia en sentido ascendente en todo el rango de trabajo. Queda por
determinar como vamos a calcular K0(w) y K1(w) con la precisión requerida.

38
3) Fracciones continuas. Para Gn(w), por la relación de recurrencia llegamos a:
1
G n (w) = ⇒ b0 = 0, bi = -2.(n-1+i)/w ∀i ≥ 1 ai = 1 ∀i ≥ 1
2.n
− + G n +1 (w)
w
La evaluaremos mediante el algoritmo de Steed, ya que, podemos demostrar que
bj + [Link]-1 es siempre un número negativo y alejado de 0.
El problema es que cuando w→ 0+ Gn(w)→2.(n-1)/w y por tanto, bi será de igual
orden que la expresión a su derecha y no se verificará la condición en la que se
basan los algoritmos de fracciones continuas.

p −w  ∞
(n, k) 
4) Expansión asintótica para u grande ≡ u >> n. Kn(w) = e 1 + ∑ k .
2.w  k =1 (2.u) 
La precisión alcanzable depende de lo grande que sea u respecto de n. A mayor
u/n, mayor precisión. Como u y n están en el rango de 0 a 20 la relación u/n no
va a ser muy elevada, por lo que esta expresión sólo será útil en precisión simple
y para calcular K0(w) y K1(w) en la zona w > wo, con wo > 1.

5) Utilizando las funciones hipergeométricas confluentes.

Sea ν fijo. Definimos la secuencia de funciones hipergeométricas confluentes:
zm(w) = U(ν + ½ + m, 2.ν +1, 2.w), con m∈{0}∪N.
En estas condiciones, utilizando la expresión que liga las funciones modificadas
de Bessel de primera especie, a las hipergeométricas confluentes, tenemos:
Kν(w) = π1/2.(2.w)ν.e-w.z0(w) y teniendo en cuenta las relaciones entre funciones
hipergeométricas confluentes contiguas llegamos a:
1  1  2 1  z 1 ( w)  .
G ν +1 (w) = ν + + w + ν −  
w 2  4  z 0 ( w) 
Además, las funciones zm(w) verifican la siguiente relación de recurrencia:
zm-1(w) = bm(w).zm(w) + am+[Link]+1(w), donde: bm(w) = 2.(m+w)
am+1 = ν2 - (m + ½ )2
Esto nos permite llegar a la expresión en fracciones continuas:
z1 (w ) 1
= , que evaluamos por el método de Steed.
z 0 (w ) a2
b 1 ( w) +
a3
b 2 ( w) +
b 3 ( w) + ⋅ ⋅ ⋅
Un problema para este método es que no converge cuando w es pequeño.
Para calcular Kν(w) usaremos la condición de normalización de Temme:
 1 
∞  1 Γ ν + + m 
- ν + 
 
m
(-1) 2
∑C
m =0
m .z m ( w) = ( 2.w)  2
, con C m =
m!  1 
Γ ν + − m 
 2 
a m +1
Los coeficientes Cm cumplen ley de recursión: Cm+1 = − C m , C0 = 1.
m +1

39
 ∞
z m ( w)
Si definimos: A(w) = ∑ C m 1 1
, entonces z 0 ( w) =
m =1 z 0 ( w) ν+
1
2
1 + A(w)
( 2.w)
Si queremos tener simultáneamente Kν(w), Kν+1(w) hay un método para calcular
A(w) al mismo tiempo que ejecutamos el algoritmo de Steed.
∞
z1 (w )
Si
z 0 (w )
= ∑ ∆h
m =0
m
, con los términos ∆hm calculados por el algoritmo de Steed.
N
z m ( w)
Sea AN la aproximación a A(w) tras sumar N términos. A N (w) = ∑C
m =1
m
z 0 ( w)
N m
Podemos poner A N (w) = ∑T
m =1
m .∆ h m , con Tm (w) = ∑C
k =1
k .t k , los coeficientes tk

verifican la siguiente ley de recurrencia tk+1 = (tk-1 –bk. tk)/ ak+1, con t0 = 0 y t1 = 1.
Otra limitación de este método es que, tanto la relación en fracciones continuas
como la condición de normalización de Temme, son válidas para ν: |ν| ≤ ½. Por
tanto, con ν∈{0}∪N y |ν| ≤ ½, sólo está ν = 0 obteniendo G1(w), K0(w) y K1(w).

6) Series para argumentos pequeños.

∞ ∞ 2.k
2 1 w
K ν (w) = ∑ c k .f k y K ν +1 (w) =
k =0 w
∑ c k .h k , con: c k =
k =0
 
k!  2 
, hk = -[Link] + pk,

fk = ([Link]-1 + pk-1 + qk-1)/(k2 - ν2), pk = pk-1/(k - ν) y qk = qk-1/(k + ν). Los valores

iniciales para estas relaciones de recurrencia son:
-ν ν

p0 =   Γ (1 + ν ) , q0 =   Γ (1 − ν )
1 w 1 w
2 2  2 2 
p.?  senh (s )  2  
f0 =  cosh (s ).G1 (? ) + ln  .G2 (? ) , con σ = ν.ln(2/w), además
sen (p.? )  s w 
1  1  1 1 
G1 (? ) = y G2 (? ) = 
1 1
 −  + .
2.?  G(1 - ? ) G(1 + ? )  2  G(1 - ? ) G(1 + ? ) 
senh (s ) sen (x )
Sabiendo que: G2 (0 ) = 1 , lim G1 (? ) = − ? , lim = 1 y lim = 1.
?→0 s →0 s x→0 x
Para ν=0, simplifica y llegar a: p0=½, pk=qk∀k≥0, f0 = -γ + ln(2/w) y f0 = -γ f0 = -
γ.
De nuevo, este método sólo será válido si ν: |ν| ≤ ½. Por tanto, con ν∈{0}∪N y
|ν| ≤ ½, sólo está ν = 0, con lo que podremos calcular G1(w), K0(w) y K1(w).

7) Métodos integrales. Complicados y lentos.

Podríamos usar fracciones continuas cuando w > wt y series ascendentes para w < wt. El
problema es saber cuando usar un sistema u otro, pues wt depende del orden . Además, esto
sólo sirve para calcular Gn(w). Por ello usaremos recurrencia ascendente en todo el rango
de trabajo y para cualquier precisión. La forma de calcular K0(w) y K1(w) dependerá de la
precisión que queramos obtener.

40
 Para precisión simple, tomaremos un valor de referencia wo, que nos divida la semi-
recta positiva en dos partes. La zona 0< w ≤ wo, donde usaremos una aproximación basada
en el desarrollo en series ascendentes. wo< w, para la cual es útil una aproximación basada
en la expansión asintótica para órdenes grandes. Para precisión doble, al no ser uniforme la
expansión asintótica, tendríamos que ir a un wo muy grande y Vmax < wo, o bien, tener un
orden polinómico muy alto en la aproximación. Si wo es grande, tendríamos que trabajar
sólo con la aproximación basada en el desarrollo en series ascendentes, lo que nos llevaría a
usar polinomios de grado elevado. Tanto en un caso como en otro acabamos usando
polinomios de grado elevado, y por tanto, sensibles al redondeo y a las cancelaciones. En
resumen, no podemos usar este método para doble precisión.

Para precisión doble usaremos las funciones hipergeométricas confluentes y las series
para argumentos pequeños.

3.3 EXPANSIONES ASINTÓTICAS PARA LA ECUACIÓN

DE DISPERSIÓN EN LA ZONA CERCANA AL CORTE
Cerca del corte w→0, y Kl(w) → +∞, ∀l ≥ 0. Esto hace que la evaluación de la
ecuación de dispersión sea poco precisa, con resultado poco fiables y perjudicando la
robustez del algoritmo para la obtención de las raíces. Usando las expansiones en series
ascendentes de las funciones de Bessel modificadas de primera especie cerca de 0 podemos
obviar estos problemas, evitar cancelaciones y llegar a una ecuación de dispersión más
sencilla y precisa en el rango cercano al corte. A partir de esta ecuación, podemos llegar a
expresiones de la forma w = w(u,s) que, para algunos casos, simplemente sustituyendo u
por V15, lograríamos poner w en función de V, con lo que obtenemos la constante
adimensional de propagación b en función de la frecuencia adimensional V y de s.

La ecuación y la expresión b(V,s) son válidas dentro de un cierto rango de frecuencias

adimensionales cercanas a la frecuencia de corte. Si este rango es suficientemente grande,
podemos utilizar directamente la función b(V,s). Si el rango es pequeño, entonces estas
expresiones nos permitirán obtener los valores en el corte. En el resto usamos la expresión
general de la ecuación de dispersión. Por último, indicar que de la ecuación en el corte, por
derivación respecto de V, obtenemos expresiones más sencillas para poder calcular el valor
de las derivadas en el corte.

1) Modos TE0,k, TM0,k ∀k ≥ 1: sea K=1 modos TE

K = s-2 modos TM
La frecuencia adimensional de corte es la solución de la ecuación16:
J0(x) = 0 equivalente a H1(x) = ∞.

15
V2 =u2 + w2, por lo que en el corte, cuando w→0, u→V.
16
En las soluciones a las ecuaciones para obtener las frecuencias de corte nunca se incluye la solución x = 0, a
no ser que se indique expresamente lo contrario.

41
 H 1 (u) K
0= − ; Ecuación de dispersión en el corte
u   w  
w 2 .ln  + ? 
 2 
K.u(V)
−?
w(V) = 2.e [V ]
2
- u(V) 2 .H1 (u(V))
que nos sirve para obtener las derivadas en el corte17.

En este caso no podemos obtener una expresión funcional del tipo b(V,s), pero
de todas formas el rango de validez de la aproximación es pequeño, por lo que no
nos es de utilidad práctica.

2) Modos EHl,k ∀l ≥ 1,∀k ≥ 1:

H (u) s 2 + 1 l
0 = − l +1 − 2 ; Ecuación de dispersión en el corte
u s w2
La frecuencia adimensional de corte es la solución de la ecuación:
Jl(x) = 0 equivalente a Hl+1(x) = ∞.
s 2 + 1 u(V)
w(V) = −l , de aquí obtenemos las derivadas en el corte
s 2 H l +1 (u(V))
s2 +1 V
w(V) = − l 2 , expresión en función de V. No la usaremos pues su
s H l +1 (V)
rango de validez es muy pequeño

3) Modos HEl,k ∀l ≥ 1,∀k ≥ 1:

−1 2.r (w)
0= + 2l ; Ecuación de dispersión en el corte.
u.H l (u) s + 1
  w 
 − ln  + ?  l =1
  2 
 1  2.A A w w
2
Donde: rl (w) =  +  2 − + +
   l=2
2  u 2. s 2 + 1
ln
2  4(
?
)
 1  A  
 1 +  2 + B .w 2  l≥3
 2.(l − 1)   u  

yA =
(
s2 −1
2
)
; B=
s2 −1 
2


( 1 ) s2 −1
+ 
2
1 
−
(
s2 +1 )2

(
2. s 2 + 1 )2.s 2 .(l − 1)  4.l.(l − 2)  s 2 + 1 

8.l  16.s 2 .(l − 2)


Así tenemos:
Para l = 1,Vc es la solución a la ecuación J1(x) = 0, equivalente a H1(x) = 0.
En este caso si se incluye la solución Vc = 0.
Para l ≥ 2, Vc es la solución a la ecuación (s2 +1).(l-1).Jl-1(x) - [Link](x) = 0, o
lo que es lo mismo, (s2 +1).(l-1) – [Link](x) = 0. Teniendo en cuenta las

17
Junto a: u´(V)= [V-w(V).w´(V)]/u(V), que se obtiene derivando respecto de V la expresión: V2 =u2 + w2.

42
 relaciones de recurrencia de las funciones de Bessel podemos ponerlo como,
[Link]-2(x) + (s2 –1). Jl-1(x) = 0, o, x + (s2 –1).Hl-1(x) = 0. Si s→1, la ecuación
tiende a Jl-2(x) = 0, o, Hl-1(x) = ∞.

• l = 1:
 s 2 +1 
− +? 
 2.u(V).H (u(V)) 
 1 
w(V) = 2.e , nos permite obtener las derivadas
 s 2 +1 
− +? 
 2.V.H (V) 
 1 
w(V) = 2.e , usamos esta forma directa en lugar de la ecuación de
dispersión, pues su rango de utilidad es suficientemente amplio. Este rango se
calcula haciendo que los términos que despreciamos en el desarrollo de la
ecuación en el corte sean menores que la precisión de cálculo. Estos términos
son del orden de [Link](w). Sea wo: -[Link](wo) = Pm, con Pm la precisión
máxima alcanzable, 0 < Pm < 1. Sea K = -ln(wo/2)-γ, entonces el valor límite
de frecuencia adimensional para el que se puede utilizar las simplificaciones
en el corte, Vo, verificará:
s 2 +1
wo = w(Vo), o sea: Vo .H 1 (Vo ) = . Tiene infinitas soluciones, tomamos
2.K
la más próxima a la frecuencia adimensional de corte del modo en cuestión.

Para k=1, se podría hacer una simplificación mayor haciendo el desarrollo

en series de H1(V); H1(V) = V/2 + O(V2). Pero tomaríamos un rango mayor
del real. La diferencia entre los dos rangos sólo sería grande cuando s >>1.

• l = 2:
2  s 2 +1  A 2.A
−1 + − −?
w(V) = 2.e [V - u(V) ] ( )

  2. s +1 u(V)
2 2 u(V).H 2 2
2 (u(V))
, de aquí podemos obtener las
derivadas en el corte.
No hay expresión explícita para b(V,s).

• l ≥ 3:
(l - 1).(s 2 + 1)
−1
u(V).H l (u(V))
w(V) = , de aquí obtenemos las derivadas en el corte.
A.u(V) − 2 + B
(l - 1).(s 2 + 1)
−1
V.H l (V)
w(V) = , no la usaremos al ser muy pequeño el rango útil.
A.V −2 + B

4) Modos LPm,k ∀m ≥ 0,∀k ≥ 1:

• m = 0:
−1  w 
0= − ln  + ?  ; Ecuación de dispersión en el corte.
u.H 1 (u)   2  

43
 La frecuencia adimensional de corte es la solución de la ecuación:
J1(x) = 0 equivalente a H1(x) = ∞.
 1 
−  + ? 
 u(V).H1 (u(V)) 
w(V) = 2.e , nos permite obtener las derivadas. De hecho este
caso es igual al de HE1,k, poniendo s = 1.
 1 
− + ? 
 V.H1 (V) 
w(V) = 2.e , usamos esta forma directa en lugar de la ecuación de
dispersión, pues su rango de utilidad es suficientemente amplio. Este rango se
calcula exactamente igual que para HE1,k, poniendo s = 1.

• m≥1
u
0= + rm (w) ; Ecuación de dispersión en el corte.
H m (u)
  w  2
− ln  + ? .w m =1

Donde: rm (w) =     2 
2
 w
m≥2
 2.(m − 1)

La frecuencia adimensional de corte es la solución de la ecuación:

Jm-1(x) = 0 equivalente a Hm(x) = ∞.

o m = 1:
u(V)
−?
w(V) = 2.e [V ]
2
- u(V) 2 .H1 (u(V))
, para obtener las derivadas en el corte.
No hay expresión explícita para b(V,s).

o m ≥ 2:
2.(m - 1).u(V)
w(V) = − , para obtener las derivadas en el corte.
H m (u(V))
2.(m - 1).V
w(V) = − , no la usaremos al ser pequeño el rango útil.
H m (V)

3.4 RAICES DE ECUACIONES ALGEBRAICAS

Tenemos que resolver la ecuación de dispersión y hallar las frecuencias adimensionales
de corte. El primer caso, se basa Gn(w) y Hn(u) y la ecuación tiene una forma compleja. En
el segundo usaremos Hn(u) y la ecuación es más sencilla. El cálculo es costoso, por ello, es
vital tener el mayor orden de convergencia. Esto nos lleva al método de Newton-Raphson,
con convergencia cuadrática. El nº de dígitos de precisión se dobla en cada evaluación.
Pero la capacidad de este método para hallar el resultado deseado depende del valor inicial
utilizado. Hemos de dar valores iniciales generales que lleven a una respuesta fiable y que
sean lo más precisos posible, pues de este modo influimos en la rapidez de la respuesta.

44
 Para la ecuación de dispersión partiremos con el valor inicial dado por la aproximación
de guía débil del apartado [Link], “Expresión analítica para la constante de propagación
axial”. Este valor inicial es válido para nuestro rango de trabajo, s ∈ [1,smax]. El algoritmo
converge de forma bastante rápida, 2 ó 3 evaluaciones en casi todo el rango. El problema
está, sobre todo, en la zona de corte y fundamentalmente para los modos LP0,1 y HE1,1. En
ese caso utilizamos los resultados del apartado anterior. Para mayor seguridad, acotamos el
cero y utilizamos un híbrido entre el Newton-Raphson (N-R) y el algoritmo de bisección,
lento pero seguro. De esta forma, la bisección sólo entra cuando el N-R se sale de la
acotación o no converge con suficiente rapidez. El criterio de control es sencillo y no
penaliza prácticamente el algoritmo.

Para la obtención de las frecuencias adimensionales de corte tenemos dos casos:

1) Modos HEn,k, n ≥ 2: hemos de obtener los ceros de [Link]-2(x) + (s2-1).(n-1).Jn-1(x) = 0.
Aplicaremos N-R unido con bisección para evitar problemas de convergencia. Los
ceros estarán entre Cn-2,k (s=1) y Vsup, que es el mínimo entre Vmax y Cn-1,k (s → +∞).
Como valor inicial, Vini, tomamos: si s ≤ 1.3: Vini = Cn-2,k + s*(Vsup – Cn-2,k)
si s > 1.3: Vini = Cn-2,k + 0.75*(Vsup – Cn-2,k)
2) El resto de modos: aquí se reduce a obtener los ceros de Jn(x) = 0. Usaremos el N-R
puro, sin ni siquiera acotar, pues basta con una adecuada elección del valor inicial.
Vamos a ver como elegimos esos valores. Vi denota el valor inicial elegido.
Cn,1: n = 0, Vi = 2,404826
n ≥ 1, Vi = n + 1.85575.n1/3 + 1.03315.n-1/3
Cn,2: n ≤ 10, Vi = Cn,1 + π
n >10, Vi = Cn,1+0,76*(Cn,1 - n)
Cn,3: n ≤ 10, Vi = Cn,2 + π
n >10, Vi = Cn,2+0,82*(Cn,2 - Cn,1)
Cn,k y k ≥ 4: n ≤ 10, Vi = Cn,k-1 + π
n >10, Vi = Cn,k-3+(Cn,k-3 - Cn,k-2)2/(Cn,k-2 - Cn,k-1)
El N-R converge si la distancia entre dos resultados consecutivos es menor que la
tolerancia. Vamos a comprobar que este criterio es válido para los casos propuestos. Ya
vimos el comportamiento de los términos relativos a las funciones de Bessel de primera
especie, Js(u). Para el resto de los términos teníamos:

Modos EH y TE (l = 0):
2
l (1 + s 2 )  l  s2 −1  − l  l2  1 1   s2 1 
t (u, w, s) = 2 +
l
 − g l
(w)  −   + g l
(w)   +  + .
2  2
+ 2 
u 2.s  w
2 2
  2.s  w
2 2
 s u
2 2
w  u w 
con tl(u,w,s) < 0 en todo el rango u,w ∈ [0,V]

Modos HE y TM (l = 0):
2
l (1 + s 2 )  l  s2 −1  l  l2  1 1   s2 1 
t l (u, w, s) = 2 +  + g l (w)  −   + g l (w)   +  + .
2  2
+ 2 
u 2.s  w
2 2
  2.s  w
2 2
 s u
2 2
w u w 
con tl(u,w,s) > 0 en todo el rango u,w ∈ [0,V].

45
 Modos LP:
- K 0 (w ) - K 0 (w) −1
m = 0, − g -1 ( w ) = = = 2 1
w.K −1 ( w ) w.K 1 ( w ) w .g ( w )
t m (u, w) = −g m −1 (w);
- K m (w)
m ≥ 1, − g m −1 ( w ) =
w.K m −1 ( w )
m
con t (u,w) < 0 en todo el rango u,w∈ [0,V].
Kl − 1(w) Kl + 1(w)
con : g l (w) ≡ ; g l (w) ≡ , V2 =u2 + w2: 0 ≤ u,w ≤ V
[Link](w) [Link] (w)

Nos limitaremos a ver su comportamiento en los extremos del rango de definición:

1) Cuando u → 0+ y w → V-.
2) Cuando w → 0+ y u → V-.

Modos TE, TM y LP:

, w ∈ [0, V ]
K 1 (w)
TE; t 0 (u, w) = −g 0 (w) = −
w.K 0 (w)

, w ∈ [0, V ], s ∈ (1,+∞ )
1 1 K 1 (w)
TM; t 0 (u, w, s) = 2 g 0 (w) = − 2
s s w.K 0 (w)
- K 0 (w) - K 0 (w ) −1
LP; m = 0, − g -1 ( w ) = = = 2 1
, w ∈ [0, V ]
w.K −1 ( w ) w.K 1 ( w ) w .g ( w )
m −1 - K m (w )
m ≥ 1, − g ( w ) =
w.K m −1 ( w )
Se puede demostrar que todas estas funciones valen -∞ cuando w = 0, tienen
signo negativo y son monótonas decrecientes en todo el rango de definición.

Para V = 10.
0 2 4 6 8 10

0.1

− K0 ( w( u) )
0.2
w( u) ⋅ K1 ( w( u) )

− K1 ( w( u) )
0.3
w( u) ⋅ K0 ( w( u) )

− Kn ( 3 , w( u) )
w( u) ⋅ Kn ( 2 , w( u) ) 0.4

0.5

0.6

46
Modos EH:
Cuando u → 0+ y w → V-:
s 2 + 1 K l +1 (V)
t l (0, V, s) = −
2.s 2 V.K l (V)
lim t l (0, V, s) = - ∞
V →0 +
Cuando w → 0+ y u → V-:
t l (V,0, s) = −∞

Representamos tl=3(u,w,s=1.001) para las frecuencias adimensionales V = 4,8 y 10.

0 2 4 6 8 10

0.2

t3 ( 10 , u) 0.4

t3 ( 8 , u)

t3 ( 4 , u)
0.6

0.8

Modos HE:
Cuando u → 0+ y w → V-:
s 2 + 1 K l −1 (V)
t l (0, V, s) =
2.s 2 V.K l (V)
 +∞ l =1
 2
lim t l (0, V, s) =  s + 1 1
V →0 + l≥2
 2.s 2 2.(l − 1)
Cuando w → 0+ y u → V-:
+∞ l =1
t (V,0, s) =
l
1
l≥2
(l − 1).(s 2 + 1)

47
 Representamos tl=2(u,w,s=1.001) para las frecuencias adimensionales V = 4,8 y 10.

0.35

0.3

t2 ( 10 , u)
0.25
t2 ( 8 , u)

t2 ( 4 , u) 0.2

0.15

0.1

0.05
0 2 4 6 8 10
u

Con estos datos basta con una simple inspección para comprobar la validez del criterio.

Para el caso de las frecuencias adimensionales de corte, utilizando las propiedades de

recurrencia y de la derivada de las funciones de Bessel de primera especie, tenemos que:
H´n(u) = H2n(u) – (2.n-1).Hn(u)/u + 1. Por tanto, cuando:
Hn(u) → 0, tengo que el cociente Hn(u)/H´n(u) → Hn(u) → 0.
Hn(u) → ±∞, tengo que el cociente Hn(u)/H´n(u) → 1/Hn(u) → 0.

Por lo que el criterio de convergencia también es válido para este caso.

Otro problema es verificar como influye la precisión en el cálculo de las funciones de

Bessel de primera especie en el resultado del N-R. El problema estaría en las zonas donde
la precisión de cálculo no es relativa, sino absoluta. Pero, ya vimos que , en ambas zonas el
cociente Hn(u)/H´n(u) → 0, por lo que el resultado obtenido tendrá la precisión exigida.

3.5 DERIVACIÓN NUMÉRICA

Sea una función f derivable en (a,b). Podemos calcular f(x) en el intervalo, no así su
derivada y queremos obtener f´(xo), con xo∈ (a,b). Denominaremos xc = (f(xo)/f´´(xo))1/2, a
la escala de curvatura de la función. El valor de xc nos indica la escala característica sobre
la que f varía del orden de sí misma.

Una familia de algoritmos se basan en la exploración del comportamiento de f en

escalas del orden de xc y en la suposición de cierta suavidad, de forma que los términos de

48
orden alto en la expansión de Taylor de f tengan significado. Estos algoritmos dan buena
precisión, pero con el coste de bastantes evaluaciones de f.

f(x + h) − f(x)
Otro método parte de la definición de derivada: f ′(x) = lim y utiliza el
h→0 h
acercamiento diferido al límite de Richardson, que consiste en extrapolar el valor del límite,
y por tanto, de la derivada, mediante el cálculo de diferencias finitas con valores de h cada
vez menores. Mediante el algoritmo de Neville, cada nuevo cálculo de diferencias finitas
nos permite obtener una extrapolación de orden superior y extrapolaciones de los ordenes
anteriores pero con h menor.

El algoritmo nos devuelve una estimación del error cometido, que intentaremos hacer
menor que una cierta precisión requerida. Hemos de suministrarle como valor de entrada
una estimación del parámetro xc. El error depende de lo bien que estimemos este parámetro.
Típicamente, el error se hace menor al ir creciendo xc, hasta que alcanzamos un punto a
partir del cual, la extrapolación da valores absurdos y, por tanto, errores muy grandes. Por
ello, para nuestra estimación, utilizaremos valores iniciales muy pequeños para ir
aumentándolos hasta alcanzar la precisión requerida, o bien, llegar a un punto en que
llegamos a un error grande, lo cual indica que no podemos alcanzar la precisión pedida.

Este algoritmo es relativamente lento, por lo que no nos va a interesar para calcular de
las derivadas de los puntos, sino sólo para obtener las condiciones de contorno para la
interpolación, (derivada en los extremos del intervalo de trabajo).

3.6 SISTEMA DE INTERPOLACIÓN

En primer lugar consideramos las aproximaciones de tipo polinómico, ya que es fácil
calcular sus derivadas y están muy documentadas. Al no estar prefijados los puntos de
interpolación, la mejor opción es la aproximación de Chebyshev, pues el polinomio
obtenido es similar al mínimax18, pero el esfuerzo para obtenerlo es mucho menor. Sin
embargo, al implementar dicho desarrollo, observamos que los resultados no eran
adecuados, con errores importantes y estos no disminuían apreciablemente al aumentar el
grado del polinomio. Como el comportamiento de las funciones a representar era bueno
dentro del intervalo de trabajo, pasamos a estudiar su comportamiento en el corte. De esta
forma, comprobamos que, para algunos modos, la derivada segunda de b(V,s) es infinita en
el corte, por lo que las aproximaciones polinómicas no sirven.
 w
2

   s =1
 V
Tenemos que b(V,s) =  2

( )
w
 1 + s − 1 .  − 1
2

 V
s >1
 s −1

18
Es, dentro de los polinomios de un mismo grado, aquel con la menor desviación máxima respecto a la
función de referencia en el intervalo de trabajo.

49
 Para s = 1, partimos de la expresión para s > 1, desarrollamos la raíz en series de Taylor
alrededor del origen, pues 0 ≤ w/V ≤ 1 y cuando s → 1, s2-1 → 0 y tomamos límite s → 1.

Derivando y tomando límite cuando la frecuencia adimensional V tiende a Vc, con Vc la

frecuencia adimensional de corte, llegamos a:
dw 
w 
d(V.b(V)) 
 = (s + 1) lim dV 
dV  V =Vc V → Vc V 

 
dw d w dw    ∀s ≥ 1
 w   
d (V.b(V))  dV   
= (s + 1) lim 2 dV + 
2
V  
2 V → Vc  dV  
dV  V = Vc V
 
 
Por las ecuaciones cerca del corte de la sección 3.7 “Expansiones asintóticas para la
ecuación de dispersión en la zona cercana al corte”, alcanzamos los siguientes resultados:

d(V.b(V))  d 2 (V.b(V)) 
 V 
dV  V = Vc dV 2  V = Vc
LP0,k 0 0
HE1,k 0 0
LP1,k 0 +∞
TE0,k 0 +∞
TM0,k 0 +∞
HE2,k 0 +∞
LPm,k m −1 (m − 1)(m + 2)
2 2
m≥2 m m2
2.( l + 1) + [1 + A ]
2
EHl,k A
(s + 1) (s + 1).A
l≥1 1+ A [1 + A ]3
HEl,k [
(s + 1). B + Vc2 ] E
l≥3 B − C + (1 − D).V c2

l s2 + 1
Donde: A= , B = (l-1)2.(s4-1), C = (l-1).(s2-1)2
2 s2
s2 + 1  (s 2 − 1) 2   s2 −1
2
 
D=   2 +  2  .( l − 2)  − (s 2 + 1) 2 
8.s 2 .( l − 2)  l   s +1  
E = E(s,l,Vc); 0 < E < +∞. No nos interesa conocer el valor de E, nos basta
con saber que es un número finito.

Entonces consideramos dos opciones: aproximaciones racionales o splines. En el primer

caso, no teníamos garantizado que el método fuera válido, ni una acotación del error en la
interpolación. El método de los splines es válido en nuestro caso, facilita el cálculo de las

50
derivadas y nos da cotas del error cometido en la interpolación. Optamos por los splines
cúbicos pues se calculan de forma eficiente19 y están bien documentados.

Una vez decidido el tipo de interpolación se nos plantea como disponemos los puntos
de interpolación. Esto es importante, ya que se puede obtener una precisión dada con el
mínimo de puntos y manejar un sistema menor y de más rápida resolución. El problema es
que para elegir adecuadamente la disposición hay que tener información de cómo se porta
la función en el rango de trabajo. Para ello utilizaríamos unos pocos puntos equiespaciados,
para tener idea del comportamiento. Con esta información hay un algoritmo que halla la
disposición óptima. Para ello son necesarios tres pasos: interpolación burda, calculo de la
disposición óptima de puntos y por último la interpolación fina con el mínimo número de
puntos. Otra opción es utilizar una partición fina de puntos equiespaciados.

Optamos por la segunda opción, pues en nuestro caso, es más rápido y sencillo de esta
forma. El número de puntos viene dictado por la cota del error de interpolación, con la
salvedad de los modos cuya derivada segunda en el corte es infinita (LP1,k, TE0,k, TM0,k y
HE2,k), o bien vale 0 (LP0,k y HE1,k) y en el modo HE3,1. Para estos, en la zona cercana al
corte, usaremos la máxima densidad de puntos, sea cual sea la precisión de trabajo.

También hemos de imponer condiciones de contorno en los dos extremos. Esto es

importante pues el valor del error de interpolación depende de ello. Las condiciones pueden
ser de tres tipos y el tipo de condición en un extremo no depende de la del otro:

1) Imponer el valor de la derivada primera.

2) Imponer el valor de la derivada segunda.
3) Condición “sin nudo”. Esta condición es útil cuando no tenemos información que
nos permita imponer una de las dos anteriores. Básicamente se trata de imponer que
en el segundo nodo (condición en el extremo izquierdo), o en el penúltimo nodo
(condición en el extremo derecho), la derivada tercera sea igual por la derecha y por
la izquierda. La desventaja es que el ajuste no es tan bueno como en los otros casos
y, por tanto, la cota de error es mucho peor. No la utilizaremos.

Usaremos la opción 1 para el cálculo de la derivada primera, y la 2 cuando

representemos la derivada segunda ya que es la opción que mejor resultados dá.

El error teórico de interpolación al imponer condiciones del tipo 1 en ambos extremos

viene dado por: ek(h) = K.h4-k, con ek(h) el error cometido al interpolar la derivada k-esima,
k = 0,1,2,3, con una distancia entre puntos h. La constante depende del tipo de función. En
nuestro caso ha sido elegida heurísticamente. La formula anterior nos permite calcular el
número de puntos equiespaciados necesarios para obtener el error exigido en la

19
La interpolación polinómica sobre n puntos implica resolver un sistema cuadrado de orden n y el algoritmo
tendrá un orden de complejidad O(n2), pero en los splines cúbicos cada punto está acoplado sólo con el
intervalo que lo contiene y los dos vecinos, lo que nos lleva a un sistema cuadrado tridiagonal, con diagonal
dominante y el algoritmo pasa a tener una complejidad de O(n). Para simplificar el algoritmo, hemos tenido
en cuenta que el sistema es de diagonal dominante, o sea, |ai,i| > |ai-1,i| + |ai,i+1|, ∀i= 1,2,...,n. Esta propiedad
garantiza que, en la sustitución hacia delante, no aparecerá ningún pivote nulo y que no habrá que pivotear.

51
interpolación de la derivada k-esima.20. Hemos usado esta fórmula incluso para el caso de
usar condiciones del tipo 2 en ambos extremos, ya que no disponíamos de más información
y estamos en el lado de la seguridad.

3.7 REPRESENTACIÓN GRÁFICA

Para la presentación de las curvas utilizamos el componente MSChart, debido a su
comodidad de uso. Los datos a representar se le pasan como una matriz de 2m columnas y
n filas, donde m es el número de curvas que el usuario quiere representar y n el número de
puntos de cada curva. El componente traza cada curva uniendo mediante rectas los n
puntos. La rapidez de presentación y actualización del componente depende en gran manera
del tamaño de la matriz. Como m viene fijado por el usuario, sólo podemos actuar sobre el
número de puntos n.

Partimos con n = 2, la colección de puntos sólo contiene los extremos del intervalo de
trabajo, comprobamos si la diferencia, en el punto medio, entre el valor real y el valor de la
recta que une los extremos es menor que la resolución gráfica;

Si no es menor, insertamos dentro de la colección el punto medio, colocándolo entre los

extremos del intervalo e incrementamos el valor de n. Ahora pasamos a verificar el
subintervalo izquierdo entre el extremo izquierdo del intervalo anterior y el punto insertado.

Si, por el contrario, es menor, entonces pasamos a verificar el siguiente subintervalo, el

comprendido entre el último punto introducido en la colección y el extremo derecho del
intervalo actual.

Así sucesivamente hasta cubrir todo el rango de trabajo. Como último detalle hemos de
resaltar que este método funciona bien excepto si hay un punto de inflexión dentro del
rango de trabajo. Para evitar este problema, en lugar de empezar sólo con dos puntos, en
realidad partimos de una colección que contiene a los extremos y a varios puntos
equiespaciados del interior. La cantidad de puntos iniciales en la colección es función de la
anchura del rango de trabajo.

En la gráfica se muestra que la idea de medir el

error en el punto medio sólo puede fallar en los puntos
de inflexión, donde el error puede ser cero pero la línea
de la gráfica puede tener errores grandes dentro del
rango de trabajo.

Con este método nos ahorramos hasta un 50% de puntos para la misma calidad de
representación y el componente MSChart se actualiza y redibuja muchísimo más rápido,
ahorrando mucho más tiempo del que se tarda en calcular los puntos de la representación.
20
En estas circunstancias, la máxima densidad de puntos corresponde a la máxima precisión alcanzable por el
programa.

Registrado por Microsoft Corporation.

52
3.8 REFERENCIA BIBLIOGRÁFICA
Adjuntamos una breve reseña bibliográfica para profundizar en las consultas, indicando
los temas que se han consultado en cada libro.

- Handbook of Mathematical Functions. With formulas, graphs and mathematical tables.

Editores: Milton Abramowitz, Irene A. Stegun. Recopilación de varios autores.
Dover Publications, New York
• Funciones de Bessel. (Autor F. W. J. Olver).

- A treatise on the theory of Bessel functions.

G. N. Watson
Cambridge University Press
• Funciones de Bessel

- Numerical Recipes in C. The Art of Scientific Computing. Second edition.

William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian [Link]
Cambridge University Press
• Funciones de Bessel
• Búsqueda de raíces de ecuaciones algebraicas
• Derivación numérica
• Interpolación
• Splines cúbicos

- A practical guide to splines.

Carl de Boor
Springer-Verlag
• Interpolación
• Splines

- Introduction to numerical analysis. Second edition.

J. Stoer, R. Bulirsch
Springer-Verlag
• Interpolación
• Splines

- ELF and GNOME: Two tiny codes to evaluate the real zeros of the Bessel functions of
the firs kind for real orders.
J. Segura, A. Gil
Computer Physics Communications 117 (1999), 250-260.
• Obtención de los ceros de las funciones de Bessel de primera clase Js(x) = 0

53