POD of time signal

import numpy as np
import scipy.linalg as splin
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams["figure.dpi"] = "100"
plt.rcParams["animation.html"] = "html5"
import time

# length, speed, time
x0 = .2 ; length = 1. ; k = 5. ; speed = 10. ; tf = .5 
n = 1000 # number of points
m = 200  # number of snapshots

def f_steadypulse(x, t):
    return 1+np.exp(-20*(x-x0)**2)*np.cos(2*np.pi*k*t)

def f_oscpulse(x, t):
    return 1+np.exp(-20*(x-length/2)**2)*np.sin(2*np.pi*k*(x-length/2+x0*np.sin(5*k*t)))

def f_convpulse(x, t):
    return 1+np.exp(-20*(x-x0-speed*t)**2)

def f_pulse_per(x, t):
    return 1+np.exp(-10*((x-x0-speed*t))**2)

def f_per(x, t):
    return 1+np.exp(np.sin(2*k/length*np.pi*(x-x0-speed*t)-1.)) 

xx = np.linspace(0., length, n)
tt = np.linspace(0., tf, m)

# build data
snapshots = np.zeros((n,m))
for j,t in enumerate(tt):
    snapshots[:,j] = f_oscpulse(xx, t)

import matplotlib.animation as animation
fig, ax = plt.subplots()
line, = ax.plot(xx, snapshots[:,0])
def animate(i):
    line.set_ydata(snapshots[:,i])  # update the data.
    return line,
ani = animation.FuncAnimation(fig, animate, interval=25, blit=True)
ani

/tmp/ipykernel_2147/2954509547.py:43: UserWarning: frames=None which we can infer the length of, did not pass an explicit *save_count* and passed cache_frame_data=True.  To avoid a possibly unbounded cache, frame data caching has been disabled. To suppress this warning either pass `cache_frame_data=False` or `save_count=MAX_FRAMES`.
  ani = animation.FuncAnimation(fig, animate, interval=25, blit=True)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
File ~/work/numerical-analysis/numerical-analysis/venv/lib/python3.11/site-packages/IPython/core/formatters.py:406, in BaseFormatter.__call__(self, obj)
    404     method = get_real_method(obj, self.print_method)
    405     if method is not None:
--> 406         return method()
    407     return None
    408 else:

File ~/work/numerical-analysis/numerical-analysis/venv/lib/python3.11/site-packages/matplotlib/animation.py:1385, in Animation._repr_html_(self)
   1383 fmt = mpl.rcParams['animation.html']
   1384 if fmt == 'html5':
-> 1385     return self.to_html5_video()
   1386 elif fmt == 'jshtml':
   1387     return self.to_jshtml()

File ~/work/numerical-analysis/numerical-analysis/venv/lib/python3.11/site-packages/matplotlib/animation.py:1302, in Animation.to_html5_video(self, embed_limit)
   1299 path = Path(tmpdir, "temp.m4v")
   1300 # We create a writer manually so that we can get the
   1301 # appropriate size for the tag
-> 1302 Writer = writers[mpl.rcParams['animation.writer']]
   1303 writer = Writer(codec='h264',
   1304                 bitrate=mpl.rcParams['animation.bitrate'],
   1305                 fps=1000. / self._interval)
   1306 self.save(str(path), writer=writer)

File ~/work/numerical-analysis/numerical-analysis/venv/lib/python3.11/site-packages/matplotlib/animation.py:121, in MovieWriterRegistry.__getitem__(self, name)
    119 if self.is_available(name):
    120     return self._registered[name]
--> 121 raise RuntimeError(f"Requested MovieWriter ({name}) not available")

RuntimeError: Requested MovieWriter (ffmpeg) not available

<matplotlib.animation.FuncAnimation at 0x7fe6dcb86f10>

def pod_svd(snapshots, maxmode=-1):
    # SVD decomposition of A = U.S.V 
    modes, sval, tcoefs = splin.svd(snapshots)  # returns U, S, V
    return modes[:,:maxmode], sval[:maxmode], tcoefs.T[:,:maxmode] # time evolution must be transposed to be columnwise

def pod_cov(snapshots, maxmode=-1):
    # SVD decomposition of A = U.S.V through eigenvalue problem AT.A
    m = snapshots.shape[1]
    cov = np.matmul(snapshots.T, snapshots)    # compute covariance matrice cov = AT.A = VT.S2.V
    sqval, rmode = splin.eigh(cov)             # eigenvalue of symmetric problem cov = L.D.R
    # ineg = np.argwhere(sqval <0.)
    # if ineg.size > 0:
    #   print(f"some bad values (negative square): {np.squeeze(ineg)}")
    idex = np.flip(np.argsort(sqval))          # compute (increasing) order and reverse 
    sval = np.sqrt(abs(sqval[idex]))           # singular values (S) are square root of D (and sort)
    rmode = rmode[:,idex]                      # R = V (and sort)
    modes = np.matmul(snapshots, rmode.T/sval) # U.S = A.VT 
    return modes[:,:maxmode], sval[:maxmode], rmode[:,:maxmode]

for thispod in [ pod_svd, pod_cov]:
    method = thispod.__name__
    tic = time.perf_counter()
    modes, sval, tcoefs = thispod(snapshots-np.average(snapshots), maxmode=40)
    toc = time.perf_counter()
    print(sval[:10])
    print(f"process time of pod {method} is {toc-tic:.3e} seconds")
    fig, ax = plt.subplots(4, 2, figsize=(10,12))
    fig.suptitle(method)
    ax[0,0].semilogy(sval, 'o',markersize=2)
    ax[0,1].plot(xx, snapshots[:,0])
    for i in range(3):
        ax[i+1,0].plot(xx, modes[:,i])
        ax[i+1,1].plot(tt, tcoefs[:,i])
    plt.show()
    #plt.hold()

[1.27279797e+02 1.08612071e+02 5.04089089e-06 4.47469002e-14
 4.22158999e-14 3.21870389e-14 2.79608020e-14 2.40106483e-14
 1.81826533e-14 1.53412180e-14]
process time of pod pod_svd is 3.958e-02 seconds

[1.27279797e+02 1.08612071e+02 4.99915635e-06 1.39084970e-06
 1.13254048e-06 7.56902547e-07 6.26975264e-07 5.72921764e-07
 5.42697881e-07 5.09152833e-07]
process time of pod pod_cov is 1.211e-01 seconds