---
jupytext:
  cell_metadata_filter: -all
  formats: md:myst,ipynb
  main_language: python
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.17.0
kernelspec:
  name: python3
  display_name: Python 3 (ipykernel)
  language: python
---

```{code-cell} ipython3
import mmf_setup; mmf_setup.nbinit()
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

```{code-cell} ipython3
%load_ext autoreload
%autoreload 2
from importlib import reload
from gpe.contexts import FPS
```

```{code-cell} ipython3
import time
import viscosity
f = viscosity.Flow(n_min=0)
n0 = np.where(f.basis.x<0, 1.0, 0.1)
n1 = np.where(f.basis.x<0, 1.0, 0.01)
y0 = f.pack(n0, 0*n0)
y1 = f.pack(n1, 0*n1)
dt = 0.01
t = 0
nfevs0 = []
nfevs1 = []
fig, axs = fig_axs = f.plot(y0)
for n in FPS(50, fig=fig, display=True):
    y0 = f.evolve(y0, t=t, dt=dt)
    nfevs0.append(f._res.nfev)
    y1 = f.evolve(y1, t=t, dt=dt)
    nfevs1.append(f._res.nfev)
    t += dt
    fig_axs = f.plot(y=y0, t=t, fig_axs=fig_axs, c='C0', label=f"$n_R=0.1$: nfev={nfevs0[-1]}")
    fig, axs = f.plot(y=y1, t=t, fig_axs=fig_axs, c='C1', label=f"$n_R=0.01$: nfev={nfevs1[-1]}", cla=False)
    axs[0].legend(loc='center left')
```

```{code-cell} ipython3
import viscosity; reload(viscosity)
f = viscosity.Flow(n_min=0)
n = np.where(f.basis.x<0, 1.0, 0.01)
y = f.pack(n, 0*n)
dt = 0.01
t = 0
nfevs = []
for n in FPS(50, fig=fig, display=True):
    #y, t = f.step_euler(y, t=t, dt=dt), t+dt
    y, t = f.evolve(y, t=t, dt=dt), t+dt
    nfevs.append(f._res.nfev)
    fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
```

# Viscosity 4: Lagrangian

```{code-cell} ipython3
Nx = 128
x0 = np.linspace(-1.0, 1.0, Nx)
dx0 = np.diff(x0).mean()

g_m = 1.0  # Coupling constant g/m
nu = 0.2  # Viscosity coefficient

w = np.sqrt(2)  # Trap frequency
V0_m = 0.5  # Bump height
sigma = 0.1  # Bump width

# Initial density provile 
X0_shift = 0.5  # Initial shift of the cloud
X0 = x0 + X0_shift
n0 = w**2/2/g_m * (1 - x0**2)
n0_x0 = -w**2 / g_m * x0

# Finite difference approximations for the derivatives
o = np.ones(Nx)
Df = (np.diag(o[1:], 1) - np.diag(o, 0))/dx0
Db = (np.diag(o, 0) - np.diag(o[:-1], -1))/dx0
Dc = (Df + Db)/2
D2 = Db @ Df

# Neumann bc's
Dc[0, :2] = [-1/dx0, 1/dx0]
Dc[-1, -2:] = [-1/dx0, 1/dx0]
D2[0, :] = 0
D2[-1, :] = 0

def pack(X, X_t):
    return np.ravel([X, X_t])

def unpack(y):
    return np.reshape(y, (2, Nx))

def unpack_xnu(y):
    X, X_t = unpack(y)
    X_x0 = Dc @ X
    n = n0 / X_x0
    u = X_t
    x = X
    return x, n, u
    
def compute_dy_dt(t, y):
    X, X_t = unpack(y)
    X_x0 = Dc @ X
    X_x0x0 = D2 @ X
    n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3
    
    u_xx = (D2 @ X_t) / X_x0**2 - Dc @ X_t * X_x0x0 / X_x0**3
    a_ext = - w**2 * X
    a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
    
    X_tt = a_ext - g_m * n_x + nu * u_xx
    return pack(X_t, X_tt)

y0 = pack(X0, 0*X0)
X, X_t = unpack(compute_dy_dt(0, y0))
plt.plot(x0, X_t)
```

```{code-cell} ipython3
%%time
from scipy.integrate import solve_ivp

Tw = 2*np.pi / w  # Trapping period
T = 40*Tw

Nt = 500
t_eval = np.linspace(0, T, Nt) 
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
                t_eval=t_eval, method="BDF")
if not res.success:
    print(res.message)

ts = res.t
Nt = len(ts)
X, X_t = res.y.reshape((2, Nx, Nt))
fig, ax = plt.subplots()
ax.plot(ts / Tw, X[::10].T)
ax.set(xlabel="$t/T_w$", ylabel="X");
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp

T = 4.0
Nt = 500
t_eval = np.linspace(0, T, Nt) 
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
                t_eval=t_eval, method="BDF")
if not res.success:
    print(res.message)

ts = res.t
Nt = len(ts)
X, X_t = res.y.reshape((2, Nx, Nt))
plt.plot(ts, X[::10].T);
```

```{code-cell} ipython3
X_x0 = Dc @ X
n = n0[:, None] / X_x0

x_ = np.linspace(X.min(), X.max(), 200)
V_ext = w**2 * X**2 / 2
V_ext += V0_m * np.exp(-X**2/2/sigma**2)
b = V_ext / g_m

b_ = w**2 * x_**2 / 2
b_ += V0_m * np.exp(-x_**2/2/sigma**2)
b_ /= g_m 

h = b + n

fig, ax = plt.subplots()
for i in FPS(np.arange(Nt)[::1], fig=fig, display=True):
    t = ts[i]
    ax.cla()
    #ax.plot(X[:, i], n[:, i])
    ax.plot(x_, b_, 'k-')
    ax.fill_between(X[:, i], b[:, i], h[:, i], color='cyan')
    ax.set(title=f"{t=:.4f}", 
           xlim=(X.min(), X.max()), 
           ylim=(0, h.max()))
    
```

```{code-cell} ipython3
# Spectral version
def s():
    import scipy.fft
    sp = scipy

    Nx = 128
    th0 = (2*np.arange(Nx) + 1) * np.pi / 2 / Nx
    k = np.arange(Nx)
    x0 = np.cos(th0)

    X0_shift = 0.5  # Initial shift of the cloud

    X0 = x0 + X0_shift

    g_m = 1.0  # Coupling constant g/m
    nu = 0.2  # Viscosity coefficient

    w = np.sqrt(2)  # Trap frequency
    V0_m = 0.5  # Bump height
    sigma = 0.1  # Bump width

    # Initial density provile 
    n0 = w**2/2/g_m * (1 - x0**2)
    n0_x0 = -w**2 / g_m * x0

    def diff(f, d=1, _k=k, _th=th0):
        ft = sp.fft.dct(f, type=2, norm='forward')
        # Shift the coefficients and pad with zero.
        df_dth = sp.fft.dst(np.concatenate([-(ft*_k)[1:], [0]]), type=3)
        s = np.sin(_th)
        df_dx = -df_dth / s
        if d == 1:
            return df_dx 
        elif d == 2:
            d2f_dth2 = sp.fft.idct(-ft*_k**2, type=2, norm='forward')
            d2f_dx2 = (d2f_dth2 - df_dth / np.tan(_th)) / s**2
            return d2f_dx2
        return diff(df_dx, d=d-1)

    def pack(X, X_t):
        return np.ravel([X, X_t])

    def unpack(y):
        return np.reshape(y, (2, Nx))

    def unpack_xnu(y):
        X, X_t = unpack(y)
        X_x0 = diff(X)
        n = n0 / X_x0
        u = X_t
        x = X
        return x, n, u

    def compute_dy_dt(t, y):
        X, X_t = unpack(y)
        X_x0 = diff(X)
        X_x0x0 = diff(X, d=2)
        n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3

        u_xx = (diff(X_t, d=2)) / X_x0**2 - diff(X_t) * X_x0x0 / X_x0**3
        a_ext = - w**2 * X
        a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
        X_tt = a_ext - g_m * n_x + nu * u_xx
        return pack(X_t, X_tt)

    y0 = pack(X0, 0*X0)

    return x0, y0, compute_dy_dt, n0_x0
```

```{code-cell} ipython3
# Spectral version
import scipy.fft
sp = scipy

Nx = 128
th0 = (2*np.arange(Nx) + 1) * np.pi / 2 / Nx
k = np.arange(Nx)
x0 = np.cos(th0)

X0_shift = 0.5  # Initial shift of the cloud

X0 = x0 + X0_shift

g_m = 1.0  # Coupling constant g/m
nu = 0.2  # Viscosity coefficient

w = np.sqrt(2)  # Trap frequency
V0_m = 0.5  # Bump height
sigma = 0.1  # Bump width

# Initial density provile 
n0 = w**2/2/g_m * (1 - x0**2)
n0_x0 = -w**2 / g_m * x0

def diff(f, d=1, _k=k, _th=th0):
    ft = sp.fft.dct(f, type=2, norm='forward')
    # Shift the coefficients and pad with zero.
    df_dth = sp.fft.dst(np.concatenate([-(ft*_k)[1:], [0]]), type=3)
    s = np.sin(_th)
    df_dx = -df_dth / s
    if d == 1:
        return df_dx 
    elif d == 2:
        d2f_dth2 = sp.fft.idct(-ft*_k**2, type=2, norm='forward')
        d2f_dx2 = (d2f_dth2 - df_dth / np.tan(_th)) / s**2
        return d2f_dx2
    return diff(df_dx, d=d-1)

def pack(X, X_t):
    return np.ravel([X, X_t])

def unpack(y):
    return np.reshape(y, (2, Nx))

def unpack_xnu(y):
    X, X_t = unpack(y)
    X_x0 = diff(X)
    n = n0 / X_x0
    u = X_t
    x = X
    return x, n, u

def compute_dy_dt(t, y):
    X, X_t = unpack(y)
    X_x0 = diff(X)
    X_x0x0 = diff(X, d=2)
    n_x = n0_x0 / X_x0**2 - n0 * X_x0x0 / X_x0**3

    u_xx = (diff(X_t, d=2)) / X_x0**2 - diff(X_t) * X_x0x0 / X_x0**3
    a_ext = - w**2 * X
    a_ext += V0_m * np.exp(-X**2/2/sigma**2) * X/sigma**2
    X_tt = a_ext - g_m * n_x + nu * u_xx
    return pack(X_t, X_tt)

y0 = pack(X0, 0*X0)
```

```{code-cell} ipython3
x0_, y0_, compute_dy_dt_, n0_x0_ = s()

dX, dX_t = unpack(compute_dy_dt(0, y0))
X, X_t = unpack(y0)
plt.plot(x0, X)
plt.plot(x0, dX_t)

dX_, dX_t_ = unpack(compute_dy_dt_(0, y0_))
X_, X_t_ = unpack(y0_)
plt.plot(x0_, X_)
plt.plot(x0_, dX_t_)
```

```{code-cell} ipython3
res_ = solve_ivp(compute_dy_dt_, y0=y0_, t_span=(0, 4), method="BDF") 
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 4), method="BDF")
plt.plot(x0, unpack(res.y[:, -1])[0])
plt.plot(x0_, unpack(res_.y[:, -1])[0])
```

```{code-cell} ipython3
%%time
from scipy.integrate import solve_ivp

Tw = 2*np.pi / w  # Trapping period
T = Tw

Nt = 500
t_eval = np.linspace(0, T, Nt)
ress = []
for V0_m in [0.0, 0.5]:
    res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T),
                    t_eval=t_eval, method="BDF")
    if not res.success:
        print(res.message)
    ress.append((V0_m, res))

ts = res.t
Nt = len(ts)
fig, axs = plt.subplots(2, 1, figsize=(6, 3), sharex=True, gridspec_kw=dict(hspace=0))
for ax, (V0_m, res) in zip(axs, ress):
    X, X_t = res.y.reshape((2, Nx, Nt))
    ax.plot(ts / Tw, X[::10].T)
    ax.set(ylabel=f"$X$ ($V_0/m = {V0_m:.1f})$")
axs[1].set(xlabel="$t/T_w$");
```

```{code-cell} ipython3
y = f.get_ground_state(x0=0.0)
fig_axs = fig, axs = f.plot(y)

dt = 0.01
t = 0
for n in FPS(200, fig=fig, display=True):
    #y, t = f.step_euler(y, t=t, dt=dt), t+dt
    y, t = f.evolve(y, t=t, dt=dt), t+dt
    fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
```

```{code-cell} ipython3
import viscosity; reload(viscosity)
f = viscosity.FlowTrap(n_min=0.1)
y = f.get_ground_state(x0=0.0)
fig_axs = fig, axs = f.plot(y)

dt = 0.01
t = 0
for n in FPS(200, fig=fig, display=True):
    #y, t = f.step_euler(y, t=t, dt=dt), t+dt
    y, t = f.evolve(y, t=t, dt=dt), t+dt
    fig_axs = f.plot(y=y, t=t, fig_axs=fig_axs)
```

```{code-cell} ipython3
%load_ext autoreload
%autoreload 2
from gpe.Examples.hydro_1d import Experiment, u

res0 = []

for n in [6, 7, 8, 9, 10]:
    e = Experiment(V0p=0.5, Nx=2**n)
    s = e.get_state()
    #s.plot()
    hist = s.evolve(steps=2000, t_final=10*u.ms, hist=True, show_plot=True);
    res0.append((np.array([s.t for s in hist]),
                 np.array([(s.get_density()*s.x).mean() / s.get_density().mean() for s in hist])))
```

```{code-cell} ipython3
plt.subplots(figsize=(20,3))
[plt.plot(t/s.experiment.t_unit, x, label=str(n)) for n, (t, x) in enumerate(res0)];plt.legend()
```

```{code-cell} ipython3
plt.plot(ts0/s.experiment.t_unit, x_cms0)
plt.plot(ts1/s.experiment.t_unit, x_cms1)
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp
import scipy as sp

class Flow:
    m = 1.2
    g = 1.3
    nu = 0.1  # Viscosity
    V0 = 0.2  # Potential barrier height
    L = 0.2   # Length of barrier

    n0 = 1.3   # Steady-state density on the left.
    v0 = 0.14  # Steady-state velocity on the left.
    dV = 0.15  # Potential drop (magnitude)
    
    potential = 'gaussian'
    potential_k = 1
    
    # Boundary conditions on the left
    def __init__(self, **kw):
        for key in kw:
            if not hasattr(self, key):
                raise ValueError(f"Unknown {key=}")
            setattr(self, key, kw[key])
        self.init()
        
    def init(self):
        if self.potential == 'gaussian':
            f = 10.0
            s = self.L/2
            self._V = lambda x: self.V0*np.exp(-f*(x/s-1)**2)
            self._dV = lambda x: -2*f*(x/s-1)/s*self.V0*np.exp(-f*(x/s-1)**2)
        else:
            xp, fp = [0, self.L/2, self.L], [self.dV, self.V0, 0]
            if self.potential_k > 1:
                xp = self.L*np.linspace(0, 1, 7)
                fp = [self.dV]*3 + [self.V0] + [0]*3
                xp = self.L*np.array([0, 0.01, 0.5, 0.99, 1])
                fp = [self.dV]*2 + [self.V0] + [0]*2
                
            self._V = sp.interpolate.InterpolatedUnivariateSpline(
                xp, fp, k=self.potential_k, ext='const')
            self._dV = self._V.derivative()

        # Compute the nu=0 bounds assuming the GPE
        j = self.v0*self.n0
        m = self.m
        nc = self.nc = (m * j**2/self.g)**(1/3)  # Relies on GPE EoS.
        self.mu0 = self.get_Eint(n=self.n0, d=1)
        self.E0 = m*self.v0**2/2 + self.mu0
        self.Vmax = self.E0 - m*(j/nc)**2/2 - self.get_Eint(nc, d=1)
        
    def get_Eint(self, n, d=0):
        if d == 0:
            return self.g * n**2 /2
        elif d == 1:
            return self.g * n
        elif d == 2:
            return self.g
        else:
            return 0

    def get_P(self, n):
        E, mu = [self.get_Eint(n, d=d) for d in [0, 1]]
        P = n*mu - E
        return P
            
    def compute_dy_dx(self, x, y):
        n, u, du, F = y
        dn = -n*du/u
        ddu = n * (u*du + (self._dV(x) + self.get_Eint(n, d=2)*dn)/self.m) / self.nu
        dF = n * self._dV(x)
        return (dn, du, ddu, dF)

    def solve(self, dv0=0.0, x0_L=-0.1, x_L=2.0, 
              method='BDF', max_step=0.1, x_eval=None, **kw):
        x_span = (x0_L*self.L, x_L*self.L)
        y0 = (self.n0, self.v0, dv0, 0.0)
        res = solve_ivp(self.compute_dy_dx,  y0=y0, t_span=x_span, t_eval=x_eval,
                        method=method, max_step=max_step, **kw)
        if not res.success:
            raise Exception(res.message)
        res.x = res.t
        res.n, res.u, res.du, res.F = res.y
        res.dn, du, res.ddu, res.dF = self.compute_dy_dx(res.x, res.y)
        self.res = res
        return res

