3.4.2. Beispiel zur Fourierentwicklung

import numpy as np
from sympy import symbols, integrate, lambdify, exp, sin, cos, pi
import scipy.integrate as scint
import matplotlib.pyplot as plt
from IPython.display import Math, display

Für das Beispiel benutzen wir zwei verschiedene orthonormale Basen für den \(L_2[-1,1]\). Wir definieren das Skalarprodukt im \(L_2\) und die Norm:

t = symbols('t')

def dot(x,y):
    return integrate(x*y,(t,-1,1))
def norm(x):
    return dot(x,x)**(1/2)
def dotN(x,y):
    xy = lambdify(t, x*y,'numpy')
    return scint.quad(xy, -1, 1)[0]
def normN(x):
    return dotN(x,x)**(1/2)

Legendre’sche Polynome

Die Legendre’schen Polynome folgenden aus der Kontruktion eines ONS mit Hilfe des Verfahren nach Schmidt (vgl. vorangehendes Beispiel):

# Orthonormalbasis nach Schmidt
N = 9
yi = [t**i for i in range(N)]
xi = [yi[0]/norm(yi[0])]
for i in range(1,N):
    s = 0
    for j in range(i):
        s += dot(yi[i],xi[j])*xi[j]
    zi=yi[i]-s
    xi.append(zi/norm(zi))

for xii in xi:
    display(xii)

\[\displaystyle 0.707106781186547\]

\[\displaystyle 1.22474487139159 t\]

\[\displaystyle 2.37170824512628 t^{2} - 0.790569415042095\]

\[\displaystyle 4.67707173346743 t^{3} - 2.80624304008046 t\]

\[\displaystyle 9.2807765030735 t^{4} - 7.95495128834872 t^{2} + 0.795495128834872\]

\[\displaystyle 18.4685120543048 t^{5} - 20.5205689492276 t^{3} + 4.39726477483448 t\]

\[\displaystyle 36.8085471137468 t^{6} - 50.1934733369281 t^{4} + 16.7311577789763 t^{2} - 0.796721798998898\]

\[\displaystyle 73.4290553655431 t^{7} - 118.616166359723 t^{5} + 53.9164392544196 t^{3} - 5.99071547271328 t\]

\[\displaystyle 146.570997825634 t^{8} - 273.599195941123 t^{6} + 157.845689965986 t^{4} - 28.6992163574405 t^{2} + 0.797200454372895\]

Mit Hilfe dieses ONS soll die Fourierentwicklung der Funktion

\[f(t) = e^{-5 t^2}\]

berechnet werden.

f = exp(-5*t**2)
fn = lambdify(t, f, 'numpy')

tp = np.linspace(-1,1,400)
plt.plot(tp, fn(tp))
plt.grid()
plt.show()

Die Fourierkoeffizienten sind gegeben durch

\[\alpha_k = (f,x_i)_{L_2}\quad \forall i = 1, \ldots, n.\]

Die Berechnung nutzt das Skalarprodukt mit numerischer Integration.

alpha = [dotN(f, xii) for xii in xi]
alpha

[0.5596217150491691,
 0.0,
 -0.4411693576778274,
 0.0,
 0.21481237948395457,
 0.0,
 -0.07762200243541718,
 0.0,
 0.02214213441643933]

Damit erhalten wir die Fourierentwicklung von \(f(t)\)

s = 0
for alphai, xii in zip(alpha,xi):
    s += alphai*xii
sn = lambdify(t, s, 'numpy')
s

\[\displaystyle 3.24539473540683 t^{8} - 8.91522330646548 t^{6} + 9.38478407796754 t^{4} - 4.68931489413068 t^{2} + 0.994864373177097\]

tp = np.linspace(-1,1,400)
plt.plot(tp, fn(tp),label='f(t)') 
plt.plot(tp, sn(tp),label='Legendre Basis')
plt.legend()
plt.grid()
plt.show()

Wir überprüfen die Parsevallsche Gleichung (Vollständigkeitsrelation):

\[\sum_{k=1}^N |(x,x_k)|^2 \le \sum_{k=1}^\infty |(x,x_k)|^2 = \|x\|^2\]

np.sum(np.array(alpha)**2)-normN(f)**2

np.float64(-2.810714668699532e-05)

Trigonometrische Funktionen

Als zweites Set von Basisfunktionen benutzen wir \(\sin(\pi i t)\) und \(\cos(\pi i t)\).

yi = []
N = 5
yi.append(1/2)
for i in range(1,N):
    yi.append(cos(pi*i*t))
    yi.append(sin(pi*i*t))
