<jupyter_start><jupyter_text>This notebook explores the instability in the filter.I use an exact (spectral) discrete Laplacian on a 1D periodic domain with 256 points.The filter is not generated from the code; instead, I just set it up by specifying the roots `s[k]`.I'll show a range of results before getting to the point.<jupyter_code>import numpy as np
from scipy import fft
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>The block below sets up wavenumbers `k`, the Laplacian in Fourier space `L`, and the data that I will apply the filter to `f0`. The data is constructed so that the DFT of the data is 1 for all wavenumbers. Once we've applied the filter, if we compute and plot the absolute value of the FFT of `f` then what we see is the effective filter shape.<jupyter_code>%matplotlib inline
k = np.concatenate([np.arange(0,129), np.arange(-127,0)])
L = - k**2
fHat = np.exp(1j*np.random.uniform(0,2*np.pi,256))
f0 = fft.ifft(fHat)<jupyter_output><empty_output><jupyter_text>First show that the problem does not depend on `n_steps` by setting up a filter with a large number of steps that has no problem.<jupyter_code>n_steps = 128
s = (np.arange(64,n_steps+64))**2
f = f0
for i in range(n_steps):
    fHat = fft.fft(f)
    f = fft.ifft(fHat + (1/s[i])*L*fHat)

plt.plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>Next show that the problem begins to arise as you include more small values of `s[i]`. The plot below shows we're still OK at min `s[i]` = 54<jupyter_code>n_steps = 128
s = (np.arange(54,n_steps+54))**2
f = f0
for i in range(n_steps):
    fHat = fft.fft(f)
    f = fft.ifft(fHat + (1/s[i])*L*fHat)

plt.plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>The plot below shows we have instability at min `s[i]=44`.<jupyter_code>n_steps = 128
s = (np.arange(44,n_steps+44))**2
f = f0
for i in range(n_steps):
    fHat = fft.fft(f)
    f = fft.ifft(fHat + (1/s[i])*L*fHat)

plt.plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>In the plot above there's no problem at small scales, but the large scales are wrong. The plot below shows that if you reverse the order of the roots `s[i]` then it still blows up, but now at the small scales:<jupyter_code>n_steps = 128
s = -np.sort(-(np.arange(44,n_steps+44))**2)
f = f0
for i in range(n_steps):
    fHat = fft.fft(f)
    f = fft.ifft(fHat + (1/s[i])*L*fHat)

plt.plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>The sequence of figures below shows how the solution evolves through the internal steps of the filter. When the roots are sorted in descending order we see good behavior until right at the end.<jupyter_code>s = -np.sort(-(np.arange(44,n_steps+44))**2)
fig, axarr = plt.subplots(ncols=1, nrows=5, figsize=[11, 7])
for j in range(1,6):
    f = f0
    for i in range(25*j):
        fHat = fft.fft(f)
        f = fft.ifft(fHat + (1/s[i])*L*fHat)
    axarr[j-1].plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>The plot below also shows the evolution of the solution through the internal filter stages, with the roots sorted ascending. It blows up extremely quickly, starting at the small scales.<jupyter_code>s = (np.arange(44,n_steps+44))**2
fig, axarr = plt.subplots(ncols=1, nrows=5, figsize=[11, 7])
for j in range(1,6):
    f = f0
    for i in range(2*j):
        fHat = fft.fft(f)
        f = fft.ifft(fHat + (1/s[i])*L*fHat)
    axarr[j-1].plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output><jupyter_text>At each stage we multiply in Fourier space by $(1-k^2/s_i)$. If $s_i$ is small then this will amplify small scales. Conversely, for large $s_i$ we are damping small scales. We can test whether any individual filter stage will amplify scales that are present in the solution as follows. If $$\left|1 - \frac{k_{\text{max}}^2}{s_i}\right|\le 1$$then there is no amplification, otherwise there is amplification.The problem with the above filter is that the first several stages all amplify the small scales to the point where by the time the damping stages begin we've already started to lose numerical precision.(When the roots are reversed, the small scales are immediately damped almost to zero, but not quite because of roundoff. When the later stages arrive they amplify any roundoff present at small scales.)This suggests that we can postpone the instability by interleaving steps that amplify the small scales and steps that damp them. The following code bears this out.I start by sorting the roots from most-damping to most-amplifying. Then I split the list of roots into two lists: one damping and one amplifying.Then I interleave the two lists (using some code from [here](https://stackoverflow.com/questions/48199961/how-to-interleave-two-lists-of-different-length/48200036)) so that the most-damping stage is first, then the most-amplifying, etc.<jupyter_code>from itertools import chain, zip_longest
def twolists(l1, l2):
    return [x for x in chain(*zip_longest(l1, l2)) if x is not None]
s = (np.arange(44,n_steps+44))**2
ind = np.argsort(np.abs(1 - 128**2/s))
s = s[ind] # sorted from most damping to most amplifying
indD = np.argwhere(np.abs(1 - 128**2/s)<=1)
indA = np.argwhere(np.abs(1 - 128**2/s)>1)
sD = s[indD].tolist() # Damping roots, sorted most to least damping
sA = np.flip(s[indA]).tolist() # Amplifying roots, sorted most to least amplifying
s = twolists(sD,sA)
s = np.array([y for x in s for y in x])
fig, axarr = plt.subplots(ncols=1, nrows=7, figsize=[11, 7])
for j in range(1,7):
    f = f0
    for i in range(8*j):
        fHat = fft.fft(f)
        f = fft.ifft(fHat + (1/s[i])*L*fHat)
    axarr[j-1].plot(k,np.abs(fft.fft(f)))

f = f0
for i in range(n_steps):
    fHat = fft.fft(f)
    f = fft.ifft(fHat + (1/s[i])*L*fHat)
axarr[6].plot(k,np.abs(fft.fft(f)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12437981180496151, 'fraction_numerical': 0.02737382378100941, 'mean_word_length': 4.562321598477641, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '236801', 'license': 'MIT', 'n_tokens_mistral': 2283, 'n_tokens_neox': 2156, 'n_words': 703, 'path': 'scratch_roundoff.ipynb', 'repo_name': 'iangrooms/CPT-Snippets', 'size': '5845'}