f = Flow()
res = f.solve()
```

```{code-cell} ipython3
class Flow1(Flow):  
    def plot(self, x0_L=0.0, x_L=1.2, ax=None):
        """Plot the density and velocity of flow over a barrier."""
        res = self.solve(x0_L=x0_L, x_L=x_L)
        x, n = res.x, res.n
        Vx = self._V(x)
        n_TF = (self.mu0 - Vx)/self.g
        if ax is None:
            fig, ax = plt.subplots(1, 1, figsize=(3, 1.5))
        ax.plot(x, n, label="$n$")
        ax.plot(x, n_TF, '--', label='$n_{TF}$')
        ax.axhline(self.nc, color='y', ls=':')
        ax.set(ylim=(0, n.max()*1.05), 
               xlabel="$x$", ylabel="$n$")
        ax1 = ax.twinx()
        ax1.plot(x, Vx/self.V0*self.nc, 'C2-.', label='$V(x)$')
        ax.legend(loc='upper right')
        plt.sca(ax)
        res.ax = ax
        res.ax1 = ax1
        return res
        
Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.3).plot()
Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.425).plot()
Flow1(nu=0.003, n0=2.0, V0=1.0, v0=0.455).plot();
```

```{code-cell} ipython3
f = Flow1(nu=0.0001, n0=2.0, V0=1.0, v0=0.3)
nc = 1.0
jc = nc*np.sqrt(nc*f.get_Eint(nc, d=2)/f.m)
```

```{code-cell} ipython3
nus = np.linspace(0.001, 0.2)
fs = []
Fs = []
for nu in nus:
    f = Flow1(n0=2.0, v0=1.0, V0=0.19, nu=nu)
    res = f.solve(atol=1e-12, rtol=1e-12)
    fs.append(f)
    Fs.append(res.F[-1]-res.F[0])
plt.plot(nus, Fs)
```

```{code-cell} ipython3
nu = 0.05
V0s = np.linspace(0, 0.19, )
fs = []
Fs = []
for V0 in V0s:
    f = Flow1(n0=2.0, v0=1.0, V0=V0, nu=nu)
    res = f.solve(x_L=5.0, atol=1e-12, rtol=1e-12)
    fs.append(f)
    Fs.append(res.F[-1]-res.F[0])
fig, ax = plt.subplots(1, 1)
ax.plot(V0s, Fs)
ax.set(xlabel=r"$V_0$", ylabel="$F$")

ind = 25
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
for f, ax in zip([fs[0], fs[ind], fs[-1]], axs):
    f.plot(ax=ax)
```

```{code-cell} ipython3
def plot_flow(x0_L=0.0, x_L=2.0, **kw):
    f = Flow(**kw)
    res = f.solve(x0_L=x0_L, x_L=x_L, max_step=0.001, rtol=1e-12, atol=1e-12)
    
    m, nu = f.m, f.nu
    x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
    e, de, dde = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
    V, dV = f._V(x), f._dV(x)
    j = f.v0*f.n0
    n0 = f.n0
    K0 = m*f.v0**2/2
    mu0 = f.get_Eint(f.n0, d=1)
    E0 = K0 + mu0
    nc = (f.m * j**2 / f.g)**(1/3)
    Vmax = E0 - 3/2*(f.m*f.g**2*j**2)**(1/3)
    print(f.nc, f.Vmax, f.v0, Vmax)
    nc0 = 2*(E0 - f.V0)/3/f.g
    n_TF = (mu0 - V)/f.g
    n_min = (mu0 - f.V0)/f.g
    ns = [n_TF]
    for _i in range(10):
        _n = ns[-1]
        dn = -m*j**2/2/f.g * (1/_n**2 - 1/n0**2)
        ns.append(n_TF+dn)
    
    α_min = (mu0 - f.V0)/(m*j**2/2*(1/n_min**2-1/n0**2))
    print(n_min, n.min(), nc, nc0, n_min**2*(E0-f.V0)/(m*j**2/2))
    #(de/(m*u**2/2)).min()
    print(f"{α_min=:.2g}")
    print(E0, f.V0)
    fig, ax = plt.subplots(1, 1, figsize=(3, 2))
    ax.plot(x, n)
    for _n in ns:
        ax.plot(x, _n, ':')
    ax.set(ylabel="n")
    ax1 = ax.twinx()
    ax1.plot(x, f._V(x), 'C8--')
    res.E0, res.j = E0, j
    return f, res
    
f, res = plot_flow(dV=0.0, V0=1.0, nu=0.0001, n0=10.0, v0=1.795, #v0=1.795, 
          g=1.0, x_L=0.6)
```

```{code-cell} ipython3
Pn = [f.g, f.V0-res.E0, 0, f.m/2*res.j**2]
n = np.linspace(0, f.n0)
plt.plot(n, np.polyval(Pn, n))
```

```{code-cell} ipython3
Pu = [f.m/2, 0, f.V0-res.E0, f.g*res.j]
u = np.linspace(f.v0, 2*f.v0)
plt.plot(u, np.polyval(Pu, u))
```

## Flow Explorer

```{code-cell} ipython3
from ipywidgets import interact


kw = dict(L=0.1, x0_L=-100.0, x_L=100.2)
#plot_flow(dV=0.0, V0=8.46, nu=0.1, n0=10.0, v0=0.1, g=1.0);
plot_flow(dV=-0.1, V0=0.001, nu=0.1, n0=1.0, v0=1.0, g=1.0, potential='step', **kw);
```

```{code-cell} ipython3

```

```{code-cell} ipython3
plot_flow(dV=0.0, V0=1.0, nu=0.1, n0=10.0, v0=0.1, g=1.0)
plt.title("Slow flow over a barrier");

plot_flow(dV=0.0, V0=1.0, nu=0.1, n0=10.0, v0=1.9, g=1.0)
plt.title("Fast flow over a barrier");
```

```{code-cell} ipython3
plot_flow(dV=0.0, V0=1.0, nu=0.01, n0=10.0, v0=10.9, g=1.0, x_L=0.13)
plt.title("Fast flow over a barrier");
```

```{code-cell} ipython3
fig, ax = plt.subplots()
f = Flow(dV=0.05, V0=0.1, nu=0.1, n0=10.0, v0=2.743, g=1.0)
res = f.solve(x_L=2.0, max_step=0.001, rtol=1e-12, atol=1e-12)

m, nu = f.m, f.nu
x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
E, dE, ddE = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
V, dV = f._V(x), f._dV(x)
#ax.plot(x, n)

if False:
    ax.plot(x, u*du)
    ax.plot(x, -(dV + ddE*dn)/m)
    ax.plot(x, nu*ddu/n)
else:
    ax.plot(x, n)

#ax.set(ylim=(-0.1, 0.1))

ax1 = ax.twinx()
ax1.plot(x, f._V(x), 'C8--')
```

```{code-cell} ipython3
fig, ax = plt.subplots()
f = Flow(dV=0, V0=0.1, nu=0.1, n0=1.0, v0=0.15, g=0.0)
res = f.solve(x_L=2.0, max_step=0.001, rtol=1e-12, atol=1e-12)

m, nu = f.m, f.nu
x, u, du, ddu, n, dn = res.x, res.u, res.du, res.ddu, res.n, res.dn
E, dE, ddE = [f.get_Eint(n, d=d) for d in [0, 1, 2]]
V, dV = f._V(x), f._dV(x)
#ax.plot(x, n)

if False:
    ax.plot(x, u*du)
    ax.plot(x, -(dV + ddE*dn)/m)
    ax.plot(x, nu*ddu/n)
else:
    ax.plot(x, n)

#ax.set(ylim=(-0.1, 0.1))

ax1 = ax.twinx()
ax1.plot(x, f._V(x), 'C8--')
```

```{code-cell} ipython3
ec = (res.n*res.u*(
    f.m*res.u**2/2 + f._V(res.x) + f.get_Eint(res.n, d=1))
    -m*f.nu*res.u*res.du
)
plt.plot(x, ec)
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline


v0 = f.v0
L = f.L
T = L/v0

print(v0**2/2, f.V0)

def V_m(x, f=10.0, L=L, V0=f.V0, d=0):
    s = L/2
    if d == 0:
        return V0 * np.exp(-f*(x/s-1)**2)
    elif d == 1:
        return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
    
def compute_dy_dt(t, y):
    x, v = y
    dx = v
    dv = -V_m(x, d=1)
    return (dx, dv)


