4.1.11. Series infinitas

En la secciones anteriores hemos considerado las funciones como entidades abstractas. En todos los casos asumimos que que el valor de esas funciones puede obtenerse con algún procedimiento que no hemos especificado. Considere por ejemplo el caso de la función trigonométrica seno. Sabemos como se define, conocemos su relación con otras funciones, su antiderivada y todas las derivadas de orden arbitrario de la función. Pero ¿sabemos cómo calcular el valor de $\sin t$?

Para ser estrictos las únicas funciones cuyos valores se pueden calcularse de manera elemental, es decir ejecutando un número finito de sumas (o restas), multiplicaciones y divisiones son las funciones polinómicas y las funciones racionales.

Definimos una función polinómica como aquella que puede escribirse de la forma: $$ P_k(t)=a_0+a_1 t+a_2 t^2+\ldots a_k t^k $$ decimos que $P_k(t)$ es un polinomio de orden $k$.

Las funciones racionales son aquellas que pueden escribirse de la forma:

$$ R(t)=\frac{P_k(t)}{Q_m(t)} $$

done $P_k(t)$ y $Q_l(t)$ son polinomios de orden $k$ y $m$ respectivamente.

4.1.11.1. Series de potencias

El primer resultado de interés para el cálculo del valor de funciones no polinómicas ni racionales depende del siguiente teorema:

Proposición: Teorema de Taylor. Si una función $f(t)$ tiene derivadas hasta orden $k$ en el punto $t=t_0$ entonces el polinomio:

$$P_k(x)=\sum_{n=0}^{k} \frac{f^{(n)}(t_0)}{n!}(t-t_0)^n$$

ofrece una aproximación al valor de la función para valores de $t\rightarrow t_0$. En otras palabras:

$$\lim_{t\rightarrow t_0}\;[f(t)-P_k(t)]=0$$

En términos prácticos, el teorema de Taylor se puede escribir diciendo que para valores de $t$ cercano a $t_0$, la función $f(t)$ se puede escribir como:

$$ \mathrm{Para}\;t\rightarrow t_0:\;f(t)\approx\sum_{n=0}^{k} \frac{f^{(n)}(t_0)}{n!}(t-t_0)^n $$

A la suma de la forma del lado izquierdo de la ecuación anterior, $s_k\equiv\sum_{n=0}^{k} a_n$ la llamamos en general una suma parcial y a la secuencia $\{s_0,s_1,s_2,...\}$ se la conoce como una serie infinita o en breve una serie.

En particular a la serie formada por las sumas parciales en el teorema de Taylor la llamamos una serie de potencias o serie de Taylor.

Definición 11.3. Serie convergente. Decimos que una serie $\{s_0,s_1,s_2,...\}$ converge a $S$ si:

$$\lim_{k\rightarrow\infty} s_k=S$$

en otros términos:

$$\sum_{n=0}^{\infty}a_n\equiv\lim_{k\rightarrow\infty} \sum_{n=0}^{k} a_n=S$$

En la definción anterior hemos introducido la notación $\sum_{n=0}^{\infty}a_n$ para representar el valor al que converge una serie (si lo hace.) Esta notación es corriente en el cálculo y será muy utilizada en este libro.

En términos de esta notación y del concepto de convergencia, el teorema de Taylor se puede escribir como:

\begin{equation} \label{eq:serie_potencias_limite} f(t)=\lim_{t\rightarrow t_0}\sum_{n=0}^{\infty} \frac{f^{(n)}(t_0)}{n!}(t-t_0)^n \end{equation}

Ahora bien, en algunos casos la serie de potencias asociada a una función puede converger incluso en puntos a una distancia finita de $t_0$. Considere por ejemplo la serie de potencias de la función exponencial $\exp(t)$ alrededor de $t_0=0$ (ver problemas al final del capítulo):

$$ \exp(t)\approx s_k=\sum_{n=0}^{k} \frac{t^n}{n!} $$

Es posible mostrar que esta serie converge para cualquier valor real de $t$. Por la misma razón, usando la notación en la Def. (funcional_integral) podemos escribir:

\begin{equation} \label{eq:exp_taylor} \exp(t)=\sum_{n=0}^{\infty} \frac{t^n}{n!} \end{equation}

e introducir una nueva definición:

Definición 11.4. Convergencia uniforme y función analítica. Decimos que una serie $s_k$, cuyos coeficientes son función de una variable independiente $t$, converge uniformemente a $f(t)$ si lo hace para todo valor de $t$ que cumple $|t-t_0|<r$, donde $r>0$ se conoce como el radio de convergencia. En otras palabras:

$$f(t)=\sum_{n=0}^{\infty} a_n(t)$$

A la función $f(t)$ que cumple esta condición la llamamos una función analítica.

Muchas de las series de potencias definidas con la Ec. (serie_potencias_limite) convergen uniformemente al valor de la función que expanden. En ese caso escribimos:

\begin{equation} \label{eq:serie_taylor} f(t)=\sum_{n=0}^{\infty} \frac{f^{(n)}(t_0)}{n!}(t-t_0)^n \end{equation}

Las series de taylor de algunas funciones analíticas conocidas que serán de gran utilidad en la mecánica celeste, así como su radio de convergencia uniforme, son provistas en la lista a continuación:

  • Funciones seno y coseno. Para todo $t\in {\rm I\!R}$: \begin{equation} \label{eq:sin_taylor} \sin t=\sum_{n=0}^{\infty} \frac{(-1)^n t^{2n+1}}{(2n+1)!}=t-\frac{t^3}{3!}+\frac{t^5}{5!}-\frac{t^7}{7!}+\ldots \end{equation}

    \begin{equation} \label{eq:cos_taylor} \cos t=\sum_{n=0}^{\infty} \frac{(-1)^n t^{2n}}{2n!}=1-t^2+\frac{t^4}{4!}-\frac{t^6}{6!}+\ldots \end{equation}

  • Funciones seno y coseno hiperbólicas. Para todo $t\in {\rm I\!R}$: \begin{equation} \label{eq:sinh_taylor} \sinh t=\sum_{n=1}^{\infty} \frac{t^{2n+1}}{(2n+1)!}=t+\frac{t^3}{3!}+\frac{t^5}{5!}+\frac{t^7}{7!}+\ldots \end{equation}

    \begin{equation} \label{eq:cosh_taylor} \cosh t=\sum_{n=0}^{\infty} \frac{t^{2n}}{2n!}=1+t^2+\frac{t^4}{4!}+\frac{t^6}{6!}+\ldots \end{equation}

Cálculo numérico de series infinitas

Las series de potencias juegan un papel central en la computación porque permiten reducir el cálculo de funciones que pueden ser muy complicadas a la ejecución de un conjuntos finito de operaciones numéricas simples: sumas (o restas), multiplicaciones y divisiones. Sin embargo usar un computador para calcular eficientemente el valor de convergencia de una serie infinita no es una tarea trivial.

Tres problemas comunes se enfrentan al calcular series infinitas en un computador. Para empezar la precisión con la que podemos representar cantidades numéricas en un computador esta limitada a un cierto número de dígitos significativos. Para ilustrarlo, considere esta operación:

In [1]:
x=1+1e-16
1.00000000000000000000

Llamamos a la cantidad epsilon ($\approx 10^{-16}$) que cumple la condición 1+epsilon=1 la precisión de la máquina. Esta limitación implica que para un número representado en notación científica:

$$ \mathrm{mantisa}\times 10^{\mathrm{exponente}} $$

donde $0\leq\mathrm{mantisa}<1$, tan solo los primeros 16 dígitos decimales son realmente significativos. Para hacer esto más evidente considere esta operación:

In [3]:
x=0.123456789012345678901234567890
0.123456789012345677369886232100

Nótese que aunque hemos querido almacenar 30 dígitos significativos, tan solo los primeros 16 fueron efectivamente recuperados de la memoria. A partir del decimal 17 lo que obtenemos son cifras al azar.

Al manipular algebraicamente cantidades que tienen esta limitación, los resultados pueden ser aún menos significativas. En un algortimo (en especial aquellos con los que se calculan el valor de convergencia de series) después de muchas operaciones entre cantidades con precisión limitada el resultado puede llegar a acumular importantes errores de redondeo.

