## effective increment

### created by Yuying Liu, 04/29/2020

For autonomous system 
\begin{equation}
    \bf{\dot{x} = f(x)}
\end{equation}

we visualize 
\begin{equation}
    F_{\Delta}(x(t)):= x(t+\Delta)-x(t)
\end{equation}

This script can be used to visualize Figure 10.


```python
%matplotlib inline

import numpy as np
from scipy import integrate

from tqdm.notebook import tqdm
from matplotlib import pyplot as plt
from matplotlib import gridspec
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import cnames
from matplotlib import animation
from ipywidgets import interact, interactive
from IPython.display import clear_output, display, HTML
```


```python
# shared parameters
max_time = 10.24
dt = 0.16
t = np.linspace(0, max_time, int(max_time/dt)+1)
```

### Hyperbolic Fixed Point

\begin{split}
    \dot{x} &= \mu x \\
    \dot{y} &= \lambda(y-x^2)     
\end{split}


```python
mu = -0.05
lam = -1.0
def hyperbolic_deriv(x_y):
        """Compute the time-derivative of the hyperbolic example."""
        x, y = x_y
        return [mu*x, lam*(y-x**2)]
    
# simulate
xvalues, yvalues = np.meshgrid(np.arange(-1.0,1.0,0.02), np.arange(-1.0,1.0,0.02))
inits = np.stack([xvalues, yvalues], 2)
sols = np.zeros((inits.shape[0], inits.shape[1], len(t), 2))
for i in tqdm(range(inits.shape[0])):
    for j in range(inits.shape[1]):
        init = inits[i, j]
        sol = integrate.solve_ivp(lambda _, x: hyperbolic_deriv(x), [0, max_time], init, t_eval=t)
        sols[i, j, :, :] = sol.y.T

# compute increments
vmin0 = [float('inf')]*2
vmax0 = [float('-inf')]*2
incre0 = np.zeros((inits.shape[0], inits.shape[1], len(t)-1, 2))
for i in tqdm(range(1, len(t))):
    for j in range(2):
        vmin0[j] = min(np.min(sols[:, :, i, j] - sols[:, :, 0, j]), vmin0[j])
        vmax0[j] = max(np.max(sols[:, :, i, j] - sols[:, :, 0, j]), vmax0[j])
    incre0[:, :, i-1, :] = sols[:, :, i, :] - sols[:, :, 0, :]

def viz_cubic_increment(time=25):
    fig = plt.figure(figsize=(10, 16))
    gs = gridspec.GridSpec(nrows=2, ncols=1, hspace=0.2)
    ax0 = fig.add_subplot(gs[0, :])
    ax1 = fig.add_subplot(gs[1, :])

    # prepare the axes limits
    ax0.set_xlim((-2, 2))
    ax0.set_ylim((-2, 2))
    ax1.set_xlim((-2, 2))
    ax1.set_ylim((-2, 2))

    # plot
    idx = int(time/dt)-1
    im0 = ax0.imshow(incre0[:, :, idx, 0], extent=[-2,2,-2,2], vmin=vmin0[0], vmax=vmax0[1])
    im1 = ax1.imshow(incre0[:, :, idx, 1], extent=[-2,2,-2,2], vmin=vmin0[0], vmax=vmax0[1])
    
    # colorbar
    cbar0 = fig.colorbar(im0, ax=ax0, aspect=5)
    cbar1 = fig.colorbar(im1, ax=ax1, aspect=5)
    cbar0.ax.tick_params(labelsize=60)
    cbar1.ax.tick_params(labelsize=60)

    # title
    ax0.set_title('x', fontsize=50)
    ax1.set_title('y', fontsize=50)
    
    # ticks
    ax0.tick_params(axis='both', which='major', labelsize=60)
    ax1.tick_params(axis='both', which='major', labelsize=60)
    ax0.axis('off')
    ax1.axis('off')
    
    return t, fig

# viz
w0 = interactive(viz_cubic_increment, time=(dt, max_time, dt))
display(w0)
```


    HBox(children=(FloatProgress(value=0.0), HTML(value='')))


    



    HBox(children=(FloatProgress(value=0.0, max=64.0), HTML(value='')))


    



    interactive(children=(FloatSlider(value=10.24, description='time', max=10.24, min=0.16, step=0.16), Output()),…


### Cubic Oscillator

\begin{split}
    \dot{x} &= -0.1x^3 + 2y^3 \\
    \dot{y} &= -2x^3 - 0.1y^3
\end{split}


```python
def cubic_deriv(x_y):
        """Compute the time-derivative of the cubic oscillator."""
        x, y = x_y
        return [-0.1*x**3+2*y**3, -2*x**3-0.1*y**3]
    
# simulate
xvalues, yvalues = np.meshgrid(np.arange(-1.0,1.0,0.02), np.arange(-1.0,1.0,0.02))
inits = np.stack([xvalues, yvalues], 2)
sols = np.zeros((inits.shape[0], inits.shape[1], len(t), 2))
for i in tqdm(range(inits.shape[0])):
    for j in range(inits.shape[1]):
        init = inits[i, j]
        sol = integrate.solve_ivp(lambda _, x: cubic_deriv(x), [0, max_time], init, t_eval=t)
        sols[i, j, :, :] = sol.y.T

# compute increments
vmin1 = [float('inf')]*2
vmax1 = [float('-inf')]*2
incre1 = np.zeros((inits.shape[0], inits.shape[1], len(t)-1, 2))
for i in tqdm(range(1, len(t))):
    for j in range(2):
        vmin1[j] = min(np.min(sols[:, :, i, j] - sols[:, :, 0, j]), vmin1[j])
        vmax1[j] = max(np.max(sols[:, :, i, j] - sols[:, :, 0, j]), vmax1[j])
    incre1[:, :, i-1, :] = sols[:, :, i, :] - sols[:, :, 0, :]

def viz_cubic_increment(time=25):
    fig = plt.figure(figsize=(10, 16))
    gs = gridspec.GridSpec(nrows=2, ncols=1, hspace=0.2)
    ax0 = fig.add_subplot(gs[0, :])
    ax1 = fig.add_subplot(gs[1, :])
    
    # prepare the axes limits
    ax0.set_xlim((-1, 1))
    ax0.set_ylim((-1, 1))
    ax1.set_xlim((-1, 1))
    ax1.set_ylim((-1, 1))

    # plot
    idx = int(time/dt)-1
    im0 = ax0.imshow(incre1[:, :, idx, 0], extent=[-1,1,-1,1], vmin=vmin1[0], vmax=vmax1[1])
    im1 = ax1.imshow(incre1[:, :, idx, 1], extent=[-1,1,-1,1], vmin=vmin1[0], vmax=vmax1[1])
    
    # colorbar
    cbar0 = fig.colorbar(im0, ax=ax0, aspect=5)
    cbar1 = fig.colorbar(im1, ax=ax1, aspect=5)
    cbar0.ax.tick_params(labelsize=60)
    cbar1.ax.tick_params(labelsize=60)

    # title
    ax0.set_title('x', fontsize=50)
    ax1.set_title('y', fontsize=50)
    
    # ticks
    ax0.tick_params(axis='both', which='major', labelsize=60)
    ax1.tick_params(axis='both', which='major', labelsize=60)
    ax0.axis('off')
    ax1.axis('off')
    
    return t, fig

# viz
w1 = interactive(viz_cubic_increment, time=(dt, max_time, dt))
display(w1)
```

### Van der Pol Oscillator

\begin{split}
    \dot{x} &= y \\
    \dot{y} &= \mu(1-x^2)y - x   
\end{split}

where $\mu=2.0$


```python
mu = 2.0
def vdp_deriv(x_y):
        """Compute the time-derivative of the vdp."""
        x, y = x_y
        return [y, mu*(1-x**2)*y-x]
    
# simulate
xvalues, yvalues = np.meshgrid(np.arange(-2.0,2.0,0.02), np.arange(-1.0,1.0,0.01))
inits = np.stack([xvalues, yvalues], 2)
sols = np.zeros((inits.shape[0], inits.shape[1], len(t), 2))
for i in tqdm(range(inits.shape[0])):
    for j in range(inits.shape[1]):
        init = inits[i, j]
        sol = integrate.solve_ivp(lambda _, x: vdp_deriv(x), [0, max_time], init, t_eval=t)
        sols[i, j, :, :] = sol.y.T

# compute increments
vmin2 = [float('inf')]*2
vmax2 = [float('-inf')]*2
incre2 = np.zeros((inits.shape[0], inits.shape[1], len(t)-1, 2))
for i in tqdm(range(1, len(t))):
    for j in range(2):
        vmin2[j] = min(np.min(sols[:, :, i, j] - sols[:, :, 0, j]), vmin2[j])
        vmax2[j] = max(np.max(sols[:, :, i, j] - sols[:, :, 0, j]), vmax2[j])
    incre2[:, :, i-1, :] = sols[:, :, i, :] - sols[:, :, 0, :]

def viz_vdp_increment(time=25):
    fig = plt.figure(figsize=(10, 16))
    gs = gridspec.GridSpec(nrows=2, ncols=1, hspace=0.2)
    ax0 = fig.add_subplot(gs[0, :])
    ax1 = fig.add_subplot(gs[1, :])

    # prepare the axes limits
    ax0.set_xlim((-2, 2))
    ax0.set_ylim((-1, 1))
    ax1.set_xlim((-2, 2))
    ax1.set_ylim((-1, 1))

    # plot
    idx = int(time/dt)-1
    im0 = ax0.imshow(incre2[:, :, idx, 0], extent=[-2,2,-1,1], vmin=vmin2[0], vmax=vmax2[0], aspect=2)
    im1 = ax1.imshow(incre2[:, :, idx, 1], extent=[-2,2,-1,1], vmin=vmin2[1], vmax=vmax2[1], aspect=2)
    
    # colorbar
    cbar0 = fig.colorbar(im0, ax=ax0, aspect=5)
    cbar1 = fig.colorbar(im1, ax=ax1, aspect=5)
    cbar0.ax.tick_params(labelsize=60)
    cbar1.ax.tick_params(labelsize=60)
    
    # title
    ax0.set_title('x', fontsize=50)
    ax1.set_title('y', fontsize=50)
    
    # ticks
    ax0.tick_params(axis='both', which='major', labelsize=60)
    ax1.tick_params(axis='both', which='major', labelsize=60)
    ax0.axis('off')
    ax1.axis('off')

    return t, fig

# viz
w2 = interactive(viz_vdp_increment, time=(dt, max_time, dt))
display(w2)
```

### Hopf bifurcation

\begin{split}
    \dot{\mu} &= 0 \\
    \dot{x} &= \mu x + y -x(x^2+y^2) \\
    \dot{y} &= \mu y - x -y(x^2+y^2)
\end{split}


```python
def hopf_deriv(mu_x_y):
    """Compute the time-derivative of the hopf"""
    mu, x, y = mu_x_y
    return np.array([0, mu*x+y-x*(x**2+y**2), -x+mu*y-y*(x**2+y**2)])
    
# simulate
muvalues, xvalues, yvalues = np.meshgrid(np.arange(-0.2,0.6,0.1), np.arange(-1,2,0.02), np.arange(-1,1,0.02), indexing='ij')
inits = np.stack([muvalues, xvalues, yvalues], 3)
sols = np.zeros((inits.shape[0], inits.shape[1], inits.shape[2], len(t), 3))
for i in tqdm(range(inits.shape[0])):
    for j in range(inits.shape[1]):
        for k in range(inits.shape[2]):
            init = inits[i, j, k, :]
            sol = integrate.solve_ivp(lambda _, x: hopf_deriv(x), [0, max_time], init, t_eval=t)
            sols[i, j, k, :, :] = sol.y.T

# compute increments
vmin3 = [float('inf')]*3
vmax3 = [float('-inf')]*3
incre3 = np.zeros((inits.shape[0], inits.shape[1], inits.shape[2], len(t)-1, 3))
for i in tqdm(range(1, len(t))):
    for j in range(3):
        vmin3[j] = min(np.min(sols[:, :, :, i, j] - sols[:, :, :, 0, j]), vmin3[j])
        vmax3[j] = max(np.max(sols[:, :, :, i, j] - sols[:, :, :, 0, j]), vmax3[j])
    incre3[:, :, :, i-1, :] = sols[:, :, :, i, :] - sols[:, :, :, 0, :]

def viz_hopf_increment(time=25, mu=0):
    fig = plt.figure(figsize=(10, 24))
    gs = gridspec.GridSpec(nrows=3, ncols=1, hspace=0.2)
    ax0 = fig.add_subplot(gs[0, :])
    ax1 = fig.add_subplot(gs[1, :])
    ax2 = fig.add_subplot(gs[2, :])

    # prepare the axes limits
    ax0.set_xlim((-1, 2))
    ax0.set_ylim((-1, 1))
    ax1.set_xlim((-1, 2))
    ax1.set_ylim((-1, 1))
    ax2.set_xlim((-1, 2))
    ax2.set_ylim((-1, 1))

    # plot
    idx_t = int(time/dt)-1
    idx_mu = min(int((mu + 0.2) / 0.1), len(np.arange(-0.2,0.6,0.1))-1)
    im0 = ax0.imshow(incre3[:, :, idx_mu, idx_t, 0], extent=[-1,2,-1,1], vmin=vmin3[0], vmax=vmax3[0], aspect=1.5)
    im1 = ax1.imshow(incre3[:, :, idx_mu, idx_t, 1], extent=[-1,2,-1,1], vmin=vmin3[1], vmax=vmax3[1], aspect=1.5)
    im2 = ax2.imshow(incre3[:, :, idx_mu, idx_t, 2], extent=[-1,2,-1,1], vmin=vmin3[2], vmax=vmax3[2], aspect=1.5)
    
    # colorbar
    cbar0 = fig.colorbar(im0, ax=ax0, aspect=5)
    cbar1 = fig.colorbar(im1, ax=ax1, aspect=5)
    cbar2 = fig.colorbar(im2, ax=ax2, aspect=5)
    cbar0.ax.tick_params(labelsize=60)
    cbar1.ax.tick_params(labelsize=60)
    cbar2.ax.tick_params(labelsize=60)
    
    # title
    ax0.set_title('$\mu$', fontsize=50)
    ax1.set_title('x', fontsize=50)
    ax2.set_title('y', fontsize=50)
    
    # ticks
    ax0.tick_params(axis='both', which='major', labelsize=60)
    ax1.tick_params(axis='both', which='major', labelsize=60)
    ax2.tick_params(axis='both', which='major', labelsize=60)
    ax0.axis('off')
    ax1.axis('off')
    ax2.axis('off')

    return t, fig

# viz
w3 = interactive(viz_hopf_increment, time=(dt, max_time, dt), mu=(-0.2,0.6,0.1))
display(w3)
```

### Lorenz

\begin{split}
    \dot{x} &= \sigma(y-x) \\
    \dot{y} &= x(\rho-z)-y \\
    \dot{z} &= xy - \beta z    
\end{split}

where $\sigma=10, \rho=28, \beta=8/3$


```python
sigma = 10
rho = 28
beta = 8/3
def lorenz_deriv(x_y_z):
        """Compute the time-derivative of the lorenz system."""
        x, y, z = x_y_z
        return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]
    
# simulate
xvalues, yvalues, zvalues = np.meshgrid(np.arange(-9,-7,0.02), np.arange(6,8,0.02), np.arange(26,28,0.5), indexing='ij')
inits = np.stack([xvalues, yvalues, zvalues], 3)
sols = np.zeros((inits.shape[0], inits.shape[1], inits.shape[2], len(t), 3))
for i in tqdm(range(inits.shape[0])):
    for j in range(inits.shape[1]):
        for k in range(inits.shape[2]):
            init = inits[i, j, k, :]
            sol = integrate.solve_ivp(lambda _, x: lorenz_deriv(x), [0, max_time], init, t_eval=t)
            sols[i, j, k, :, :] = sol.y.T

# compute increments
vmin4 = [float('inf')]*3
vmax4 = [float('-inf')]*3
incre4 = np.zeros((inits.shape[0], inits.shape[1], inits.shape[2], len(t)-1, 3))
for i in tqdm(range(1, len(t))):
    for j in range(3):
        vmin4[j] = min(np.min(sols[:, :, :, i, j] - sols[:, :, :, 0, j]), vmin4[j])
        vmax4[j] = max(np.max(sols[:, :, :, i, j] - sols[:, :, :, 0, j]), vmax4[j])
    incre4[:, :, :, i-1, :] = sols[:, :, :, i, :] - sols[:, :, :, 0, :]

def viz_lorenz_increment(time=5, z=0):
    fig = plt.figure(figsize=(10, 24))
    gs = gridspec.GridSpec(nrows=3, ncols=1, hspace=0.2)
    ax0 = fig.add_subplot(gs[0, :])
    ax1 = fig.add_subplot(gs[1, :])
    ax2 = fig.add_subplot(gs[2, :])

    # prepare the axes limits
    ax0.set_xlim((-9, -7))
    ax0.set_ylim((6, 8))
    ax1.set_xlim((-9, -7))
    ax1.set_ylim((6, 8))
    ax2.set_xlim((-9, -7))
    ax2.set_ylim((6, 8))

    # plot
    idx_t = int(time/dt)-1
    idx_z = min(int((z - 26)/0.5), len(np.arange(26,28,0.1))-1)
    im0 = ax0.imshow(incre4[:, :, idx_z, idx_t, 0], extent=[-9,-7,6,8], vmin=vmin4[0], vmax=vmax4[0])
    im1 = ax1.imshow(incre4[:, :, idx_z, idx_t, 1], extent=[-9,-7,6,8], vmin=vmin4[1], vmax=vmax4[1])
    im2 = ax2.imshow(incre4[:, :, idx_z, idx_t, 2], extent=[-9,-7,6,8], vmin=vmin4[2], vmax=vmax4[2])
    
    # colorbar
    cbar0 = fig.colorbar(im0, ax=ax0, aspect=5)
    cbar1 = fig.colorbar(im1, ax=ax1, aspect=5)
    cbar2 = fig.colorbar(im2, ax=ax2, aspect=5)
    cbar0.ax.tick_params(labelsize=60)
    cbar1.ax.tick_params(labelsize=60)
    cbar2.ax.tick_params(labelsize=60)
    
    # title
    ax0.set_title('x', fontsize=50)
    ax1.set_title('y', fontsize=50)
    ax2.set_title('z', fontsize=50)
    
    # ticks
    ax0.tick_params(axis='both', which='major', labelsize=60)
    ax1.tick_params(axis='both', which='major', labelsize=60)
    ax2.tick_params(axis='both', which='major', labelsize=60)
    ax0.axis('off')
    ax1.axis('off')
    ax2.axis('off')

    return t, fig

# viz
w4 = interactive(viz_lorenz_increment, time=(dt, max_time, dt), z=(26, 28, 0.1))
display(w4)
```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```


```python

```

b9cce0ff4682f181b4828f2229f14204cb42f096
ipynb
Jupyter Notebook
scripts/others/visualize_increment.ipynb
luckystarufo/multiscale_HiTS
d1bb69114b384994a4493bede0d38deb9f4740ec
MIT
2020-09-09T05:38:21.000Z
2022-03-12T04:07:42.000Z
scripts/others/visualize_increment.ipynb
luckystarufo/multiscale_HiTS
d1bb69114b384994a4493bede0d38deb9f4740ec
MIT
scripts/others/visualize_increment.ipynb
luckystarufo/multiscale_HiTS
d1bb69114b384994a4493bede0d38deb9f4740ec
MIT
2021-01-19T08:52:02.000Z
2021-11-15T19:30:36.000Z
Qwen/Qwen-72B
1. YES
2. YES