y0 = [0, v0]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, T), 
                max_step=0.001,
                t_eval=np.linspace(0, T, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline

def V_m(x, f=10.0, L=1.0, V0=0.2, d=0):
    s = L/2
    if d == 0:
        return V0 * np.exp(-f*(x/s-1)**2)
    elif d == 1:
        return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
    
def compute_dy_dt(t, y):
    x, v = y
    dx = v
    dv = -V_m(x, d=1)
    return (dx, dv)

y0 = [0, 1]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 1), 
                max_step=0.001,
                t_eval=np.linspace(0, 1, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
```

```{code-cell} ipython3
f.V0
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp
from scipy.interpolate import InterpolatedUnivariateSpline as Spline

def V_m(x, f=10.0, L=1.0, V0=0.2, d=0):
    s = L/2
    if d == 0:
        return V0 * np.exp(-f*(x/s-1)**2)
    elif d == 1:
        return -2*f*(x/s-1)/s*V0*np.exp(-f*(x/s-1)**2)
    
def compute_dy_dt(t, y):
    x, v = y
    dx = v
    dv = -V_m(x, d=1)
    return (dx, dv)

y0 = [0, 1]
res = solve_ivp(compute_dy_dt, y0=y0, t_span=(0, 1), 
                max_step=0.001,
                t_eval=np.linspace(0, 1, 300))
t = res.t
x, v = res.y
n = f.n0/np.diff(x)*(x[1]-x[0])
x, v = [(v[1:]+v[:-1])/2 for v in [x, v]]
plt.plot(x, n)
```

```{code-cell} ipython3
m = 1
v0 = 1
t = np.linspace(0, 2)
x = v0*t - 0.1*np.exp(-20*(t-1)**2)
v = v0 + 0.1*np.exp(-20*(t-1)**2)*40*(t-1)
V = m*(v0**2 - v**2)/2
a = 0.1*np.exp(-20*(t-1)**2)*40*(1 - 40*(t-1)**2)
plt.plot(x, a)
```

```{code-cell} ipython3
fig, ax = plt.subplots()
axn = ax.twinx()

js = np.linspace(0.001, 10)
nus = [0.025, 0.05, 0.1, 0.2]
nL = 1.0
g = 1.0
F_unit = g*nL**2/2
for nu in nus:
    Fs = []
    nRs = []
    for j in js:
        f = Flow(n0=nL, g=g, v0=j/nL, dV=0.0, nu=nu)
        res = f.solve(method='LSODA', max_step=0.1, rtol=1e-8)
        Fs.append(res.F[-1])
        nRs.append(res.n[-1])
    Fs, nRs = map(np.asarray, [Fs, nRs])
    l, = ax.plot(js/nu, Fs/F_unit, label=f"{nu=}")
    axn.plot(js/nu, nRs/nL, ls='--', c=l.get_c())
ax.set(xlabel=r'$j/\nu$', ylabel='$F$ [$gn_L^2/2$] (solid)',
       xlim=(0, 100))
axn.set(ylabel='$n_R/n_L$ (dashed)')
ax.legend(loc='center left')
ax.grid('on');
```

```{code-cell} ipython3
fig, ax = plt.subplots()
axn = ax.twinx()

js = np.linspace(0.001, 3)
nus = [0.025, 0.05, 0.1, 0.2]
nL = 1.0
g = 1.0
F_unit = g*nL**2/2
for nu in nus:
    Fs = []
    nRs = []
    for j in js:
        f = Flow(n0=nL, g=g, v0=j/nL, dV=0.0, nu=nu)
        res = f.solve(method='LSODA', max_step=0.1, rtol=1e-8)
        Fs.append(res.F[-1])
        nRs.append(res.n[-1])
    Fs, nRs = map(np.asarray, [Fs, nRs])
    l, = ax.plot(js, Fs/F_unit, label=f"{nu=}")
    axn.plot(js, nRs/nL, ls='--', c=l.get_c())
ax.set(xlabel='$j$', ylabel='$F$ [$gn_L^2/2$] (solid)')
axn.set(ylabel='$n_R/n_L$ (dashed)')
ax.legend(loc='center left')
ax.grid('on');
```

In the previous plot, we look at static flow across a barrier as a function of the
flow-rate (current) $j$.  We plot the total force on the barrier (solid), and the
density drop (dashed).



We start with the two conservation laws for density and momentum far to the left and
right of the potential, assuming that $u_{,x}=0$ in these regions, and setting $F(x_L) =
0$ and $F(x_R) = F$:

\begin{gather*}
  n_Lu_L = n_Ru_R = j = \text{const.}\\
  mn_Lu_L^2 + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) = 
  mn_Ru_R^2 + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
  \frac{m j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) = 
  \frac{m j^2}{n_R} + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\end{gather*}





\begin{gather*}
   nu^2 + \frac{P}{m} - ν u_{,x} 
   = nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
   u_{,x}
   = \text{const.},\\
  P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
  F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}


Assuming that $u_{,x} = 0$ far from the potential (static flow), and setiing $V(x_R) =
0$, we have the following after integrating the energy current:
\begin{gather*}
\left.
  u\left(
    \frac{m}{2}nu^2 + nV + n\mathcal{E}'(n) - m\nu u_{,x}
  \right)
\right|^{x_R}_{x_L}
= - \int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x} \leq 0,\\
  \tfrac{m}{2}(u_{R}^2 - u_{L}^2) + V_L + \mathcal{E}'(n_R) -
  \mathcal{E}'(n_L) =
  -\frac{1}{j}\int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x}
\end{gather*}
































## Correction Needed Below

The following analysis is incorrect.  It is missing the factor of $1/n$ discussed above,
and neglects to include a difference in the potential from one side to the other (DC
offset) required to maintain a persistent flow in the presence of dissipation.

## Steady State (Incorrect)

Now consider steady\-state flow $n(x, t) = n(x)$, $u(x, t) = u(x)$.  These become:
\begin{gather*}
  (nu)_{,x} = 0, \qquad
  muu_{,x} + V_{,x} + \mu'(n)n_{,x} = mνu_{,xx}.
\end{gather*}
The first equation establishes that the flux $Φ = nu$ is the same at every point, allowing us to exchange velocity $u$ for density $n$:
\begin{gather*}
  n = \frac{Φ}{u}, \qquad
  n_{,x} = -\frac{Φ}{u^2}u_{,x}.
\end{gather*}
Let's assume that $V(x)$ has compact support over $[0, L]$, then we can give boundary conditions $u(0) = u_0$, $u_{,x}(0) = 0$.  The second equation and be integrated once to give the non\-linear first\-order equation
\begin{gather*}
  mνu_{,x} = \frac{mu^2}{2} + V(x) + \overbrace{\mu\left(\frac{Φ}{u}\right)}^{\mu(n)}.
\end{gather*}

For general $V(x)$, this has no analytic solution, but we can compute a perturbative solution if we let $V(x) \rightarrow \lambda V(x)$ where $\lambda \ll 1$.  We note that $u(x) = u_0$ is a solution for $\lambda = 0$ and write
$$u(x) = u_0 + \lambda u_1(x) + \lambda^2 u_2(x) + O(\lambda^3).$$_Note: In what
follows, we revert back to the notation_ $u_1' = [u_1]_{,x}$_...

I originally wrote it this way, and have not changed it yet._  
As shown below, we originally tried expanding to first order in $\lambda$, but found this was insufficient, so we expand to second order in $\lambda$ obtaining  the following linear equations for $u_1$ and $u_2$:
\begin{align*}
  mνu_1' + \overbrace{\left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)}^{-mνc}u_1
  &= V(x),\\
  mνu_2' + \left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)u_2
  &= \underbrace{\left(\frac{m}{2}
    + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2u_0^2}
  \right)}_{mνb}u_1^2,\\
  u_1' - c u_1
  &= V(x)/mν,\\
  u_2' - cu_2
  &= bu_1^2,
\end{align*}
where we have introduced the constants
\begin{gather*}
  c = \frac{u_0}{ν}-\frac{n_0\mathcal{E}''(n_0)}{mνu_0}
  = 
  \frac{u_0}{ν}-\frac{n_0\mu'(n_0)}{mνu_0},\\
  b = \frac{1}{2ν} + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2mνu_0^2}
    = \frac{1}{2ν} + \frac{\bigl(nP'(n)\bigr)'|_{n=n_0}}{2mνu_0^2}.
\end{gather*}

## Do it!

To solve this problem, one needs to find a solution to the following inhomogeneous differential equation:
\begin{gather*}
  au_1'(x) - u_1(x) = f(x),
\end{gather*}
where $a>0$.  We are particularly interested in the limit where $a\rightarrow 0$, which will correspond to zero viscosity.

You should find that
\begin{gather*}
  u_1(x) = e^{x/a}\int_{x_0}^{x} d{y}\; f(y)e^{-y/a}
\end{gather*}

Notice that this has a non-analytic dependence $e^{x/a}$ on the parameter $a \propto ν$.  Show that, when expanding from the appropriate limit, all coefficients of the Taylor series for $e^{x/a}$ in powers of $a$ are zero. This makes it somewhat tricky to expand for small viscosities since all terms with such a factor will vanish.

Argue, however, that there may be an analytic *particular* solution $u_1(x, a)$ whose form can be found by expanding in powers of $a$.  We can then write the general solution as
\begin{gather*}
  u_1(x) = u_1(x, a) + Ce^{x/a}.
\end{gather*}
The non-analytic piece, thus describes transients (exponentially decaying in the appropriate direction) induced by the initial (boundary) conditions.

Show that this particular solution is:
\begin{gather*}
  u_1(x, a) = - \sum_{n=0}^{\infty} a^{n} f^{(n)}(x).
\end{gather*}
For example, this allows us to easily solve for $f(x) = \sin(kx)$:
\begin{gather*}
  -u_1(x, a) = \sin(x) + ak \cos(x) - (ak)^2 \sin(x) - (ak)^3 \cos(x) + \cdots\\
             = \Bigl(1 - (ak)^2 + (ak)^4 - \cdots\Bigr)\Bigl(\sin(kx) + ak \cos(kx)\Bigr)\\
             = \frac{\sin(kx) + ak \cos(kx)}{1 + (ak)^2},\\
   u_1(x) = Ce^{x/a} - \frac{\sin(kx) + ak \cos(kx)}{1+(ak)^2}.
\end{gather*}
We can easily check by differentiating that this satisfies the original equation.

## Force

The \(buoyant\) force acting on the potential $V(x)$ is
\begin{gather*}
  F = -\int_{0}^{L}d{x}\;nV'
    = -Φ\int_{0}^{L}d{x}\;\frac{V'}{u}
\end{gather*}
Returning to our original equation, we have
\begin{gather*}
  -nV'(x) = m\overbrace{nu}^{Φ}u_{,x} + \overbrace{nμ_{,x}}^{[P(n)]_{,x}} - mνnu_{,xx},\\
  F = \Bigl.(mΦu + P)\Bigr|_{0}^{L} - mν\int_0^{L}d{x}\;nu_{,xx}
\end{gather*}


# A Numerical Example

To make sure what we are doing makes sense, we give a concrete example.  To minimize the
chance for error, we implement the original equations in the form $y=(n, u, u', F)$,
allowing for a density\-dependent viscosity: 

\begin{gather*}
  (nu)' = 0, \qquad
  muu' + \Bigl(V + \mathcal{E}'(n)\Bigr)' = \frac{(mν(n)u')'}{n}
  = m\frac{ν'(n)n'u'+ν(n)u''}{n}\\
  \begin{pmatrix}
    n\\
    u\\
    u'\\
    F
  \end{pmatrix}'
  =
  \begin{pmatrix}
    -\frac{nu'}{u}\\
    u'\\
    \frac{n\Bigl(uu' + \frac{V'(x) + \mathcal{E}''(n)n'}{m}\Bigr)-ν'n'u'}{ν}\\
    nV'(x)
  \end{pmatrix}, \qquad
  y(0) = \begin{pmatrix}
    n_0\\
    u_0\\
    0\\
    0
  \end{pmatrix}.
\end{gather*}

```{code-cell} ipython3
assert False
%matplotlib inline
import matplotlib
from functools import partial
import numpy as np, matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (4,3)

m = 1.8
nu0 = 0.032
L = 1.3
n0 = 0.85
u0 = 0.773/2
Φ = n0*u0
g_m = 1.5
lam = 0.02

def mu(n, d=0):
    """Return the GPE chemical potential derivatives."""
    if d == 0:
        return m*g_m*n
    elif d == 1:
        return m*g_m+0*n
    elif d == 2:
        return 0*n

def mu(n, d=0):
    """Return the UFG chemical potential derivatives."""
    if d == 0:
        return m*g_m*n**(5/3)
    elif d == 1:
        return 5/3*m*g_m*n**(2/3)
    elif d == 2:
        return 10/9*m*g_m*n**(-1/3)

def mu_(n, d=0):
    """Return the GPE chemical potential derivatives."""
    if d == 0:
        return m*g_m*n**3
    elif d == 1:
        return 3*m*g_m*n**2
    elif d == 2:
        return 6*m*g_m*n

def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1, m*x*(x-L), 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1, m*(2*x-L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1, m*2, 0)

def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        m*np.sin(np.pi*(x-L)/L)**2, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*m*np.sin(2*np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        2*np.pi**2/L**2*m*np.cos(2*np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -4*np.pi**3/L**3*m*np.sin(2*np.pi*(x-L)/L), 0)

def V(x, d=0):
    """Return the external potential or derivative."""
    s = np.sin(np.pi*(x-L)/L)
    ds = np.cos(np.pi*(x-L)/L)*np.pi/L
    ds2 = -np.sin(np.pi*(x-L)/L)*(np.pi/L)**2
    ds3 = -np.cos(np.pi*(x-L)/L)*(np.pi/L)**3

    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        m*s**4, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        m*4*s**3*ds, 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        m*(12*s**2*ds**2 + 4*s**3*ds2), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        m*(24*s*ds**3 + 36*s**2*ds*ds2 + 4*s**3*ds3), 0)


def nu(n, nu0=1.0, gamma=0.0, d=0):
    """Return the density-dependent viscosity."""
    if d == 0:
        return nu0*n**gamma
    elif d == 1:
        return nu0*gamma*n**(gamma-1)

def compute_dy_dt(x, y, lam=lam, nu0=nu0, 
                  u0=u0, n0=n0, m=m):
    n, u, du, F = y
    dn = -n*du/u
    dV = lam*V(x, d=1)
    #ddu = (m*u*du + dV + mu(n, d=1)*dn)/m/nu
    nu_, dnu_ = [nu(n, nu0=nu0, d=d) for d in [0, 1]]
    ddu = (n*(m*u*du + dV + mu(n, d=1)*dn)/m - du*dnu_*dn)/nu_
    dF = n*dV
    return (dn, du, ddu, dF)

y0 = (n0, u0, 0, 0)

res = solve_ivp(
    partial(compute_dy_dt, nu0=0.2, lam=0.5),
    y0=y0,
    t_span=(-L, 5*L), max_step=0.001)
x = res.t
n, u, du, F = res.y
assert np.allclose(n*u, n0*u0)
fig, ax = plt.subplots()
ax.plot(x/L, n/n0, label="$n$")
ax.plot(x/L, u/u0, label="$u$")
ax.axvspan(0, 1, color='y', alpha=0.5)
ax.set(xlabel="$x/L$", ylabel="$n/n_0$", title=f"{n[0]=}, {n[-1]=}", );
plt.legend(loc="upper left")
ax = plt.twinx()
ax.plot(x/L, F, '--C1', label="$F$")
ax.legend(loc="lower right")
ax.set(ylabel="$F$");
```

Notice that the potential causes the density to change, then the viscosity allows this to relax outside of the shaded region where the potential is active.  We now check some relationships.

```{code-cell} ipython3
# Check that the potential and mu derivatives are correct
fig, axs = plt.subplots(1, 2)
ax = axs[0]
ax.plot(x/L, V(x), 'C0')
ax.plot(x/L, V(x, d=1), 'C1-')
ax.plot(x/L, np.gradient(V(x), x), 'C1:')
ax.plot(x/L, V(x, d=2), 'C2-')
ax.plot(x/L, np.gradient(V(x, d=1), x), 'C2:')
ax.plot(x/L, V(x, d=3), 'C3-')
ylim = ax.get_ylim()
ax.plot(x/L, np.gradient(V(x, d=2), x), 'C3:')
ax.set(ylim=ylim)

ax = axs[1]
ns = np.linspace(n0/2, 2*n0)
ax.plot(ns/n0, mu(ns), 'C0')
ax.plot(ns/n0, mu(ns, d=1), 'C1-')
ax.plot(ns/n0, np.gradient(mu(ns), ns), 'C1:')
ax.plot(ns/n0, mu(ns, d=2), 'C2-')
ax.plot(ns/n0, np.gradient(mu(ns, d=1), ns), 'C2:')
```

```{code-cell} ipython3
from functools import partial

def get_F(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
    y0 = (n0, u0, 0, 0)
    dy_dt = partial(compute_dy_dt, lam=lam, nu=nu,
                    u0=u0, n0=n0, m=m)
    res = solve_ivp(dy_dt,
                    y0=y0, t_span=(0, L), 
                    method='BDF', 
                    rtol=tol, atol=tol)
    x = res.t
    n, u, du, F = res.y
    return F[-1]

lams = np.linspace(0, 0.5)[1:-1]
nus = np.linspace(0, 0.5)[1:-1]
Fs_lam = [get_F(lam=lam) for lam in lams]
Fs_nu = [get_F(nu=nu) for nu in nus]

fig, axs = plt.subplots(1, 2, figsize=(5, 2))
axs[0].plot(lams**2, Fs_lam)
axs[0].set(xlabel=r"$\lambda^2$", ylabel="$F$")
axs[1].plot(nus, Fs_nu)
axs[1].set(xlabel=r"$\nu$", ylabel="$F$");
```

Now let's check the small $λ$ approximations:

```{code-cell} ipython3
from functools import partial
from scipy.integrate import cumtrapz

c = u0/nu - n0*mu(n0, d=1)/m/nu/u0
b = (1 + (2*n0*mu(n0, d=1) + n0**2*mu(n0, d=2))/m/u0**2)/2/nu
a = c*nu
print(f"{a=}, {nu*b}")


def f(x, us):
    u1, u2 = us
    du1 = V(x)/m/nu + c*u1
    du2 = b*u1**2 + c*u2
    return (du1, du2)
res1 = solve_ivp(f, y0=[0, 0], t_span=(-L, 2*L), t_eval=x,
                 method='BDF', atol=1e-8, rtol=1e-8)
u1, u2 = res1.y
n1 = Φ/(u0+lam*u1)
F1 = cumtrapz(lam*n1*V(x, d=1), x, initial=0)

# Small nu approximations
u10 = -V(x)/m/a
u11 = -V(x, d=1)/m/a**2
u12 = -V(x, d=2)/m/a**3

fig, axs = plt.subplots(2, 2, figsize=(10, 5))
ax = axs[0,0]
ax.plot(x/L, u1/u0, label='$u_1/u_0$')
ax.plot(x/L, u10/u0, label=r"small $\nu$ approx")
ax.plot(x/L, (u-u0)/u0/lam, label="full")

ax.set(xlabel=r"$x/L$", ylabel="$u_1/u_0$")
ax.legend()

ax = axs[0, 1]
ax.plot(x/L, F1, label='$F_1$')
ax.plot(x/L, F, label="full")
ax.legend()

ax = axs[1, 0]
ax.plot(x/L, u1/u0, 'C0-')
ax.plot(x/L, u10/u0, 'C0:')
ax.plot(x/L, (u1-u10)/nu/u0, 'C1-')
ax.plot(x/L, u11/u0, 'C1:')
ax.plot(x/L, u12/u0, 'C2:')
ylim = ax.get_ylim()
ax.plot(x/L, (u1-u10-nu*u11)/nu**2/u0, 'C2-', scaley=False)
ax.set(ylim=ylim)

ax = axs[1, 1]
ax.plot(x/L, (u-u0)/lam/u0, 'C0-')
ax.plot(x/L, u1/u0, 'C0:')
ax.plot(x/L, (u-u0-lam*u1)/lam**2/u0, 'C1-')
ax.plot(x/L, u2/u0, 'C1:')
```

The following code demonstrates that the order $λ$ solution $u_1(x)$ gives the correct force when integrated. *(There is no need to compute $u_2(x)$.)*

```{code-cell} ipython3
def get_F1(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
    a = u0 - n0*mu(n0, d=1)/m/u0
    c = a/nu
    def f(x, y):
        u1, F1 = y
        du1 = V(x)/m/nu + c*u1
        dF = -lam**2 * n0 / u0 * u1 * V(x, d=1)
        return (du1, dF)

    res1 = solve_ivp(f,
        y0=[0, 0], t_span=(0, L),
        method='BDF', atol=tol, rtol=tol)
    u1, F1 = res1.y
    return F1[-1]


lams = np.linspace(0.01, 0.3, 20)
qs = [get_F1(lam=lam, tol=lam**2 * 1e-8) /
      get_F(lam=lam, tol=lam**2 * 1e-8) for lam in lams]
plt.plot(lams, qs, '-+')
plt.plot([0], [1], 'o')
plt.gca().set(xlabel=r'$\lambda$', ylabel='$F_1/F$');
```

Thus, the complete solution to order $λ$ is to solve the following problem:
\begin{gather*}
  u_1'(x) + cu_1(x) = V(x)/mν,\qquad
  c = \frac{u_0 - P'_0/mu_0}{ν}\\
  F %= λ\int_{0}^{L}d(x)\;\frac{Φ V'(x)}{u_0+λu_1}
    = \frac{-λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1(x) V'(x)
    = \frac{λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1'(x) V(x)
\end{gather*}
Note that the coefficient $c$ here is generally less than zero.  We can see this by noting that the speed of sound in a material is $c_s$ is:
\begin{gather*}
  c_s = \sqrt{\frac{n_0\mu'(n_0)}{m}} = \sqrt{\frac{P'(n_0)}{m}},\\
  c = \frac{u_0-c_s^2/u_0}{\nu} = \frac{u_0}{\nu}\left(1 - \frac{c_s^2}{u_0^2}\right).
\end{gather*}
Thus, for subsonic flow $u_0<c_s$, $c<0$.  We keep this sign convention to avoid many other signs below.  The homogeneous solution is $u_1(x) = Ae^{cx}$.  To find the general solution, there are several options:

- One is the method of [variation of parameters](https://en.wikipedia.org/wiki/Variation_of_parameters) can be used, promoting $A → A(x)$ to find
  \begin{gather*}
    u_1(x) = A(x)e^{cx}, \qquad
    A'(x) e^{cx} = V(x)/{mν}, \qquad
    A(x) = \frac{1}{mν}\int_0^xd{y}\; e^{-cy}V(y), \\
    u_1(x) = \frac{1}{mν}\int_0^xd{y}\; e^{c(x-y)}V(y),\\
    F = \frac{-λ^2n_0}{mu_0}\int_{0}^{L}d{x}\;\int_0^xd{y}\; e^{c(x-y)}V(y)V'(x).
  \end{gather*}
  The integral is over the lower triangle $y<x$ of the box $(x, y) ∈[0, L]×[0,L]$. To better express the dimensionality of the solution, we write the potential as $V(x) = V_0U(x/L)$ so that the shape is described by the dimensionless function $U(\xi)$, and we introduce the dimensionless quantity $\kappa \propto \nu$:
  \begin{gather*}
    λV(x) = λV_0U(x/L), \qquad
    \kappa = \frac{1}{-cL} = \frac{\nu}{-aL},\\
    \begin{aligned}
    F &= \frac{-λ^2n_0V_0^2L}{mu_0}
    \underbrace{\int_{0}^{1}d{x}\;\int_0^xd{y}\; e^{(y-x)/\kappa}U(y)U'(x)}_{-\kappa^2 f(\kappa, U)},\\
     &= \frac{Lu_0\kappa^2}{\nu}\frac{λ^2V_0^2n_0}{mu_0^2}f(\kappa, U).
     \end{aligned}
  \end{gather*}
  While formally correct, this solution obscures the viscosity dependence.  We will see below that the normalization chosen here for $f(\kappa, U)$ makes sense in the limit of small viscosity:
  \begin{gather*}
    f(\kappa, U) = 
    -\int_{0}^{1}d{x}\;\int_0^xd{y}\; \frac{e^{(y-x)/\kappa}}{\kappa^2}U(y)U'(x)
  \end{gather*}

- Another is formal manipulation.  We first introduce the parameter $a=cν$ which is independent of the viscosity, then solve algebraically treating the derivative operator as a matrix:
  \begin{align*}
    νu_1'(x) &= au_1(x) + V(x)/m, \\
    u_1(x) &= \frac{1}{\left(1-\frac{ν}{a}\frac{d}{d{x}}\right)}\frac{-V(x)}{am},\\
    &= -\frac{V(x)}{am} - \frac{νV'(x)}{a^2m} - \frac{ν^2V''(x)}{a^3m} + O(ν^3).
  \end{align*}
  This is nice because it immediately allows us to expand in powers of $ν$ by expanding $1/(1-x) =1+x+x^2+x^3+\cdots$:
  \begin{gather*}
    F = \frac{λ^2n_0 V_0^2}{mau_0} \int_{0}^{1}d{x}\;\left(
      U(x)U'(x) - \kappa[U'(x)]^2 + \kappa^2U''(x)U'(x) + O(ν^3)
    \right)
  \end{gather*}
  Since $U(0)=U(1)=0$, the first term, which is a total derivative $(U^2)'/2$, vanishes. This gives:
  \begin{align*}
    F &= \frac{mλ^2n_0}{au_0}\frac{V_0^2}{m^2}(-\kappa)\overbrace{
    \Biggl(
      \int_{0}^{1}d{x}\;[U'(x)]^2
      -\kappa\left.\frac{[U'(x)]^2}{2}\right|_{0}^{1}
      +\kappa^2\int_{0}^{1}d{x}\;U'(x)U'''(x)
      + O(ν^3)
    \Biggr)}^{f(\kappa, U)}\\
    &= \frac{mλ^2n_0}{au_0}\frac{ν}{a} \frac{V_0^2}{m^2L}
       f(\kappa, U).
  \end{align*}
  This justifies the form of $f(\kappa, U)$ chosen above.
  If we choose a symmetric potential so that $|U'(0)| = |U'(1)|$, then the second term vanishes, allowing us to precisely define what a small viscosity means:
  \begin{gather*}
      \frac{ν^2}{a^2}
      \frac{\int_{0}^{L}d{x}\;V'(x)V'''(x)}
      {\int_{0}^{L}d{x}\;[V'(x)]^2}
      =
      \kappa^2
      \frac{\int_{0}^{1}d{x}\;U'(x)U'''(x)}
      {\int_{0}^{1}d{x}\;[U'(x)]^2}
      \ll 1.
  \end{gather*}

```{code-cell} ipython3
u10 = -V(x)/a/m
u11 = -1/a * V(x, d=1)/a/m
u12 = -1/a**2 * V(x, d=2)/a/m
u13 = -1/a**3 * V(x, d=3)/a/m

plt.plot(x/L, (u1 - u10)/nu)
plt.plot(x/L, (u1 - u10 - nu*u11)/nu**2)
plt.plot(x/L, (u1 - u10 - nu*u11 - nu**2*u12)/nu**3)

coeff = m*lam**2*n0/a/u0*nu/a
t1 = np.trapz(V(x, d=1)**2, x)/m**2
t2 = np.trapz(V(x, d=1)*V(x, d=2), x)/m**2/a
t3 = np.trapz(V(x, d=1)*V(x, d=3), x)/m**2/a**2

print(f"{t1=}, {t2=}, {t3=}")

F_ = get_F(tol=1e-13)/coeff
F1_ = get_F1(tol=1e-13)/coeff
print(F_, (F_-F1_)/lam)
print(F1_, t1)
print((F1_ - t1)/nu,
      (F1_ - t1 - nu*t2)/nu**2,
      (F1_ - t1 - nu*t2 - nu**2*t3)/nu**3,
      )
```

We now replace the parameter $a$ with the speed of sound $c_s$.  The final result, to leading order in $λ$ is:

\begin{gather*}
  F = \overbrace{\frac{ν}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)^2}\frac{λ^2V_0^2n_0}{mu_0^2}
    f(\kappa, U), \qquad
    \mathrm{Re} = \frac{u_0L}{ν}, \qquad
    \kappa = \overbrace{\frac{\nu}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}, \qquad
    c_s = \sqrt{\frac{P'(n_0)}{m}},\\
  F = \frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}\frac{λ^2V_0^2n_0}{mu_0^2}
    \Bigl(\kappa f(\kappa, U)\Bigr)
\end{gather*}

Here we have expressed the results in terms of the dimensionless Reynold's number $\mathrm{Re}$ and the ratio of the speed of sound $c_s$ to the flow $u_0$.  Note that the force diverges as the flow speed approaches the speed of sound.  The dependence on the shape appears with the dimensionless parameter $\kappa = ν/(-aL)$:
\begin{align*}
  f(\kappa, U) &= \frac{-1}{\kappa^2}\int_0^{L} dx\int_0^{x}d{y}\;e^{(y-x)/\kappa}
  U(y)U'(x),\\
  &\approx \sum_{n=1}^{\infty}(-\kappa)^{n-1}
  \int_0^1d{x}\;U^{(n)}(x)U'(x).
\end{align*}

**Caution:** The expansion in $\kappa \propto \nu$ here justifies our choice for the form of $f(\kappa, U)$ and the corresponding separation of dimensionful parameters, but it turns out that this series does not converge to the solution.  We now demonstrate this.  Consider for example, the dimensionless potential
\begin{align*}
  U(x) &= \sin(\pi(x-1)), &
  f(\kappa, U) &= \frac{\pi^2}{2}
    \frac{1+\pi^2\kappa^2 - 2\kappa(1+e^{-1/\kappa})}
         {(1+\pi^2\kappa^2)^2},\\
  U(x) &= \sin^2(\pi(x-1)), &
  f(\kappa, U) &= \frac{\pi^2}{2}
  \frac{\Bigl(
    1 + 4\pi^2\kappa^2 
    + 8\pi^2\kappa^3(1 - e^{-1/\kappa})\Big)}
    {(1+ 4\pi^2\kappa^2)^2}
 \end{align*}
This [solution](https://www.wolframalpha.com/input?i=-integrate%28integrate%28e%5E%7B%28y-x%29%2Fkappa%7D*sin%28pi*%28y-1%29%29*cos%28pi*%28x-1%29%29*pi%2Fkappa%5E2%2C+y%3D0..x%29%2C+x%3D0..1%29) contains the non\-analytic term $e^{-1/\kappa}$ whose Taylor series for positive $\kappa$ is zero.  Thus, the series expansion converges to the incorrect form of:
\begin{gather*}
  \sum_{n=1}^{\infty}(-\kappa)^{n-1}
  \int_0^1d{x}\;U^{(n)}(x)U'(x) =
  \frac{\pi^2}{2}\left(1 - (\pi\kappa)^2 + (\pi\kappa)^4 - \cdots\right)\\
  = \frac{\pi^2}{2}\frac{1}{1+\pi^2\kappa^2}
\end{gather*}

fails to converge.

```{code-cell} ipython3
def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        np.sin(np.pi*(x-L)/L), 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*np.cos(np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        -np.pi**2/L**2*np.sin(np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -np.pi**3/L**3*np.cos(np.pi*(x-L)/L), 0)

def get_f(nu, u0=u0, n0=n0, V0=1, L=L, m=m, lam=lam):
    cs2 = n0*mu(n0, d=1)/m
    coeff = nu/u0/L / (cs2/u0**2-1)**2 * lam**2 * V0**2 *n0/m/u0**2
    F1 = get_F1(nu=nu, u0=u0, n0=n0, L=L, m=m, lam=lam, tol=1e-8)
    kappa = nu/u0/L/(cs2/u0**2-1)
    return kappa, F1/coeff

nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])

xs = (np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
         * (1 + xs - 2*kappas*(1+np.exp(-1/kappas)))
         / (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
         * (1 + xs - 2*kappas*(1+0*np.exp(-1/kappas)))
         / (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
```

This figure shows that the integral gives the correct formula, and that for large $\kappa$ the non\-analytic piece is correct.  **It also shows that we have made a mistake in our series approximation.  Please fix!**

Here is $U(x) = \sin^2\bigl(\pi(x-1)\bigr)$:

```{code-cell} ipython3
def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        np.sin(np.pi*(x-L)/L)**2, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*np.sin(2*np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        2*np.pi**2/L**2*np.cos(2*np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -4*np.pi**3/L**3*np.sin(2*np.pi*(x-L)/L), 0)

nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])

xs = (2*np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
         * (1 + xs + 2*xs*kappas*(1-np.exp(-1/kappas)))
         / (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
         * (1 + xs + 2*xs*kappas*(1 - 0*np.exp(-1/kappas)))
         / (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
```

```{code-cell} ipython3
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
    return sympy.sin(sympy.pi*(x-1))

integrand = -1/k**2 * sympy.exp((y-x)/k)*U(x).diff(x)*U(y)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1)).simplify()
display(f.factor())
display(sympy.series((1+(sympy.pi*k)**2 - 2*k)/(1+(sympy.pi*k)**2), k, 0, 10))

f_series = sum(sympy.simplify(
        (U(x).diff(*(x,)*n)*U(x).diff(x))
        .integrate((x, 0, 1))*(-k)**(n-1))
    for n in range(1, 9)[::2])
display(f_series)
```

If we choose $U(\xi) = U(1-\xi)$ to be symmetric, then the even terms vanish, and we can write:
\begin{gather*}
  f(\kappa, U) = \sum_{n=0}^{\infty}c_{n}\kappa^{2n}, \qquad
  c_n = \int_0^1d{\xi}\;U^{(2n+1)}(\xi)U'(\xi).
\end{gather*}
Here are a couple of special cases.  Note: For numerical work, it might be better to use a potential that is sufficiently flat at the boundaries, hence the higher power.  This is what I have used in the code.  (Please check these.)
\begin{gather*}
  f(\kappa, U) = \pi^2g(\pi^2\kappa^2),\qquad
  \kappa = \frac{\nu/u_0L}{\Bigl(1-\frac{c_s^2}{u_0^2}\Bigr)},\\
  \begin{aligned}
  U(\xi) = \sin\bigl(\pi(\xi-1)\bigr)&:& 
  c_n &= \frac{\pi^2}{2}(-\pi^2)^{n}, &
  g(x) &= \frac{1}{2(1+x)},\\
  U(\xi) = \sin^4\bigl(\pi(\xi-1)\bigr)&:& 
  c_n &= \frac{\pi^2}{2}(4^{2n-1} + 4^{n})(-\pi^2)^{n}, &
  g(x) &= \frac{68x + 5}{8(1+4x)(1+16x)}.
  \end{aligned}
\end{gather*}
Thus, for a sinusoidal potential (please check):
\begin{gather*}
  \begin{aligned}
  \frac{F}{L} &= 
  \frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa}{2(1+\pi^2\kappa^2)}, &
  U(\xi) &= \sin\bigl(\pi(\xi-1)\bigr),\\
  \frac{F}{L} &= 
  \frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa(68\pi^2\kappa^2+5)}{8(1+4\pi^2\kappa^2)(1+16\pi^2\kappa^2)}, &
  U(\xi) &= \sin^4\bigl(\pi(\xi-1)\bigr),\\
  \end{aligned}
\end{gather*}

These can also be obtained by :

\begin{gather*}
  %-integrate(integrate(e^{(x-y)/kappa}*sin(pi*(y-1))*cos(pi*(x-1))*pi/kappa^2, y=0..x), x=0..1)
  f(\kappa, U) = \frac{-1}{\kappa}\int_0^L d{x}\int_0^x d{y}\; e^{(x-y)/\kappa} U'(x)U(y),\\
  \frac{\pi^2}{2}\frac{(1 + 2 (1 + e^{1/\kappa}) κ + π^2 κ^2)}{(1 + \pi^2 \kappa^2)^2}
\end{gather*}

```{code-cell} ipython3
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
    return sympy.sin(sympy.pi*(x-1))**2

integrand = -1/k*sympy.exp((x-y)/k)*U(x).diff(x)*U(y)
display(integrand)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1))
f.expand().factor()
```

```{code-cell} ipython3
sympy.series(sympy.exp(1/k), k, 0, 5)
```

*(I found these from the expansion https://oeis.org/A092812)*:

\begin{gather*}
  f(\kappa, U) = \pi^2 \sum_{n=0}^{\infty}
  \left(\frac{16^{n}}{8} + \frac{4^{n}}{2}\right)(-\pi^2\kappa^2)^{n}
  = \pi^2 g(\pi^2\kappa^2)
  , \\
  g(x) %= \frac{1 + 16x + 24x^2}{(1+4x)(1+16x)} - \frac{3}{8}
       = \frac{68x + 5}{8(1+4x)(1+16x)}
\end{gather*}

```{code-cell} ipython3
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))

def get_c(n):
    return k**(2*n)*(V.diff(*(x,)*(2*n+1))*V.diff(x)).integrate((x, 0, 1))

[get_c(n) for n in [0, 1, 2, 3, 4, 5]]
```

```{code-cell} ipython3
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))**4
def get_c(n):
    f = sympy.lambdify([x], V.diff(*(x,)*(n+1))*V.diff(x))
    return int(np.round(quad(f, 0, 1)[0]) / np.pi**(n+2)) * k**(n)*sympy.pi**(n+2)