print('Anzahl Basisfunktionen: '+str(len(yi)))

Anzahl Basisfunktionen: 9

Diese Funktionen sind orthogonal, entsprechend erhalten wir dieses Set unverändert aus dem Schmidtschen Orthogonalisierungsverfahren:

# Orthonormalbasis nach Schmidt
xi = [yi[0]/normN(yi[0])]
for i in range(1,2*N-1):
    s = 0
    for j in range(i):
        s += dotN(yi[i],xi[j])*xi[j]
    zi=yi[i]-s
    xi.append(zi/normN(zi))
    display(xi[-1])

\[\displaystyle 1.0 \cos{\left(\pi t \right)} - 4.6478198642589 \cdot 10^{-17}\]

\[\displaystyle 1.0 \sin{\left(\pi t \right)}\]

\[\displaystyle - 3.35746622125652 \cdot 10^{-17} \cos{\left(\pi t \right)} + 1.0 \cos{\left(2 \pi t \right)} + 9.0918731340168 \cdot 10^{-17}\]

\[\displaystyle - 1.33059369296921 \cdot 10^{-17} \sin{\left(\pi t \right)} + 1.0 \sin{\left(2 \pi t \right)}\]

\[\displaystyle 2.70104702138758 \cdot 10^{-17} \cos{\left(\pi t \right)} - 9.66870421923116 \cdot 10^{-17} \cos{\left(2 \pi t \right)} + 1.0 \cos{\left(3 \pi t \right)} - 1.00259879875121 \cdot 10^{-17}\]

\[\displaystyle - 2.36046319725668 \cdot 10^{-17} \sin{\left(\pi t \right)} - 1.8867319136135 \cdot 10^{-17} \sin{\left(2 \pi t \right)} + 1.0 \sin{\left(3 \pi t \right)}\]

\[\displaystyle - 4.93289800545064 \cdot 10^{-17} \cos{\left(\pi t \right)} + 1.24818653922218 \cdot 10^{-16} \cos{\left(2 \pi t \right)} + 9.03788555451802 \cdot 10^{-17} \cos{\left(3 \pi t \right)} + 1.0 \cos{\left(4 \pi t \right)} - 1.00476584178745 \cdot 10^{-17}\]

\[\displaystyle 8.81033073617638 \cdot 10^{-17} \sin{\left(\pi t \right)} - 1.77734401643418 \cdot 10^{-17} \sin{\left(2 \pi t \right)} - 6.14528970199674 \cdot 10^{-17} \sin{\left(3 \pi t \right)} + 1.0 \sin{\left(4 \pi t \right)}\]

Auf den ersten Blick sieht das Resultat nicht unverändert aus. Man beachte jedoch, dass die Koeffizienten numerisch null sind!

Für die Fourier-Koeffizienten erhalten wir:

alpha2 = [dotN(f, xii) for xii in xi]
alpha2

[0.5596217150491691,
 0.48508058459128284,
 0.0,
 0.10914699708164805,
 0.0,
 0.01007774067484598,
 0.0,
 -0.00025517128989833646,
 0.0]

und damit die Fourierreihe

s2 = 0
for alphai, xii in zip(alpha2,xi):
    s2 += alphai*xii
s2n = lambdify(t, s2, 'numpy')
s2

\[\displaystyle 0.485080584591283 \cos{\left(\pi t \right)} + 0.109146997081648 \cos{\left(2 \pi t \right)} + 0.010077740674846 \cos{\left(3 \pi t \right)} - 0.000255171289898337 \cos{\left(4 \pi t \right)} + 0.395712309610513\]

tp = np.linspace(-1,1,400)
plt.plot(tp, fn(tp),label='f(t)') 
plt.plot(tp, sn(tp),label='Legendre Basis')
plt.plot(tp, s2n(tp),label='trigonometrische Basis')
plt.legend()
plt.grid()
plt.show()

Parsevallsche Gleichung (Vollständigkeitsrelation):

\[\sum_{k=1}^N |(x,x_k)|^2 \le \sum_{k=1}^\infty |(x,x_k)|^2 = \|x\|^2\]

np.sum(np.array(alpha2)**2)-normN(f)**2

np.float64(-4.505698534273961e-07)

Vergleich

plt.plot(tp, sn(tp)-fn(tp),label='Legendre Basis')
plt.plot(tp, s2n(tp)-fn(tp),label='trigonometrische Basis')
plt.legend()
plt.grid()
plt.show()