Para ilustrar la manera como los errores de redondeo se propagan considere el cálculo de esta sencilla cantidad:

$$ (Nx+y)-Nx $$

donde $N\gg 1$. Si bien matemáticamente el resultado es obviamente $y$, en el computador obtenemos:

In [5]:
N=1e12
x=1.0
y=0.1
z=(N*x+y)-(N*x)
x = +1.000000000000000e+00
y = +1.000000000000000e-01
z = +9.997558593750000e-02

A este tipo de errores, que se producen normalmente cuando se restan o dividen cantidades muy parecidas, se los conoce como errores de cancelación substractiva y pueden llegar a ser muy importantes, por ejemplo, al calcular series alternantes.

Para ilutstrar como los errores de cancelación substractiva pueden volverse significatiovos en la estimación del valor de convergencia de una serie alternante consideren el siguiente algoritmo en el que se calcula la suma parcial de la serie de Taylor para la función seno:

In [7]:
from math import factorial
from numpy import pi

#Numero de términos de la serie
k=18
#Valor del ángulo
t=3*pi
#Suma parcial
sk=0
for n in range(k):
    sk+=(-1)**n*t**(2*n+1)/factorial(2*n+1)
s_18(9.42478) = -7.6538436107548528e-08
sin(9.42478) = 3.6739403974420594e-16

Si bien el resultado obtenido con la suma parcial es correcto hasta la séptima cifra significativa ($\sim 10^{-7}$), hay una manera más eficiente de calcular este resultado evitando la propagación de errores de redondeo (cancelación substractiva). Podemos aprovechar la propiedad de la función:

$$ \sin(t)=\sin(t+2n\pi) $$

para que reconocer que $\sin(3\pi)=\sin(\pi)$ y calcular la suma parcial en el valor más pequeño de $t$:

In [9]:
#Valor de t mapeado en el intervalo [0,2 pi)
from numpy import mod
tn=mod(t,2*pi)

#Suma parcial
sk=0
for n in range(k):
    sk+=(-1)**n*mod(tn,2*pi)**(2*n+1)/factorial(2*n+1)
s_18(9.42478) = 3.3280566951896521e-16
sin(9.42478) = 3.6739403974420594e-16

La diferencia en los dos resultados se hace notar notar inmediatamente.

Otro truco que puede ser útil, se ilustra en el cálculo de la función:

\begin{equation} \label{eq:exp_taylor_general} \exp(-t)=\sum_{n=0}^{\infty} \frac{(-1)^n t^n}{n!} \end{equation}

In [11]:
#Número de términos en la serie
k=10
#Valor de t
t=5.0
#Suma pacial de exp(-t)
sk=0
for n in range(k):
    sk+=(-1)**n*t**n/factorial(n)
s_10(1) = -1.8271053791887111e+00
exp(1) = 2.7182818284590451e+00

Si en lugar de la serie alternante usamos la identidad $\exp(-t)=1/\exp(t)$ y calculamos la suma parcial correspondiente a $\exp(t)$ (Ec. exp_taylor) obtenemos:

In [13]:
#Suma parcial de exp(t)
sk=0
for n in range(k):
    sk+=t**n/factorial(n)
1/s_10(5) = 6.9594528636495370e-03
exp(5) = 6.7379469990854670e-03

El tercer problema es la convergencia lenta, que puede hacer que se requieran decenas o cientos de términos de una serie para conseguir una precisión aceptable en el valor estimado de convergencia. La convergencia lenta no solo aumenta el tiempo requerido para calcular una serie, sino que también puede hacer que los errores de redondeo y cancelación sustractiva se acumulen al punto en el que el resultado sea más error que otra cosa.