get_c(4)
4**(0-1)/2 + 2**(0-1), 5/8

sympy.simplify((1+16*x+24*x**2)/(1+4*x)/(1+16*x) - sympy.S(3)/8).factor()
```

```{code-cell} ipython3
from scipy.integrate import quad
def f(x):
    return V(x, d=1)**2

cs = np.sqrt(n0*mu(n0, d=1)/m)
Re = u0*L/nu
tmp = lam**2 * quad(f, 0, L)[0]
print(a/u0, (1 - (cs/u0)**2))
print(
  (1/Re) * (1/(1-(cs/u0)**2)**2) * (n0*L/m/u0**2) * tmp,
  get_F())
```

```{code-cell} ipython3
from scipy.integrate import quad
def f(x):
    return V(x)*V(x,d=1)**2

tmp = quad(f, 0, L)[0]
lam**2*n0*nu/m/a**2/u0 * tmp, get_F()
```

```
>>>>>>> variant B
---
jupytext:
  cell_metadata_filter: -all
  formats: md:myst,ipynb
  main_language: python
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.13.7
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
import mmf_setup; mmf_setup.nbinit()
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

[Rankine-Hugoniot conditions]: <https://en.wikipedia.org/wiki/Rankine%E2%80%93Hugoniot_conditions>
[Shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations>
[conservative form of the shallow water equations]: <https://en.wikipedia.org/wiki/Shallow_water_equations#Conservative_form>





----








# Effective Viscosity (Old)

## Conservation Laws

### Particle Number

The first hydrodynamic equation demonstrates the conservation of particle number
\begin{gather*}
  N = \int\d{x}\; n(x),\\
  \dot{N} = \int\d{x}\; \dot{n}(x)
          = \int\d{x}\; (-nu)_{,x}
          = -\Bigl.nu\Bigr|_{x_L}^{x_R}
          = \Phi_{n}(x_L) - \Phi_n(x_R).
\end{gather*}
This says that that total rate of change of particles in a region between $x_L$ and
$x_R$ is equal to the flux of particles $\Phi_{n}(x_L)$ entering the region from the
left at $x_L$ minus the flux of particles $\Phi_{n}(x_R)$ leaving the region at $x_R$ on
the right.  Notice that a positive flux corresponds to a flow to the right, consistent
with the signs in this equation.


:::::{admonition} To Do.  What is the meaning of $D_t n = -nu_{,x}$?
:class: dropdown

What is the physical meaning of the convective derivative of the density?
\begin{gather*}
  D_t n = \dot{n} + u n_{,x} = -nu_{,x}.
\end{gather*}
Why is this not zero?

:::{admonition} Solution
:class: dropdown

The convective derivative characterizes how a quantity changes as you follow the
particle flow.  If you are following a particle, how does the density change?  Well,
this depends on how close the neighbouring particles are.  The average distance to the
neighbours changes iff there is a gradient in the velocity field.  To see this, consider
two particles at $x_0 < x_1$ so that
\begin{gather*}
  n \propto \frac{1}{x_1-x_0} = \frac{1}{\delta},\\
  \dot{n} \propto -\frac{\dot{x}_1-\dot{x}_0}{(x_1-x_0)^2}
  = -\frac{\dot{x}_1-\dot{x}_0}{\delta^2} \approx -\frac{u_{,x}}{\delta} = -n u_{,x}.
\end{gather*}
I.e. the density only changes if the velocity of the relative particles is different --
this is exactly what a gradient in the velocity indicates.
:::
:::::

### Momentum Density

Multiply the continuity equation by $u$ and add it to $n$ times the Bernoulli equation:
\begin{gather*}
  u\dot{n} + u(nu)_{,x} + 
  n\dot{u} + nuu_{,x} + \frac{n}{m}\Bigl(
    V(x) + \mathcal{E}'(n)\Bigr)_{,x} 
    - (ν u_{,x})_{,x} = 0\\
  (nu)_{,t} + \left(nu^2 + \frac{P}{m} - ν u_{,x}\right)_{,x} = 0,\\
\end{gather*}
This means that, if we have a static solution, then the following is constant
everywhere: 
\begin{gather*}
   nu^2 + \frac{P}{m} - ν u_{,x} 
   = nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
   u_{,x}
   = \text{const.},\\
  P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
  F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}
Unfortunately, as mentioned above, the force depends on the solution through the
density, so we can use this to check and interpret results, but it does not help us
solve anything.

### Energy

We can similarly consider the energy flux.  The energy density is the sum of the
kinetic, potential, and internal energies.
\begin{gather*}
  \varepsilon = \frac{m}{2}nu^2 + nV + \mathcal{E}(n),\\
  \dot{\varepsilon} + \left[
    u\left(\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n)\right)
  \right]_{,x} = n \dot{V} + mu(\nu u_{,x})_{,x},\\
  \dot{\varepsilon} + \left[
    u\left(\frac{m}{2}nu^2 + nV + n\mathcal{E}'(n)
    - m\nu u_{,x}
    \right)
\right]_{,x} = n \dot{V} - m\nu u_{,x}^2.
\end{gather*}
Here we see explicitly that energy conservation is violated by time-dependent
potentials, and by a viscosity **if** the velocity field has a gradient $u_{,x} \neq
0$. Thus, there will be no dissipation from this viscosity if we have simple sloshing
in a harmonic trap where the velocity is constant.  The second form gives a current
density whose gradient is strictly decreasing if viscosity is active and $\dot{V}=0$.

# An Analytic Example

Consider static flow: $\dot{n} = \dot{u} = 0$.  The density and momentum continuity
equations give the following in terms of the constants $j$ and $p$:
\begin{gather*}
  un = j, \qquad
  mnu^2 + F(x) + n \mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} = mp
  = m \frac{j^2}{n_L} + P(n_L).
\end{gather*}
Consider a step-like potential with small $\epsilon\rightarrow 0$:
\begin{gather*}
  V(x) = \begin{cases}
    V_L & x < 0,\\
    0 & \epsilon < x,
  \end{cases} \qquad
  V_{,x}(x) \approx -V_L\delta(x).
\end{gather*}
Working backwards, let $u(x)$ be linear in the region $[0, \epsilon]$
\begin{gather*}
  u'(x) = \begin{cases}
    \frac{u_R - u_L}{\epsilon} & 0 < x < \epsilon\\
    0 & \text{otherwise}
  \end{cases}.
\end{gather*}
In this interval, we have
\begin{gather*}
  u(x) = u_L + x\frac{u_R-u_L}{\epsilon}, \qquad
  n(x) = \frac{j}{u(x)} = \frac{u_L n_L}{u_L + x\frac{u_R-u_L}{\epsilon}},\\
\end{gather*}





