Spectral discretization¶
A spectral discretization uses a basis of functions which dimension is the maximum order of representation of a solution. In an spectral collocation method, the variables are the values at N collocation points. The expected error is then proportional to $(1/N)^N$ unlike classical p-order methods whose error is proportional to $(1/N)^p$.
To analyze the convergence, a function f is discretized on collocation points. The derivated function is computed thanks to the differentiation operator and compared to the theoretical function.
In [1]:
Copied!
import aerokit.common.numspectral as sp
import numpy as np
def f(x):
return np.exp(-((x - 0.1) ** 2))
def df(x):
return -2 * (x - 0.1) * f(x)
def df2(x):
return -2 * (1 - 2.*(x - 0.1)**2) * f(x)
import aerokit.common.numspectral as sp
import numpy as np
def f(x):
return np.exp(-((x - 0.1) ** 2))
def df(x):
return -2 * (x - 0.1) * f(x)
def df2(x):
return -2 * (1 - 2.*(x - 0.1)**2) * f(x)
Validation of Chebyshev derivation operator¶
In [2]:
Copied!
import matplotlib.pyplot as plt
#
alln = range(5,31) # mesh sizes
err = np.zeros(len(alln))
#
for i, n in enumerate(alln):
SpOp = sp.ChebCollocation(n)
x = SpOp.x # collocation points
df_th = df(x)
df_num = SpOp.matder(1) @ f(x)
err[i] = np.sqrt(np.sum((df_th-df_num)**2)/n )
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,6))
ax1.plot(x, df_th, '-', x, df_num, 'o')
ax1.set_title('Gaussian function first derivative')
ax1.set_xlabel('$x$')
ax1.set_ylabel('$\dfrac{{d} f}{{d} x}$',rotation='horizontal', ha='right')
ax2.set_title('Error')
ax2.semilogy(alln, err)
ax2.set_xlabel('$N$')
ax2.set_ylabel('$\epsilon$',rotation='horizontal', ha='right');
import matplotlib.pyplot as plt
#
alln = range(5,31) # mesh sizes
err = np.zeros(len(alln))
#
for i, n in enumerate(alln):
SpOp = sp.ChebCollocation(n)
x = SpOp.x # collocation points
df_th = df(x)
df_num = SpOp.matder(1) @ f(x)
err[i] = np.sqrt(np.sum((df_th-df_num)**2)/n )
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,6))
ax1.plot(x, df_th, '-', x, df_num, 'o')
ax1.set_title('Gaussian function first derivative')
ax1.set_xlabel('$x$')
ax1.set_ylabel('$\dfrac{{d} f}{{d} x}$',rotation='horizontal', ha='right')
ax2.set_title('Error')
ax2.semilogy(alln, err)
ax2.set_xlabel('$N$')
ax2.set_ylabel('$\epsilon$',rotation='horizontal', ha='right');