La siguiente rutina implementa el cálculo de la serie de potencias de la función exponencial pero en lugar de escoger un número fijo de términos de la suma parcial, como lo hicimos en los Algs. (exp_negativo y (exp_positivo), la rutina va agregando nuevos términos hasta que consigue un nivel de precisión o tolerancia deseado $\delta\ll 1$. Para esto, por cada término que agrega a la serie calcula el error relativo el valor de la suma como:

$$ \Delta_n\equiv\left|\frac{s_{n+1}-s_{n}}{\bar{s}}\right| $$

donde $\bar{s}=(s_n+s_{n+1})/2$ es el promedio entre dos sumas parciales sucesivas.

El cálculo de la serie converge cuando se cumple:

$$ \Delta_n<\delta $$

In [15]:
def exp_taylor(t,delta=1e-15):
    #Valor inicial del error relativo (número arbitrario)
    Dn=1
    #Primer valor de la suma parcial
    sn=0
    #Inicializa contador de términos
    n=0
    #El ciclo continúa mientras no se alcance la tolerancia
    while Dn>delta:
        #Almacena el último valor de la suma
        s=sn
        #Calcula el término n-esimo de la suma
        from math import factorial
        an=t**n/factorial(n)
        #Agrega el término a la suma n-esima
        sn+=an
        #Promedio de dos sumas parciales sucesivas
        smed=(sn+s)/2
        #Error relativo
        Dn=abs(an/smed)
        #Incremento en el contador n
        n+=1
    return sn,Dn,n

Aunque esta no es la forma más compacta y efectiva de calcular la función exponencial, la mayoría de las rutinas que usan iteraciones en este libro tendrán una estructura y notación similares a esta. Además, esta forma de escribirla nos permitirá ilustrar algunos de las limitaciones mencionados anteriormente.

Considere por ejemplo el cálculo de la función $\exp(-t)$:

In [16]:
t=20
x,error_x,n_x=exp_taylor(-t)
y,error_y,n_y=exp_taylor(+t)
x = exp(-20) = 5.4781029165292094e-10 (error 7.00e-16, n = 96)
1/y = 1/exp(+20) = 2.0611536224385591e-09 (error 9.22e-16, n = 66)
exp(-20) = 2.0611536224385579e-09

Nótese que todos los efectos descritos hasta aquí son otra vez evidentes: (1) el error en la serie alternada es grande y debido principalmente a cancelaciones substractivas, (2) la serie alternante requiere un número mayor de términos para conseguir el mismo nivel de precisión de la serie no alternante (convergencia lenta.)

4.1.11.2. Otras series

A continuación se enumeran otras series que serán de utilidad en el texto:

  • Serie geométrica: para $t\neq 1$,

    \begin{equation} \label{eq:serie_geometrica} \frac{1}{1-t}=\sum_{n=0}^{\infty} t^n=1+t+t^2+t^3+\ldots \end{equation}

  • Serie binomial: para todos los valores de $a\in{\rm I\!R}$ y $|t|<1$,

    \begin{equation} \label{eq:serie_binomio} (1+t)^a=\sum_{n=0}^{\infty} \left(\begin{array}{c}a\\k\end{array}\right) t^n=1+a t+\frac{a(a-1)}{2!}t^2+\frac{a(a-1)(a-2)}{3!}t^3+\cdots \end{equation}

    Los coeficientes de esta serie tienen la forma general:

    $$ \left(\begin{array}{c}a\\k\end{array}\right)\equiv \frac{a (a-1) (a-2) \cdots (a-k+1)}{k!} $$

    Cuando para aproximar la función $(1+t)^a$ se usa una de las sumas parciales de la serie binomial, decimos que se hace una aproximación binomial de $(1+t)^a$.

4.1.11.3. Series de Stumpff

Una familia de funciones analíticas de gran interés en la mecánica celeste son las funciones o series de Stumpff:

$$ c_k(t)=\sum_{n=0}^{\infty} \frac{(-1)^n t^{n}}{(2n+k)!}=\frac{1}{k!}-\frac{t}{(k+2)!}+\frac{t^2}{(k+4)!}-\frac{t^3}{(k+6)!}+\cdots $$

que convergen para todo $t\in{\rm I\!R}$

Estas series tienen una serie de propiedades muy interesantes que serán explotadas en el Capítulo El problema de los dos cuerpos y que enumeramos a continuación:

  • Relación de recurrencia: Si se conoce el valor de la función de orden $k+2$ es posible calcular el valor de la función $k$ usando:

    \begin{equation} \label{eq:stumpff_recurrencia} c_{k}(t)=\frac{1}{k!}-t c_{k+2}(t) \end{equation}

  • Relación de recurrencia para las derivadas: Las derivadas de las series de Stumpff satisfacen las siguientes relaciones de recurrencia:

    \begin{equation} \label{eq:stumpff_derivadas_recurrencia} \begin{array}{rcll} \frac{\mathrm{d}}{\mathrm{d}t}\left[c_0(\alpha s^2)\right] & = & -\alpha sc_{1}(\alpha s^2) & \\ \frac{\mathrm{d}}{\mathrm{d}t}\left[t^{k+1}c_{k+1}(\alpha s^2)\right] & = & s^{k}c_{k}(\alpha s^2) & \mathrm{para}\;k=0,1,2,\ldots\\ \end{array} \end{equation} donde $\alpha$ es cualquier número real. La parametrización de la variable independiente de las series de Stumpff en la forma $\alpha s^2$ será muy común en mecánica celeste.

  • Equivalencia con las funciones trigonométricas e hiperbólicas: Todas series de Stumpff se pueden escribir en términos de las funciones trigonométricas e hiperbólicas. Por esta razón se consideran una generalización de estas últimas. Se puede probar, usando las series de Taylor en las Ecs. (sin_taylor)-(cos_taylor) que las primeras dos series de Stumpff se pueden escribir como:

    \begin{equation} \label{eq:c0_trigo} c_0(\alpha s^2)=\left\{ \begin{array}{ll} \cos(\alpha^{1/2}s) & \mathrm{Si}\;\alpha\geq 0\\ \cosh(|\alpha|^{1/2}s) & \mathrm{Si}\;\alpha<0\\ \end{array} \right. \end{equation} y

    \begin{equation} \label{eq:c1_trigo} c_1(\alpha s^2)=\left\{ \begin{array}{ll} \sin(\alpha^{1/2}s)/(\alpha^{1/2}s) & \mathrm{Si}\;\alpha\geq 0\\ \sinh(|\alpha|^{1/2}s)/(|\alpha|^{1/2}s) & \mathrm{Si}\;\alpha<0\\ \end{array} \right. \end{equation}

  • Relaciones de escalado: Si bien las series de Stumpff convergen para todos los valores reales de la variable independiente, el cálculo numérico de las mismas para valores muy grandes de $|t|$ puede sufrir de errores de redondeo o ser muy ineficiente. Como vimos en el caso del cálculo numérico de la serie de Taylor del seno, una solución posible para este inconveniente es la de mapear un valor arbitrario de $t$ en un intervalo más pequeño. En el caso de la función coseno, por ejemplo, se puede calcular el coseno de cualquier ángulo grande, por ejemplo $2t$ usando del coseno del ángulo más pequeño $t$ y la identidad:

    $$ \cos(2t)=2\cos(t)^2-1 $$

    Las siguientes juegan un papel análogo en el caso del cálculo de las funciones de Stumpff.

    \begin{eqnarray} c_0(4t) & = & 2c_0(t)^2-1\\ c_1(4t) & = & c_0(t)c_1(t)\\ c_2(4t) & = & \frac{1}{2}c_1(t)^2\\ c_3(4t) & = & \frac{1}{4}c_2(t)+\frac{1}{4}c_0(t)c_3(t)\\ \end{eqnarray}

  • Regla iterativa: Las series de Stumpff se pueden calcular también usando la regla iterativa (ver problemas al final del capítulo):

    \begin{equation} \label{eq:stumpff_iteracion} p_n=\frac{t}{(2n+k+1)(2n+k+2)}(1-p_{n+1}) \end{equation}

    Si se comienza con un valor $p_{N}=0$ para un valor de $N$ suficientemente grande, al descender con la regla iterativa hasta $p_0$ se obtiene una aproximación para la serie de Stumpff de orden $k$:

    $$ c_k(t)\approx\frac{1}{k!}(1-p_0) $$

    El error de esta aproximación es aproximadamente igual a:

    \begin{equation} \label{eq:stumpff_iteracion_error} \epsilon_{N}=\frac{t^N}{(2N+k+2)!} \end{equation}

En el Apéndice Algoritmos y rutinas útiles hemos incluido una corta rutina, serie_stumpff(t,k) que usa la regla iterativa para calcular el valor de las series de Stumpff de forma muy eficiente. Para ilustrar el uso de esta rutina, en el algoritmo a continuación ponemos a prueba algunas de las propiedades mencionadas en esta sección.

Así por ejemplo, podemos verificar la relación de recurrencia (Ec. stumpff_recurrencia):

In [18]:
#Verifica las reglas de recurrencia
t=2
k=4

from pymcel.export import serie_stumpff
ck=serie_stumpff(t,k)
ckp2=serie_stumpff(t,k+2)
Regla de recurrencia:
c_4(2) = 0.03898592369134362
1/4!-t c_6(2) = 0.03898592369134361
In [20]:
#Verifica las reglas de escalado
t=2
k=4

from pymcel.export import serie_stumpff
c0_4t=serie_stumpff(4*t,0)
c0_t=serie_stumpff(t,0)
c1_4t=serie_stumpff(4*t,1)
c1_t=serie_stumpff(t,1)
c2_4t=serie_stumpff(4*t,2)
c2_t=serie_stumpff(t,2)
c3_4t=serie_stumpff(4*t,3)
c3_t=serie_stumpff(t,3)
c_0(8) = -0.9513631281258474, con escalado = -0.9513631281258474
c_1(8) = 0.10891980905843202, con escalado = 0.10891980905843208
c_2(8) = 0.24392039101573093, con escalado = 0.24392039101573093
c_3(8) = 0.111385023867696, con escalado = 0.11138502386769598

4.1.11.4. Series de Fourier

A lo largo de esta sección hemos visto como es posible aproximar algunas funciones analíticas usando series de potencias o series de polinomios. Pero los polinomios no son las únicas funciones capaces de realizar esta tarea. Así por ejemplo, cualquier familia contable de funciones de variable real $f_n(t)$, que en un determinado dominio, obedezca la condición:

$$ \int_{\cal D} g_n(t)g_m(t)\;\mathrm{d}t=0,\mathrm{Si}\;n\neq m $$

puede aproximar funciones definidas en el dominio ${\cal D}$ con series del tipo:

$$ f(t)\approx\sum_{n=1}^{k} A_n g_n(t) $$

donde:

\begin{equation} \label{eq:An_serie_funciones} A_n=\frac{\int_D g_n(t)f(t)\;\mathrm{d}t}{\int_D g_n(t)^2\;\mathrm{d}t} \end{equation}

Uno de los conjuntos contables de funciones más notables usados en la historia para aproximar otras funciones es el de las funciones trigonométricas seno y coseno definidas en el dominio finito ${\cal D}=[0,T]$:

$$ \begin{array}{rcl} g_{0}(t) & = & 1/2\\ g_{2n-1}(t) & = & \sin(n\omega t)\\ g_{2n}(t) & = & \cos(n\omega t)\\ \end{array} $$

donde $n=1,2,\cdots$ y $\omega\equiv2\pi/T$. Es posible mostrar que:

\begin{eqnarray} \nonumber \int_0^{T} g_n(t) g_m(t) & = & 0\\ \nonumber \int_0^{T} g_n(t)^2 & = & \frac{T}{2}\\ \end{eqnarray}

Las series resultantes, que se puede escribir como

\begin{equation} \label{eq:serie_fourier} s_k(t)=\frac{1}{2}A_0+\sum_{n=1}^{k}A_n\cos(n\omega t)+\sum_{n=1}^{k}B_n\sin(n\omega t) \end{equation}

y que permiten aproximar funciones $f(t)$ en el dominio $[0,T]$ se conocen como Series de Fourier.

Si una función $f(t)$ es aproximada con una serie de Fourier, los coeficientes correspondientes se calculan usando la Ec. (An_series_funciones):

\begin{eqnarray} \label{eq:A0_fourier} A_0 & = & \frac{2}{T}\int_0^{T} f(t)\;\mathrm{d}t\\ \label{eq:An_fourier} A_n & = & \frac{2}{T}\int_0^{T} \cos(n\omega t) f(t)\;\mathrm{d}t\\ \label{eq:Bn_fourier} B_n & = & \frac{2}{T}\int_0^{T} \sin(n\omega t) f(t)\;\mathrm{d}t\\ \end{eqnarray}

Como sucede con las series de potencias, una cosa es aproximar una función en un dominio usando series de funciones y otra muy distinta es asegurar la convergencia uniforme de la serie infinita. En el caso de las series de Fourier el teorema de Dirichlet ofrece las condiciones para dicha convergencia:

Proposición: Teorema de Dirichlet para series de Fourier. Dada una función $f(t)$ definida en el dominio $[0,T]$ y que cumple las siguientes condiciones dentro de este intervalo:

  1. Es integrable.
  2. Esta acotada.
  3. Tiene un número finito de discontinuidades.

La serie de Fourier (Ec. serie_fourier) converge uniformemente a $f(t)$ para cualquier valor de $t\in[0,T]$

Para ilustrar lo visto aquí sobre las series de Fourier, a continuación presentamos el algoritmo de una rutina que calcula el valor de los coeficientes $A_0$, $A_n$ y $B_n$ para cualquier función provista por el usuario. El algoritmo se vale de la rutina quad para calcular numéricamente las integrales en las Ecs. (A0_fourier)-(Bn_fourier):

In [22]:
def coeficientes_fourier(funcion,T,k,args=()):
    #Funciones externas
    from scipy.integrate import quad
    from numpy import sin,cos

    #Parametro omega
    w=2*pi/T
    
    #Determina los coeficientes en t:
    f=lambda t:funcion(t,*args)
    As=[2*quad(f,0,T,args=args)[0]/T]
    Bs=[0]
    for n in range(1,k+1):
        f_cos_n=lambda t:funcion(t,*args)*cos(n*w*t)
        As+=[2*quad(f_cos_n,0,T)[0]/T]
        f_sin_n=lambda t:funcion(t,*args)*sin(n*w*t)
        Bs+=[2*quad(f_sin_n,0,T)[0]/T]
    
    return As,Bs

Nótese que el usuario debe proveer la lista de las opciones de la función funcion en la forma de la tupla args (ver Sección Funciones). Esta rutina solo permite calcular el valor de los coeficientes de la serie (almacenados en las listas As y Bs). Para calcular la suma parcial se debe usar la siguiente rutina:

In [23]:
def serie_fourier(t,As,Bs,T,k):
    from numpy import sin,cos

    #Parametro omega
    w=2*pi/T

    #Calcula la suma
    sk=As[0]/2
    for n in range(1,k+1):
        sk+=As[n]*cos(n*w*t)+Bs[n]*sin(n*w*t)
    
    return sk

Como un ejemplo, calculemos los primeros 30 coeficientes de la serie de Fourier de una línea recta:

In [24]:
def f(t):
    y=t
    return y
#Dominio de aproximación
T=2*pi
#Número de coeficientes
k=30
#Coeficientes
As,Bs=coeficientes_fourier(f,T,k)
As = [ 6.283  0.    -0.    ... -0.    -0.     0.   ]
Bs = [ 0.    -2.    -1.    ... -0.071 -0.069 -0.067]

Para comprobar el resultado podemos hacer un gráfico de la rutina serie_fourier con algunos o todos los coeficientes calculados:

In [26]:
#Valores de la variable independiente
import matplotlib.pyplot as plt
fig=plt.figure();
ax=fig.gca();

#Aproximaciones
from numpy import linspace
ts=linspace(0,T,200)
k=1
ax.plot(ts,serie_fourier(ts,As,Bs,T,k),label=f"k = {k}");
k=5
ax.plot(ts,serie_fourier(ts,As,Bs,T,k),label=f"k = {k}");
k=10
ax.plot(ts,serie_fourier(ts,As,Bs,T,k),label=f"k = {k}");
#Gráfico de la función
ts=linspace(0,T,30)
ax.plot(ts,f(ts),'k+',label="Función");

#Decoración
ax.legend();
ax.set_xlabel("$t$");
ax.set_ylabel("$f(t)$");

Figura 4.13.

En las libretas provistas con la versión electrónica del libro encontrarán una versión interactiva de este gráfico.

Figura 4.14.