For $x \gg \epsilon$ we have
\begin{gather*}
  n(x \gg \epsilon) = n_L + \delta n, \qquad
  u(x \gg \epsilon) = u_L + \delta u = \frac{j}{n_L + \delta n}, \quad
  F(x \gg \epsilon) = \delta F,
\end{gather*}
where
\begin{gather*}
  m\frac{j^2}{n_L + \delta n} + \delta F + P(n_L + \delta n) = m \frac{j^2}{n_L} + P(n_L),\\
  \delta F = \int_{0}^{\epsilon} nV'(x)\d{x} 
           = -V_L(n_L + \alpha\delta n)
           \in [-n_LV_L, -(n_L + \delta n)V_L],\\
           0 \leq \alpha \leq 1.
\end{gather*}
For small $\delta n$ we can write
\begin{gather*}
  P(n_L + \delta n) - P(n_L) \approx \delta n P'(n_L) + \frac{(\delta n)^2}{2}P''(n_L),
\end{gather*}
which is exact for the GPE $P(n) = gn^2/2$.  This gives:
\begin{gather*}
  \delta n = \frac{V_L n_L}{
  P'(n_L) + \frac{\delta n}{2} P''(n_L)
  - \frac{mj^2}{(n_L + \delta n)n_L} 
  - V_L\alpha}.
\end{gather*}
This can be iterated a couple of times to converge.
If $\delta n$ is small, we can expand:
\begin{gather*}
  \frac{\delta n}{n_L} \approx \frac{1}{
  \frac{n_L\mathcal{E}''(n_L) - mu_L^2}{V_L} - \alpha}
\end{gather*}







Let's start with a perfect fluid without viscosity $\nu=0$ flowing over a step
The momentum continuity equation gives
\begin{gather*}
  mnu^2 + F(x) + n\mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} = \text{const.}
\end{gather*}
Starting with $n(x<0) = n_L$, $u(x<0) = u_L$, and $F(x<0) = 0$ with $j=u_Ln_L$ constant,
we have
\begin{gather*}
  m\frac{j^2}{n} + F(x) + n\mathcal{E}'(n) - \mathcal{E}(n) - m\nu u_{,x} = 
  m\frac{j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L).
\end{gather*}



