---
execution:
  timeout: 60
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.16.6
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
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

(sec:EffectiveViscosity)=
# Effective Viscosity

Here we consider viscous hydrodynamics in 1D of a fluid with number density $n(x, t)$,
velocity $u(x, t)$, and energy density $\mathcal{E}(n)$.  Thermodynamics gives the
"pressure" $P(n)$ (just a force in 1D) as follows:

\begin{gather*}
  \mu(n) = \mathcal{E}'(n),\\
  P(n) = \mu n - \mathcal{E} 
  = n\mathcal{E}'(n) - \mathcal{E}(n)
  = n^2
  \left(\frac{\mathcal{E}}{n}\right)',\\
  P'(n) = n\mathcal{E}''(n) = n \mu'(n).
\end{gather*}
*Here $ϵ(n) = \mathcal{E}/n = E/N$ is the energy per particle, which sometimes appears.
We will not use this here.*

Conservation of particle number and momentum give the hydrodynamic equations:
\begin{gather*}
  \dot{n} + (nu)_{,x} = 0\\
  D_tu = \dot{u} + uu_{,x} = -\frac{1}{m}\Bigl(
    V(x) + \overbrace{\mathcal{E}'(n)}^{\mu}
  \Bigr)_{,x} + \frac{\bigl(ν(n) u_{,x}\bigr)_{,x}}{n},\\
  mnD_tu = mn(\dot{u} + uu_{,x}) = \overbrace{-
    n\Bigl(
      V(x) + \mathcal{E}'(n)\Bigr)_{,x}}^{-\mathcal{P}_{,x}}
    + m\bigl(ν(n) u_{,x}\bigl)_{,x},
\end{gather*}
where $ν$ is a viscosity coefficient, and $V(x, t)$ is an external potential.

:::::{admonition} A Comment on the Form of these Equations

In the applied mathematics literature, these equations are often expressed in different
forms.  Here we reconcile these.  A standard form is:
\begin{align*}
  \dot{\rho} + (\rho u)_{,x} &= 0\\
  (\rho u)_{,t} + \bigl(\rho u^2 + P(\rho) \bigr)_{,x} 
  &= \varepsilon\bigl(ν(\rho) u_{,x}\bigr)_{,x}.
\end{align*}
Note that this is equivalent to our second form above after using the continuity
equation:
\begin{gather*}
  (\rho u)_{,t} + (\rho u^2)_{,x} = 
  \underbrace{\dot{\rho} u}_{\mathclap{-u(\rho u)_{,x}}} 
  + \rho \dot{u} + \rho_{,x} u^2 + 2\rho u u_{,x}
  = \underbrace{\rho(\dot{u} + u u_{,x})}_{mn D_{t} u}
\end{gather*}
These equations here are expressed in terms of the mass-density $\rho = m n$ and
velocity $u$, and are usually simplified with power-law forms for $P(\rho) =
\rho^\gamma$ with $\gamma > 1$ and $ν(\rho) = \rho^{\alpha}$ with $\alpha > 0$. The
latter is sometimes called a **Lamé viscosity**.

With $\alpha = 1$ and $\gamma = 2$, these are equivalent to the Saint-Venant system when
the pressure is quadratic $P(\rho) = \tfrac{1}{2}\kappa \rho^2$ which is the case with
the standard GPE with $m^2 \kappa = g$:
\begin{gather*}
  \mathcal{E}_{\text{GPE}}(n) = \frac{g}{2} n^2, \qquad
   P_{\text{GPE}}(n) = n\mathcal{E}'_{\text{GPE}}(n) - \mathcal{E}_{\text{GPE}}(n)
                     = \frac{g}{2}n^2.
\end{gather*}

Some references:

* {cite}`Strani:2022` uses these with the variable $w = u$ and looks at stability
  properties on bounded intervals.
* {cite}`Joseph:2011` uses the same form $\alpha=1$ with $\gamma = 2/5$ to describe
  shockwaves in an elongated unitary Fermi gas, which has equation of state:
  \begin{gather*}
    u_{,t} = -\left(\frac{u^2}{2} + Cn^{2/5} + V(x)\right)_{,x} 
           + ν \frac{\bigl(n u_{,x}\bigr)_{,x}}{n} \tag{6}\\
    n(u_{,t} + u u_{,x}) = - n\Bigl(\underbrace{Cn^{2/5}}_{\mathcal{E}'(n)} + V(x)\Bigr)_{,x}
    + ν \bigl(n u_{,x}\bigr)_{,x}.
  \end{gather*}
  Note: the equation of state for the UFG is
  \begin{gather*}
    \mathcal{E}_{\mathrm{UFG}}(n) \propto n^{5/3}
  \end{gather*}
  where $[n] = 1/L^3$ is the 3D density.  They argue via the Thomas-Fermi approximation
  that
  \begin{gather*}
    n_{3D}(r, z) = \tilde{n}\left(1 - \frac{r^2}{R_{\perp}^2} - \frac{V(z)}{\mu}\right)^{3/2}
                 = \frac{\tilde{n}}{R_\perp^{3}}\Bigl(r^2 - R^2(z)\Bigr)^{3/2}, \\
    R(z) = R_{\perp}\sqrt{1-\frac{V(z)}{\mu}},\\
    n_{1D}(z) %= \int_0^{R(z)}\!\!\!\!\!\!\!\!\!\!\d{r}\; 2\pi r n_{3D}(r, z)
              = \pi \frac{\tilde{n}}{R_\perp^3}
              \int_0^{R^2(z)}\!\!\!\!\!\!\!\!\!\!\d{r^2}\;
              \left(r^2 - R^2(z)\right)^{3/2}
              = \frac{2\pi \tilde{n}}{5 R_\perp^3}R^{5}(z)\\
              = \frac{2\pi \tilde{n}}{5 R_\perp^3}
              \left(1 - \frac{V(z)}{\mu}\right)^{5/2}.
\end{gather*}
Hence, $\mu  = \mathcal{E}'(n_{1D}) \propto n_{1D}^{2/5}$.  However, see also
{ref}`sec:NPSEQ`. 
:::::

:::::{admonition} A Comment on the Pressure
As discussed below, the momentum current is best expressed in terms of a pressure
$\mathcal{P}$.  To do this, we are looking for a function $\mathcal{P}$ whose gradient
is $n$ times the gradient of the potential
\begin{gather*}
    \mathcal{P}_{,x} = n \Bigl(V+\mathcal{E}'(n)\Bigr)_{,x},\\
    \mathcal{P} = F(x) + P(n) = F(x) + n\mathcal{E}'(n) - \mathcal{E}(n),\\
    P(n) = n\mathcal{E}'(n) - \mathcal{E}(n), \qquad
    F(x) = \int_{x_0}^{x}n(x)V'(x)\d{x}.
\end{gather*}
This contains both the intrinsic pressure $P(n)$ from the equation of state of the
fluid, and a contribution from the external potential $F(x)$.
The force $F(x)$ here is the integrated force on the potential such that
$F(\infty)$ would be the force exerted by the fluid on the potential.  I.e., if this
potential were created by you sticking an object into the flow, then this would be the
total force on the object.  This is effectively a **buoyant force** but under
appropriate boundary conditions -- i.e. no potential drop -- represents the **drag
force**.  This should not be confused with the local force $-V_{,x}$ of the potential on
the fluid.

It is a peculiarity of 1D that force and pressure have the same dimensions -- more
generally, one would obtain similar relationships after integrating a section.

:::{admonition} Do It! Show that $(n\mathcal{E}'(n) - \mathcal{E})_{,x} = \bigl(\mathcal{E}'(n)\bigr)_{,x}$.
:class: dropdown

Just apply the chain rule and cancel terms
\begin{gather*}
  \Bigl(n\mathcal{E}'(n) - \mathcal{E}(n)\Bigr)_{,x} = 
  n_{,x}\mathcal{E}'(n) + n\mathcal{E}''(n) n_{,x} - \mathcal{E}'(n)n_{,x} = 
  n\mathcal{E}''(n) n_{,x}
  =
  n\Bigl(\mathcal{E}'(n)\Bigr){,x}.
\end{gather*}
:::
Unfortunately, since $n(x)$ is a property of the solution, there is no simple way of
expressing the force, so we cannot use this as a method of solution.
:::::

:::{admonition} Interface with Vacuum

It is known (see e.g. {cite}`Lian:2010` and references therein) that these equations,
with a constant viscosity $ν(n)$ leads to singular behavior at the interface with the
vacuum (zero density).  This can be alleviated by allowing $ν(n)$ to depend on density.
We do not consider the vacuum interface here, assuming that $n>0$ everywhere in the
problem domain.
:::

:::{margin}
The term **shallow water** here refers to the approximation that treat each column of
water as an individual "compressible" element that moves horizontally.  The shallow
approximation neglects the complicated rolling motion that can develop in deeper water.
Instead, we separate the vertical and horizontal motion such that the vertical motion is
expressed solely in terms of the height of the column (non-inertial).
:::
:::::{admonition} A Physical Analogy: Shallow Water

To develop intuition, it is useful to have a familiar physical analogy.  Consider an
incompressible fluid of mass density $\rho$ like water flowing along a channel of width $w$
whose bottom height $b(x)$ varies, but not enough to excite appreciable vertical dynamics.  If
the fluid height is $\eta(x)$, then the internal energy density is given by the
gravitational potential energy (here $a_g$ is the acceleration due to gravity)
\begin{gather*}
  \int_{b(x)}^{b(x) + \eta(x)} a_g \rho w z\;\d{z}
              = a_g w\rho \frac{(b+\eta)^2 - b^2}{2}
              = a_g w\rho \frac{\eta^2}{2} + a_g w b \eta
\end{gather*}
This suggests identifying the following 1D number density and coupling constant to match the GPE.
\begin{gather*}
    n_{1}(x) \equiv \frac{\rho w \eta}{m}, \qquad
    g \equiv \frac{m^2 a_g}{w \rho^2}, \qquad
    V(x) \equiv \frac{m a_g b(x)}{\rho}.
\end{gather*}

:::{admonition} Comparison with [Wikipedia][Shallow water equations]
:class: dropdown

To compare with the [conservative form of the shallow water equations][], note that, in
2D, the following are equivalent, after application of the continuity equation:
\begin{gather*}
   \rho \eta D_t u = \pdiff{(\rho \eta u)}{t} 
   + \pdiff{\rho\eta u^2}{x} + \pdiff{\rho\eta u v}{y}.
\end{gather*}
Their form does not include the potential $V(x)$ (they consider a flat bed).  Explicitly
\begin{gather*}
  \dot{n} = - \sum_{i}\nabla_{i}(nv_{i}), \qquad
  D_t v_i = \dot{v}_{i} + \sum_{j} v_j \nabla_{i}v_i,\\
  \pdiff{n v_i}{t} + \sum_{j}\nabla_j(nv_iv_j) 
  = \underbrace{n \dot{v}_i + n\sum_{j}v_j\nabla_j(v_i)}_{nD_tv_i} 
  + \underbrace{\dot{n}v_i + v_i\sum_{j}\nabla_j(nv_j)}_{0}.
\end{gather*}
:::
:::::

[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>
[weir]: https://en.wikipedia.org/wiki/Weir
[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>

## Conservation Laws

Local conservation laws can be expressed in terms of a **density** $a$ and a
**current** $c$ such that, under strict conservation, the integrated quantity $A$ does
not change with time:
\begin{gather*}
  A(t) = \int a(x, t)\d{x}, \qquad
  \dot{A}(t) = 0.
\end{gather*}
The local expression of this conservation is
\begin{gather*}
  \dot{a}(x, t) + c_{,x}(x, t) = 0.
\end{gather*}
I.e. if the current vanishes far away, or if we work in a periodic box
\begin{gather*}
  \dot{A}(t) = \int \dot{a}(x, t)\d{x}
             = -\int c_{,x}(x, t)\d{x}
             = -\left.c(x, t)\right|^{\infty}_{-\infty} = 0.
\end{gather*}

:::{margin}
\begin{gather*}
  a = n, \qquad
  c = j = un.
\end{gather*}
:::
### 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.
:::
:::::

:::{margin}
\begin{gather*}
  a = nu, \qquad
  c = nu^2 + \frac{\mathcal{P}}{m} - ν(n) u_{,x}.
\end{gather*}
:::
### 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} 
    - (ν(n) u_{,x})_{,x} = 0\\
  (nu)_{,t} + \Biggl(
    \underbrace{nu^2 + \frac{\mathcal{P}}{m} - ν(n) u_{,x}}_{p}
  \Biggr)_{,x} = 0.
\end{gather*}

:::{margin}
\begin{gather*}
  a = \varepsilon = \frac{mnu^2}{2} + nV + \mathcal{E}(n),\\
  c = j\left(\frac{mu^2}{2} + V + n\mathcal{E}'(n)\right).
\end{gather*}
:::
### Energy Density

We can similarly consider the energy flux.  The energy density is the sum of the
kinetic, potential, and internal energies.
:::{margin}
These have only been check for constant $ν(n)$ but are probably true.
:::
\begin{align*}
  \varepsilon = \frac{m}{2}nu^2 + nV + \mathcal{E}(n),&\\
  \dot{\varepsilon} + \Bigl[
    \underbrace{un}_{j}\bigl(\tfrac{m}{2}u^2 + V + \mathcal{E}'(n)\bigr)
  \Bigr]_{,x} &= n \dot{V} + mu\bigl(ν(n) u_{,x}\bigr)_{,x},\tag{checked}\\
  \dot{\varepsilon} + \Bigl[
    j\bigl(\tfrac{m}{2}u^2 + V + \mathcal{E}'(n)\bigr) - mν(n) uu_{,x}
  \Bigr]_{,x} &= n \dot{V} - mν(n) u_{,x}^2.\tag{checked}
\end{align*}
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 shows explicitly
that positive $ν(n)$ will dissipate energy if $\dot{V} = 0$ and also gives an energy
density whose gradient is strictly decreasing for static solutions, which will be useful
in our analysis below.

:::::{admonition} Numerical Checks

To make sure what we are doing makes sense, we give a concrete numerical example.  We
formulate the PDE as an ODE in time:
\begin{align*}
  \dot{n} &= - (nu)_{,x}\\
  \dot{u} &= - uu_{,x} - \frac{V'(x) + \mathcal{E}''(n)n_{,x}}{m} 
             + \frac{\bigl(\nu(n) u_{,x}\bigr)_{,x}}{n}.
\end{align*}
For details, see {ref}`sec:Viscosity1`
:::::




:::{margin}
Here we use the static momentum and energy conservation, and
\begin{gather*}
  F_{,x} = -nV_{,x},\\
  P(n) = n \mathcal{E}'(n) - \mathcal{E}(n),\\
  \mathcal{P}_{,x} = n \Bigl(V+\mathcal{E}'(n)\Bigr)_{,x}.
\end{gather*}
:::
## Some Analytic Results for Static Flow

For static problems, the current and momentum currents are constants everywhere:
\begin{gather*}
  j = nu, \qquad
  p = nu^2 + \frac{\mathcal{P}}{m} - ν(n) u_{,x}.
\end{gather*}

### Rankine-Hugoniot Conditions

We start by considering solutions that are stationary far to the left $j_L=u_Ln_L$ and
far to the right $j_R = u_Rn_R$ such that $j_{,t} = (nu)_{,t} = 0$ and $u_{,x} = 0$ in
these asymptotic regions.  The momentum conservation equation thus implies that $mnu^2
+\mathcal{P}$ is constant in these regions, but the constant might differ on the left
and right of some non-stationary region.  Continuity implies that, if $j_L \neq j_R$,
then the transition region must move with velocity $v$ to accommodate the change in
particle number.  This is one of the **[Rankine-Hugoniot conditions][]**:

:::{margin}
  To check the sign, consider $n_L = j_L = 0$.  The trailing edge must have positive
  velocity $v>0$.
:::
\begin{gather*}
  v(n_R - n_L) = j_R - j_L.
\end{gather*}
The usual application of this is to consider the motion of a shockwave or soliton which
maintains its shape.  Boosting to the co-moving frame, the problem becomes static
everywhere and, in this frame, $j_L = j_R = j$ is constant everywhere.  Additionally,
the momentum density is constant everywhere:
\begin{gather*}
   mnu^2 + \mathcal{P} - mν(n) u_{,x} 
   = nu^2 + F(x) + P(n) - mν(n)u_{,x}
   = \text{const.},\\
  \mathcal{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*}

Considering again regions far from the potential where the flow is homogeneous, this
implies that the total force on whatever potential exists in the intermediate region is
\begin{gather*}
   F_{\text{tot}} = F_R - F_L 
   = \left.\left(\frac{m j^2}{n} + P(n)\right)\right|^{n_L}_{n=n_R}.
\end{gather*}




To be explicit, consider static flow: $\dot{n} = \dot{u} = 0$ starting with constant
density $n_L$ and velocity $u_L$ on the left with no potential $V(x) = 0$:
\begin{gather*}
  nu = j = n_Lu_L = \text{const.},\\
  nu^2 + \frac{F + P(n)}{m} - ν u_{,x} = n_Lu_L^2 + \frac{P(n_L)}{m} = \text{const.},\\
  %\left[u\Bigl(\frac{m}{2}n u^2 + nV + n\mathcal{E}'(n) - m ν
  %v_{,x}\Bigr)\right]_{,x} + m ν u_{,x}^2 = 0,\\
  %\Bigl(\frac{m}{2} u^2 + V + \mathcal{E}'(n)\Bigr)_{,x} = \frac{m ν}{j} (u u_{,xx}),\\
  %nV_{,x} + m ju_{,x} + nn_{,x}\mathcal{E}''(n) = m ν u_{,xx},\\
\end{gather*}
From this, we can solve for the force:
\begin{gather*}
  F(x) = m j^2\left(\frac{1}{n_L} - \frac{1}{n(x)}\right) + P(n_L) - P\bigl(n(x)\bigr)
  + mν u_{,x}.
  \tag{checked}
\end{gather*}
If the flow far to the right or the potential is also constant $u_{,x} \rightarrow 0$,
then the total force on the potential is
\begin{gather*}
  F(\infty) = F_{\text{total}} = \left.\left(\frac{m j^2}{n} + P(n)\right)\right|^{n_L}_{n=n_R}.
\end{gather*}
If there is no potential $V(x) = 0$, then $F_{\text{total}} = 0$ and we find the
**Michelson–Rayleigh line**:
\begin{gather*}
  \frac{P(n_R) - P(n_L)}{1/n_R - 1/n_L} = -m j^2.
\end{gather*}
*(The final condition comes from energy conservation, but this is only independent if we
consider finite temperature effects.)*


:::::{admonition} Numerical Checks

To make sure what we are doing makes sense, we give a concrete example.  To minimize the
chance for error, we first implement the original equations in the form $y=(n, ν(n)u,
\bigl(ν(n)u\bigr)', F)$, looking for stationary solutions $\dot{n} = \dot{u} = 0$.  For
an initial condition, we specify the density and velocity to the left of the barrier and
directly integrate.
\begin{gather*}
  (nu)_{,x} = 0, \qquad
  muu_{,x} = -(V+\mathcal{E}')_{,x} + m\bigl(ν(n) u_{,x}\bigr)_{,x}/n\\
  n_{,x} = -\frac{nu_{,x}}{u}, \qquad
  u_{,xx} = n\frac{uu_{,x} + \frac{(V+\mathcal{E}')_{,x}}{m}}{ν}
          % = \frac{mnuu_{,x} + n(V'(x) + \mathcal{E}''(n)n_{,x})}{mν(n)}
  \\
  \begin{pmatrix}
    n\\
    ν(n)u\\
    \bigl(ν(n)u\bigr)_{,x}\\
    F
  \end{pmatrix}_{,x}
  =
  \begin{pmatrix}
    -\frac{nu_{,x}}{u}\\
    \bigl(ν(n)u\bigr)_{,x}\\
    n\Bigl(uu_{,x} + \frac{(V+\mathcal{E}')_{,x}}{m}\Bigr)\\
    nV_{,x}
  \end{pmatrix}, \qquad
  y(0) = \begin{pmatrix}
    n_0\\
    ν(n_0)u_0\\
    0\\
    0
  \end{pmatrix}.
\end{gather*}
:::::

```{code-cell}
:tags: [hide-input]
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
    
    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='LSODA', 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)
        assert(res.success)
        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}
:tags: [margin, hide-input]
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')

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.453).plot()
```
:::{margin}
Low velocity flow (top) and high-velocity flow (bottom).  Notice that at low flow
velocities, the $n(x)\approx n_{TF}(x)$.  This approximation breaks down as the flow
velocity increases to the point that flow at the top of the barrier is almost
supersonic.  In this limit, the density (and velocity) develops a kink.  In the last
frame, we increase the viscosity, allowing for hypersonic flow.
:::
(sec:Viscosity_NoViscosity)=
### No Viscosity.
In the limit of zero viscosity, particle number and energy conservation gives immediately
\begin{gather*}
  \frac{mu^2}{2} + V(x) + \mathcal{E}'(n) 
  = \underbrace{
    \overbrace{K_0}^{mu_0^2/2} + \overbrace{\mu_0}^{\mathcal{E}'(n_0)}}_{E_0} 
  = \text{const.},\qquad
  V(x) = E_0 - \frac{m j^2}{2n^2} - \mathcal{E}'(n).
\end{gather*}
This directly relates the density to the potential and for $j=0$, corresponds to the
Thomas-Fermi (TF) approximation.   For example, with the GPE
$\mathcal{E}'(n) = gn$, so
\begin{multline*}
  n(x) = (\mathcal{E}')^{-1}\left(
    \mu_0 - V(x) + K_0 - \frac{m j^2}{2n(x)^2}
  \right)\\
  =
  \underbrace{\frac{\mu_0 - V(x)}{g}}_{n_{TF}(x)} 
  - \frac{m j^2}{2g}\left(\frac{1}{n^2(x)} - \frac{1}{n_0^2}\right).
\end{multline*}
The second term corrects the TF density, and becomes significant at high rates
of flow, or low densities.  While the TF approximation gives some insight, the
non-linear nature of the equation of state causes this to break down in an interesting
way.  This can be deduced from this expression: the correction term reduces the density,
but this increases the significance of the correction.  The feedback causes things to
break earlier than expected.

```{code-cell}
:tags: [margin, hide-input]
n = np.linspace(0.01, 5)
fig, ax = plt.subplots(figsize=(4, 3))
ax.plot(n, -1/n**2-n)
ax.set(ylim=(-6, 0), xlabel="$n$", ylabel="$-1/n^2 - n$");
```
To see this, we consider this from a different perspective: as a limit on the maximum
magnitude of the potential $V(x)$ that will allow flow:
\begin{gather*}
  V(x) \leq E_0 - \frac{m j^2}{2n^2} - \mathcal{E}'(n).
\end{gather*}
Maximizing the right-hand side for fixed $E_0$ and $j$, we find a critical density $n_c$
\begin{gather*}
   \underbrace{\frac{mj^2}{n_c^3}}_{mu_c^2/n_c} = \mathcal{E}''\bigl(n_c\bigr),\qquad
   u_c = c_s = \sqrt{\frac{n_c\mathcal{E}''\bigl(n_c\bigr)}{m}},\\
   V(x) \leq V_c = E - \frac{m j^2}{2n_c^2} - \mathcal{E}'(n_c).
\end{gather*}
The velocity $u_c$ associated with this critical density is the speed of sound $c_s$ in
the fluid at density $n_c$.  The interpretation of this was not obvious (at least to me)
but very interesting (see the discussion below).  For slow flow, $u \ll u_c$, the TF
approximation is valid, and if the potential is smooth, then the density will smoothly
follow this $n(x) \approx [\mu_0 - V(x)]/g$ in the GPE.  For a barrier, the density
always decreases, meaning that the velocity $u = j/n$ increases.  If one increases the
current or increases the barrier height, however, then the fluid velocity over the
barrier increases until it reaches the speed of sound $u_c=c_s$.  When this happens,
$n(x)$ -- and consequently $u(x)$ -- develops a kink as the to independent solutions to
the non-linear equation merge and vanish.

:::{important}
We conclude that the static hydrodynamic equations cannot support hypersonic solutions
-- and if there is any viscosity, no matter how small, this cannot be neglected at the
kink where $u_{,xx}$ diverges!

This implies that the local velocity $u$ can never exceed the local speed of sound
implying a **global** critical density $n_c$ which is related to the conserved current
$j$ by the non-linear equation
\begin{gather*}
   j = n_c\sqrt{\frac{n_c\mathcal{E}''\bigl(n_c\bigr)}{m}}.
\end{gather*}
This in turn can be used to find the maximum allowable height of the potential
supporting a given rate of flow (which depends on both $j$ and $n_0$).  Finite
viscosity will affect these bounds in a non-trivial way.
:::

What about the other direction?  Suppose we have a barrier with maximum height
$V_{\max}$ and an initial density $n_0$: what is the maximum flow $j$ that can be
supported?
\begin{gather*}
  \frac{m j^2}{2} \leq \frac{V_{\max} + \mathcal{E}'(n_c) - \mathcal{E}'(n_0)}
                            {\frac{1}{n_0^2} - \frac{1}{n_c^2}}.
\end{gather*}
Here, $n_c$ depends on $j$, giving a non-linear condition.



:::::{admonition} Nonlinearities.
:class: dropdown

When developing this, I ran into issues because of the nonlinear nature of the
equations.  I expected that the TF approximation would break down when
\begin{gather*}
  n_{\min} = \frac{\mu_0 - V_{\max}}{g}
\end{gather*}
was such that the second term started to become significant.  The equations broke down,
however, well before this, developing a kink with a sharp drop in density for
significantly lower $V_0$ than expected.  The reason turns out to be related to the
non-linear nature of the solution.

Starting with the GPE, we have a cubic:
\begin{gather*}
  P_n(n) = gn^3 - (E_0 - V)n^2 + \frac{m j^2}{2} = 0.
\end{gather*}
The problem is that the relevant solution becomes degenerate (and then vanishes) when
$P_n'(n) = 0$ at the root
\begin{gather*}
  n_c = \frac{2(E_0 - V)}{3g}, \\
  P_n(n_c) = -\frac{g}{2}n_c^3 + \frac{m j^2}{2} = 0,\\
  n_c = \left(\frac{m j^2}{g}\right)^{1/3}.
\end{gather*}
This gives a constraint on $V$ to prevent the development of a kink in the density:
\begin{gather*}
  V(x) < E_0 - \frac{3g}{2}\left(\frac{m j^2}{g}\right)^{1/3}
  = E_0 - \frac{3}{2}(m g^2 j^2)^{1/3}.
  \tag{checked}
\end{gather*}

Generalizing, we consider the function
\begin{gather*}
  P(n) = n^2\mathcal{E}'(n) - (E_0 - V)n^2 + \frac{m j^2}{2}
\end{gather*}
and solve for $P(n_c) = P'(n_c) = 0$:
\begin{gather*}
  2\mathcal{E}'(n_c) + n_c\mathcal{E}''(n_c) = 2(E_0 - V)\\
  \mathcal{E}''(n_c) n_c^3 = m j^2 = m n_c^2 u_c^2.
\end{gather*}
This defines a critical flow velocity $u_c$ which is the speed of sound $c_s$ in the fluid:
\begin{gather*}
  u_c = c_s = \sqrt{\frac{\mathcal{E}''(n_c) n_c}{m}}.
\end{gather*}
:::::


(sec:BarrierFlow)=
### Flow Past a Barrier

:::{margin}
![Picture of the Calgary Weir](https://www.harviepassage.ca/wp-content/uploads/weir.jpg)
:::
It is interesting to consider flow past a weak barrier.  In the shallow-water
hydrodynamic picture, this is what happens when a river flows past a [weir][].  When
water flows from a slow-moving river over a weir, a thin layer of water moves rapidly
over the weir as required by continuity.  From this intuition, we might expect the
density $n$ over the barrier to decrease as the velocity $u$ increases keeping $j=un$
constant.

:::{margin}
It might be interesting to explore this with time-dependence.  Presumably this will
allow hypersonic flow, but producing excitations.  I am not sure if something in the
model breaks down in this case.
:::
This intuition is in contrast with the picture of free particles, which **lose** speed
as they move up and over the potential, gaining potential energy as they lose kinetic
energy.  In this case, continuity requires that the density increase, in contrast with
the typical behavior of a weir.  This conflict suggests a transition between a
hydrodynamic regime where interactions dominate (weir) and a ballistic or kinetic regime
where kinetic energy dominates, however, as discussed above, this ballistic regime
cannot be realized in the static problem as it requires hypersonic flow.

We thus limit our analysis here to subsonic flow, asking, what force does a fluid with
viscosity exert as it flows over a barrier.  The density and momentum continuity
equations give the following in terms of the constants $j$ and $p$:
\begin{gather*}
  j = nu, \qquad
  mp = mnu^2 + \mathcal{P} - mν u_{,x}\\
  =  mnu^2 + F(x) + n \mathcal{E}'(n) - \mathcal{E}(n) - mν u_{,x}
  = m \frac{j^2}{n_L} + P(n_L).
\end{gather*}
We assume that, far from the potential, $u_{,x} = 0$ so we can drop the viscosity term
in this equation.  The net force due to the potential is thus
\begin{gather*}
  F = m j(u_L - u_R) + P(n_L) - P(n_R) 
    = mj^2 \left(\frac{1}{n_L} - \frac{1}{n_R}\right) + P(n_L) - P(n_R).
\end{gather*}
:::{admonition} GPE
Specifically, for the GPE, we have
\begin{gather*}
  F = \left(\frac{g}{2} - \frac{mj^2}{n_Ln_R(n_L + n_R)}\right)(n_L^2 - n_R^2)\\
    = \frac{gn_L^2}{2}\left(1 - \frac{2mj^2/gn_L^3}{n_R/n_L(1 + n_R/n_L)}\right)
                      \left(1 - \frac{n_R^2}{n_L^2}\right)
\end{gather*}
:::
From this we see several features:
* If we have no flow $j=0$, then there is a force in the direction of the lower
  density. This makes sense: we have a dam between two reservoirs, and the higher
  pressure wins.
* If there is flow, then the force decreases in magnitude until, at a critical 
  current, the direction of the force reverses.
* There is **no force** on the potential if $n_L = n_R$, regardless of the flow velocity.

Can we independently adjust $n_L$, $n_R$, and $j$?  This is not obvious.  We can
certainly start with $n_L = n_R$ and an arbitrary $j$ if we have no barrier.

We can also consider the energy density

#### Dimensional Analysis

We start with a naïve dimensional analysis of the 1D problem in the GPE with water
flowing over a weir of height $V_0$, width $\sigma$, and potential drop $\delta$.  Here we
consider steady state flow with (conserved) current $j$ and initial density $n_L = n_0$.
\begin{align*}
  %[\hbar] &= \frac{MD^2}{T}, &
  [m] &= M, & 
  [ν] &= \frac{D}{T}, &
  [g] &= \frac{MD^3}{T^2}, \\
  [n_0] &= \frac{1}{D}, &
  [j] &= \frac{1}{T}, \\
  [V_0] &= \frac{MD^2}{T^2}, &
  [\delta] &= \frac{MD^2}{T^2}, &
  [\sigma] &= D.
\end{align*}
From these, we can choose three independent quantities to set our units.  We start with
$m = n_0 = j = 1$:
\begin{gather*}
  1 = \underbrace{\frac{1}{n_0}}_{\text{distance}}
    = \underbrace{\frac{1}{j}}_{\text{time}}
    = \underbrace{m}_{\text{mass}}
    = \underbrace{\frac{j}{n_0}}_{\text{speed}}
    = \underbrace{\frac{m j^2}{n_0^2}}_{\text{energy}}.
\end{gather*}
This gives the following dimensionless quantities
\begin{gather*}
  %\tilde{\hbar} = \frac{\hbar n_0^2}{m j}, \qquad
  \tilde{ν} = \frac{ν n_0}{j}, \qquad
  \tilde{g} = \frac{g n_0^3}{m j^2}, \qquad
  \tilde{\delta} = \frac{\delta n_0^2}{m j^2}, \qquad
  \tilde{V}_0 = \frac{V_0 n_0^2}{m j^2}, \qquad
  \tilde{\sigma}_0 = \sigma n_0.
\end{gather*}

:::{admonition} About Naïve Dimensional Analysis
:class: dropdown

Since $m$ is intrinsic to the problem, and $j$ is conserved, we are pretty
sure that these are good parameter choices.  However, without further analysis, it is
not clear which of the other parameters will be of most relevance for an analysis.  For
example, instead of $g$, the speed of sound might be relevant, but should one use the
density $n_0$ or the minimum density at the top of the barrier?  Most likely, the
parameters discussed above in section {ref}`sec:Viscosity_NoViscosity` will play an
important role, but we need to start somewhere, so we start here.

We would also like to vary the viscosity $ν$, so we want to keep it as a
parameter.  Similarly with $V_0$ and $\sigma$.  Finally, we have $g$, but it might be
reasonable to vary this, so we choose to fix our units with $j$, $n_0$, and $m$.

Later refinements might include setting $2m = 1$ so that $\tilde{V}_0 = 1$ corresponds
to particles with just enough kinetic energy to overcome the barrier.  A better
parameterization might include the effects of pressure.
:::

To simplify our analysis, we start by considering the limit of a step-potential: $\sigma
\rightarrow 0$ and $V_0 = 0$, reducing our exploration to a 3-dimensional parameter
space characterized by the viscosity $ν$, interactions $g$ (compressibility), and step
size $\delta$.

```{code-cell}
nus = np.linspace(0.001, 0.2)
fs = []
Fs = []
for nu in nus:
    f = Flow1(n0=2.0, v0=1.0, V0=0.1, 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, figsize=(4,3))
ax.plot(nus, Fs)
ax.set(xlabel=r"$ν$", ylabel="$F$")

ind = np.argmax(Fs)
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}
nu = 0.05
V0s = np.linspace(0.001, 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, figsize=(4,3))
ax.plot(V0s, Fs)
ax.set(xlabel=r"$V_0$", ylabel="$F$")

ind = len(V0s)//2
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)
```




:::{admonition} Misc. Work (Incomplete)
:class: dropdown

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ν u_{,x} = mp
  = m \frac{j^2}{n_L} + P(n_L).
\end{gather*}
Far from the potential, $u_{,x} = 0$ so we can drop the viscosity term.  The net force
due to the potential is thus
\begin{gather*}
  F = m j(u_L - u_R) + P(n_L) - P(n_R) 
    = mj^2 \left(\frac{1}{n_L} - \frac{1}{n_R}\right) + P(n_L) - P(n_R).
\end{gather*}
Specifically, for the GPE, we have
\begin{gather*}
  F = \left(\frac{g}{2} - \frac{mj^2}{n_Ln_R(n_L + n_R)}\right)(n_L^2 - n_R^2)
    = \frac{gn_L^2}{2}\left(1 - \frac{2mj^2/gn_L^3}{n_R/n_L(1 + n_R/n_L)}\right)
                      \left(1 - \frac{n_R^2}{n_L^2}\right)
\end{gather*}
If we have no flow $j=0$, then this tells us that there is a force in the direction of
the lower density, but this decreases in magnitude as a velocity flows, until at a
critical current, the direction increases.








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 $ν=0$ flowing over a step
The momentum continuity equation gives
\begin{gather*}
  mnu^2 + F(x) + n\mathcal{E}'(n) - \mathcal{E}(n) - mν 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ν 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 first implement
the original equations in the form $y=(n, u, u', F)$, looking for stationary solutions
$\dot{n} = \dot{u} = 0$.  For an initial condition, we specify the density and velocity
to the left of the barrier and directly integrate.

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

These equations are not very easy to reason with, however: e.g., they look singular in the
limit of vanishing viscosity.  To deal with this, we replace the velocity with a
variable $q = ν u_{,x}$, and use the continuity equation to remove the velocity dependence:
\begin{gather*}
 u_{,x} = -u\frac{n_{,x}}{n}, \qquad
  j n_{,x} = n^2\frac{V'(x) + n_{,x}\mathcal{E}''(n)}{m} - n\underbrace{ν u_{,xx}}_{q_{,x}}\\
\end{gather*}

Instead, we use the continuity $j=nu$ and positivity of $n$
which suggest introducing the variable $q = \ln n$:
\begin{gather*}
  q = \ln n, \qquad n = e^{q}, \qquad
  q_{,x} = \frac{n_{,x}}{n}, \\
  q_{,xx} = -\frac{u_{,xx}}{u} + \frac{u_{,x}^2}{u^2} = -\frac{u_{,xx}}{u} + q_{,x}^2,\\
  \frac{u_{,x}}{u} = -q_{,x}, \qquad
  \frac{u_{,xx}}{u} = q_{,x}^2 - q_{,xx},\\
  -m j q_{,x} = -\frac{e^{2q}}{j}(V'(x)+\mathcal{E}''(e^{q})e^{q}q_{,x}) + mν (q_{,x}^2 - q_{,xx})\\
  m j q_{,x} = 
  \frac{e^{2q}V'(x)}{j}
  +\frac{e^{3q}}{j}\mathcal{E}''(e^{q})q_{,x}
  - mν (q_{,x}^2 - q_{,xx})\\
\end{gather*}

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

```{code-cell}
# Check some of the formulae
f = Flow()
res = f.solve(atol=1e-10, rtol=1e-10)

m, j, nL, nu = f.m, f.n0*f.v0, f.n0, f.nu

# Continuity
assert np.allclose(j, res.n*res.u)

# Force
P, PL = map(f.get_P, [res.n, nL])
x, V = res.x, f._V(res.x)
dE, dEL, ddE = f.get_Eint(res.n, d=1), f.get_Eint(nL, d=1), f.get_Eint(res.n, d=2)
F1 = m*j**2 * (1/nL - 1/res.n) + PL - P + m*nu*res.du
assert np.allclose(res.F, F1, atol=1e-7)

# This has some numerical artefacts at the kinks.  We smooth these:
res.j = res.u*res.n
e = (m*res.j*res.u**2/2 + res.j*(V + dE) - m*nu*res.u*res.du)
de = sp.interpolate.UnivariateSpline(x, e, s=1e-13).derivative()(x)
assert np.allclose(de, -m*nu*res.du**2, atol=3e-4)
```






































:::{admonition} Details.
:class: dropdown

The original equations should be solved to isolate $n_{,x}$ and $u_{,x}$:
\begin{gather*}
  \begin{pmatrix}
    u & n\\
    \mathcal{E}''(n) & u
  \end{pmatrix}
  \begin{pmatrix}
    n_{,x}\\
    u_{,x}
  \end{pmatrix}
  =
  \begin{pmatrix}
  0\\
  -V'(x)
  \end{pmatrix}\\
  \begin{pmatrix}
    n_{,x}\\
    u_{,x}
  \end{pmatrix}
  =
  \overbrace{\frac{1}{u^2+n\mathcal{E}''(n)}
  \begin{pmatrix}
    u & -n\\
    -\mathcal{E}''(n) & u
  \end{pmatrix}
  }^{
  \begin{pmatrix}
    u & n\\
    \mathcal{E}''(n) & u
  \end{pmatrix}^{-1}}
  \begin{pmatrix}
  0\\
  -V'(x)
  \end{pmatrix}
  =
  \frac{V'(x)}{u^2+n\mathcal{E}''(n)}
  \begin{pmatrix}
    n\\
    -u
  \end{pmatrix}.
\end{gather*}
Dividing by $n$ and $u$, we obtain
\begin{gather*}
  (\ln n)_{,x} = -(\ln u)_{,x} = \frac{-V'(x)}{u^2+n\mathcal{E}''(n)}.
\end{gather*}
This follows from the continuity equation which directly implies
\begin{gather*}
  \frac{n_{,x}}{n} = -\frac{u_{,x}}{n}.
\end{gather*}
:::

\begin{gather*}
         \\
  n_{,x} = \frac{n(V+\mathcal{E}')_{,x}}{u^2}, \qquad
  u_{,x} = -\frac{(V+\mathcal{E}')_{,x}}{u}\\
  \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}}{ν}\\
    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}}{ν}.
\end{gather*}
The momentum current is
\begin{gather*}
  j_{p} = \frac{j_0^2}{n} + \frac{P}{m} + ν\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 with height $V_0$, width $L$, and potential difference
$\delta V$ dropping from left to right.

```{code-cell}
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()
```

This gives us the following intrinsic parameters
\begin{gather*}
  [m] = M, \qquad
  [g] = \frac{MD^3}{T^2}, \qquad
  [L] = D, \qquad
  [V_0] = [\delta V] = \frac{MD^2}{T^2}, \qquad
  [ν] = \frac{1}{T}.
\end{gather*}
In addition, we have the extrinsic parameters related to the initial state
\begin{gather*}
  [n_0] = \frac{1}{D}, \qquad
  [u_0] = \frac{D}{T}.
\end{gather*}
We now choose three constants to define the units of mass, distance, and time.  Without
performing a full analysis, it is hard to tell what will ultimately be most appropriate,
so we start with an educated guess.  The mass scale is naturally set by $m=1$.

Some other quantities of relevance include the speed of sound $c_s = \sqrt{gn_0/m}$ and the
conserved current $j = n_0u_0$, $[j] = 1/T$.  More complicated parameters might be the
minimum speed required to overcome the barrier.

To start, we will focus on the initial state


### Analysis

```{code-cell}
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/ν$', 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}
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ν u_{,x}
  \right)
\right|^{x_R}_{x_L}
= - \int_{x_L}^{x_R} mν 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ν 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}
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}
# 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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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
```

----








# 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}
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}
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}
%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}
# 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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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}
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()
```