Show code cell content
import mmf_setup;mmf_setup.nbinit()
from pathlib import Path
import os
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
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Droplets#
Here we consider the problem of finding the solution for a spherical drop of some material with a saturating equation of state \(\mathcal{E}(n)\) that has a similar structure to that in nuclei. This problem has several interesting features and challenges we would like to explore:
Surface tension breaks the traditional convexity of thermodynamic relations for small systems. We start by formulating a model inspired by nuclear physics in Model: A Saturating Equation of State, which we then manipulate to simplify the mathematical structure.
Efficient numerical solutions requires solving boundary-value problems, which are somewhat finicky. Here we explore some techniques for taming this problem.
The numerical work solving the BVP starts from Asymptotic Behaviour. If you are most interested in how to solve this, please skip to that section.
Model: A Saturating Equation of State#
As a physical model, we choose a mean-field model for a quantum system like a Bose-Einstein condensate, but with a saturating equation of state, whose energy density (energy per unit volume) is
As we shall see shortly, homogeneous matter in such a system will minimize the energy per particle at the saturation density \(n_0\) with saturation energy \(-\epsilon_0\) per particle:
Our model will be to minimize the following energy functional for a “condensate wavefunction” \(\psi(\vect{x})\). This gives rise to Gross-Pitaevskii–like equations (GPE) with particle number density \(n = \abs{\psi}^2\).
Minimizing with respect to fixed particle number, we obtain the following GPE-like equation
Thermodynamic Limit#
First we note some thermodynamic properties of this system. Here is the equation of state:
Show code cell source
class EoSBase:
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):
pass
def E(self, n, d=0):
"""Return the dth derivative of the energy density."""
raise NotImplementedError
def __call__(self, n):
return self.E(n)
def E_N(self, n):
"""Return the energy per particle."""
return self.E(n)/n
def mu(self, n):
"""Return the chemical potential."""
return self.E(n, d=1)
def P(self, n):
"""Return the pressure."""
return self.mu(n) * n - self.E(n)
class EoS(EoSBase):
e0 = 1.0
n0 = 1.0
kappa = 1.0
def init(self):
k = self.kappa
self.E_poly = self.e0*self.n0*np.poly1d([k, -2*k, k-1, 0])
self.n_poly = np.poly1d([self.n0, 0])
self.mu_poly = self.E_poly.deriv() / self.n0
self.P_poly = self.mu_poly * self.n_poly - self.E_poly
self.mu_min = self.mu(self.mu_poly.deriv().roots[-1])
super().init()
def E(self, n, d=0):
"""Return the dth derivative of the energy density."""
n_ = n / self.n0
return self.E_poly.deriv(d)(n_) / self.n0**d
e0 = n0 = 1
fig, axs = plt.subplots(1, 4, figsize=(12, 2.5))
for i, k in enumerate([0.5, 1, 1.5]):
E = EoS(kappa=k, e0=e0, n0=n0)
E_N = E.E_N
mu = dE = E.mu
mu0 = mu(n0)
n_unstable = np.linspace(0, n0)[1:]
n_stable = np.linspace(n0, 2)
n = n0*np.poly1d([1, 0])
ax = axs[0]
fmt_stable = f'-C{i}'
fmt_unstable = f':C{i}'
fmt_mixed = 'k--.'
ax.plot(n_stable/n0, E(n_stable)/n0/e0, fmt_stable)
ax.plot([0, 0], [0, 0.5], fmt_stable)
ax.plot(n_unstable/n0, E(n_unstable)/n0/e0, fmt_unstable)
ax.plot(n_unstable[[0, -1]]/n0, mu0*n_unstable[[0, -1]]/n0/e0, fmt_mixed)
ax.set(ylim=(-1.5, 0.5),
xlabel="$n$ [$n_0$]", ylabel=r"$\mathcal{E}=E/V$ [$\epsilon_0 n_0$]")
ax = axs[1]
ax.plot(n_stable/n0, E(n_stable)/n_stable/e0, fmt_stable)
ax.plot(n_unstable/n0, E(n_unstable)/n_unstable/e0, fmt_unstable)
ax.plot([0, 0], [mu0/e0, 0.1], fmt_stable)
ax.plot([0, 1], [mu0/e0, mu0/e0], fmt_mixed)
ax.set(ylim=(-1.1, 0.0),
xlabel="$n$ [$n_0$]", ylabel=r"$E/N = \mathcal{E}/n$ [$\epsilon_0$]")
ax = axs[2]
ax.plot(n_stable/n0, mu(n_stable)/e0, fmt_stable)
ax.plot([0, 0], [mu0/e0, -1.5], fmt_stable)
ax.plot([0, 1], [mu0/e0, mu0/e0], fmt_mixed)
ax.plot(n_unstable/n0, mu(n_unstable)/e0, fmt_unstable)
assert np.allclose(E.mu_min, -e0*(1+k/3)) # Check mu_min formula
ax.set(ylim=(-1.5, 0.5),
xlabel="$n$ [$n_0$]", ylabel=r"$\mu = \mathcal{E}'(n)$ [$\epsilon_0$]")
ax = axs[3]
ax.plot(mu(n_stable)/e0, E.P(n_stable)/n0/e0, fmt_stable, label=f"$\kappa={k:.1f}$")
ax.plot([-1.5, mu0/e0], [0, 0], fmt_stable)
ax.plot([mu0/e0], [0], fmt_mixed)
ax.plot(mu(n_unstable)/e0, E.P(n_unstable)/n0/e0, fmt_unstable)
ax.set(xlim=(-1.5, 0.5), ylim=(-0.6, 1.1),
xlabel="$\mu$ [$\epsilon_0$]", ylabel=r"$P=\mu n - \mathcal{E}$ [$\epsilon_0 n_0$]")
axs[3].legend()
plt.tight_layout()
<>:91: SyntaxWarning: invalid escape sequence '\k'
<>:96: SyntaxWarning: invalid escape sequence '\m'
<>:91: SyntaxWarning: invalid escape sequence '\k'
<>:96: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_5133/1817750480.py:91: SyntaxWarning: invalid escape sequence '\k'
ax.plot(mu(n_stable)/e0, E.P(n_stable)/n0/e0, fmt_stable, label=f"$\kappa={k:.1f}$")
/tmp/ipykernel_5133/1817750480.py:96: SyntaxWarning: invalid escape sequence '\m'
xlabel="$\mu$ [$\epsilon_0$]", ylabel=r"$P=\mu n - \mathcal{E}$ [$\epsilon_0 n_0$]")
This describes a saturating system with an unstable branch between \(n_0 \in [0, 1]\). In this region, one will have a mixed phase through the Maxwell construction (dashed black line) between the two phases shown: one at \(n=0\) corresponding to the vacuum, and the other at \(n=n_0\). Note that the Maxwell construction reduces to a point in the \((\mu, P)\) plot (grand-canonical ensemble) at \(\mu_0 = -\epsilon_0\).
In infinite matter, or an infinitely large droplet, the surface will have no curvature, so the pressures of the two phases must match. In this case \(P=0\) in both the vacuum and the saturating homogeneous phase with density \(n_0\) and chemical potential \(\mu_0 < 0\). Note: until \(\mu_0 < -mc^2\), the vacuum is stable (no anti-particles are produced), and so can be in chemical equilibrium for any negative \(\mu < 0\).
Finite Droplets#
In a finite droplet, the surface is curved, and the surface tension will increase the pressure inside the droplet, causing \(n>n_0\) and increasing \(\mu > \mu_0\). Larger and larger surface tension will correspond to greater curvature – i.e. smaller droplets. Thus, as \(\mu\) increases towards zero, the size of the droplet will shrink. This is slightly paradoxical since one typically associates increasing the chemical potential with increasing the overall particle number. In infinite matter, this arises from the convexity of the energy \(E(N)\):
Finite droplets break this convexity.
To get an idea of how this works, consider a spherical droplet with a thin surface of tension \(\sigma\). The surface has an energy \(E_{\sigma}(r) = 4\pi \sigma r^2\) proportional to the area of the sphere. This leads to a compressional force \(F = E_{\sigma}'(r) = 8 \pi \sigma r\) and a pressure (force per unit area) \(P = F/4\pi r^2 = 2\sigma /r\). This should equal the pressure in the droplet:
This can also be derived using the compressible liquid drop model (CLDM). Let \(n = n_i > n_0\) be the density inside the droplet, and \(n_0 = 0\) be the density outside. The energy \(E\) and particle number \(N\) of a droplet of radius \(r\) are:
The size of the droplet will minimize the energy for fixed \(N\):
This gives the same relationship between the surface tension and the pressure we found above.
Show code cell source
fig, ax = plt.subplots(figsize=(5,3))
E = EoS(kappa=1)
sigma = 0.2
for i, R in enumerate([1, 2, 3]):
rs = np.linspace(0, R, 1000)[1:]
Vs = 4/3*np.pi * rs**3
ns = np.array([(E.P_poly - _P).roots[0] for _P in 2*sigma / rs]).real
Es = E(ns)*Vs + sigma*4*np.pi * rs**2
Ns = ns*Vs
V = 4/3*np.pi*R**3
Nmax = Ns[-1]
ns_ = np.concatenate([Ns/V, [Nmax/V]])
Es_ = np.concatenate([Es/V, [E(Nmax/V)]])
ax.plot(ns_, Es_, f"-C{i+1}", label=f"${R=}$")
ns = np.linspace(E.n0, 2.1)
ax.plot(ns, E.E(ns), '-C0')
ns = np.linspace(0, E.n0)
ax.plot(ns, E.mu(E.n0)*ns, '--k')
ax.set(xlabel="n", ylabel="$\mathcal{E}(n)$",
xlim=(-0.1, 2.1), ylim=(-1.2, 1.0),
title=rf"Droplet with $\sigma={sigma:.4g}$")
ax.legend()
plt.tight_layout()
<>:23: SyntaxWarning: invalid escape sequence '\m'
<>:23: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_5133/1395581838.py:23: SyntaxWarning: invalid escape sequence '\m'
ax.set(xlabel="n", ylabel="$\mathcal{E}(n)$",
What about a bubble? First we consider a fixed total volume of radius \(R\). Now the volume of matter is that excluded by the bubble:
Thus, unless our equation of state supports negative pressure OR the surface tension is negative, then we cannot support bubbles.
Note
We might worry that by restricting our attention to spherical drops, we are missing some sort of mixed phase. To generalize the problem, we note that we can perform a constrained minimization over all possible states \(\psi\), which we do with a Legendre transform:
If there are multiple states that satisfy the stationarity condition, the one which minimize the energy \(E[\psi]\) should be chosen. This guarantees that \(P(\mu)\) is convex, but not that \(E(N)\) is convex as the following example demonstrates.
Consider the following artificial problem with \(\psi = (n, \theta)\):
The stationarity condition gives
Of the possible solutions \(\theta \in [0, \pi]\), only \(\theta = 0\) minimizes \(E\), so we have \(\mu = - \sin 2 n\)
N = np.linspace(0, 6*np.pi, 1000)
mu = - np.sin(2*N)
E = - np.cos(N)**2
P = mu*N - E
plt.plot(mu, P)
[<matplotlib.lines.Line2D at 0x7fe042f2ffb0>]
Units#
To simplify our formula and numerical work, we choose a convenient set of units. We simplify the differential equation to the following form
The EoS suggests using the following units \(\epsilon_0 = n_0 = 2m = 1\). To restore results, use the appropriate factors of unity
The ground state can be chosen to be real, so we shall only consider real functions \(\psi\) and write \(n = \psi^2\) for simplicity.
Homogeneous States#
First we consider the homogeneous states where \(\nabla^2 \psi = 0\). This will be the asymptotic form of our solutions far from the center of a droplet or bubble and requires
For our equation of states,
For a smaller \(\mu\), the only asymptotic solution is the vacuum \(n=0\). Conversely, the vacuum \(n=0\) can only exist if \(\mu \leq \mathcal{E}'(0) = \epsilon_0(\kappa - 1)\). Thus, droplets or bubbles require
Spherical Solutions#
Switching to spherical coordinates and our dimensionless units, we have the following equations:
Although there is an apparent singularity at \(r=0\), from physical considerations, nothing special happens at the origin, so we expect \(\psi(n)\) to be smooth for small \(n\).
Spherical Bubbles#
Another possibility is a spherical bubble with \(\psi(\infty) = \sqrt{n_0}\):
We are looking for solutions that are finite at \(r=0\), but which fall to zero as \(r\rightarrow \infty\).
To study spherical droplets, we need \(\psi(\infty) \rightarrow 0\) so that \(n(\infty) = 0\). To elucidate the asymptotic The asymptotic , which requires
It will turn out that droplets require \(\mu < \mathcal{E}'(0)\), and so we will rewrite the equation for numerical work as follows:
We can then choose the coordinate \(x \propto r\) to absorb the constant:
This gives us our working equation:
The qualitative properties of the system are described by the following dimensionless constants:
Scaling everything out, we have the following differential equation:
First we note a couple of trivial solutions: first, we have homogenous solutions \(\psi(r) = 0\) (corresponding to the vaccum), and \(\psi = \sqrt{n_0}\) where \(\mathcal{E}'(n_0) = \mu\) is also a solution. For our form, \(n_0 = (1 \pm \sqrt{1+4\mu})/2\)
As discussed above, droplets require \(\mu < \mathcal{E}'(0)\), so we can write
where the coefficient on the left is positive and can be absorbed into a redefinition of \(x \propto r\), and we can scale out the dependence on \(\mu\) into a single parameter \(\alpha\):
Important
This gives us our working equation, where we use \(y(x) = \psi(r)\) to emphasize this change:
where the bounds and form of \(V(n)\) are for our EoS.
Numerical Attempt 1#
We start with a brute-force solution of the problem as written:
Show code cell source
from scipy.integrate import solve_bvp
class BVP1(EoS):
hbar2_2m = 1.0
alpha = 2.0
def V(self, n):
return 3/4*n*(3*n-4)
def rhs(self, x, q):
"""Return dq/dx. Note q = (y, dy)."""
y, dy = q
ddy = (1 + self.alpha * self.V(y**2))*y - 2*dy/x
dq = np.asarray([dy, ddy])
return dq
def bc(self, qa, qb):
"""Return the boundary conditions."""
ya, dya = qa
yb, dyb = qb
bc = [dya, yb]
return np.asarray(bc)
def get_initial_guess(self, sol=None, k=1):
"""Return `(xs, ys, dys)` for an initial guess.
Note: this function needs self.alpha, self._x0 and self._x1.
"""
alpha = self.alpha
x0, x1 = self._x0, self._x1
xs = np.linspace(x0, x1)
if sol is None:
ys = 1/np.cosh(k*xs)
dys = k*np.sinh(k*xs)/np.cosh(k*xs)**2
else:
(ys, dys) = sol.sol(xs)
return map(np.asarray, (xs, ys, dys))
def solve(self, sol0=None, x0=0.01, x1=20.0, alpha=None, k=1, **kw):
"""Return the solution `sol` to the BVP.
Arguments
---------
sol0 : Solution, None
A previous solution to use as a starting point.
x0, x1 : float
Inner and out radii for boundary conditions.
mu : float
Chemical potetial.
"""
self._x0, self._x1 = x0, x1
if alpha is not None:
self.alpha = alpha
xs, ys, dys = self.get_initial_guess(sol=sol0, k=k)
y0s = np.asarray([ys, dys])
sol = solve_bvp(self.rhs, self.bc, x=xs, y=y0s, **kw)
mu = self.mu(0) - (self.mu(0) - self.mu_min) / self.alpha
x, (y, dy) = sol.x, sol.y
sol.alpha = self.alpha
sol.mu = mu
sol.r = x/np.sqrt((self.mu(0) - mu)/self.hbar2_2m)
sol.n = y**2
return sol
def find_sols(self, alpha0=2, dalphas=[0.2, -0.1], n_min=0.1, N=10000):
"""Starting from alpha0, iterate to try to find as many sols as possible."""
alpha = alpha0
sol_ = self.solve(alpha=alpha)
sols = [sol_]
alphas = [alpha]
for dalpha in dalphas:
sol0 = sol_
alpha = alpha0
for _N in range(N):
while abs(dalpha) > 1e-8 and 4/3 < alpha:
with np.errstate(all="ignore"):
# Suppress floating point warnings.
sol = b.solve(sol0=sol0, alpha=alpha+dalpha)
if sol.success and sol.n[0] > n_min:
break
dalpha /= 2
if not sol.success or sol.n[0] < n_min:
break
sol0 = sol
alpha = b.alpha
sols.append(sol)
alphas.append(alpha)
inds = np.argsort(alphas)
alphas = np.array(alphas)[inds]
sols = np.array(sols)[inds]
return alphas, sols
b = BVP1()
sol = b.solve()
r, n = sol.r, sol.n
fig, ax = plt.subplots(figsize=(3, 2))
ax.plot(r, n)
ax.set(xlabel="$r$", ylabel="$n$", title=rf"$\alpha={sol.alpha:.2g}$, $\mu={sol.mu:.2g}$");
We see that this works, but it is very sensitive to the initial guess. Here we start with this solution, and then slowly increase and decrease \(\alpha\). Note: we must take many small steps, otherwise the solution does not converge, however, the general procedure is sound.
Show code cell source
b = BVP1()
alphas, sols = b.find_sols()
fig, ax = plt.subplots(figsize=(3, 2))
for alpha, sol in zip(alphas, sols):
r, n = sol.r, sol.n
ax.plot(r, n)
ax.set(xlim=(0, 20),
xlabel="$r$", ylabel="$n$",
title=rf"${alphas[0]:.2f} \leq \alpha \leq {alphas[-1]:.2f}$");
From these solutions, we can estimate the surface tension \(\sigma\) as follows from the formula discussed above:
In the following, we estimate the radius of the surface as the point where \(n = n_0/2\), then plot \(\sigma \approx P(n_0) r/2\)
Show code cell source
from scipy.interpolate import InterpolatedUnivariateSpline
rs = np.array([sol.r[np.argmin(abs(sol.n - 0.5))] for sol in sols])
mus = np.array([sol.mu for sol in sols])
rs = []
for sol in sols:
if sol.n.max() > 0.5:
_r = InterpolatedUnivariateSpline(sol.r, sol.n-0.5).roots()[0]
else:
_r = 0
rs.append(_r)
rs = np.asarray(rs)
fig, ax = plt.subplots(figsize=(3, 2))
n0 = b.n0
E0 = b.E(n0)
ax.plot(rs, (mus*n0 - E0)*rs/2)
#plt.axhline([0.5], ls=":", c="y")
ax.set(xlabel="$r$", ylabel=r"$\sigma \approx P(n_0)r/2$");
Asymptotic Behaviour#
A good starting point is to consider the asymptotic behavior of droplets. For large \(r\rightarrow
\infty\), we expect that we can neglect the second term, and are looking for solutions
that vanish \(\psi \rightarrow 0\).
Noting that \(V(y^2) \propto y^2\), we can probably neglect the right-hand-side as \(x \rightarrow
\infty\) since it scales as \(y^3\). If we also neglect the \(2y'/x\) term, then we have
simply \(y'' = y\), which implies
We can do a bit better by nothing that
satisfies \(y'' + 2y'/x - y = 0\). Neglecting the \(\alpha V(y^2)y\) thus requires
The asymptotic forms give better boundary conditions at large \(x\), which will reduce truncation errors.
At the origin, we must have that \(y' \rightarrow 0\) at least as fast a \(x\), otherwise the second term will diverge. Since we physically expect the solution to be smooth at the origin, we expect we can expand it as a series in \(x\). Expanding to second lowest order, we have:
This allows us to refine the boundary conditions at finite \(x\).
Numerical Attempts 2 and 3#
We now include these refined boundary conditions:
BVP2: \(y'(x_1) + y(x_1) = 0\).BVP3: \(y'(x_1) + (1+x_1^{-1})y(x_1) = 0\).
Show code cell source
from scipy.integrate import solve_bvp
class BVP2(BVP1):
def bc(self, qa, qb):
"""Return the boundary conditions."""
ya, dya = qa
yb, dyb = qb
xa, xb = self._x0, self._x1
bc = [
3*dya - xa*ya * (1 + self.alpha * self.V(ya**2)),
yb + dyb
]
return np.asarray(bc)
class BVP3(BVP1):
def bc(self, qa, qb):
"""Return the boundary conditions."""
ya, dya = qa
yb, dyb = qb
xa, xb = self._x0, self._x1
bc = [
3*dya - xa*ya * (1 + self.alpha * self.V(ya**2)),
dyb + (1+1/xb)*yb
]
return np.asarray(bc)
b = BVP2()
sol = b.solve()
r, n = sol.r, sol.n
fig, ax = plt.subplots(figsize=(3, 2))
ax.plot(r, n)
ax.set(xlabel="$r$", ylabel="$n$", title=rf"$\alpha={sol.alpha:.2g}$, $\mu={sol.mu:.2g}$");
Show code cell source
b = BVP3()
alphas, sols = b.find_sols()
fig, ax = plt.subplots(figsize=(3, 2))
for alpha, sol in zip(alphas, sols):
r, n = sol.r, sol.n
ax.plot(r, n)
ax.set(xlim=(0, 20),
xlabel="$r$", ylabel="$n$",
title=rf"${alphas[0]:.2f} \leq \alpha \leq {alphas[-1]:.2f}$");
This shows that refining the boundary conditions does not improve the overall convergence if the guess is not good, however, it improves the accuracy, as we shall now see.
Error Analysis#
The main improvement is with the boundary conditions, which improve the rate of convergence, allowing us to solve the problem on a restricted domain.x
Important
You should check these relationships quantitatively. The errors should scale as expected: this provides a very strong and invaluable check on aspects of your code.
In our case, we expect the following three sources of errors due to the boundary conditions at \(x_1\), which we obtain by keeping the first neglected term:
For estimates, we use the asymptotic solution, though dropping the denominator would also give approximately correct results:
If the coefficient \(A\) cannot be easily estimated (i.e. from a solution), then you can choose it manually to roughly match the actual errors seen. The important thing is that the slope of the scaling is correct:
Show code cell source
from matplotlib.gridspec import GridSpec
gs = GridSpec(2, 2)
fig = plt.figure(figsize=(7, 6))
ax = fig.add_subplot(gs[0, :])
ax0 = fig.add_subplot(gs[1, 0])
ax1 = fig.add_subplot(gs[1, 1])
b0 = BVP2(alpha=1.6)
sol0 = b0.solve(k=0.33)
assert sol0.success
x0s = 10**np.linspace(0.5, -3)
x1s = np.linspace(8, 30)
x_test = np.linspace(0, x1s.min())
ax.plot(sol0.x, sol0.n)
ax.plot(x_test, sol0.sol(x_test)[0]**2, lw=2)
ax.set(xlabel='$x$', ylabel='$n$')
kw = dict(tol=1e-10, max_nodes=10000)
# Estimate coefficient A
x1 = 20.0
A = sol0.sol(x1)[0]/(np.exp(-x1)/x1)
y1 = A*np.exp(-x1s)/x1s
err1_est = [
abs(y1),
abs(y1)**2,
abs(y1)**3
]
for i, _BVP in enumerate([BVP1, BVP2, BVP3]):
b = _BVP(alpha=1.6)
sol0 = b.solve(k=0.33)
assert sol0.success
sols0 = [b.solve(x0=x0, sol0=sol0, **kw) for x0 in x0s]
sols1 = [b.solve(x1=x1, sol0=sol0, **kw) for x1 in x1s]
assert all(sol.success for sol in sols0)
assert all(sol.success for sol in sols1)
errs0 = [abs(sol.sol(x_test) - sols0[-1].sol(x_test)).max() for sol in sols0]
errs1 = [abs(sol.sol(x_test) - sols1[-1].sol(x_test)).max() for sol in sols1]
ax0.loglog(x0s[:-1], errs0[:-1], f'+-C{i}', label=f"BVP{i+1}")
ax1.semilogy(x1s[:-1], errs1[:-1], f'+-C{i}')
ax1.plot(x1s, err1_est[i], f'--C{i}', label=f"$y(x_1)^{{{i+1}}}$")
ax0.set(xlabel="$x_0$", ylabel="max abs err");
ax1.set(ylim=(1e-12, 1), xlabel="$x_1$", ylabel="max abs err");
ax0.legend()
ax1.legend()
plt.tight_layout()
From the scaling, we see that our error analysis is not completely correct. In particular, both
BPV1 and BVP2 seem to have the same error scaling, but as \(y_1^2\) rather than
linear in \(y_1 = y(x_1)\). I don’t fully understand this yet, but has to do with the
fact that the maximum error is incurred somewhere for \(x < x_1\) less than \(x_1\).
Likely, the error in the boundary condition is enforced by the \(e^{x}\) solution which
dampens extremely quickly before the location where the maximum error is incurred.
The scaling for BVP3 is consistent with, but slightly better than \(y_1^3\).
Additional Ideas (Incomplete)#
Use a reasonable ansatz to minimize \(E[\psi]\) subject to fixed \(N[\psi]\). See below. This can be used to find a good initial guess.
Use the asymptotic form to construct a better solution which factors off the singular terms.
Use variables like \(y = e^{q(x)}\) which will make manifest the positivity of the solution.
Initial Guess#
If the problem is sufficiently unstable, then one might be able to tame it by introducing a constrained form of the solution. For example, the following functional form has the correct asymptotic properties for a bubble with radius \(X\) and width \(\lambda\):
One can use these to directly try to minimize the energy functional as a function of \(X\), \(y_0\), and \(\lambda\). Note that the asymptotic behaviours constrain, i.e. \(\lambda = 2\), but we are looking for global minimization, so it is best not to try to use the properties of the “tails” to “wag the dog”.
One can proceed analytically, but we will write a simple numerical routine:
To find a solution, we return to the original definition.
Show code cell source
class BVPGuess(BVP1):
def get_initial_guess(self, sol=None):
x = np.linspace(self._x0, self._x1, 1000)
def y(self, x, p):
y0, X, lam = p
return y0/4*(1 + np.tanh((x+X)/lam))*(1 - np.tanh((x-X)/lam))
def _get_N(self, x, y):
return np.trapz(x**2 * y**2, x)
def _get_E(self, x, y):
dy = np.gradient(y, x, edge_order=2)
n = y**2
Vint = self.alpha * 3/4*(n**3 - 2*n**2)
return np.trapz(x**2 * (dy**2 + Vint + n), x)
def go(self, p, plot=True):
x = np.linspace(0.001, 20, 1000)
#X, lam = p
y = self.y(x, p)
#N_ = self._get_N(x, y)
#y0 = np.sqrt(N/N_)
#y *= y0
#assert np.allclose(self._get_N(x, y), N)
E = self._get_E(x, y)
if plot:
plt.plot(x, y**2)
plt.title(f"{N=}, {E=}")
return E
g = BVPGuess()
g.go((10, 2))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[11], line 32
30 return E
31 g = BVPGuess()
---> 32 g.go((10, 2))
Cell In[11], line 21, in BVPGuess.go(self, p, plot)
19 x = np.linspace(0.001, 20, 1000)
20 #X, lam = p
---> 21 y = self.y(x, p)
22 #N_ = self._get_N(x, y)
23 #y0 = np.sqrt(N/N_)
24 #y *= y0
25 #assert np.allclose(self._get_N(x, y), N)
26 E = self._get_E(x, y)
Cell In[11], line 6, in BVPGuess.y(self, x, p)
5 def y(self, x, p):
----> 6 y0, X, lam = p
7 return y0/4*(1 + np.tanh((x+X)/lam))*(1 - np.tanh((x-X)/lam))
ValueError: not enough values to unpack (expected 3, got 2)
Show code cell source
from scipy.interpolate import InterpolatedUnivariateSpline
n0s = np.array([sol.sol(0)[0]**2 for sol in sols])
Xs = []
for sol in sols:
if sol.n.max() > 0.5:
_X = InterpolatedUnivariateSpline(sol.x, sol.n-0.5).roots()[0]
else:
_X = 0
Xs.append(_X)
Xs = np.asarray(Xs)
fig, ax = plt.subplots(figsize=(3, 2))
plt.plot(Xs, 1+b.alpha*b.V(n0s))
plt.plot(Xs, -6*np.exp(3*Xs)/(1+np.exp(Xs))**4)
[<matplotlib.lines.Line2D at 0x7fe03b426990>]
Numerical Attempt 4#
Now we note that \(A(x) \geq 0\) everywhere. This can be enforced by letting \(A(x) = e^{B(x)}\):
b1 = BVP1()
b = BVP2()
sol1 = b1.solve(R=5)
sol = b.solve(R=5)
x1, (y, dy) = sol1.x, sol1.y
x, (A, dA) = sol.x, sol.y
dA, ddA = b.rhs(x, sol.y)
B = np.log(A)
n = np.exp(B-x)**2
dB = dA/A
ddB = ddA/A - dB**2
assert np.allclose(dB, np.gradient(B, x, edge_order=2), rtol=0.02, atol=0.01)
assert np.allclose(ddB, np.gradient(dB, x, edge_order=2), rtol=0.05, atol=0.05)
assert np.allclose(ddB, -dB**2 + 2*(1-1/x)*dB + b.V(n) + 2/x, rtol=0.05, atol=0.05)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[13], line 3
1 b1 = BVP1()
2 b = BVP2()
----> 3 sol1 = b1.solve(R=5)
4 sol = b.solve(R=5)
5 x1, (y, dy) = sol1.x, sol1.y
Cell In[5], line 58, in BVP1.solve(self, sol0, x0, x1, alpha, k, **kw)
56 xs, ys, dys = self.get_initial_guess(sol=sol0, k=k)
57 y0s = np.asarray([ys, dys])
---> 58 sol = solve_bvp(self.rhs, self.bc, x=xs, y=y0s, **kw)
60 mu = self.mu(0) - (self.mu(0) - self.mu_min) / self.alpha
61 x, (y, dy) = sol.x, sol.y
TypeError: solve_bvp() got an unexpected keyword argument 'R'
class BVP4(BVP1):
def rhs(self, x, q):
"""Return dq/dr. Note q = (B, dB)."""
B, dB = q
n = np.exp(B-x)**2
ddB = -dB**2 + 2*(1-1/x)*dB + self.V(n) + 2/x
dq = np.asarray([dB, ddB])
return dq
def bc(self, qa, qb):
"""Return the boundary conditions."""
B0, dB0 = qa
Bb, dBb = qb
x0 = self._r0
bc = [
dB0 - (1+x0/3*(1 + self.V(np.exp(2*B0)))),
dBb
]
return np.asarray(bc)
def solve(self, sol0=None, r0=0.01, R=20.0, mu=-0.3, k=1, **kw):
"""Return the solution `sol` to the BVP.
Arguments
---------
sol0 : Solution, None
A previous solution to use as a starting point.
r0, R : float
Inner and out radii for boundary conditions.
mu : float
Chemical potetial.
"""
self._r0 = r0
self._mu = mu
if sol0 is None:
xs = np.linspace(r0, R)
ys = 1/np.cosh(k*xs)
dys = k*np.sinh(k*xs)/np.cosh(k*xs)**2
As = ys*np.exp(xs)
dAs = As + dys*np.exp(xs)
Bs = np.log(As)
dBs = dAs/As
else:
xs, (Bs, dBs) = sol0.x, sol0.y
y0s = np.asarray([Bs, dBs])
sol = solve_bvp(self.rhs, self.bc, x=xs, y=y0s, **kw)
x, (B, dB) = sol.x, sol.y
A = np.exp(B)
y = A*np.exp(-x)
sol.r = x/np.sqrt(-mu)
sol.n = y**2
return sol
b = BVP4()
sol = b.solve(mu=-0.3, R=5)
r, n = sol.r, sol.n
fig, ax = plt.subplots()
ax.plot(r, n)
ax.set(xlabel="$r$", ylabel="$n$", title=f"$\mu={b._mu:.2f}$");
<>:59: SyntaxWarning: invalid escape sequence '\m'
<>:59: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_5133/2650077862.py:59: SyntaxWarning: invalid escape sequence '\m'
ax.set(xlabel="$r$", ylabel="$n$", title=f"$\mu={b._mu:.2f}$");
/tmp/ipykernel_5133/2650077862.py:5: RuntimeWarning: overflow encountered in exp
n = np.exp(B-x)**2
/home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/scipy/integrate/_bvp.py:309: RuntimeWarning: invalid value encountered in subtract
0.125 * h * (f[:, 1:] - f[:, :-1]))
/tmp/ipykernel_5133/2650077862.py:16: RuntimeWarning: overflow encountered in exp
dB0 - (1+x0/3*(1 + self.V(np.exp(2*B0)))),
/tmp/ipykernel_5133/2650077862.py:5: RuntimeWarning: overflow encountered in square
n = np.exp(B-x)**2
/tmp/ipykernel_5133/1176432483.py:8: RuntimeWarning: overflow encountered in multiply
return 3/4*n*(3*n-4)
/tmp/ipykernel_5133/1176432483.py:8: RuntimeWarning: overflow encountered in scalar multiply
return 3/4*n*(3*n-4)
/tmp/ipykernel_5133/2650077862.py:6: RuntimeWarning: overflow encountered in square
ddB = -dB**2 + 2*(1-1/x)*dB + self.V(n) + 2/x
/tmp/ipykernel_5133/2650077862.py:6: RuntimeWarning: invalid value encountered in add
ddB = -dB**2 + 2*(1-1/x)*dB + self.V(n) + 2/x
/home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/scipy/integrate/_bvp.py:40: RuntimeWarning: invalid value encountered in subtract
df_dy[:, i, :] = (f_new - f0) / hi
/home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/scipy/integrate/_bvp.py:566: RuntimeWarning: invalid value encountered in divide
r_middle /= 1 + np.abs(f_middle)
/home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/scipy/integrate/_bvp.py:567: RuntimeWarning: invalid value encountered in divide
r1 /= 1 + np.abs(f1)
/home/docs/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/scipy/integrate/_bvp.py:568: RuntimeWarning: invalid value encountered in divide
r2 /= 1 + np.abs(f2)
Bessel Functions#
One approach is to try to use some known functions to construct approximate solutions, then numerically find corrections. A similar equation is solved by the spherical Bessel functions, specifically for \(n=0\)
Thus, we can consider a solution of the form \(\psi(r) = Aj_0(r\sqrt{\mu})\):
#a = 1/\sqrt{\mu}
from sympy import var, sinh, sqrt
r, A, a, mu = var('r, A, a, mu')
a = 1/sqrt(-mu)
psi = A*sinh(r/a)/(r/a)
display((psi.diff(r, r) + 2/r*psi.diff(r) + mu*psi).simplify())
r = np.linspace(0, 10)[1:]
plt.plot(r, np.sinh(r)/r)
[<matplotlib.lines.Line2D at 0x7fe03902da30>]