Integrating the momentum conservation equation across the potential we have:
\begin{gather*}
  mj_0^2 n_{,x} = n^3  (V+\mathcal{E}'')_{,x}, \qquad
  mj_0^2\delta n = n_0^2 (-V_L + \mathcal{E}''(n)\delta
\end{gather*}




# A Numerical Example: Steady-State Flow

We start with steady-state flow over a barrier.  To make sure what we are doing makes
sense, we give a concrete example.  To minimize the chance for error, we implement the
original equations in the form $y=(n, u, u', F)$, looking for stationary solutions
$\dot{n} = \dot{u} = 0$:

\begin{gather*}
  (nu)_{,x} = 0, \qquad
  muu_{,x} = -(V+\mathcal{E}')_{,x} + m(\nu u_{,x})_{,x}/n\\
  n_{,x} = -\frac{nu_{,x}}{u}, \qquad
  u_{,xx} = \frac{mnuu_{,x} + n(V+\mathcal{E}')_{,x}}{m\nu}\\
  \begin{pmatrix}
    n\\
    u\\
    u_{,x}\\
    F
  \end{pmatrix}_{,x}
  =
  \begin{pmatrix}
    -\frac{nu_{,x}}{u}\\
    u_{,x}\\
    n\frac{uu_{,x} + \frac{(V+\mathcal{E}')_{,x}}{m}}{\nu}\\
    nV_{,x}
  \end{pmatrix}, \qquad
  y(0) = \begin{pmatrix}
    n_0\\
    u_0\\
    0\\
    0
  \end{pmatrix}.
\end{gather*}

:::{admonition} Removing the Velocity

Alternatively, we can use the conservation equation to remove the velocity from our
equations since $nu = n_0u_0 = j_0 = \text{const.}$ everywhere:
\begin{gather*}
  u = \frac{j_0}{n}, \quad
  u_{,x} = \frac{-j_0 n_{,x}}{n^2}, \quad
  u_{,xx} = \frac{-j_0 n_{,xx}}{n^2}+\frac{2j_0 (n_{,x})^2}{n^3},\\
  n_{,xx} = \frac{2(n_{,x})^2}{n} - \frac{-j_0 n_{,x} + n^3\frac{(V+\mathcal{E}')_{,x}}{mj_0}}{\nu}.
\end{gather*}
The momentum current is
\begin{gather*}
  j_{p} = \frac{j_0^2}{n} + \frac{P}{m} + \nu\frac{j_0 n_{,x}}{n^2}
\end{gather*}
:::

We will use the GPE equation of state:
\begin{gather*}
  \mathcal{E} = \frac{gn^2}{2}, \qquad  
  \mathcal{E}' = gn, \qquad  
  \mathcal{E}'' = g,
\end{gather*}
and a piecewise linear potential.

```{code-cell} ipython3
import scipy as sp
x = np.linspace(-0.1, 1.1)
V0 = 1.2
dV = 0.2
xp, fp = [0, 0.5, 1.0], [dV, V0, 0]
V = sp.interpolate.InterpolatedUnivariateSpline(xp, fp, k=1, ext='const')
plt.plot(x, V(x), label='$V(x)$')
plt.plot(x, V.derivative()(x), label="$V'(x)$")
plt.legend()
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp
import scipy as sp

class Flow:
    m = 1.0
    g = 1.0
    nu = 0.1  # Viscosity
    V0 = 0.2  # Potential barrier height
    L = 0.2   # Length of barrier

    n0 = 1.0  # Steady-state density on the left.
    v0 = 0.1  # Steady-state velocity on the left.
    dV = 0.1  # Potential drop (magnitude)
    
    # Boundary conditions on the left
    def __init__(self, **kw):
        for key in kw:
            if not hasattr(self, key):
                raise ValueError(f"Unknown {key=}")
            setattr(self, key, kw[key])
        self.init()
        
    def init(self):
        xp, fp = [0, self.L/2, self.L], [self.dV, self.V0, 0]
        self._V = sp.interpolate.InterpolatedUnivariateSpline(xp, fp, k=1, ext='const')
        self._dV = self._V.derivative()
        
    def get_Eint(self, n, d=0):
        if d == 0:
            return self.g * n**2 /2
        elif d == 1:
            return self.g * n
        elif d == 2:
            return self.g
        else:
            return 0
            
    def compute_dy_dx(self, x, y):
        n, u, du, F = y
        dn = -n*du/u
        ddu = n * (u*du + (self._dV(x) + self.get_Eint(n, d=2)*dn)/self.m) / self.nu
        dF = n * self._dV(x)
        return (dn, du, ddu, dF)

    def solve(self, dv0=0.0, x0_L=-0.1, x_L=2.0):
        x_span = (x0_L*self.L, x_L*self.L)
        y0 = (self.n0, self.v0, dv0, 0.0)
        res = solve_ivp(self.compute_dy_dx,  y0=y0, t_span=x_span, max_step=0.001)
        assert(res.success)
        res.x = res.t
        res.n, res.v, res.dv, res.F = res.y
        res.dn, dv, res.ddv, res.dF = self.compute_dy_dx(res.x, res.y)
        self.res = res
        return res

f = Flow()
res = f.solve()
```

### Analysis

We start with the two conservation laws for density and momentum far to the left and
right of the potential, assuming that $u_{,x}=0$ in these regions, and setting $F(x_L) =
0$ and $F(x_R) = F$:

\begin{gather*}
  n_Lu_L = n_Ru_R = j = \text{const.}\\
  mn_Lu_L^2 + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) = 
  mn_Ru_R^2 + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
  \frac{m j^2}{n_L} + n_L\mathcal{E}'(n_L) - \mathcal{E}(n_L) = 
  \frac{m j^2}{n_R} + F + n_R\mathcal{E}'(n_R) - \mathcal{E}(n_R)\\
\end{gather*}





\begin{gather*}
   nu^2 + \frac{P}{m} - ν u_{,x} 
   = nu^2 + \frac{F(x) + n\mathcal{E}' - \mathcal{E}}{m} - ν
   u_{,x}
   = \text{const.},\\
  P = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
  F(x) = \int_{x_0}^{x} n(x)V'(x)\d{x}.
\end{gather*}


Assuming that $u_{,x} = 0$ far from the potential (static flow), and setiing $V(x_R) =
0$, we have the following after integrating the energy current:
\begin{gather*}
\left.
  u\left(
    \frac{m}{2}nu^2 + nV + n\mathcal{E}'(n) - m\nu u_{,x}
  \right)
\right|^{x_R}_{x_L}
= - \int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x} \leq 0,\\
  \tfrac{m}{2}(u_{R}^2 - u_{L}^2) + V_L + \mathcal{E}'(n_R) -
  \mathcal{E}'(n_L) =
  -\frac{1}{j}\int_{x_L}^{x_R} m\nu u_{,x}^2 \d{x}
\end{gather*}
































## Correction Needed Below

The following analysis is incorrect.  It is missing the factor of $1/n$ discussed above,
and neglects to include a difference in the potential from one side to the other (DC
offset) required to maintain a persistent flow in the presence of dissipation.

## Steady State (Incorrect)

Now consider steady\-state flow $n(x, t) = n(x)$, $u(x, t) = u(x)$.  These become:
\begin{gather*}
  (nu)_{,x} = 0, \qquad
  muu_{,x} + V_{,x} + \mu'(n)n_{,x} = mνu_{,xx}.
\end{gather*}
The first equation establishes that the flux $Φ = nu$ is the same at every point, allowing us to exchange velocity $u$ for density $n$:
\begin{gather*}
  n = \frac{Φ}{u}, \qquad
  n_{,x} = -\frac{Φ}{u^2}u_{,x}.
\end{gather*}
Let's assume that $V(x)$ has compact support over $[0, L]$, then we can give boundary conditions $u(0) = u_0$, $u_{,x}(0) = 0$.  The second equation and be integrated once to give the non\-linear first\-order equation
\begin{gather*}
  mνu_{,x} = \frac{mu^2}{2} + V(x) + \overbrace{\mu\left(\frac{Φ}{u}\right)}^{\mu(n)}.
\end{gather*}

For general $V(x)$, this has no analytic solution, but we can compute a perturbative solution if we let $V(x) \rightarrow \lambda V(x)$ where $\lambda \ll 1$.  We note that $u(x) = u_0$ is a solution for $\lambda = 0$ and write
$$u(x) = u_0 + \lambda u_1(x) + \lambda^2 u_2(x) + O(\lambda^3).$$_Note: In what
follows, we revert back to the notation_ $u_1' = [u_1]_{,x}$_...

I originally wrote it this way, and have not changed it yet._  
As shown below, we originally tried expanding to first order in $\lambda$, but found this was insufficient, so we expand to second order in $\lambda$ obtaining  the following linear equations for $u_1$ and $u_2$:
\begin{align*}
  mνu_1' + \overbrace{\left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)}^{-mνc}u_1
  &= V(x),\\
  mνu_2' + \left(\frac{n_0}{u_0}\mu'(n_0) - mu_0\right)u_2
  &= \underbrace{\left(\frac{m}{2}
    + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2u_0^2}
  \right)}_{mνb}u_1^2,\\
  u_1' - c u_1
  &= V(x)/mν,\\
  u_2' - cu_2
  &= bu_1^2,
\end{align*}
where we have introduced the constants
\begin{gather*}
  c = \frac{u_0}{ν}-\frac{n_0\mathcal{E}''(n_0)}{mνu_0}
  = 
  \frac{u_0}{ν}-\frac{n_0\mu'(n_0)}{mνu_0},\\
  b = \frac{1}{2ν} + \frac{2n_0\mu'(n_0)+n_0^2\mu''(n_0)}{2mνu_0^2}
    = \frac{1}{2ν} + \frac{\bigl(nP'(n)\bigr)'|_{n=n_0}}{2mνu_0^2}.
\end{gather*}

## Do it!

To solve this problem, one needs to find a solution to the following inhomogeneous differential equation:
\begin{gather*}
  au_1'(x) - u_1(x) = f(x),
\end{gather*}
where $a>0$.  We are particularly interested in the limit where $a\rightarrow 0$, which will correspond to zero viscosity.

You should find that
\begin{gather*}
  u_1(x) = e^{x/a}\int_{x_0}^{x} d{y}\; f(y)e^{-y/a}
\end{gather*}

Notice that this has a non-analytic dependence $e^{x/a}$ on the parameter $a \propto ν$.  Show that, when expanding from the appropriate limit, all coefficients of the Taylor series for $e^{x/a}$ in powers of $a$ are zero. This makes it somewhat tricky to expand for small viscosities since all terms with such a factor will vanish.

Argue, however, that there may be an analytic *particular* solution $u_1(x, a)$ whose form can be found by expanding in powers of $a$.  We can then write the general solution as
\begin{gather*}
  u_1(x) = u_1(x, a) + Ce^{x/a}.
\end{gather*}
The non-analytic piece, thus describes transients (exponentially decaying in the appropriate direction) induced by the initial (boundary) conditions.

Show that this particular solution is:
\begin{gather*}
  u_1(x, a) = - \sum_{n=0}^{\infty} a^{n} f^{(n)}(x).
\end{gather*}
For example, this allows us to easily solve for $f(x) = \sin(kx)$:
\begin{gather*}
  -u_1(x, a) = \sin(x) + ak \cos(x) - (ak)^2 \sin(x) - (ak)^3 \cos(x) + \cdots\\
             = \Bigl(1 - (ak)^2 + (ak)^4 - \cdots\Bigr)\Bigl(\sin(kx) + ak \cos(kx)\Bigr)\\
             = \frac{\sin(kx) + ak \cos(kx)}{1 + (ak)^2},\\
   u_1(x) = Ce^{x/a} - \frac{\sin(kx) + ak \cos(kx)}{1+(ak)^2}.
\end{gather*}
We can easily check by differentiating that this satisfies the original equation.

## Force

The \(buoyant\) force acting on the potential $V(x)$ is
\begin{gather*}
  F = -\int_{0}^{L}d{x}\;nV'
    = -Φ\int_{0}^{L}d{x}\;\frac{V'}{u}
\end{gather*}
Returning to our original equation, we have
\begin{gather*}
  -nV'(x) = m\overbrace{nu}^{Φ}u_{,x} + \overbrace{nμ_{,x}}^{[P(n)]_{,x}} - mνnu_{,xx},\\
  F = \Bigl.(mΦu + P)\Bigr|_{0}^{L} - mν\int_0^{L}d{x}\;nu_{,xx}
\end{gather*}


# A Numerical Example

To make sure what we are doing makes sense, we give a concrete example.  To minimize the
chance for error, we implement the original equations in the form $y=(n, u, u', F)$,
allowing for a density\-dependent viscosity: 

\begin{gather*}
  (nu)' = 0, \qquad
  muu' + \Bigl(V + \mathcal{E}'(n)\Bigr)' = \frac{(mν(n)u')'}{n}
  = m\frac{ν'(n)n'u'+ν(n)u''}{n}\\
  \begin{pmatrix}
    n\\
    u\\
    u'\\
    F
  \end{pmatrix}'
  =
  \begin{pmatrix}
    -\frac{nu'}{u}\\
    u'\\
    \frac{n\Bigl(uu' + \frac{V'(x) + \mathcal{E}''(n)n'}{m}\Bigr)-ν'n'u'}{ν}\\
    nV'(x)
  \end{pmatrix}, \qquad
  y(0) = \begin{pmatrix}
    n_0\\
    u_0\\
    0\\
    0
  \end{pmatrix}.
\end{gather*}

```{code-cell} ipython3
%matplotlib inline
import matplotlib
from functools import partial
import numpy as np, matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (4,3)

m = 1.8
nu0 = 0.032
L = 1.3
n0 = 0.85
u0 = 0.773/2
Φ = n0*u0
g_m = 1.5
lam = 0.02

def mu(n, d=0):
    """Return the GPE chemical potential derivatives."""
    if d == 0:
        return m*g_m*n
    elif d == 1:
        return m*g_m+0*n
    elif d == 2:
        return 0*n

def mu(n, d=0):
    """Return the UFG chemical potential derivatives."""
    if d == 0:
        return m*g_m*n**(5/3)
    elif d == 1:
        return 5/3*m*g_m*n**(2/3)
    elif d == 2:
        return 10/9*m*g_m*n**(-1/3)

def mu_(n, d=0):
    """Return the GPE chemical potential derivatives."""
    if d == 0:
        return m*g_m*n**3
    elif d == 1:
        return 3*m*g_m*n**2
    elif d == 2:
        return 6*m*g_m*n

def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1, m*x*(x-L), 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1, m*(2*x-L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1, m*2, 0)

def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        m*np.sin(np.pi*(x-L)/L)**2, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*m*np.sin(2*np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        2*np.pi**2/L**2*m*np.cos(2*np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -4*np.pi**3/L**3*m*np.sin(2*np.pi*(x-L)/L), 0)

def V(x, d=0):
    """Return the external potential or derivative."""
    s = np.sin(np.pi*(x-L)/L)
    ds = np.cos(np.pi*(x-L)/L)*np.pi/L
    ds2 = -np.sin(np.pi*(x-L)/L)*(np.pi/L)**2
    ds3 = -np.cos(np.pi*(x-L)/L)*(np.pi/L)**3

    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        m*s**4, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        m*4*s**3*ds, 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        m*(12*s**2*ds**2 + 4*s**3*ds2), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        m*(24*s*ds**3 + 36*s**2*ds*ds2 + 4*s**3*ds3), 0)


def nu(n, nu0=1.0, gamma=0.0, d=0):
    """Return the density-dependent viscosity."""
    if d == 0:
        return nu0*n**gamma
    elif d == 1:
        return nu0*gamma*n**(gamma-1)

def compute_dy_dt(x, y, lam=lam, nu0=nu0, 
                  u0=u0, n0=n0, m=m):
    n, u, du, F = y
    dn = -n*du/u
    dV = lam*V(x, d=1)
    #ddu = (m*u*du + dV + mu(n, d=1)*dn)/m/nu
    nu_, dnu_ = [nu(n, nu0=nu0, d=d) for d in [0, 1]]
    ddu = (n*(m*u*du + dV + mu(n, d=1)*dn)/m - du*dnu_*dn)/nu_
    dF = n*dV
    return (dn, du, ddu, dF)

y0 = (n0, u0, 0, 0)

res = solve_ivp(
    partial(compute_dy_dt, nu0=0.2, lam=0.5),
    y0=y0,
    t_span=(-L, 5*L), max_step=0.001)
x = res.t
n, u, du, F = res.y
assert np.allclose(n*u, n0*u0)
fig, ax = plt.subplots()
ax.plot(x/L, n/n0, label="$n$")
ax.plot(x/L, u/u0, label="$u$")
ax.axvspan(0, 1, color='y', alpha=0.5)
ax.set(xlabel="$x/L$", ylabel="$n/n_0$", title=f"{n[0]=}, {n[-1]=}", );
plt.legend(loc="upper left")
ax = plt.twinx()
ax.plot(x/L, F, '--C1', label="$F$")
ax.legend(loc="lower right")
ax.set(ylabel="$F$");
```

Notice that the potential causes the density to change, then the viscosity allows this to relax outside of the shaded region where the potential is active.  We now check some relationships.

```{code-cell} ipython3
# Check that the potential and mu derivatives are correct
fig, axs = plt.subplots(1, 2)
ax = axs[0]
ax.plot(x/L, V(x), 'C0')
ax.plot(x/L, V(x, d=1), 'C1-')
ax.plot(x/L, np.gradient(V(x), x), 'C1:')
ax.plot(x/L, V(x, d=2), 'C2-')
ax.plot(x/L, np.gradient(V(x, d=1), x), 'C2:')
ax.plot(x/L, V(x, d=3), 'C3-')
ylim = ax.get_ylim()
ax.plot(x/L, np.gradient(V(x, d=2), x), 'C3:')
ax.set(ylim=ylim)

ax = axs[1]
ns = np.linspace(n0/2, 2*n0)
ax.plot(ns/n0, mu(ns), 'C0')
ax.plot(ns/n0, mu(ns, d=1), 'C1-')
ax.plot(ns/n0, np.gradient(mu(ns), ns), 'C1:')
ax.plot(ns/n0, mu(ns, d=2), 'C2-')
ax.plot(ns/n0, np.gradient(mu(ns, d=1), ns), 'C2:')
```

```{code-cell} ipython3
from functools import partial

def get_F(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
    y0 = (n0, u0, 0, 0)
    dy_dt = partial(compute_dy_dt, lam=lam, nu=nu,
                    u0=u0, n0=n0, m=m)
    res = solve_ivp(dy_dt,
                    y0=y0, t_span=(0, L), 
                    method='BDF', 
                    rtol=tol, atol=tol)
    x = res.t
    n, u, du, F = res.y
    return F[-1]

lams = np.linspace(0, 0.5)[1:-1]
nus = np.linspace(0, 0.5)[1:-1]
Fs_lam = [get_F(lam=lam) for lam in lams]
Fs_nu = [get_F(nu=nu) for nu in nus]

fig, axs = plt.subplots(1, 2, figsize=(5, 2))
axs[0].plot(lams**2, Fs_lam)
axs[0].set(xlabel=r"$\lambda^2$", ylabel="$F$")
axs[1].plot(nus, Fs_nu)
axs[1].set(xlabel=r"$\nu$", ylabel="$F$");
```

Now let's check the small $λ$ approximations:

```{code-cell} ipython3
from functools import partial
from scipy.integrate import cumtrapz

c = u0/nu - n0*mu(n0, d=1)/m/nu/u0
b = (1 + (2*n0*mu(n0, d=1) + n0**2*mu(n0, d=2))/m/u0**2)/2/nu
a = c*nu
print(f"{a=}, {nu*b}")


def f(x, us):
    u1, u2 = us
    du1 = V(x)/m/nu + c*u1
    du2 = b*u1**2 + c*u2
    return (du1, du2)
res1 = solve_ivp(f, y0=[0, 0], t_span=(-L, 2*L), t_eval=x,
                 method='BDF', atol=1e-8, rtol=1e-8)
u1, u2 = res1.y
n1 = Φ/(u0+lam*u1)
F1 = cumtrapz(lam*n1*V(x, d=1), x, initial=0)

# Small nu approximations
u10 = -V(x)/m/a
u11 = -V(x, d=1)/m/a**2
u12 = -V(x, d=2)/m/a**3

fig, axs = plt.subplots(2, 2, figsize=(10, 5))
ax = axs[0,0]
ax.plot(x/L, u1/u0, label='$u_1/u_0$')
ax.plot(x/L, u10/u0, label=r"small $\nu$ approx")
ax.plot(x/L, (u-u0)/u0/lam, label="full")

ax.set(xlabel=r"$x/L$", ylabel="$u_1/u_0$")
ax.legend()

ax = axs[0, 1]
ax.plot(x/L, F1, label='$F_1$')
ax.plot(x/L, F, label="full")
ax.legend()

ax = axs[1, 0]
ax.plot(x/L, u1/u0, 'C0-')
ax.plot(x/L, u10/u0, 'C0:')
ax.plot(x/L, (u1-u10)/nu/u0, 'C1-')
ax.plot(x/L, u11/u0, 'C1:')
ax.plot(x/L, u12/u0, 'C2:')
ylim = ax.get_ylim()
ax.plot(x/L, (u1-u10-nu*u11)/nu**2/u0, 'C2-', scaley=False)
ax.set(ylim=ylim)

ax = axs[1, 1]
ax.plot(x/L, (u-u0)/lam/u0, 'C0-')
ax.plot(x/L, u1/u0, 'C0:')
ax.plot(x/L, (u-u0-lam*u1)/lam**2/u0, 'C1-')
ax.plot(x/L, u2/u0, 'C1:')
```

The following code demonstrates that the order $λ$ solution $u_1(x)$ gives the correct force when integrated. *(There is no need to compute $u_2(x)$.)*

```{code-cell} ipython3
def get_F1(lam=lam, nu=nu, u0=u0, n0=n0, m=m, L=L, tol=1e-7):
    a = u0 - n0*mu(n0, d=1)/m/u0
    c = a/nu
    def f(x, y):
        u1, F1 = y
        du1 = V(x)/m/nu + c*u1
        dF = -lam**2 * n0 / u0 * u1 * V(x, d=1)
        return (du1, dF)

    res1 = solve_ivp(f,
        y0=[0, 0], t_span=(0, L),
        method='BDF', atol=tol, rtol=tol)
    u1, F1 = res1.y
    return F1[-1]


lams = np.linspace(0.01, 0.3, 20)
qs = [get_F1(lam=lam, tol=lam**2 * 1e-8) /
      get_F(lam=lam, tol=lam**2 * 1e-8) for lam in lams]
plt.plot(lams, qs, '-+')
plt.plot([0], [1], 'o')
plt.gca().set(xlabel=r'$\lambda$', ylabel='$F_1/F$');
```

Thus, the complete solution to order $λ$ is to solve the following problem:
\begin{gather*}
  u_1'(x) + cu_1(x) = V(x)/mν,\qquad
  c = \frac{u_0 - P'_0/mu_0}{ν}\\
  F %= λ\int_{0}^{L}d(x)\;\frac{Φ V'(x)}{u_0+λu_1}
    = \frac{-λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1(x) V'(x)
    = \frac{λ^2n_0}{u_0} \int_{0}^{L}d{x}\;u_1'(x) V(x)
\end{gather*}
Note that the coefficient $c$ here is generally less than zero.  We can see this by noting that the speed of sound in a material is $c_s$ is:
\begin{gather*}
  c_s = \sqrt{\frac{n_0\mu'(n_0)}{m}} = \sqrt{\frac{P'(n_0)}{m}},\\
  c = \frac{u_0-c_s^2/u_0}{\nu} = \frac{u_0}{\nu}\left(1 - \frac{c_s^2}{u_0^2}\right).
\end{gather*}
Thus, for subsonic flow $u_0<c_s$, $c<0$.  We keep this sign convention to avoid many other signs below.  The homogeneous solution is $u_1(x) = Ae^{cx}$.  To find the general solution, there are several options:

- One is the method of [variation of parameters](https://en.wikipedia.org/wiki/Variation_of_parameters) can be used, promoting $A → A(x)$ to find
  \begin{gather*}
    u_1(x) = A(x)e^{cx}, \qquad
    A'(x) e^{cx} = V(x)/{mν}, \qquad
    A(x) = \frac{1}{mν}\int_0^xd{y}\; e^{-cy}V(y), \\
    u_1(x) = \frac{1}{mν}\int_0^xd{y}\; e^{c(x-y)}V(y),\\
    F = \frac{-λ^2n_0}{mu_0}\int_{0}^{L}d{x}\;\int_0^xd{y}\; e^{c(x-y)}V(y)V'(x).
  \end{gather*}
  The integral is over the lower triangle $y<x$ of the box $(x, y) ∈[0, L]×[0,L]$. To better express the dimensionality of the solution, we write the potential as $V(x) = V_0U(x/L)$ so that the shape is described by the dimensionless function $U(\xi)$, and we introduce the dimensionless quantity $\kappa \propto \nu$:
  \begin{gather*}
    λV(x) = λV_0U(x/L), \qquad
    \kappa = \frac{1}{-cL} = \frac{\nu}{-aL},\\
    \begin{aligned}
    F &= \frac{-λ^2n_0V_0^2L}{mu_0}
    \underbrace{\int_{0}^{1}d{x}\;\int_0^xd{y}\; e^{(y-x)/\kappa}U(y)U'(x)}_{-\kappa^2 f(\kappa, U)},\\
     &= \frac{Lu_0\kappa^2}{\nu}\frac{λ^2V_0^2n_0}{mu_0^2}f(\kappa, U).
     \end{aligned}
  \end{gather*}
  While formally correct, this solution obscures the viscosity dependence.  We will see below that the normalization chosen here for $f(\kappa, U)$ makes sense in the limit of small viscosity:
  \begin{gather*}
    f(\kappa, U) = 
    -\int_{0}^{1}d{x}\;\int_0^xd{y}\; \frac{e^{(y-x)/\kappa}}{\kappa^2}U(y)U'(x)
  \end{gather*}

- Another is formal manipulation.  We first introduce the parameter $a=cν$ which is independent of the viscosity, then solve algebraically treating the derivative operator as a matrix:
  \begin{align*}
    νu_1'(x) &= au_1(x) + V(x)/m, \\
    u_1(x) &= \frac{1}{\left(1-\frac{ν}{a}\frac{d}{d{x}}\right)}\frac{-V(x)}{am},\\
    &= -\frac{V(x)}{am} - \frac{νV'(x)}{a^2m} - \frac{ν^2V''(x)}{a^3m} + O(ν^3).
  \end{align*}
  This is nice because it immediately allows us to expand in powers of $ν$ by expanding $1/(1-x) =1+x+x^2+x^3+\cdots$:
  \begin{gather*}
    F = \frac{λ^2n_0 V_0^2}{mau_0} \int_{0}^{1}d{x}\;\left(
      U(x)U'(x) - \kappa[U'(x)]^2 + \kappa^2U''(x)U'(x) + O(ν^3)
    \right)
  \end{gather*}
  Since $U(0)=U(1)=0$, the first term, which is a total derivative $(U^2)'/2$, vanishes. This gives:
  \begin{align*}
    F &= \frac{mλ^2n_0}{au_0}\frac{V_0^2}{m^2}(-\kappa)\overbrace{
    \Biggl(
      \int_{0}^{1}d{x}\;[U'(x)]^2
      -\kappa\left.\frac{[U'(x)]^2}{2}\right|_{0}^{1}
      +\kappa^2\int_{0}^{1}d{x}\;U'(x)U'''(x)
      + O(ν^3)
    \Biggr)}^{f(\kappa, U)}\\
    &= \frac{mλ^2n_0}{au_0}\frac{ν}{a} \frac{V_0^2}{m^2L}
       f(\kappa, U).
  \end{align*}
  This justifies the form of $f(\kappa, U)$ chosen above.
  If we choose a symmetric potential so that $|U'(0)| = |U'(1)|$, then the second term vanishes, allowing us to precisely define what a small viscosity means:
  \begin{gather*}
      \frac{ν^2}{a^2}
      \frac{\int_{0}^{L}d{x}\;V'(x)V'''(x)}
      {\int_{0}^{L}d{x}\;[V'(x)]^2}
      =
      \kappa^2
      \frac{\int_{0}^{1}d{x}\;U'(x)U'''(x)}
      {\int_{0}^{1}d{x}\;[U'(x)]^2}
      \ll 1.
  \end{gather*}

```{code-cell} ipython3
u10 = -V(x)/a/m
u11 = -1/a * V(x, d=1)/a/m
u12 = -1/a**2 * V(x, d=2)/a/m
u13 = -1/a**3 * V(x, d=3)/a/m

plt.plot(x/L, (u1 - u10)/nu)
plt.plot(x/L, (u1 - u10 - nu*u11)/nu**2)
plt.plot(x/L, (u1 - u10 - nu*u11 - nu**2*u12)/nu**3)

coeff = m*lam**2*n0/a/u0*nu/a
t1 = np.trapz(V(x, d=1)**2, x)/m**2
t2 = np.trapz(V(x, d=1)*V(x, d=2), x)/m**2/a
t3 = np.trapz(V(x, d=1)*V(x, d=3), x)/m**2/a**2

print(f"{t1=}, {t2=}, {t3=}")

F_ = get_F(tol=1e-13)/coeff
F1_ = get_F1(tol=1e-13)/coeff
print(F_, (F_-F1_)/lam)
print(F1_, t1)
print((F1_ - t1)/nu,
      (F1_ - t1 - nu*t2)/nu**2,
      (F1_ - t1 - nu*t2 - nu**2*t3)/nu**3,
      )
```

We now replace the parameter $a$ with the speed of sound $c_s$.  The final result, to leading order in $λ$ is:

\begin{gather*}
  F = \overbrace{\frac{ν}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)^2}\frac{λ^2V_0^2n_0}{mu_0^2}
    f(\kappa, U), \qquad
    \mathrm{Re} = \frac{u_0L}{ν}, \qquad
    \kappa = \overbrace{\frac{\nu}{u_0L}}^{1/\mathrm{Re}}\frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}, \qquad
    c_s = \sqrt{\frac{P'(n_0)}{m}},\\
  F = \frac{1}{\left(\frac{c_s^2}{u_0^2}-1\right)}\frac{λ^2V_0^2n_0}{mu_0^2}
    \Bigl(\kappa f(\kappa, U)\Bigr)
\end{gather*}

Here we have expressed the results in terms of the dimensionless Reynold's number $\mathrm{Re}$ and the ratio of the speed of sound $c_s$ to the flow $u_0$.  Note that the force diverges as the flow speed approaches the speed of sound.  The dependence on the shape appears with the dimensionless parameter $\kappa = ν/(-aL)$:
\begin{align*}
  f(\kappa, U) &= \frac{-1}{\kappa^2}\int_0^{L} dx\int_0^{x}d{y}\;e^{(y-x)/\kappa}
  U(y)U'(x),\\
  &\approx \sum_{n=1}^{\infty}(-\kappa)^{n-1}
  \int_0^1d{x}\;U^{(n)}(x)U'(x).
\end{align*}

**Caution:** The expansion in $\kappa \propto \nu$ here justifies our choice for the form of $f(\kappa, U)$ and the corresponding separation of dimensionful parameters, but it turns out that this series does not converge to the solution.  We now demonstrate this.  Consider for example, the dimensionless potential
\begin{align*}
  U(x) &= \sin(\pi(x-1)), &
  f(\kappa, U) &= \frac{\pi^2}{2}
    \frac{1+\pi^2\kappa^2 - 2\kappa(1+e^{-1/\kappa})}
         {(1+\pi^2\kappa^2)^2},\\
  U(x) &= \sin^2(\pi(x-1)), &
  f(\kappa, U) &= \frac{\pi^2}{2}
  \frac{\Bigl(
    1 + 4\pi^2\kappa^2 
    + 8\pi^2\kappa^3(1 - e^{-1/\kappa})\Big)}
    {(1+ 4\pi^2\kappa^2)^2}
 \end{align*}
This [solution](https://www.wolframalpha.com/input?i=-integrate%28integrate%28e%5E%7B%28y-x%29%2Fkappa%7D*sin%28pi*%28y-1%29%29*cos%28pi*%28x-1%29%29*pi%2Fkappa%5E2%2C+y%3D0..x%29%2C+x%3D0..1%29) contains the non\-analytic term $e^{-1/\kappa}$ whose Taylor series for positive $\kappa$ is zero.  Thus, the series expansion converges to the incorrect form of:
\begin{gather*}
  \sum_{n=1}^{\infty}(-\kappa)^{n-1}
  \int_0^1d{x}\;U^{(n)}(x)U'(x) =
  \frac{\pi^2}{2}\left(1 - (\pi\kappa)^2 + (\pi\kappa)^4 - \cdots\right)\\
  = \frac{\pi^2}{2}\frac{1}{1+\pi^2\kappa^2}
\end{gather*}

fails to converge.

```{code-cell} ipython3
def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        np.sin(np.pi*(x-L)/L), 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*np.cos(np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        -np.pi**2/L**2*np.sin(np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -np.pi**3/L**3*np.cos(np.pi*(x-L)/L), 0)

def get_f(nu, u0=u0, n0=n0, V0=1, L=L, m=m, lam=lam):
    cs2 = n0*mu(n0, d=1)/m
    coeff = nu/u0/L / (cs2/u0**2-1)**2 * lam**2 * V0**2 *n0/m/u0**2
    F1 = get_F1(nu=nu, u0=u0, n0=n0, L=L, m=m, lam=lam, tol=1e-8)
    kappa = nu/u0/L/(cs2/u0**2-1)
    return kappa, F1/coeff

nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])

xs = (np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
         * (1 + xs - 2*kappas*(1+np.exp(-1/kappas)))
         / (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
         * (1 + xs - 2*kappas*(1+0*np.exp(-1/kappas)))
         / (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
```

This figure shows that the integral gives the correct formula, and that for large $\kappa$ the non\-analytic piece is correct.  **It also shows that we have made a mistake in our series approximation.  Please fix!**

Here is $U(x) = \sin^2\bigl(\pi(x-1)\bigr)$:

```{code-cell} ipython3
def V(x, d=0):
    """Return the external potential or derivative."""
    if d == 0:
        return np.where(abs(2*x/L-1) < 1,
        np.sin(np.pi*(x-L)/L)**2, 0)
    elif d == 1:
        return np.where(abs(2*x/L-1) < 1,
        np.pi/L*np.sin(2*np.pi*(x-L)/L), 0)
    elif d == 2:
        return np.where(abs(2*x/L-1) < 1,
        2*np.pi**2/L**2*np.cos(2*np.pi*(x-L)/L), 0)
    elif d == 3:
        return np.where(abs(2*x/L-1) < 1,
        -4*np.pi**3/L**3*np.sin(2*np.pi*(x-L)/L), 0)

nus = np.linspace(0.1, 10.0)
kappas, fs = np.transpose([get_f(nu=nu) for nu in nus])

xs = (2*np.pi*kappas)**2
plt.plot(xs, fs, '.')
plt.plot(xs, np.pi**2/2/(1+xs), label='Series')
plt.plot(xs, np.pi**2/2
         * (1 + xs + 2*xs*kappas*(1-np.exp(-1/kappas)))
         / (1+xs)**2, label='Full integral')
plt.semilogy(xs, np.pi**2/2
         * (1 + xs + 2*xs*kappas*(1 - 0*np.exp(-1/kappas)))
         / (1+xs)**2, label='Integral w/o non-analytic piece')
plt.legend()
plt.gca().set(xlabel=r"$\pi^2\kappa^2$", ylabel="$f(\kappa, U)$");
```

```{code-cell} ipython3
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
    return sympy.sin(sympy.pi*(x-1))

integrand = -1/k**2 * sympy.exp((y-x)/k)*U(x).diff(x)*U(y)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1)).simplify()
display(f.factor())
display(sympy.series((1+(sympy.pi*k)**2 - 2*k)/(1+(sympy.pi*k)**2), k, 0, 10))

f_series = sum(sympy.simplify(
        (U(x).diff(*(x,)*n)*U(x).diff(x))
        .integrate((x, 0, 1))*(-k)**(n-1))
    for n in range(1, 9)[::2])
display(f_series)
```

If we choose $U(\xi) = U(1-\xi)$ to be symmetric, then the even terms vanish, and we can write:
\begin{gather*}
  f(\kappa, U) = \sum_{n=0}^{\infty}c_{n}\kappa^{2n}, \qquad
  c_n = \int_0^1d{\xi}\;U^{(2n+1)}(\xi)U'(\xi).
\end{gather*}
Here are a couple of special cases.  Note: For numerical work, it might be better to use a potential that is sufficiently flat at the boundaries, hence the higher power.  This is what I have used in the code.  (Please check these.)
\begin{gather*}
  f(\kappa, U) = \pi^2g(\pi^2\kappa^2),\qquad
  \kappa = \frac{\nu/u_0L}{\Bigl(1-\frac{c_s^2}{u_0^2}\Bigr)},\\
  \begin{aligned}
  U(\xi) = \sin\bigl(\pi(\xi-1)\bigr)&:& 
  c_n &= \frac{\pi^2}{2}(-\pi^2)^{n}, &
  g(x) &= \frac{1}{2(1+x)},\\
  U(\xi) = \sin^4\bigl(\pi(\xi-1)\bigr)&:& 
  c_n &= \frac{\pi^2}{2}(4^{2n-1} + 4^{n})(-\pi^2)^{n}, &
  g(x) &= \frac{68x + 5}{8(1+4x)(1+16x)}.
  \end{aligned}
\end{gather*}
Thus, for a sinusoidal potential (please check):
\begin{gather*}
  \begin{aligned}
  \frac{F}{L} &= 
  \frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa}{2(1+\pi^2\kappa^2)}, &
  U(\xi) &= \sin\bigl(\pi(\xi-1)\bigr),\\
  \frac{F}{L} &= 
  \frac{n_0 λ^2V_0^2}{mu_0^2}\frac{1}{\Bigl(1 - \frac{c_s^2}{u_0^2}\Bigr)}\frac{\pi^2\kappa(68\pi^2\kappa^2+5)}{8(1+4\pi^2\kappa^2)(1+16\pi^2\kappa^2)}, &
  U(\xi) &= \sin^4\bigl(\pi(\xi-1)\bigr),\\
  \end{aligned}
\end{gather*}

These can also be obtained by :

\begin{gather*}
  %-integrate(integrate(e^{(x-y)/kappa}*sin(pi*(y-1))*cos(pi*(x-1))*pi/kappa^2, y=0..x), x=0..1)
  f(\kappa, U) = \frac{-1}{\kappa}\int_0^L d{x}\int_0^x d{y}\; e^{(x-y)/\kappa} U'(x)U(y),\\
  \frac{\pi^2}{2}\frac{(1 + 2 (1 + e^{1/\kappa}) κ + π^2 κ^2)}{(1 + \pi^2 \kappa^2)^2}
\end{gather*}

```{code-cell} ipython3
import sympy
x, y, k = sympy.var('x,y,kappa', positive=True)
def U(x):
    return sympy.sin(sympy.pi*(x-1))**2

integrand = -1/k*sympy.exp((x-y)/k)*U(x).diff(x)*U(y)
display(integrand)
f = integrand.integrate((y, 0, x)).integrate((x, 0, 1))
f.expand().factor()
```

```{code-cell} ipython3
sympy.series(sympy.exp(1/k), k, 0, 5)
```

*(I found these from the expansion https://oeis.org/A092812)*:

\begin{gather*}
  f(\kappa, U) = \pi^2 \sum_{n=0}^{\infty}
  \left(\frac{16^{n}}{8} + \frac{4^{n}}{2}\right)(-\pi^2\kappa^2)^{n}
  = \pi^2 g(\pi^2\kappa^2)
  , \\
  g(x) %= \frac{1 + 16x + 24x^2}{(1+4x)(1+16x)} - \frac{3}{8}
       = \frac{68x + 5}{8(1+4x)(1+16x)}
\end{gather*}

```{code-cell} ipython3
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))

def get_c(n):
    return k**(2*n)*(V.diff(*(x,)*(2*n+1))*V.diff(x)).integrate((x, 0, 1))

[get_c(n) for n in [0, 1, 2, 3, 4, 5]]
```

```{code-cell} ipython3
import numpy as np
from scipy.integrate import quad
import sympy
x, k = sympy.var('x, kappa', real=True)
V = sympy.sin(sympy.pi*(x-1))**4
def get_c(n):
    f = sympy.lambdify([x], V.diff(*(x,)*(n+1))*V.diff(x))
    return int(np.round(quad(f, 0, 1)[0]) / np.pi**(n+2)) * k**(n)*sympy.pi**(n+2)
get_c(4)
4**(0-1)/2 + 2**(0-1), 5/8

sympy.simplify((1+16*x+24*x**2)/(1+4*x)/(1+16*x) - sympy.S(3)/8).factor()
```

```{code-cell} ipython3
from scipy.integrate import quad
def f(x):
    return V(x, d=1)**2

cs = np.sqrt(n0*mu(n0, d=1)/m)
Re = u0*L/nu
tmp = lam**2 * quad(f, 0, L)[0]
print(a/u0, (1 - (cs/u0)**2))
print(
  (1/Re) * (1/(1-(cs/u0)**2)**2) * (n0*L/m/u0**2) * tmp,
  get_F())
```

```{code-cell} ipython3
from scipy.integrate import quad
def f(x):
    return V(x)*V(x,d=1)**2

tmp = quad(f, 0, L)[0]
lam**2*n0*nu/m/a**2/u0 * tmp, get_F()
```

# PDE

```{code-cell} ipython3
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
%load_ext autoreload
from gpe.contexts import FPS
from IPython.display import clear_output
%autoreload 2
```

```{code-cell} ipython3
import viscosity
import itertools
for (periodic, fft) in itertools.product(*([[False, True]]*2)):
    f = viscosity.FlowStep(V_dV=0.2, V_sigma=1.0, V0=1.0, n0=1.0, v0=0.0,
                           nu=1.0, x0=0, x1=2*np.pi, Nx=200,
                           periodic=periodic, fft=fft, 
                           ghost_points=1)
    y0 = f.get_y0()

    x = f.x
    y = np.sin(x)
    dy = np.cos(x)
    ddy = -np.sin(x)

    #assert np.allclose(f.get_derivative(y)[1:-1], dy[1:-1], rtol=1.7e-4)
    #assert np.allclose(f.get_derivative(y, d=2)[1:-1], ddy[1:-1], rtol=9e-5)
#plt.plot((f.get_derivative(y, d=1) / dy - 1)[1:-1])
#plt.plot((f.get_derivative(y, d=2) / ddy - 1)[1:-1])

#plt.plot(f.x, f.get_derivative(y, d=1))
#plt.plot(f.x, dy)
#plt.ylim(-1, 1)
```

### Dam Breaking

```{code-cell} ipython3
from scipy.integrate import solve_ivp

f = viscosity.Flow(x0=-1.0, x1=1.0, Nx=100)
assert np.allclose(f._D2x, -f._Dx.T @ f._Dx)
x = f.x

n0 = np.where(x < 0, 1.0, 0.00001)
v0 = 0*n0
y0 = f.pack(n0, v0)

t_span = (0, 0.1)

res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF", rtol=1e-12, atol=1e-12)
if not res.success:
    print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])

plt.plot(x, nF)
```

### Burger's Equation

```{code-cell} ipython3
from importlib import reload
reload(viscosity)

f = viscosity.FlowBurgers(x0=0, x1=2*np.pi, Nx=256, nu=0.07)#, fft=True)
n0, u0 = f.get_exact_sol(f.x, t=0)
y0 = f.pack(n0*0, u0)
dx = 2*np.pi/99
dt = f.nu * dx
T = dt * 100
print(T)
t_span = (0, T)

res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF")
if not res.success:
    print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])

plt.plot(f.x, vF)
plt.plot(f.x, nF)
n_, v_ = f.get_exact_sol(f.x, t=T)
plt.plot(f.x, v_)
```

```{code-cell} ipython3
fig, ax = plt.subplots()
ax1 = ax.twinx()
axs = [ax, ax1]

n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
    n, u = f.unpack(y)
    n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
    u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]

skip = 1
fps = FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig, embed=True, filename="tmp.m4v")
for t, y in fps:
    n, u = f.unpack(y)
    [ax.cla() for ax in axs]
    ax.plot(f.x, n)
    ax1.plot(f.x, u, '--C1')
    ax.set(title=f"{t=:.4f}", ylim=n_range)
    ax1.set(ylim=u_range)
```

```{code-cell} ipython3
def step0(u, dt=dx * f.nu):
    nx = len(u)
    un = u
    u = u.copy()
    for i in range(1, nx-1):
        u[i] = (un[i] + dt * ((-un[i] * (un[i] - un[i-1]) / dx
             + f.nu * (un[i+1] - 2 * un[i] + un[i-1]) / dx**2)))
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + f.nu * dt / dx**2 *\
           (un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]
    return u

def step0(u, dt=dx * f.nu):
    nx = len(u)
    un = u
    u = u.copy()
    du = u.copy()
    ddu = u.copy()
    for i in range(1, nx-1):
        du[i] = (un[i+1] - un[i]) / dx
        ddu[i] = (un[i+1] - 2 * un[i] + un[i-1]) / dx**2
    for i in range(1, nx-1):
        du[i] = (un[i] - un[i-1]) / dx
        #ddu[i] = (un[i+1] - 2 * un[i] + un[i-1]) / dx**2
        u[i] = (un[i] + dt * (-un[i] * du[i] + f.nu * ddu[i]))
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + f.nu * dt / dx**2 *\
           (un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]
    return u


def step1(u, dt=dx * f.nu):
    du = f.get_derivative(u)
    ddu = f.get_derivative(u, d=2)

    du_dt = (-u * du  + f.nu * ddu)
    return u + du_dt * dt
```

```{code-cell} ipython3
def step(u, dt=dx * f.nu):
    du_dt = f.unpack(f.compute_dy_dt(0, f.pack(u,u)))[1]
    return u + du_dt*dt

us = [u0]
t = 0
for n in range(100):
    t += dt
    us.append(step0(us[-1]))
plt.plot(f.x, us[-1])
plt.plot(f.x, f.get_exact_sol(f.x, t=t)[1])
```

```{code-cell} ipython3
from scipy.integrate import solve_ivp

f = viscosity.FlowBurgers(x0=-1.0, x1=1.0, Nx=100)
x = f.x

n0 = np.where(x < 0, 1.0, 0.5)
v0 = 0*n0
y0 = f.pack(n0, v0)

t_span = (0, 1)

res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=t_span, method="BDF", rtol=1e-12, atol=1e-12)
if not res.success:
    print(res.message)
print(res.t.max())
nF, vF = f.unpack(res.y[:, -1])

plt.plot(x, nF)
```

```{code-cell} ipython3
fig, ax = plt.subplots()
ax1 = ax.twinx()
axs = [ax, ax1]

n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
    n, u = f.unpack(y)
    n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
    u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]

skip = 1
fps = FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig, embed=True, filename="tmp.m4v")
for t, y in fps:
    n, u = f.unpack(y)
    [ax.cla() for ax in axs]
    ax.plot(f.x, n)
    ax1.plot(f.x, u, '--C1')
    ax.set(title=f"{t=:.4f}", ylim=n_range)
    ax1.set(ylim=u_range)
```

```{code-cell} ipython3
fig, axs = plt.subplots(2, 1, height_ratios=(3, 1), gridspec_kw=dict(hspace=0))
ax = axs[0]
axV = axs[1]
ax1 = ax.twinx()
axs = [ax, axV, ax1]

n_range = [0, 0]
u_range = [0, 0]
for t, y in zip(res.t, res.y.T):
    n, u = f.unpack(y)
    n_range = [min(n.min(), n_range[0]), max(n.max(), n_range[1])]
    u_range = [min(u.min(), u_range[0]), max(u.max(), u_range[1])]

skip = 1
for t, y in FPS(zip(res.t[::skip], res.y.T[::skip]), timeout=1000, display=fig):
    n, u = f.unpack(y)
    [ax.cla() for ax in axs]
    ax.plot(f.x, n)
    axV.plot(f.x, f.get_Vext(f.x, t))
    ax1.plot(f.x, u, '--C1')
    ax.set(title=f"{t=:.4f}", ylim=n_range)
    ax1.set(ylim=u_range)
```

```{code-cell} ipython3
import viscosity
f = viscosity.Flow(V_dV=-0.2, V_sigma=1.0, V0=1.0, n0=2.0, v0=0.0)
f = viscosity.Flow(V_dV=0, V_sigma=0.1, V0=0, n0=2.0, v0=0.0, V_w=1, nu=0.1, nu_alpha=1)
y0 = f.solve_ivp(x0=-3, x1=3, Nx=200, ghost=1)
x = f.x
n0 = np.maximum(0, 1 - f.get_Vext(x, t=0))/f.g + 1e-5
v0 = 0*n0
y0 = f.pack(n0, v0)

from scipy.integrate import solve_ivp

t0 = 0
t1 = 0.5

res = solve_ivp(f.compute_dy_dt, y0=y0, t_span=(t0, t1), method='BDF')
if not res.success:
    print(res.message)
print(res.t.max())
```

$$
  f(x \pm h) = f \pm h f' + \frac{h^2}{2} f'' \pm \frac{h^3}{6} f'''\\
  8\frac{f(x + h) - f(x-h)}{12h} - \frac{f(x + 2h) - f(x-2h)}{12h} = f'\\
$$

```{code-cell} ipython3
%pylab inline
from scipy.integrate import solve_ivp
from gpe.contexts import FPS

Lx = 10.0
Nx = 256//2
dx = Lx/Nx
x = np.arange(Nx) * dx - Lx/2
k = 2*np.pi * np.fft.fftfreq(Nx, dx)

# Fourier-space representation of derivatives
D = {1: None, 2: None}
D[1] = 1j*k
assert Nx % 2 == 0
D[1][Nx//2] = 0  # Kill highest unmatched frequency
D[2] = -k**2

def diff(f, d=1):
    """Compute the derivative of f."""
    global D
    return np.fft.ifft(D[d] * np.fft.fft(f)).real

def pack(n, u):
    return np.concatenate([n, u])
    
def unpack(y):
    return np.reshape(y, (2, Nx))

def get_V_m(t, x, d=0):
    """Return the dth derivative of V/m."""
    global V0_m, sigma, dV_m_dx
    V_m = V0_m * np.exp(-(x/sigma)**2/2)
    if d == 0:
        dV = dV_m_dx*x
        return V_m + dV - np.min(dV)  # Constant to make V_min = 0
    elif d == 1:
        return - x * V_m / sigma**2 + dV_m_dx

def get_E_m(n, d=0):
    """Return the dth derivative of E/m where E is the energy density."""
    global g_m
    if d == 0:
        return g_m * n**2 / 2 
    elif d == 1:
        return g_m * n
    elif d == 2:
        return g_m

def get_nu_n(n, d=0):
    """Return the dth derivative of the vicosity diff(nu(n), n, d)/n, divided by n."""
    if d == 0:
        return nu # * n / n
    elif d == 1:
        return nu / n
    else:
        return 0
        
def compute_dy_dt(t, y):
    global get_V_m, get_E_m, get_nu_n
    n, u = unpack(y)
    j = u*n    
    dn_dt = -diff(j)
    n_x = diff(n)
    u_x = diff(u)
    u_xx = diff(u, d=2)
    du_dt = -u*u_x
    du_dt -= get_V_m(t, x, d=1) + get_E_m(n, d=2)*n_x
    du_dt += get_nu_n(n, d=1)*n_x*u_x + get_nu_n(n) * u_xx / n
    
    # Correction needed for constant j = j_avg
    du = j_avg/n - u    
    du_dt += j_corr_factor * du
    return pack(dn_dt, du_dt)
```

```{code-cell} ipython3
k = 2*np.pi / Lx
s, c = np.sin(k*x), np.cos(k*x)
f = np.exp(s)
df = k*c*f
ddf = k**2*f*(c**2 - s)
assert np.allclose(diff(f), df)
assert np.allclose(diff(f, 2), ddf)
```

```{code-cell} ipython3
V0_m = 1.0
dV_m_dx = -0.05
sigma = 0.5
g_m = 1.0
n0 = 1.5
u0 = 0.3
j_avg = n0*u0
nu = 0.8

t = 0

b = get_V_m(t, x) / g_m
h = n0
eta = h - b

n = eta
u = 0*x + u0
y0 = pack(n, u)

dt = dx / nu * 0.002
j_corr_factor = 0

fig, axs = plt.subplots(2, 1, sharex=True)
ax = axs[0]
ax_c = ax.twinx()
ax_stream = axs[1]
ax_j = ax_stream.twinx()
y = y0.copy()
for n in FPS(80, fig=fig, display=True, embed=True, fps=20):
    for _ in range(500):
        dy_dt = compute_dy_dt(t, y)
        y += dy_dt * dt
        t += dt
    n, u = unpack(y)
    ax.cla()
    ax.plot(x, n)
    ax.set(ylim=(0, 2.0), title=f"{t=:.4g}, N={sum(n)*dx:.4g}",
           ylabel="$n$")
    
    # Speed of sound
    u_c = np.sqrt(n*get_E_m(n, d=2))
    ax_c.cla()
    ax_c.yaxis.set_label_position("right")
    ax_c.plot(x, u/u_c, '--C1')
    ax_c.axhline(1.0, color='k', ls=':')
    ax_c.set(ylim=(-0.1, 1.3), ylabel="Mach number $u/u_c$")

    ax_stream.cla()
    b = get_V_m(t, x) / g_m
    eta = n
    h = b + eta
    ax_stream.plot(x, b, '-k')
    ax_stream.fill_between(x, b, h, color='cyan')
    ax_stream.set(ylim=(0, 2.0), xlim=(x.min(), x.max()),
                  xlabel="$x$", ylabel="Shallow-water analog.")

    ax_j.cla()
    ax_j.plot(x, n*u)
    ax_j.plot(x, u)
    ax_j.yaxis.set_label_position("right")
    ax_j.set(ylim=(0, 1.0), 
             ylabel="$j$ (blue) and $u$ (orange)")
```

```{code-cell} ipython3
dt = dx / nu * 0.002
j_corr_factor = 0

fig, axs = plt.subplots(2, 1, sharex=True)
ax = axs[0]
ax_c = ax.twinx()
ax_stream = axs[1]
ax_j = ax_stream.twinx()
for n in FPS(80, fig=fig, display=True, embed=True, fps=20):
    for _ in range(500):
        dy_dt = compute_dy_dt(t, y)
        y += dy_dt * dt
        t += dt
    n, u = unpack(y)
    ax.cla()
    ax.plot(x, n)
    ax.set(ylim=(0, 2.0), title=f"{t=:.4g}, N={sum(n)*dx:.4g}",
           ylabel="$n$")
    
    # Speed of sound
    u_c = np.sqrt(n*get_E_m(n, d=2))
    ax_c.cla()
    ax_c.yaxis.set_label_position("right")
    ax_c.plot(x, u/u_c, '--C1')
    ax_c.axhline(1.0, color='k', ls=':')
    ax_c.set(ylim=(-0.1, 1.3), ylabel="Mach number $u/u_c$")

    ax_stream.cla()
    b = get_V_m(t, x) / g_m
    eta = n
    h = b + eta
    ax_stream.plot(x, b, '-k')
    ax_stream.fill_between(x, b, h, color='cyan')
    ax_stream.set(ylim=(0, 2.0), xlim=(x.min(), x.max()),
                  xlabel="$x$", ylabel="Shallow-water analog.")

    ax_j.cla()
    ax_j.plot(x, n*u)
    ax_j.plot(x, u)
    ax_j.yaxis.set_label_position("right")
    ax_j.set(ylim=(0, 1.0), 
             ylabel="$j$ (blue) and $u$ (orange)")
```

```{code-cell} ipython3
import sympy
x = sympy.symbols('x', real=True)
u_exact = sympy.exp((x**4-1)/4)

def deq(u):
    return (u.diff(x, x) - (x**6 + 3*x**2)*u)

assert deq(u_exact).simplify() == 0
assert u_exact.subs([(x, -1)]) == 1
assert u_exact.subs([(x, 1)]) == 1
```

```{code-cell} ipython3
N = 2
an = sympy.symbols(f'a:{N}', real=True)
kn = tuple(sympy.pi * (2*sympy.S(n) + 1)/2 for n in range(N))
u = 1 + sum(a_ * sympy.cos(k_ * x) for a_, k_ in zip(an, kn))
```

```{code-cell} ipython3
import numpy as np
from scipy.optimize import curve_fit
def model(x, *a):
    n_ = np.arange(len(a), float)[None, :]
    k_ = np.pi * (2*n_ - 1)/2
    a_ = np.array(a)[None, :]
    x_ = np.asarray(x)[:, None]
    return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)

xdata = np.linspace(-1, 1, 200)
ydata = sympy.lambdify(x, u_exact, "numpy")(xdata)
curve_fit(model, xdata, ydata, p0=[1,1])
```

```{code-cell} ipython3
import numpy as np
from scipy.optimize import curve_fit
from uncertainties import correlated_values, unumpy as unp
N = 2

def model(x, *a):
    n_ = np.arange(len(a), dtype=float)[None, :]
    k_ = np.pi * (2*n_ + 1)/2
    a_ = np.array(a)[None, :]
    x_ = np.asarray(x)[:, None]
    return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)

xdata = np.linspace(-1, 1, 200)
ydata = sympy.lambdify(x, u_exact, "numpy")(xdata)
p0, pcov = curve_fit(model, xdata, ydata, p0=np.ones(N))
p = correlated_values(p0, pcov)
print(p)
y_ = model(xdata, *p)
y, dy = unp.nominal_values(y_), unp.std_devs(y_)

f = 3
plt.plot(xdata*f, ydata, '-')
plt.plot(xdata*f, y, ':')
plt.fill_between(xdata*f, y-dy, y+dy, color='C1', alpha=0.5)
plt.grid('on')
plt.xticks([0, 1, 2, 3, 4, 5])
```

```{code-cell} ipython3
import numpy as np
from scipy.optimize import curve_fit

def model(x, *a):
    n_ = np.arange(len(a), dtype=float)[None, :]
    k_ = np.pi * (2*n_ + 1)/2
    a_ = np.array(a)[None, :]
    x_ = np.asarray(x)[:, None]
    return 1 + (a_ * np.cos(k_*x_)).sum(axis=-1)

exact = sympy.lambdify(x, u_exact, "numpy")

fig, ax = plt.subplots(figsize=(4,3), constrained_layout=True)
xdata = np.linspace(-1, 1, 200)
ydata = exact(xdata)
ax.plot(xdata, ydata, '-k')
for N in range(2, 5):
    p, pcov = curve_fit(model, xdata, ydata, p0=np.ones(N))
    ax.plot(xdata, model(xdata, *p), '--', label=f"{N=}")
fig.legend(loc='upper center');
```

```{code-cell} ipython3
plt.plot(xdata, model(xdata, *get_a(np.array([0, 0.999])).ravel()))
plt.plot(xdata, ydata, '-k')
```

```{code-cell} ipython3
import sympy
def get_spectral_approx(N, xi=None, show=False):
    an = sympy.symbols(f'a:{N}', real=True)
    kn = tuple(sympy.pi * (2*n + 1)/2 for n in sympy.S(np.arange(N)))
    if xi is None:
        xi = tuple((2*n+1)/(2*N+1) for n in sympy.S(np.arange(N)))
        xi = tuple((2*n+1)/(2*N+1) for n in sympy.S(np.arange(N)))
    #print(xi)
    #xi = tuple(sympy.cos(sympy.pi * (N-n)/2/N) for n in sympy.S(np.arange(N)))
    x = sympy.symbols('x')
    u = 1 + sum(a_ * sympy.cos(k_ * x) for a_, k_ in zip(an, kn))
    Eqs = [(u.diff(x,x) - (x**6 + 3*x**2)*u).subs([(x, x_)]).evalf() for x_ in xi]
    sol = sympy.solve(Eqs, an)
    u_x = u.subs(sol)
    if show:
        display(sol)
        display(u.evalf())
    u_ = sympy.lambdify(x, u_x, "numpy")
    return u_

fig, ax = plt.subplots()
x = np.linspace(-1, 1)
u_exact = np.exp((x**4-1)/4)
ax.plot(x, u_exact, label="Exact: $u(x) = e^{(x^4-1)/4}$")
for N in [2, 3, 4, 5, 6]:
    u = get_spectral_approx(N=N)
    ax.plot(x, u(x), label=f"Spectral: ${N=}$")
u = get_spectral_approx(N=2, xi=[0.1,0.2])
ax.plot(x, u(x), label=f"Spectral: ${N=}$")
fig.legend()
```

```{code-cell} ipython3
from scipy.optimize import minimize

res = minimize(get_err, [0.3, 0.75])
res = minimize(get_err2, [0.3, 0.75])
u = get_spectral_approx(N=2, xi=res.x)
X = np.linspace(-1, 1, 200)
plt.plot(X, u(X))
plt.plot(X, exact(X))
```

```{code-cell} ipython3
from tqdm.auto import tqdm
from itertools import product

def get_a(x, exact=exact):
    n_ = np.arange(len(x), dtype=float)[None, :]
    k_ = np.pi * (2*n_ + 1)/2
    x_ = np.asarray(x)[:, None]
    A = np.cos(k_*x_)
    a = np.linalg.solve(A, exact(x) - 1)
    return a
    
def get_err(x, weight=False):
    x = np.asarray(x)
    a = get_a(x)
    X = np.linspace(-1, 1, 201)[1:-1]
    w = 1/np.sqrt(1-X**2) if weight else 1
    return (w*(model(X, *a) - exact(X))**2).mean()

def get_err2(x, weight=False):
    u = get_spectral_approx(N=2, xi=x)
    X = np.linspace(-1, 1, 201)[1:-1]
    w = 1/np.sqrt(1-X**2) if weight else 1
    return (w*(u(X) - exact(X))**2).mean()

fig, axs = plt.subplots(1, 2, figsize=(9, 3))
x0, x1 = np.meshgrid(np.linspace(0.001, 0.99), np.linspace(0.002, 0.999), indexing='ij')
err = [get_err([x0, x1]) for (x0, x1) in tqdm(list(zip(x0.ravel(), x1.ravel())))]
err = np.reshape(err, x0.shape)
mesh = axs[0].pcolormesh(x0, x1, 1/err.T)
fig.colorbar(mesh, ax=axs[0])

x0, x1 = np.meshgrid(np.linspace(0.001, 0.99, 30), np.linspace(0.002, 0.999, 30), indexing='ij')
err = [get_err2([x0, x1]) for (x0, x1) in tqdm(list(zip(x0.ravel(), x1.ravel())))]
err = np.reshape(err, x0.shape)
mesh = axs[1].pcolormesh(x0, x1, 1/err.T)
fig.colorbar(mesh, ax=axs[1])
```

```{code-cell} ipython3
def model(x, *a):
    n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
    a_ = np.array(a)[None, :]
    x_ = np.asarray(x)[:, None]
    return 1 + (a_ * (np.cos(2*n_*np.arccos(x_))-1)).sum(axis=-1)

def get_xi(N):
    """Return the N collocation points."""
    i = np.arange(N)
    xi = np.cos((i + 0.5) * np.pi / 2 / N)
    return xi
    
def get_a(N, exact=exact):
    x_ = get_xi(N)[:, None]
    n_ = np.arange(1, N+1, dtype=float)[None, :]
    k_ = 2*n_
    A = np.cos(k_*np.arccos(x_)) - 1
    a = np.linalg.solve(A, exact(x_.ravel()) - 1)
    return a

xi = get_xi(2)
a = get_a(2)
plt.plot(xdata, model(xdata, *a))
plt.plot(xdata, exact(xdata))
list(map(plt.axvline, xi))
```

```{code-cell} ipython3
x_ = get_xi(N)[:, None]
n_ = np.arange(1, N+1, dtype=float)[None, :]
k_ = 2*n_
A = np.cos(k_*np.arccos(x_)) - 1
a = np.linalg.solve(A, exact(x_.ravel()) - 1)
A @ a, exact(x_.ravel()) - 1
def model(x, *a):
    n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
    a_ = np.array(a)[None, :]
    x_ = np.asarray(x)[:, None]
    return 1 + (a_ * (np.cos(2*n_*np.arccos(x_))-1)).sum(axis=-1)
model(x_, *a)
```

```{code-cell} ipython3
model(x_, *a)
n_ = np.arange(1, len(a)+1, dtype=float)[None, :]
a_ = np.array(a)[None, :]
(a_ * (np.cos(k_*np.arccos(x_))-1)).sum(axis=-1), exact(x_) - 1
A @ np.ravel(a_),  (a_ * (np.cos(k_*np.arccos(x_))-1)).sum(axis=-1), exact(x_) - 1
```

```{code-cell} ipython3
def phi(x, n, d=0):
    """Return the dth derivative of phi_n(x)"""
    th = np.arccos(x)
    cn = np.cos(2*n*th)

    if d == 0:
        return cn - 1
    
    sn = np.sin(2*n*th)
    s = np.sin(th)

    if d == 1:
        return (2*n) * sn / s

    c = np.cos(th)
    if d == 2:
        return -((2*n)**2 * cn * s - (2*n) * sn * c) / s**3

x0 = 0.3
h = 1e-8
n = 3
assert np.allclose((phi(x0+h, n) - phi(x0-h, n))/2/h, phi(x0, n, d=1))
assert np.allclose((phi(x0+h, n, d=1) - phi(x0-h, n, d=1))/2/h, phi(x0, n, d=2))

def get_a(N):
    i = np.arange(N)
    xi = np.cos((i + 0.5) * np.pi / 2 / N)
    phis = np.array([phi(xi, n=n) for n in 1+np.arange(N)])
    d2phis = np.array([phi(xi, n=n, d=2) for n in 1+np.arange(N)])
    rhs = xi**6 + 3*xi**2
    A = np.transpose(d2phis - rhs[None, :]*phis)
    print(A.shape)
    a = np.linalg.solve(A, rhs)
    return a
    
plt.plot(xdata, ydata, ':')
plt.plot(xdata, model(xdata, *get_a(6)))
```

# Stiff Example

https://www.mathworks.com/company/technical-articles/stiff-differential-equations.html

```{code-cell} ipython3
%pylab inline
from ipywidgets import interact
from scipy.special import lambertw

from scipy.integrate import solve_ivp

def F(t, y):
    return y**2*(1-y)

def exact(t, delta):
    a = 1/delta - 1
    return 1/(lambertw(a*np.exp(a-t)) + 1)
    

@interact(method=['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA'],
          log_delta=(-6, -1))
def go(log_delta=-2.0, method="RK45"):
    delta = 10**(log_delta)
    res = %time solve_ivp(F, y0=[delta], t_span=[0, 2/delta], atol=1e-12, rtol=1.e-12, method=method)
    plt.plot(res.t*delta, exact(res.t, delta=delta))
    plt.plot(res.t*delta, res.y.T)
    plt.title(f"{res.nfev=}")
    display(plt.gcf())
    plt.close('all')
```

```{code-cell} ipython3
solve_ivp??
```

```{code-cell} ipython3

```