 |
VOLUMEN 34, NÚMERO 2 | JULIO-DICIEMBRE 2022 | PP. 97-108
ISSN: 2250-6101
DOI: https://doi.org/10.55767/2451.6007.v34.n2.39487
Resolución de ecuaciones trascendentes de la física mediante el método de Chebyshev
Solving transcendental equations of physics with Chebyshev's method
1 Departamento de Física y Química, IES José Vilaplana, Avinguda Gil d'Atrosillo, 26, 12500 Vinaròs, Castelló, España.
*E-mail: d.marquesvillarroy@edu.gva.es
Recibido el 9 de agosto de 2022 | Aceptado el 20 de octubre de 2022
En este trabajo se proporciona a los docentes de física un método optimizado para resolver ecuaciones trascendentes que no pueden resolverse de manera algebraica y aparecen muy a menudo en cursos de enseñanza superior de física e ingeniería. El método se basa en una interpolación con polinomios de Chebyshev y está optimizado en tiempo computacional, facilidad de uso y precisión. El método se ha aplicado en problemas particulares de la física donde aparecen ecuaciones trascendentes como la compresión de un resorte real; la ecuación de un diodo; la solución de la ecuación de Schrödinger en un pozo de potencial; y el cálculo de números de onda de corte de un cable coaxial. Se compara el método con otros de la bibliografía para comprobar su correcto funcionamiento y las mejoras que presenta. También se proporcionan los códigos fuente de MATLAB para implementar el método y los ejemplos particulares.
In this work, physics teachers are provided with an optimized method to solve transcendental equations that cannot be solved alge-braically and appear very often in higher education physics and engineering courses. The method is based on an interpolation with Chebyshev’s polynomials, and it is optimized for computational time, manageability, and accuracy. The method has been applied in specific problems of physics where transcendental equations appear, such as the compression of a real spring; the equation of a diode; the solution of the Schrödinger equation in a potential well; and the computing of cutoff wave numbers of a coaxial wire. The method is compared with others of the literature to check its correct behavior and the improvements that it presents. MATLAB source codes to implement the method and particular examples are also provided.
Las ecuaciones trascendentes son aquellas en las que, al menos un término de la igualdad es una función trascendente, que se caracteriza por no satisfacer una ecuación polinómica, por ejemplo, las funciones logarítmicas, exponenciales o trigonométricas. En el campo de la física aplicada aparecen muchas ecuaciones de este tipo a la hora de resolver problemas particulares. Este tipo de ecuaciones no se pueden resolver de forma analítica, ya que no se puede aislar la incógnita de las mismas y, por lo tanto, deben resolverse de manera numérica (Boyd, 2014b).
Para resolver este tipo de ecuaciones se puede reescribir la ecuación trascendente en forma de f(x)=0, y encontrar las raíces de esta función, que coincidirán con las de la ecuación original.
www.revistas.unc.edu.ar/index.php/revistaEF
Existen muchos métodos para resolver ecuaciones no lineales expresadas de la forma f(x)=0. Una manera es mediante métodos iterativos como los de Newton, Traub, Ostrowski u otros de orden mayor como se presentan en Quinga (2021). No obstante, estos métodos requieren de una estimación inicial lo suficientemente próxima a la solución para converger correctamente y son bastante sensibles a funciones abruptas o con valores extremos (valores muy grandes o pequeños de f(x) o de x). Además, estos métodos están diseñados para encontrar una sola raíz, lo cual puede ser una limitación si el problema a estudiar admite varias soluciones, como veremos en este artículo.
Otra manera de resolver ecuaciones del tipo f(x)=0 es mediante una aproximación polinómica (Boyd, 2014b). De esta manera, la función f(x) se aproxima mediante un polinomio y se calculan las raíces del mismo (los métodos actuales de resolución de polinomios son muy eficientes, tanto en precisión como en velocidad computacional). Para que este método funcione correctamente, la aproximación polinómica de la función debe ser buena, por lo que el orden del polinomio será habitualmente alto, lo cual no supone un problema sustancial a la hora de encontrar raíces, ya que los métodos de resolución polinómica están muy optimizados hoy en día (aunque a mayor orden, mayor precisión, pero menor velocidad computacional). Este método tiene la ventaja de encontrar todas las raíces de un cierto intervalo en un solo paso a diferencia de los métodos iterativos, siempre y cuando la función esté bien aproximada por un polinomio en dicho intervalo (Boyd, 2014b).
En este artículo nos centraremos en el segundo método de aproximación de funciones trascendentes mediante polinomios, concretamente polinomios de Chebyshev, ya que aproximan muy fielmente funciones de todo tipo con un orden no excesivamente alto. A continuación, se muestra el método de Chebyshev para efectuar la búsqueda de raíces y se utiliza para resolver problemas reales de física aplicada. Hemos utilizado funciones trascendentes de todo tipo para que se pueda ver la versatilidad y utilidad del método.
Tanto el método de Chebyshev como la comprobación experimental del mismo han sido implementados y desarrollados con el entorno de programación MATLAB versión 7.10.0 en un ordenador con un procesador Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz 2.59 GHz y una RAM de 16 GB.
El método de Chebyshev presentado en Boyd (2002; 2014a) para obtener las raíces de una función f(x) consiste, básicamente, en aproximar la función en un intervalo [a,b], donde la función ha de ser continua, por un polinomio de Chebyshev, y obtener las raíces del mismo, que coinciden con las raíces de la función definida sobre ese intervalo.
El primer paso del método consiste en aproximar la función f(x) mediante un polinomio de Chebyshev, obteniendo los coeficientes del mismo aj. Una vez obtenidos los coeficientes, se calcula la Companion Matrix Ajk, cuyos valores propios son las raíces de la función f(x). Este procedimiento se puede implementar mediante los siguientes pasos:
Crear puntos de interpolación la función f(x) mediante un mallado Chebyshev-Lobato en el intervalo [a,b].
𝑥 = 𝑏−𝑎
𝜋𝑘
𝑏+𝑎
cos ( ) +
𝑘
2 𝑁
, 𝑘 = 0,1,2,3 … 𝑁 (1)
2
donde N es el grado del polinomio por el que se quiere aproximar la función.
Evaluar la función en los puntos del mallado.
𝑓𝑘 ≡ 𝑓(𝑥𝑘), 𝑘 = 0,1,2,3 … 𝑁 (2)
Calcular matriz de interpolación (N+1)x(N+1).
𝜏𝑗𝑘
= 2
𝑝𝑗𝑝𝑘𝑁
cos (𝑗𝜋𝑘) (3)
𝑁
donde, 𝑝𝑗 = {
2, 𝑗 = 0; 𝑗 = 𝑁
1, 1 ≤ 𝑗 ≤ 𝑁 − 1
Calcular los coeficientes del polinomio.
𝑎𝑗
𝑁
= ∑
𝑘=0
𝜏𝑗𝑘𝑓𝑘 , 𝑗 = 0,1,2, … 𝑁 (4)
Calcular la Companion Matrix.
𝛿2,𝑘, 𝑗 = 1; 𝑘 = 1,2, … 𝑁
1 (𝛿
+ 𝛿
), 𝑗 = 2, … , (𝑁 − 1); 𝑘 = 1,2, … , 𝑁
𝐴𝑗𝑘 =
2 𝑗,𝑘+1
𝑗,𝑘−1
(5)
(−1)𝑎𝑗−1 1
𝑁
{ 2𝑎 + 2 𝛿𝑘,𝑁−1, 𝑗 = 𝑁; 𝑘 = 1,2,3, … 𝑁
Calcular los valores propios de la Companion Matrix, que coinciden con las raíces del polinomio en el intervalo [a,b]
𝑟𝑜𝑜𝑡𝑠(𝑓(𝑥)|𝑥=[𝑎,𝑏]) = 𝑒𝑖𝑔(𝐴𝑗𝑘) (5)
La implementación de este método está programada en MATLAB y es de código abierto. Se puede obtener en Morris (2007). La llamada a este programa es:
raices=FindRealRoots(funhandle,a,b,n)
Se deben pasar como parámetros de entrada la función a evaluar (fundhandle), el intervalo [a,b] de estudio y el orden del polinomio n. Las funciones a evaluar se muestran en los siguientes apartados; el intervalo de estudio en cada caso se determinará graficando la función trascendente y observado las cercanías de la intersección de esta función con el eje de abscisas; y finalmente, el orden del polinomio, se recomienda utilizar un orden mínimo de n=15, ya que con este orden suele aproximar cualquier función bastante bien. No obstante, si se aumenta el orden aumentará la precisión del método (y de las raíces) aunque empeorará la velocidad computacional. Se trata, en última instancia, de encontrar un punto de equilibrio entre velocidad y precisión. Para cada problema se puede hallar por prueba-error.
A continuación, se muestran algunas ecuaciones trascendentes de la física aplicada en varios campos de estudio: mecánica, electrónica, cuántica y electromagnetismo. No se realiza una deducción teórica de esas ecuaciones, ya que el artículo se basa en su resolución. No obstante, se referencia el origen de todas ellas para que el lector pueda consultarlo en caso de necesidad.
El primer paso consiste en definir las funciones trascendentes a partir de las ecuaciones físicas de forma que tengamos funciones del tipo f(x)=0. A continuación, para obtener una primera estimación de las raíces, un procedimiento bastante sencillo es graficar la función y observar los cortes con el eje de abscisas. De esta manera, determinaremos gráficamente el intervalo de estudio, que en ningún caso ha de ser tan preciso como las estimaciones iniciales de los métodos iterativos, es suficiente con determinar un intervalo donde se encuentre la raíz.
Seguidamente aplicaremos el método de Chebyshev descrito anteriormente para calcular las raíces de cada ecuación en el intervalo de estudio. Finalmente, mostraremos el tiempo computacional que emplea el método en nuestro sistema informático y la precisión que ofrece, la que determinaremos evaluando las raíces en la propia función, de manera que nos muestre cuan alejado está la solución obtenida del cero ideal de una raíz, mostrando así el error absoluto que presenta. En los anexos se proporciona el código de MATLAB empleado para obtener los resultados mostrados en las siguientes subsecciones.
La fuerza elástica de los resortes reales no presenta una dependencia lineal con la elongación, por lo que al aplicar las leyes de Newton obtenemos un sistema de ecuaciones no lineal.
FIGURA 1. Compresión de un resorte: al caer sobre este un objeto desde una altura h se comprime una distancia d (Quinga, 2021).
Por ejemplo, en Quinga (2021) se obtiene la siguiente ecuación trascendente de este problema (no se pude considerar una ecuación polinómica porque tiene exponentes fraccionarios):
2
1
𝑚𝑔ℎ + 𝑚𝑔𝑑 = 2 𝑘1𝑑
2
+ 5 𝑘2𝑑
5/2
(6)
Por lo tanto, definimos la siguiente función trascendente f(d) cuyas raíces se obtendrán igualando f(d)=0:
𝑓(𝑑)
1
2
= 𝑚𝑔ℎ + 𝑚𝑔𝑑 − 2 𝑘1𝑑
2
− 5 𝑘2𝑑
5/2
(7)
Donde m es la masa del objeto; g es la gravedad; k1 y k2 son constantes de proporcionalidad; h es la longitud en reposo del muelle; y d es la elongación que sufre el muelle (la incógnita). En este problema en particular tomaremos
los siguientes valores: 𝑚 = 95 𝑔; 𝑔 = 9.8𝑚/𝑠2 ; 𝑘1
= 40000𝑔/𝑠2; 𝑘2
𝑔
= 40
𝑠2𝑚0.5; ℎ = 0.43𝑚.
Si representamos la función trascendente (7) sobre un intervalo adecuado en función de d, obtenemos la siguiente
gráfica:
1
2
= 𝑚𝑔ℎ + 𝑚𝑔𝑑 − 2 𝑘1𝑑
2
− 5 𝑘2𝑑
5/2
para los valo-
res: 𝑚 = 95 𝑔; 𝑔 = 9.8𝑚/𝑠2 ; 𝑘1
= 40000𝑔/𝑠2; 𝑘2
𝑔
= 40
𝑠2𝑚0.5; ℎ = 0.43𝑚.
De la imagen se deduce que la raíz estará en torno a 0.16 m. Elegimos, por lo tanto, el intervalo para aplicar el método entre 0 y 0.3 m; y el orden del polinomio de Chebyshev n=15 (si quisiéramos más precisión solo hay que aumentar este parámetro para que la aproximación sea mejor y, por tanto, las raíces también). En el anexo I se muestran las líneas de código de MATLAB para resolver este problema mediante el método propuesto de Chebyshev. Con estos parámetros, se obtiene un resultado de d=0.1667 m, que concuerda con Quinga (2021).
El tiempo empleado para encontrar la raíz de la ecuación trascendente en nuestro sistema informático ha sido de 0.0006 segundos. Si utilizamos un método iterativo convencional, como el de Traub propuesto en Quinga (2021), que en teoría es el más rápido de los métodos iterativos que se propone en ese artículo, con una estimación inicial de x0=0.16 m (prácticamente la estimación inicial ya es la solución) tarda 0.0013 segundos en encontrar la raíz con la misma precisión que Chebyshev, aproximadamente el doble de tiempo que el método propuesto en este trabajo.
Al evaluar el resultado (d) en la función (f(d)) se obtiene un valor del orden de magnitud de 10-15, como la función trascendente toma valores del orden de 102 al ser evaluada en el intervalo de estudio (se puede ver en la gráfica), el error cometido es muy bajo.
La tensión en los bornes de un diodo como el mostrado en la figura se puede obtener aplicando la ley de Kirchoff (𝑉𝐵𝐵 = 𝑉𝐷 + 𝑉𝑅). A su vez, la corriente que atraviesa el diodo (y por ende, la resistencia) se determina por la ley de
𝑉𝐷
Shockley (𝐼 = 𝐼𝑠 (𝑒𝑉𝑇 − 1)) y se llega a la siguiente ecuación trascendente como demuestran Sze y Kwok (2007):
𝑉𝐷
𝑉𝐵𝐵 = 𝑉𝐷 + 𝐼𝑠 (𝑒𝑉𝑇 − 1) 𝑅𝐿 (8)
Por lo tanto, definimos la siguiente función trascendente:
𝑉𝐷
𝑓(𝑉𝐷 ) = 𝑉𝐵𝐵 − 𝑉𝐷 − 𝐼𝑠 (𝑒𝑉𝑇 − 1) 𝑅𝐿 (9)
Donde 𝑉𝐵𝐵 es la tensión del generador; 𝑅𝐿 es la resistencia de carga; 𝐼𝑠 es la corriente de saturación inversa; 𝑉𝑇 es el equivalente en tensión de la temperatura; y 𝑉𝐷 es la tensión en los bornes del diodo (la incógnita). Este problema en concreto tomaremos los siguientes valores: 𝑉𝐵𝐵 = 1𝑉; 𝐼𝑠 = 1.2 · 10−11𝐴; 𝑉𝑇 = 26𝑚𝑉; 𝑅𝐿 = 4 Ω.
Si representamos la función trascendente (9) sobre un intervalo adecuado en función de 𝑉𝐷 obtenemos la siguiente gráfica:
𝑉𝐷
FIGURA 4. Representación gráfica de la función trascendente del diodo 𝑓(𝑉𝐷) = 𝑉𝐵𝐵 − 𝑉𝐷 − 𝐼𝑠 (𝑒 𝑉𝑇 − 1) 𝑅𝐿 con los valores 𝑉𝐵𝐵 =
1𝑉; 𝐼𝑠 = 1.2 · 10−11𝐴; 𝑉𝑇 = 26𝑚𝑉; 𝑅𝐿 = 4 Ω.
De la imagen se deduce que la raíz estará en torno a 0.59 V. Elegimos, por lo tanto, el intervalo para aplicar el método entre 0.5 y 0.65 V; y el orden del polinomio de Chebyshev n=15. En el anexo II se muestran las líneas de código de MATLAB para resolver este problema. Con estos parámetros, obtenemos un resultado de 𝑉𝐷 = 0.594 𝑉, que concuerda con Aranzabal Olea (2001).
El tiempo empleado para encontrar la raíz de la ecuación trascendente en nuestro sistema informático ha sido de 0.0005 segundos. Si utilizamos el método de Traub propuesto en Quinga (2021), que en teoría es el más rápido de los métodos iterativos que propone, con una estimación inicial de x0=0.55 V tarda 0.0014 segundos en encontrar la raíz con la misma precisión que Chebyshev, aproximadamente el doble de tiempo que el método propuesto en este trabajo, como en el caso anterior (resorte real).
Al evaluar el resultado (𝑉𝐷 ) en la función (f(𝑉𝐷 )) se obtiene un valor del orden de magnitud de 10-11, como la función trascendente toma valores del orden de 101 al ser evaluada en el intervalo de estudio (se puede ver en la gráfica), el error cometido es muy bajo.
La ecuación Schrödinger proporciona, entre otras cosas, los niveles cuantizados de energía de una partícula en una cierta región. Una de las regiones más estudiadas es el pozo de potencial finito (Galindo y Pascual, 1989; Bohm, 1989).
Del estudio de esta estructura se extraen dos ecuaciones trascendentes que nos dan los niveles de energía de la partícula en el pozo de potencial para el estado simétrico (paridad par) y antisimétrico (paridad impar):
𝐶𝑎𝑠𝑜 𝑃𝑎𝑟: √2𝑚𝐸 sin (√2𝑚𝐸 𝑎) = √2𝑚(|𝑉0|−𝐸) cos (√2𝑚𝐸 𝑎) (10)
ћ ћ ћ ћ
𝐶𝑎𝑠𝑜 𝐼𝑚𝑝𝑎𝑟: √2𝑚𝐸 cos (√2𝑚𝐸 𝑎) = − √2𝑚(|𝑉0|−𝐸) sin (√2𝑚𝐸 𝑎) (11)
ћ ћ ћ ћ
Por lo tanto, definimos las siguientes funciones trascendentes:
𝑓 (𝐸) = √2𝑚𝐸 sin (√2𝑚𝐸 𝑎) − √2𝑚(|𝑉0|−𝐸) cos (√2𝑚𝐸 𝑎) (12)
𝑝𝑎𝑟 ћ ћ ћ ћ
𝑓 (𝐸) = √2𝑚𝐸 cos (√2𝑚𝐸 𝑎) + √2𝑚(|𝑉0|−𝐸) sin (√2𝑚𝐸 𝑎) (13)
𝑖𝑚𝑝𝑎𝑟 ћ ћ ћ ћ
Donde m es la masa de la partícula; ћ es constante de Planck reducida; 2a es la anchura del pozo; V0 es la profundidad del pozo de potencial; y E es la energía cuantizada de la partícula (la incógnita). En este problema en particular
tomaremos los siguientes valores: 𝑚 = 1.6 · 10−31 𝑘𝑔; ћ = ℎ
= 6,63·10−34; 𝑎 = 0.2 𝑛𝑚;;𝑉
= 75𝑒𝑉.
2𝜋
2𝜋 0
Si representamos las funciones (12) y (13) sobre un intervalo adecuado en función de E, obtenemos la siguiente
gráfica:
FIGURA 6. Representación gráfica de las ecuaciones trascendentes del pozo de potencial 𝑓 (𝐸) = √2𝑚𝐸 sin (√2𝑚𝐸 𝑎) −
𝑝𝑎𝑟 ћ ћ
√2𝑚(|𝑉0|−𝐸) cos (√2𝑚𝐸 𝑎) en línea continua y 𝑓 (𝐸) = √2𝑚𝐸 cos (√2𝑚𝐸 𝑎) + √2𝑚(|𝑉0|−𝐸) sin (√2𝑚𝐸 𝑎) en línea discontinua para
ћ ћ
−31
𝑖𝑚𝑝𝑎𝑟
ℎ
ћ ћ ћ ћ
6.63·10−34
los valores numéricos 𝑚 = 1.6 · 10
𝑘𝑔; ћ = =
2𝜋
2𝜋 ; 𝑎 = 0.2 𝑛𝑚;𝑉0 = 75𝑒𝑉.
Cabe destacar en este caso que las funciones trascendentes con las que trabajamos (𝑓𝑝𝑎𝑟(𝐸) y 𝑓𝑖𝑚𝑝𝑎𝑟(𝐸)) adquieren valores muy elevados (del orden de 1010) para unos valores de energía (E) muy bajos (del orden de 10-17). Esto es debido a que se han introducido los datos en el sistema internacional y la energía está expresada en julios (V0 también). En este tipo de problemas es habitual utilizar electronvoltios como unidad de energía debido a los niveles con los que se trabaja, pero en este caso se han mantenido las unidades en julios para mostrar que el método funciona también en condiciones extremas.
De la imagen se deduce que hay 3 raíces pares y 3 impares en el intervalo mostrado. Elegimos, por lo tanto, el intervalo para aplicar el método entre 0 y 1.2 ·10-17 J; el orden del polinomio de Chebyshev podemos seguir tomando n=15. En el anexo III se muestran las líneas de código de MATLAB para resolver este problema mediante el método propuesto de Chebyshev. Con estos parámetros, obtenemos los resultados mostrados en la tabla I, que concuerdan con Boscá (2014).
(𝐸) = √2𝑚𝐸 sin (√2𝑚𝐸 𝑎) − √2𝑚(|𝑉0|−𝐸) cos (√2𝑚𝐸 𝑎) y
( ) √2𝑚𝐸
√2𝑚𝐸
√2𝑚(|𝑉0|−𝐸)
√2𝑚𝐸
𝑝𝑎𝑟 ћ
ћ ћ
−31
ћ
ℎ 6.63·10−34
𝑓𝑖𝑚𝑝𝑎𝑟 𝐸 =
cos (
ћ
𝑎) +
ћ
sin (
ћ
𝑎) para los valores numéricos 𝑚 = 1.6 · 10
ћ
𝑘𝑔; ћ = =
2𝜋
2𝜋 ;
𝑎 = 0.2 𝑛𝑚;;𝑉0 = 75𝑒𝑉. Estas raíces se corresponden con las energías cuantizadas para el caso par e impar.
Energía en Julios (x 10-17) | Energía en eV | |
0.0304 | 1.9008 | |
Caso Par | 0.2716 | 16.9772 |
0.7384 | 46.1503 | |
0.121 | 7.5367 | |
Caso Impar | 0.479 | 29.9395 |
1.036 | 64.7201 | |
El tiempo computacional invertido por nuestro sistema informático ha sido de 0.0005 segundo y 0.0008 segundos para el caso par e impar respectivamente. Es importante remarcar que este tiempo es el que ha tardado en encontrar las 3 raíces simultáneamente en cada caso, no en encontrar cada una de ellas, ya que el método encuentra simultáneamente todas las raíces que existen en el intervalo de estudio. Con métodos iterativos como los propuestos en Quinga (2021) deberíamos aplicar el método para cada una de las raíces y el tiempo empleado sería mucho mayor que el de Chebyshev. Además, para valores tan extremos de la función a analizar f(x) y de la variable independiente x, los métodos iterativos como los de Quinga (2021) no suelen converger fácilmente y su comportamiento es inestable, como es el caso de este problema, lo que hace al método de Chebyshev óptimo para resolver ecuaciones trascendentes con solución múltiple y aquellas que tomen valores numéricos extremos.
Al evaluar el resultado (E) en las funciones (12) y (13) (𝑓𝑝𝑎𝑟(𝐸) y 𝑓𝑖𝑚𝑝𝑎𝑟(𝐸)) se obtienen valores de un orden de magnitud de 106, como la función trascendente toma valores del orden de 1010 al ser evaluada en el intervalo de estudio (se puede ver en la gráfica), el error cometido es de 4 órdenes de magnitud más pequeño, lo cual podemos considerar aceptable.
En un cable coaxial se pueden propagar muchos modos electromagnéticos además del fundamental (que es el tras-versal electromagnético, TEM), como los transversales magnéticos (TMmn) y los transversales eléctricos (TEmn). Para saber si un modo se propaga por un cable se debe calcular su número de onda de corte o frecuencia de corte, que es la mínima frecuencia a la que debe operar el cable para que este modo se propague.
FIGURA 7. Estructura geométrica de un cable coaxial con un conductor interior de radio a y uno externo de radio b.
El estudio electromagnético de estas estructuras no es trivial, ya que, al tratarse de estructuras cilíndricas involucra funciones de Bessel de primera (𝐽𝑚(𝑥)) y segunda especie (𝑌𝑚 (𝑥)). En Pozar (1998) se puede encontrar un estudio detallado de este problema, donde se obtienen dos ecuaciones trascendentes cuyas soluciones proporcionan los números de onda de corte de los modos TM y TE respectivamente. Estas ecuaciones son:
𝐶𝑎𝑠𝑜 𝑇𝑀: 𝐽𝑚(𝑘𝑐𝑛𝑎)𝑌𝑚(𝑘𝑐𝑛𝑏) = 𝐽𝑚(𝑘𝑐𝑛𝑏)𝑌𝑚 (𝑘𝑐𝑛𝑎) (14)
𝑚
𝐶𝑎𝑠𝑜 𝑇𝐸: 𝐽′ (𝑘
𝑐𝑛
𝑎)𝑌′ (𝑘
𝑐𝑛
𝑏) = 𝐽′ (𝑘
𝑐𝑛
𝑏)𝑌′ (𝑘
𝑐𝑛
𝑎) (15)
𝑚
𝑚
𝑚
Por lo tanto, definimos las siguientes funciones trascendentes:
𝑓𝑇𝑀(𝑘𝑐𝑛) = 𝐽𝑚(𝑘𝑐𝑛𝑎)𝑌𝑚(𝑘𝑐𝑛𝑏) − 𝐽𝑚(𝑘𝑐𝑛𝑏)𝑌𝑚 (𝑘𝑐𝑛𝑎) (16)
𝑚
𝑚
𝑚
𝑚
𝑓𝑇𝐸 (𝑘𝑐𝑛) = 𝐽′ (𝑘𝑐𝑛𝑎)𝑌′ (𝑘𝑐𝑛𝑏) − 𝐽′ (𝑘𝑐𝑛𝑏)𝑌′ (𝑘𝑐𝑛𝑎) (17)
Donde a es el radio menor del cable coaxial y b el mayor; los subíndices m y n definen el modo electromagnético que se propaga por el cable; 𝑘𝑐𝑛 es el número de onda de corte del modo enésimo (la incógnita). En este problema en concreto tomaremos como valores: 𝑎 = 1 𝑚𝑚; 𝑏 = 8 𝑚𝑚; m = 0.
Si representamos las funciones (16) y (17) sobre un intervalo adecuado en función de 𝑘𝑐𝑛, obtenemos la siguiente gráfica:
FIGURA 8. Representación de las ecuaciones trascendentes de un cable coaxial 𝑓𝑇𝑀(𝑘𝑐𝑛) = 𝐽𝑚(𝑘𝑐𝑛𝑎)𝑌𝑚(𝑘𝑐𝑛𝑏) −
𝐽𝑚(𝑘𝑐𝑛𝑏)𝑌𝑚(𝑘𝑐𝑛𝑎) en línea continua y 𝑓𝑇𝐸(𝑘𝑐𝑛) = 𝐽′ (𝑘𝑐𝑛𝑎)𝑌′ (𝑘𝑐𝑛𝑏) − 𝐽′ (𝑘𝑐𝑛𝑏)𝑌′ (𝑘𝑐𝑛𝑎) en línea discontinua para los valores
numéricos : 𝑎 = 1 𝑚𝑚; 𝑏 = 8 𝑚𝑚; m=0.
𝑚 𝑚
𝑚 𝑚
Vemos que en este intervalo hay cuatro raíces para los modos TM y 4 para los modos TE, por lo tanto, el intervalo de estudio será de 0 a 2000. El orden del polinomio podemos seguir utilizando n=15. En el anexo IV se muestran las líneas de código de MATLAB para resolver este problema mediante el método propuesto de Chebyshev.
TABLA II. Raíces las ecuaciones trascendentes de un cable coaxial 𝑓𝑇𝑀(𝑘𝑐𝑛) = 𝐽𝑚(𝑘𝑐𝑛𝑎)𝑌𝑚(𝑘𝑐𝑛𝑏) − 𝐽𝑚(𝑘𝑐𝑛𝑏)𝑌𝑚(𝑘𝑐𝑛𝑎) y
𝑚
𝑚
𝑚
𝑚
𝑓𝑇𝐸(𝑘𝑐𝑛) = 𝐽′ (𝑘𝑐𝑛𝑎)𝑌′ (𝑘𝑐𝑛𝑏) − 𝐽′ (𝑘𝑐𝑛𝑏)𝑌′ (𝑘𝑐𝑛𝑎) para los valores numéricos : 𝑎 = 1 𝑚𝑚; 𝑏 = 8 𝑚𝑚; m=0. Las raíces se co-
rresponden con los números de onda de corte (kcn) de los modos TM0n y TE0n. Se evalúa la raíz en la función f para observar la precisión, siendo f la función 𝑓𝑇𝑀 o 𝑓𝑇𝐸 según corresponda el modo.
Modo | kcn con Chebyshev | f(kcn) | kcn con (Peñaranda Foix, 2010) | f(kcn) |
TM01 | 429,394 | 0,17·10-6 | 429,567 | -5,6·10-4 |
TM02 | 884,533 | 2,64·10-6 | 884,88 | 5,89·10-4 |
TM03 | 1336,632 | 4,16·10-6 | 1337,15 | -5,9·10-4 |
TM04 | 1787,437 | -5,2·10-6 | 1788,24 | 6.89·10-4 |
TE01 | 499,822 | 0,48·10-6 | 499,823 | -6·10-6 |
TE02 | 935,130 | 2,31·10-6 | 935,137 | 1,5·10-5 |
TE03 | 1375,188 | 2,43·10-6 | 1375,23 | -5·10-5 |
TE04 | 1818,285 | -6,36·10-6 | 1818,38 | -6,5·10-6 |
Con los valores numéricos mostrados anteriormente, se obtienen los resultados mostrados en la tabla II, que concuerdan con los de la aplicación elaborada por Peñaranda Foix (2010). Esta herramienta digital en línea calcula los números de onda de corte de un cable coaxial mediante el método integral de Cauchy (Austin, Kravanja y Lloyd, 2014). En la tabla se han evaluado las raíces en la función f(kcn) para observar la precisión de ambos métodos, siendo f la función 𝑓𝑇𝑀(𝑘𝑐𝑛) o 𝑓𝑇𝐸 (𝑘𝑐𝑛) según corresponda con el modo.
Los tiempos de computación invertidos en nuestro sistema informático para obtener las raíces de las ecuaciones
trascendentes han sido de 0.0011 segundos para los modos TM (los 4 a la vez) y 0.0022 para los TE (los 4 a la vez). Esta diferencia de tiempos significativa (tarda el doble en modos TE que en TM) es debida a que los modos TE han de evaluar la función derivada de Bessel que, a su vez, se constituye de dos funciones Bessel ordinarias. Los tiempos de computación del programa de Peñaranda Foix (2010) están en torno a 1 segundo, por lo que se muestra en este caso la gran diferencia en tiempo computacional que ofrece el método de Chebyshev frente a otros.
En cuanto a la precisión, se puede observar en la tabla que las raíces proporcionadas por el método de Chebyshev están más próximas a cero que las de la aplicación, por lo que el método de Chebyshev es más preciso y rápido en este caso. Concretamente, obtenemos raíces que presentan un error absoluto de un orden de magnitud de 10-6, que podemos considerar suficientemente aceptable. Aunque, si quisiéramos más precisión, bastaría con aproximar la función por un orden más alto del polinomio de Chebyshev. En este caso, por ejemplo, si tomamos n=30 en lugar de n=15, el error absoluto de las raíces pasa a ser del orden 10-14, por lo que mejora bastante la precisión del método, mientras que el tiempo computacional no empeora significativamente, pues tarda 0.0017 segundos en calcular los modos TM y 0.0035 en calcular los modos TE.
En este artículo se ha presentado un método de resolución de ecuaciones trascendentes mediante polinomios de Chebyshev. Se ha comprobado la precisión, estabilidad y rapidez del método para varias funciones trascendentes que aparecen en varios campos de la física aplicada. Se ha mostrado también la versatilidad del método para varios tipos de funciones (trigonométricas, exponenciales, de Bessel) y para todo tipo de valores numéricos (valores muy elevados de la función (del orden de 1010 en el caso de la ecuación de Schrödinger) o muy pequeños en la variable independiente (del orden de 10-17 en el caso de la ecuación de Schrödinger). También se ha podido comprobar la facilidad de implementación y uso mediante el software matemático MATLAB, pues con introducir como parámetros de entrada la función trascendente, el intervalo de estudio y el orden del polinomio, el método proporciona todas las raíces que se encuentren en ese intervalo, lo que no hacen los métodos iterativos tradicionales como el de Newton, Traub o la Secante y hace a este método óptimo para usarse en estas situaciones donde las ecuaciones trascendentes presenten múltiples soluciones y que tantas veces aparecen en el campo de la física. Además, con los códigos fuente proporcionados se pretende dotar a los docentes de física, alumnos e investigadores de una herramienta útil, fácil y fiable para la resolución de ecuaciones trascendentes. También es interesante destacar que, en todos los casos estudiados, el método de Chebyshev ha presentado la mayor velocidad computacional, lo cual puede ser un parámetro determinante en algún campo de investigación donde se precise una cierta inmediatez para conocer la solución, como en simulaciones electromagnéticas o electrotérmicas.
Para líneas futuras queda, por una parte, implementar algoritmos iterativos a partir de este método que encuentren un número determinado de raíces de una función trascendente, lo cual se podría implementar realizando una segmentación del eje real en intervalos finitos y determinando las raíces en cada intervalo hasta encontrar el número solicitado. Por otra parte, queda estudiar la posibilidad de aplicar este método para la resolución de funciones trascendentes en el plano complejo y ampliar así la versatilidad del método en el campo de la física, como por ejemplo los problemas de relatividad general donde se trabaja con espacios de Minkowsky o en el estudio cuántico del espín de las partículas subatómicas.
Aranzabal Olea, A. (2001). Modos de resolución de circuitos con diodos. Recuperado el 10 de agosto de 2022, de http://www.sc.ehu.es/sbweb/electronica/elec_basica/tema3/Paginas/Pagina6.htm
Austin, A., Kravanja, P. & Lloyd N., T. (2014). Numerical algorithms based on analytic function values at roots of unity.
SIAM Journal on Numerical Analysis, 52(4), 1795-1821. Bohm, D. (1989). Quantum Theory. New York: Dover.
Boyd, J. (2002). Computing zeros on a real interval through Chebyshev expansion and polynomial rootfinding. SIAM Journal on Numerical Analysis, 40(5), 1666-1682.
Boyd, J. (2014a). Finding the zeros of a univariative equation: Proxy rootfinders, Chebyshev interpolation and the companion matrix. Philadelphia, United States: Society for Industrial and Applied Mathematics.
Boyd, J. (2014b). Solving Transcendental Equations: The Chebyshev Polynomial Proxy and Other Numerical rootfinders, Perturbation Series, and Oracles. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics (SIAM).
FCenRed. (2014, octubre 14). El pozo cuadrado finito. Física cuántica en la red. http://www.fisicacuantica.es/pozo-cuadrado-finito/
Galindo, A. y Pascual, P. (1989). Mecánica Cuántica. Madrid: Eudema.
Ma, D. (26 de Marzo de 2010). General Mathieu functions with arbitrary parameters V1.0. Recuperado el 2 de agosto de 2022, de: https://www.mathworks.com/matlabcentral/fileexchange/27101-general-mathieu-functions-with-arbitrary-parameters-v1-0
Morris, S. (29 de mayo de 2007). Real roots on interval. Recuperado el 2 de agosto de 2022, de https://www.mathworks.com/matlabcentral/fileexchange/15122- real-roots-on-interval
Peñaranda Foix, F. L. (2010). Números de onda de corte de los modos superiores TE y TM de un coaxial. https://riu-net.upv.es/handle/10251/8070
Pozar, D. (1998). Microwave Engineering. New York: Wiley & Sons, Inc.
Quinga, S. (2021). Ecuaciones no lineales en física y su resolución mediante el uso de métodos iterativos multipaso de orden alto. Revista Enseñanza de la Física, 33(3), 145-165.
Sze, S. & Kwok, K. (2007). Physics of semiconductors devices. Taiwan: John Wiley & Sons, Inc.
Se proporciona el código fuente de MATLAB para las funciones trascendentes f(x)=0 presentadas en este trabajo y las líneas de llamada para inicializar las variables, representar las funciones gráficamente y encontrar las raíces mediante el método de Chebyshev.
function [f,df]=resorte(x,m,g,h,k1,k2)
f=m.*g.*h+m.*g.*x-1/2.*k1.*x.^2-2./5.*k2.*x.^(5/2);
df=g.*m - k1.*x - k2.*x.^(3/2); % se calcula la derivada para aplicar el método de Traub
end
% Inicialización de variables
m=95; h=0.43; k1=40000; k2=40; g=9.8;
% Líneas de llamada para representar gráficamente las funciones trascendentes x=linspace(0,0.3,1000);f=resorte(x,m,g,h,k1,k2);
figure; plot(x,f,'k');hold on;plot(x,zeros(size(x)),'k');
xlabel('d (m)');ylabel('f(d)');grid
% Líneas de llamada para encontrar las raíces F=@(xx) resorte(xx,m,g,h,k1,k2);
tic;[ROOTS_Cheby]=FindRealRoots(F,0.15,0.2,15);toc % FindRealRoots: función disponible en Morris (2007)
Precision=F(ROOTS_Cheby);
tic;[x] = Traub(@(x)F(x),0.16,Precision,100);toc; % Traub: función disponible en Quinga (2021)
function [f,df]=diodo(x,Vbb,Is,Vt,RL) f=Vbb-x-Is.*(exp(x./Vt)-1).*RL;
df=- (Is.*RL.*exp(x./Vt))./Vt - 1; % se calcula la derivada para aplicar el método de Traub
end
% Inicialización de variables Vbb=1; Is=1.2e-11; Vt=26e-3; RL=4;
% Líneas de llamada para representar gráficamente las funciones trascendentes x=linspace(0.5,0.65,1000);f=diodo(x,Vbb,Is,Vt,RL);
figure; plot(x,f,'k');hold on;plot(x,zeros(size(x)),'k');
xlabel('VD');ylabel('f(VD)');grid
% Líneas de llamada para encontrar las raíces F=@(xx) diodo(xx,Vbb,Is,Vt,RL);
tic;[ROOTS_Cheby]=FindRealRoots(F,0.5,0.65,15);toc % FindRealRoots: función disponible en Morris (2007)
Precision=F(ROOTS_Cheby);
tic; [x] = Traub(@(x)F(x),0.55,Precision,100); toc; % Traub: función disponible en Quinga (2021)
function f=schrodinger_par(x,m,a,V0) h=6.63e-34/2/pi;
f=sqrt(2.*m.*x)./h.*sin(sqrt(2.*m.*x)./h.*a)-sqrt(2.*m.*(abs(V0)-x))./h.*cos(sqrt(2.*m.*x)./h.*a);
end
function f=schrodinger_impar(x,m,a,V0) h=6.63e-34/2/pi;
f=sqrt(2.*m.*x)./h.*cos(sqrt(2.*m.*x)./h.*a)+sqrt(2.*m.*(abs(V0)-x))./h.*sin(sqrt(2.*m.*x)./h.*a);
end
% Inicialización de variables
m=9.11e-31; a=0.2e-9; V0=75*1.6e-19;
% Líneas de llamada para representar gráficamente las funciones trascendentes
x=linspace(0,75*1.6e-19,1000); fpar=schrodinger_par(x,m,a,V0);fimpar=schrodinger_impar(x,m,a,V0); figure; plot(x,fpar,'k');hold on;plot(x,fimpar,'--k');plot(x,zeros(size(x)),'k');
xlabel('E (Julios)');ylabel('f(E)'); grid
% Líneas de llamada para encontrar las raíces FPAR=@(xx) schrodinger_par(xx,m,a,V0);
tic;[ROOTS_ChebyPar]=FindRealRoots(FPAR,0,1.2e-17,15);toc
PrecisionPar= FPAR (ROOTS_ChebyPar); FIMPAR=@(xx) schrodinger_impar(xx,m,a,V0);
tic;[ROOTS_ChebyImpar]=FindRealRoots(FIMPAR,0,1.2e-17,15);toc PrecisionImpar= FIMPAR (ROOTS_ChebyImpar);
NIVELES_ENERGIA=sort([ROOTS_ChebyPar;ROOTS_ChebyImpar]./1.6e-19); %pasamos de J a eV
function f=coaxial_TM(x,nu,a, b)
f=besselj(nu,x.*a).*bessely(nu,x.*b)-bessely(nu,x.*a).*besselj(nu,x.*b);
end
function f=coaxial_TE(x,nu,a, b) nu=nu.*ones(size(x));
f=dbesselj(nu,x.*a).*dbessely(nu,x.*b)-dbessely(nu,x.*a).*dbesselj(nu,x.*b);
%Las funciones dbesselj y dbessely calculan las derivadas de las funciones de Bessel de 1º y 2º especie, se puede encontrar programadas en MATLAB en (Ma, 2010)
end
% Inicialización de variables nu=0; a=1e-3; b=8e-3;
% Líneas de llamada para representar gráficamente las funciones trascendentes x=linspace(0,2000,1000);nu=nu.*ones(size(x));fTE=coaxial_TE(x,nu,a, b);fTM=coaxial_TM(x,nu,a,b); figure; plot(x,fTM,'k');hold on;plot(x,fTE,'k--');plot(x,zeros(size(x)),'k'); xlabel('kcn');ylabel('f(kcn)');grid
% Líneas de llamada para encontrar las raíces FTM=@(xx) coaxial_TM(xx,nu(1),a,b);
tic;[ROOTS_ChebyTM]=FindRealRoots(FTM,0,2000,15);toc
PrecisionTM=FTM(ROOTS_Cheby); FTE=@(xx) coaxial_TE(xx,nu(1),a,b);
tic;[ROOTS_ChebyTE]=FindRealRoots(FTE,0,2000,15);toc PrecisionTE=F(ROOTS_Cheby);
% Para aumentar precisión: aumentar orden n, lo cambiamos de 15 a 30 para que se observe la mejora en precisión tic;[ROOTS_ChebyTM]=FindRealRoots(FTM,0,2000,30);toc
PrecisionTM=FTM(ROOTS_Cheby);
tic;[ROOTS_ChebyTE]=FindRealRoots(FTE,0,2000,30);toc PrecisionTE=F(ROOTS_Cheby);
