text:
\begin{document}
\maketitle

\begin{abstract}
 A numerical solution to the Fokker-Planck equation using a two-level scheme is presented. The Fokker-Planck (FP) equation is of parabolic type equation govern the time evolution of probability density function of the stochastic processes. The FP equation also preserves the positivity and conservative of the total probability. A Chang-Cooper discretization scheme is used to ensure the positiveness and conservation of the total probability with second-order accuracy. We investigate a two-level scheme with factor-three-coarsening strategy and have a significant reduction in computations and CPU time. Numerical experiments are performed to validate the efficiency and second-order accuracy of the proposed two-level algorithm with backward time-difference schemes.
\end{abstract}
{\bf{keywords:}} Stocastic process; Fokker-Planck equation; Chang-Cooper scheme; two-level scheme; staggered grids; finite difference \\
{\bf{MSC 2010:}} 35Q84; 49K20; 65N55

\section{Introduction}

Stochastic models are used in diverse field as ecology, genetics, economics and engineering. Closed form solutions of such models are know, however, only for some of the simplest drift and diffusion functions. Therfore, there has always been a need of numerical methods that solves comples stochastic models and hence the Fokker-Planck equation. The global dynamical behaviour of a nonlinear system with noise is formally described by the probability density function (PDF) evolution along deterministic and diffusuion that satisfies the Fokker-Planck (FP) partial differential equation \cite{Risken}. The FP system that has the capability of connecting stochastic and deterministic dynamics has been applied to various applications in physics, chemistry, biology and finance \cite{Carlini, Chang, Risken, RoyBorzi2018, Spencer}.

In this article, we consider the FP equations that corresponds to the stochastic differential equations. In particular, the stochastic process defined by the following multidimensional model \cite{Mohammadi}
\begin{eqnarray}\label{stocasticmodel}

  dX_t &=& b(X_t,t)dt+\sigma(X_t,t)dW_t \\
  X(t_0) &=& X_0,
\end{eqnarray}
where $X_t\in{\mathbb{R}}^{d}$ is the state variable and $dW_t\in{\mathbb{R}}^l$ denotes the Wiener process. Moreover, $ \sigma\in{\mathbb{R}}^{d\times l}$ is a full rank dispersion matrix. Note that a statistical distribution can describes the state of the stochastic process. For this, the probability density function (PDF) distribution and the evolution of this PDF distribution can be modelled by the FP equation.

The numerical solution of the FP equations has been obtained by several researchers. One of the most popular scheme in this regard which solves the linear FP equation is the Chang-Cooper (CC) scheme introduced by Chang and Cooper in $1970$ \cite{Chang}. One of the impotent features of CC scheme is that the discrete solution preserves some intrinsic properties of the original given problem, one such as positivity and conservation of the total probability. Later, several improvements have been done \cite{Drozdov}, where we have seen high order finite difference schemes and also the nonlinear case. Finite element schemes have also been discussed, see \cite{Spencer}. It is also worth noting that some semi-analytic techniques are employed to solve the FP equation, for example, in \cite{Dehghan} the FP equation is investigated by the Adomian decomposition method. In \cite{Tatari}, variational iteration method is presented to solve the FP equation. Moreover, a finite difference scheme with cubic $C^1$-spline collocation method for solving the nonlinear Fokker-Planck equation is presented in \cite{Sepehrian}. A fast algorithm for the numerical solution of the FP equation is presented by \cite{Palleschi90, Palleschi90} and a finite difference scheme, in one-dimension, using a staggered grid to solve the Fokker-Planck equations with drift-admitting jumps is presented in \cite{Chen}. In the year 2020, the research to find the numerical solution to the stochastic models and henece the FP equation is still on; e.g., in \cite{Carlini}, a discretization scheme is developed to solve the one-dimensional nonlinear Fokker-Planck-Kolmogorov equation that preserves the nonnegativity of the solution and conserves the mass; a solution to the Fokker–Planck Equation with piecewise-constant drift is proposed in \cite{ChenChen},  a numerical method, named as information length, for measuring distances between statistical states as represented by PDF has been proposed in \cite{Anderson}. Also, there has been work on Fractional Fokker-Planck Equation as well, e.g., a space-time Petrov-Galerkin spectral method for time fractional FP equation with nonsmooth solution has been studied in \cite{Zeng} and a numerical solution of the Cauchy problem for the fractional FP equation in connection with Sinc convolution methods is proposed in \cite{Baumann&Stenger}.

In this work, we intend to solve the FP equation with linear and nonlinear drift vector and constant diffusion tensor. By doing this, the Gaussian distribution for the FP equation, which is a parabolic type differential equation that also satisfies the positiveness and conservation condition (\ref{normalizedcond_tfpe}). We present a two-level algorithm with coarsening by a factor-of-three strategy on staggered grids c.f. \cite{Butt,Butt&Borzi,Butt&Yuan} with (backward) time-difference scheme of order one and two, i.e., BDF1 and BDF2, respectively. A Chang-Cooper discretization scheme has been used to guarantee the second-order accuracy, positiveness and conservation of the total probability.

In the next Sec. {\ref{sec_tfpe}}, the Fokker-Planck equation is presented and a Chang-Cooper discretization scheme is explained in Sec. \ref{sec_disfpe}. A two-level scheme with inter-grid transfer operators is presented in Sec. \ref{sec_mg}. Numerical results are reported in Sec. \ref{sect_num} and a Sec. of conclusions is given in the last.
\section{Fokker-Planck equation \label{sec_tfpe}}

We consider the following time-dependent Fokker-Planck equation in one-dimensional computational domain $\Omega \subset \mathbb{R}$ with $Q:=\Omega \times (0,T)$ and Lipschitz boundary $\partial\Omega$:
\begin{equation}\label{tfpeq}
\frac{\partial u(x,t)}{\partial t} - \frac{\sigma^2}{2} \frac{\partial }{\partial x^2} u(x,t) + \frac{\partial}{\partial x} (f(x,t)\,u(x,t)) = 0, \qquad \mbox{ in } Q
\end{equation}
with the initial PDF distribution
\begin{equation}\label{inicond_tfpe}
u(x,0) = u_0(x), \qquad \mbox{ in }  \Omega
\end{equation}
which satisfies the positiveness and conservation of PDF distribution condition
\begin{equation}\label{normalizedcond_tfpe}
u_0 \geqslant 0, \qquad \int_{\Omega} u_0(x) \,dx = 1.
\end{equation}

The Fokker-Planck equation (\ref{tfpeq}) can be written in flux form (with non-zero source term $g$), i.e.,
\begin{equation}\label{fluxtfpe}
\frac{\partial u}{\partial t} - \nabla \cdot F(x,t) = \frac{\partial u}{\partial t} - \frac{\partial}{\partial x} F(x,t) = g(x,t)
\end{equation}
where
\begin{equation*}
\nabla  = \frac{\partial}{\partial x}, \qquad F(x,t) = B(x,t)u(x,t)+C(x,t)\frac{\partial}{\partial x} \, u(x,t)
\end{equation*}
represents the flux and the source term $g(x,t)$ has been added for the numerical investigation purposes. However, the positivity and conservation of the PDF distribution function $u(x,t)$ for the FP equation are claimed when $g(x,t)=0$. The initial condition is given by (\ref{inicond_tfpe}) and the boundary conditions are
\begin{equation}\label{bcstimefpe}
F = 0, \qquad \mbox{ on }  \partial \Omega \times (0,T).
\end{equation}

For simplicity, we choose $C(x,t)=a_{ii}(x,t)$, $B(x,t)=\partial x_i a_{ii} (x,t)$. Also, we assume that $C(x,t)$ is a positive continuous scalar function and in the case of Ornstein-Uhlenbeck process that we shall follow, $C(x,t)$ is a positive constant function and $B(x,t)$ is constant in time and linear in the spatial variable. Further, we assume that $B(x,t)$ is a function such that satisfies the Lipschitz continuity
\begin{equation*}
|B(x+h,t)-B(x,t)|\leq y \, h
\end{equation*}
where $y$ is the Lipschitz constant.
\section{Discretization on staggered grid \label{sec_disfpe}}

In this section, we discretize the FP equation on staggered grid, see Fig. \ref{staggeredgrid}. We use the Chang–Cooper (CC) scheme which is second-order accurate and guarantees the conservation of the total probability and positive solution to the numerical solution of FP equation, see \cite{Chang}.

We consider a one-dimensional computational domain, i.e., $\Omega = (-a,a)$. For discretization, we consider a sequence of uniform grids $\{\Omega_h\}_{h>0}$ with spatial mesh size $h$ and $N$ as the number of cells
\begin{equation*}
\Omega_h = \left\{x_i = -a+ih, i= 0,1,\ldots, N \right\} \cap \Omega
\end{equation*}
On a uniform staggered grid, the flux $F$ and PDF distribution function $u$ (solution points) are
\begin{align*}
&F_{i} = F(-a+i\,h),& &0 \leq i \leq N& \\
&u_{i} = u(-a+(i-1/2)h),& &1 \leq i \leq N.&
\end{align*}
We choose the spatial mesh size $h$ such that the boundary of the domain $\Omega$ coincide with the grid points.
\begin{figure}[h!]
\begin{center}
\setlength{\unitlength}{1cm}
\begin{picture}(13,4)

\put(0.5,1.5){\framebox(1.0,1.0){$F$}}
\put(1.5,2){\line(1,0){1}}
\put(3.5,2){\line(1,0){1}}
\put(4.5,1.5){\framebox(1.0,1.0){$F$}}
\put(5.5,2){\line(1,0){1}}
\put(7.5,2){\line(1,0){1}}
\put(8.5,1.5){\framebox(1.0,1.0){$F$}}
\put(9.5,2){\line(1,0){1}}
\put(11.5,2){\line(1,0){1}}
\put(12.5,1.5){\framebox(1.0,1.0){$F$}}

\put(3,2){\circle{1}{$u$}}
\put(7,2){\circle{1}{$u$}}
\put(11,2){\circle{1}{$u$}}
\end{picture}
\caption{Staggered grid for one-dimensional FP equation} \label{staggeredgrid}
\end{center}
\end{figure}

Note that the time-dependent FP equation is a parabolic type equation with an additional (positiveness and conservation of PDF distribution) condition (\ref{normalizedcond_tfpe2d}). For discretization, let ${\tau}$ and $N_t$ be the time stepsize and number of time steps, respectively. We define
\begin{equation}
Q_{h,{\tau}} = \left\{ (x_i,t_m): x_i = -a+(i-1/2)h, t_m = m {\tau}; i= 1,\ldots, N; 0 \leq m \leq N_t \right\}
\end{equation}
where on $Q_{h,{\tau}}$, we mean $u_{i}^m$ the value of the grid function in $\Omega_h$ at $x_i$ and time $t_m$. The Chang-Cooper discretization scheme is used c.f. \cite{Chang} for the spatial variable and for $\frac{\partial u}{\partial t}$ approximation we use first-order backward difference formula (BDF1)
\begin{equation}\label{bdf1}
  \frac{\partial u}{\partial t} \approx \frac{u_{i}^{m+1} - u_{i}^m}{{\tau}}.
\end{equation}

The discetize version of $\nabla \cdot F$ at time $t_m$, corresponding to the time-dependent FP equation, is given by
\begin{equation}\label{tfluxh}
  \nabla \cdot F =  \frac{F_{i+1/2}^m - F_{i-1/2}^m}{h}
\end{equation}
where $F_{i+1/2}^m$ denotes the flux in $x$-direction at the point $x_i$. The discretized flux zero boundary conditions (\ref{bcstimefpe}) are given by
\begin{eqnarray}\label{bcsfpeh}
F(-1/2,t_m) = 0, \qquad F(N+1/2,t_m) = 0, \qquad \forall \, m=0,1,\cdots,N_t
\end{eqnarray}
Moreover,
\begin{eqnarray*}
F_{i+1/2}^m &=& \left[ (1-\delta_{i+1/2}^m) f_{i+1/2}^m + \frac{\sigma^2}{2h} \right]u_{i+1}^{m+1}
- \left[ \frac{\sigma^2}{2h} - \delta_{i+1/2}^m f_{i+1/2}^m \right]u_{i}^{m+1} \\
\end{eqnarray*}
and
\begin{eqnarray}\label{deltatimefpeh}
\delta_{i+1/2}^m &=& \frac{1}{\omega_{i+1/2}^m} - \frac{1}{exp(\omega_{i+1/2}^m)-1}, \qquad \omega_{i+1/2}^m = \frac{2hf_{i+1/2}^m}{\sigma^2}.\nonumber
\end{eqnarray}
Thus the discrete time-dependent FP equation $\frac{\partial u}{\partial t} - \nabla \cdot F = 0$ becomes
\begin{equation}\label{timefpeh}
\frac{u_{i}^{m+1} - u_{i}^m}{{\tau}} - \frac{F_{i+1/2}^m - F_{i-1/2}^m}{h} = 0.
\end{equation}

Conservation of the FP equation follows from the discrete FP equation and for this we use the flux form of the FP equation:
\begin{lemma}\rm
The conservation property holds
\begin{equation*}
  \sum_{i=0}^{N}\,u_{i}^{m+1}=\sum_{i=0}^{N} \,u_{i}^{m}, \qquad m\geq0.
\end{equation*}
\proof
Denote $D_t \, u_i^m = \frac{u_i^{m+1}-u_i^m}{{\tau}}$, and consider the time-dependent FP equation in flux form, i.e.,
\begin{equation*}
D_t \, u_i^m (\thickapprox \frac{\partial u}{\partial t}) = \frac{1}{h}\left(F_{i+1/2}^m-F_{i-1/2}^m\right)
\end{equation*}
in equation (\ref{bdf1}) then taking sum over $i$ gives
\begin{equation*}
\sum_{i=0}^{N} (u_i^{m+1}-u_i^m) = \frac{{\tau}}{h}\left(F_{i+1/2}^m-F_{i-1/2}^m\right).
\end{equation*}
Note that, at the boundaries, we have a zero right hand side because it is the difference of fluxes. Thus,
$$\sum_{i=0}^{N}\,u_i^{m+1}=\sum_{i=0}^{N} \,u_i^{m}.$$\endproof
\end{lemma}

For stability and convergence of the CC scheme with first-order time approximation (BDF1), we denote
\begin{eqnarray*}
  D_+ u_i &=& \frac{u_{i+1}-u_i}{h} \\
  D_- u_i &=& \frac{u_i - u_{i-1}}{h} \\
  M_\delta u_i &=& (1-\delta_{i-1})u_i + \delta_{i-1}u_{i-1}
\end{eqnarray*}
We have the following CC discretization scheme to the FP equation with non-zero source term $g$ using BDF1 for the time variable, see \cite{Chang,Mohammadi}
\begin{eqnarray*}
\frac{u_i^{m+1}-u_i^m}{{\tau}} &=& \frac{1}{h}[((1-\delta_{i}^{m})B_{i+1/2}^{m} + \frac{1}{h}C_{i+1/2}^{m})
u_{i+1}^{m+1} \\
&-& (\frac{1}{h}(C_{i+1/2}^{m}+C_{i-1/2}^{m})+(1-\delta_{i-1}^{m})B_{i-1/2}^{m} -\delta_{i}^{m}B_{i+1/2}^{m})u_{i}^{m+1} \\
&+& (\frac{1}{h}C_{i-1/2}^{m}-\delta_{i}^{m}B_{i-1/2}^{m})u_{i-1}^{m+1}]+g_j^{m+1}, \qquad i=0,1,\ldots N
\end{eqnarray*}
where
\begin{eqnarray*}

F_{i+1/2}^{m} &=& B_{i+1/2}^{m} \left((1-\delta_{i}^{m})u_{i+1}^{m+1} + \delta_{i}^{m} u_i^{m+1}\right) + C_{i+1/2}^m
\left(\frac{u_i^{m+1}-u_i^m}{h}\right)\\
\delta_{i}^m &=& \frac{1}{\omega_{i}^m} - \frac{1}{exp(\omega_{i}^m)-1} \\
\omega_{i}^m &=& \frac{h \, B_{i+1/2}^m}{C_{i+1/2}^m}
\end{eqnarray*}
with zero-flux boundary conditions, i.e., $F_{-1/2}^{m} = 0, F_{N+1/2}^{m}=0$. Note that at equilibrium the numerical fluxes must be zero, $F_{i+1/2}=0$. Therefore,
\begin{equation*}
\frac{u_{i+1}^{m+1}}{u_{i}^{m+1}} = \frac{\frac{1}{h}C_{i+1/2}^{m}-\delta_i^m B_{i+1/2}^{m}}{(1-\delta_i^m)B_{i+1/2}^{m}+\frac{1}{h}C_{i+1/2}^{m}},
\end{equation*}
and if we solve $F(x_{i+1/2},t^{m+1})=0$, we have
\begin{eqnarray*}
\frac{u_{i+1}^{m+1}}{u_{i}^{m+1}} &=& exp\left(- \int_{x_i}^{x_{i+1}}\frac{B(x,t^{m+1})}{C(x,t^{m+1})} \, dx_i\right) \\
                                  &\approx& \frac{h\,B_{i+1/2}^{m}}{C_{i+1/2}^{m}}.
\end{eqnarray*}

With this setting, the discretized FP equation with source term $g(x,t)$ is given by
\begin{equation}\label{CCFPE1d}
\frac{u_i^{m+1}-u_i^m}{{\tau}} = D_+ C_{i-1/2}^m D_-u_{i}^{m+1} + D_+B_{i-1/2}^{m}M_{\delta} u_i^{m+1} + g_i^{m+1}
\end{equation}
and for positivity, stability and convergence results, see \cite{Mohammadi}.

\section{Two-level scheme \label{sec_mg}}

In this section, we illustrate the proposed two-grid algorithm with intergrid transfer operators in details. As we know that multigrid scheme uses grids that we obtained after discretization (finite difference or finite element) and such grids are usually obtained from a coarse grid, for example, by halving the coarsest grid, see \cite{Trottenberg}. As a result, we obtained a non-nested hierarchy of grids and need extra efforts to construct intergrid transfer operators. Therefore, we note that when a coarsening by a factor-of-three is used, we obtain a nested sequence of grids. This allows us to use bilinear interpolation and straight injection and hence the implementation of intergrid transfer operators becomes easier, which we explain in details as follows.

Let $\Omega_k$ denotes the nested grids or levels with mesh size
$h_k=1/3^{k-1}$, where $k = 1,\dots,L$, and $L$ denotes the finest level. In this way, we have a variable $X_{I}^{k-1}$ at the coarse grid point $I$ on $\Omega_{k-1}$ that has the same spatial location as the variable $X_i^{k}$ at the fine grid point $i$ on $\Omega_{k}$, see Fig. \ref{staggeredgrid}
\begin{itemize}
\item  $u_{I+1/2}^{k-1}$ corresponds to $u_{i+1/2}^{k}$ for $i=3I-1$.
\end{itemize}

In the following, we present the intergrid transfer operators, that is, the prolongation and restriction operators, that are to be used in the two-level algorithm.

For prolongation, we use a quadratic interpolation using Lagrange polynomial, e.g., using a quadratic polynomial $p(x)=a_0+a_1x+a_2x^2$ that satisfies $p(x_i)=u_i, \qquad i = 1,2,3$, we have the following quadratic (second-order) Lagrange Interpolation formula
\begin{equation*}
p(x) = \frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}u_1+\frac{(x-x_1)(x-x_2)}{(x_2-x_1)(x_2-x_3)}u_2+\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}u_3
\end{equation*}
Moreover, consider the space $\cU_k$ of $u^k: \Omega_k \to \R$, $k=1,\ldots,L$ such that for every two grids ${\Omega}_k$ and ${\Omega}_{k-1}$, a prolongation operator, $I^k_{k-1}:\cU_{k-1} \to \cU_{k}$ is defined which is consistent with each partition or subinterval of the discretization.
\begin{figure}[h!]
\begin{center}
\setlength{\unitlength}{1cm}
\begin{picture}(13,6)

\put(4.5,1.5){\framebox(1.0,1.0){$F$}}
\put(8.6,1.5){\framebox(1.0,1.0){$F$}}
\put(6.3,1.8){\framebox(.4,.4){$F$}}
\put(7.5,1.8){\framebox(.4,.4){$F$}}

\put(5.5,2){\line(1,0){.8}}
\put(6.7,2){\line(1,0){.8}}
\put(7.9,2){\line(1,0){.7}}

\put(5.7,2){\circle{.1}}
\put(5.9,2){\circle{.1}}
\put(6.1,2){\circle{.1}}

\put(6.9,2){\circle{.1}}
\put(7.1,2){\circle{.1}}
\put(7.3,2){\circle{.1}}

\put(8.1,2){\circle{.1}}
\put(8.3,2){\circle{.1}}
\put(8.5,2){\circle{.1}}

\put(0.5,4.5){\framebox(1.0,1.0){$F$}}
\put(1.5,5){\line(1,0){1}}
\put(3.5,5){\line(1,0){1}}
\put(4.5,4.5){\framebox(1.0,1.0){$F$}}
\put(5.5,5){\line(1,0){1}}
\put(7.5,5){\line(1,0){1}}
\put(8.5,4.5){\framebox(1.0,1.0){$F$}}
\put(9.5,5){\line(1,0){1}}
\put(11.5,5){\line(1,0){1}}
\put(12.5,4.5){\framebox(1.0,1.0){$F$}}

\put(3,5){\circle{1}{$u$}}
\put(7,5){\circle{1}{$u$}}
\put(11,5){\circle{1}{$u$}}
\end{picture}
\caption{Illustration of {\it{Straight Injection operator $I_k^{k-1}$}} on a single subintrval (partition). The the coarse grid $\Omega_{k-1}$ (upper line) after coarsening by a factor-of-three of the fine grid (bottom line) points on $\Omega_{k}$}\label{injection_stagrid}
\end{center}
\end{figure}

To transfer the residuals (resp. functions) from fine ${\Omega}_k$ to coarse grid ${\Omega}_{k-1}$
a straight injection, that is, $I^{k-1}_{k}: \cU_{k} \to \cU_{k-1}$ is used as a restriction operator. For example, for the flux function $F$ and PDF distribution $u$, we have the following correspondence, see Fig. \ref{injection_stagrid}
\begin{itemize}
\item  $F^{k-1}_{I} \leftarrow F^{k}_{i}$ for $i=3I-2$;
\item  $u_{I+1/2}^{k-1} \leftarrow u_{i+1/2}^{k}$ for $i=3I-1$.
\end{itemize}
Here we remark that we can use the half or full weighting as a restriction operator $I^{k-1}_{k}$. We use the straight injection because it gives a natural choice in a coarsening by a factor of three \cite{Butt&Borzi,Butt,Butt&Yuan} and hence an advantage of using the proposed multigrid scheme with factor-three coarsening.

The two-level algorithm to solve the Fokker-Planck equation is given by:
\begin{algorithm}\label{tFP_TGalgo}
{\bf TG$(m_1,m_2)$ for solving $A_k \, u_k =g_k$}.
\begin{enumerate}
\item{Set $u_k^0$ using initial condition and normalized condition (\ref{normalizedcond_tfpe});} \\

\item  Pre-smoothing:\\
$u_k^{m,(l)}=S_k(u_k^{m, (l-1)},g_k^m)$, $l=1,\dots,m_1$; \\

\item  Compute the residual $r_k^m=u_k^m - A_k (u_k^{m, (m_1)})$;
\item  Restrict the residual $r_{k-1}^m=I^{k-1}_k r_k^m$;
\item  Solve the coarse-grid (or error equation) problem on coarser level, i.e., solve
       $A_{k-1}(e_{k-1}^m)=r_{k-1}^m$ with $e_{k-1}^m :=0$ as an initial guess;
\item  Transfer the error (using interpolation operator), i.e., $I^k_{k-1}$: $e_k^m := I^k_{k-1} e_{k-1}^m$
\item  Coarse-grid correction step: $u_k^{m,(m_1 +1)}=u_k^{m,(m_1)}+e_k^m$;
\item  Apply normalized condition (\ref{normalizedcond_tfpe});
\item  {Post-smoothing on the fine grid: \\
$u_k^{m,(l)}=S_k(u_k^{m,(l-1)},g_k^m)$, $l=m_1+2,\dots,m_1+m_2+1$;}
\end{enumerate}
\end{algorithm}
\subsection{FP equation with second-order time difference \label{sec_fpe2dt2}}

In the following, we discuss the discretization of time-dependent FP equation with second-order difference scheme to the time derivative. In particular, we consider the one-dimensional time-dependent FP equation with second-order backward time difference formula (BDF2):
\begin{equation*}
\frac{\partial u(x,t)}{\partial t} \approx \frac{3u_{i}^{m+1} - 4u_{i}^{m} + u_{i}^{m-1}}{2{\tau}}
\end{equation*}
Then, we have the following discretized FP equation c.f. (Section 3.2 \cite{Mohammadi})
\begin{equation}\label{tfpeht2}
\frac{3u_{i}^{m+1} - 4u_{i}^{m} + u_{i}^{m-1}}{2{\tau}} = D_+C_{i-1/2}^{m}D_-u_i^{m+1}+D_+B_{i-1/2}^m\,M_\delta\,u_i^{m+1}+g_i^{m+1}.
\end{equation}
For conservation property, in case of one-dimensional FP equation with second-order time difference scheme (BDF2), we have by summing over $i$ and using the zero-flux boundary conditions:
\begin{equation*}
3u_{i}^{m+1} = 4u_{i}^{m} - u_{i}^{m-1}.
\end{equation*}
Then by induction and using $\sum_{i=0}^{N} \, u_i^1 = \sum_{i=0}^{N} \, u_i^0$ and $\sum_{i=0}^{N} \, u_i^m = \sum_{i=0}^{N} \, u_i^{m-1}$, we have
\begin{equation*}
\sum_{i=0}^{N} \, u_i^{m+1} = \sum_{i=0}^{N} \, u_i^m, \qquad m\geq1.
\end{equation*}
Analogously to FP equation with first-order backward time differencing BDF1 scheme, we have the conservation property for the FP equation with second-order time difference scheme (BDF2). For detailed proof about the numerical stability, convergence and positivity of the CC scheme with second-order time difference approximation scheme (BDF2), see \cite{Mohammadi}.
\section{Numerical experiments \label{sect_num}}

In this section, we present numerical examples to solve the Fokker-Plank equation with linear and nonlinear drift function to demostarte the efficiency and second-order accuract for the proposed two-level algorithm with BDF1 and BDF2, respectively. We use Matlab $2016$ on laptop $i7$, $1.86GHz$ with$4GB$ RAM, for the numerical simulations.
\subsection{Stationary FP equation}

First, we consider a stationary FP equation on $\Omega=[-6,6]$:
\begin{equation}\label{sec_fpe1d}
\frac{d}{dx} \left[-\frac{\sigma^2}{2} \frac{d}{dx}u(x) + f(x)u(x)\right] = 0.
\end{equation}
We take the diffusion coefficient $\sigma = 1$ and the linear drift function $f(x)=-x$ so that (for comparing the numerical and analyticla solution) we have an analytic solution given by $u_e = 1/exp(x^2)$. We employ the two-level Algorithm \ref{tFP_TGalgo} with $3-pre$ and $3-post$ smoothing (i.e., $m_1 = m_2 = 3$) steps. The solution error, at the discrtized level $k$, is presented in Table \ref{tab1_errfpe_exp1} based on the following discrete $L^1-norm$
\begin{equation*}
  \||u|\|_1 = h \sum_{i=1}^{N} |u_{i}|,
\end{equation*}
and discrete $L^2-norm$
\begin{equation*}
  \|u\|_2 = h^2 \sum_{i=1}^{N} u^2_{i}.
\end{equation*}
We stop the iterations when the difference of discrete $L^2-norm$ of errors of the new and old numerical approximation to $u$, i.e., when $\| u^{new} \|_2 - \| u^{old} \|_2 < tol=10^{-8}$. Number of two-grid cycles required to reach a desired tolerance with CPU time (seconds) are also reported in Table \ref{tab1_errfpe_exp1}.

Further, the numerical and analytical solution for this FP model is depicted in Fig. \ref{figs_fpe1d} on $N=81$. This shows that the proposed algoritm track the desired PDF, that is, for the stationary case, we have the match of numerical and analytical PDF distributions.
\begin{table}[h!]
\caption{Error history for stationary FP equation}
\begin{center}
\begin{tabular}{lcccc}
\hline
$N$ & $\||u - u_e\||_1$ & $\|u - u_e\|_2$ & \#TG & CPU\\
\hline
$27$ & $1.7766e-09 $ & $1.0398e-09$ & $10$ & $0.07$\\
$81$  & $3.3659e-10$  & $1.6667e-10$ & $09$ & $0.09$\\
\hline
\end{tabular}
\end{center}\label{tab1_errfpe_exp1}
\end{table}
\begin{figure}[ht]
\centering
\includegraphics[height=5in,width=5in]{Figs_FPE1d_u_ue_N81.eps}
\caption{Stationary FP equation: Numerical solution (.) and Analytical solution (-) on $N = 81$ mesh.}\label{figs_fpe1d}
\end{figure}

\subsection{Time-dependant FP equation with linear drift}

In this section, we show the second-order convergence of the proposed two-level scheme with BDF1 (resp. BDF2). For this, we consider a time-dependent Fokker-Planck equation (\ref{tfpeq}) given by \cite{Roy&Borzi2017} with initial and boundary conditions on $Q:=\Omega \times (0,T)$ with $\Omega=[-6,6]$ and $T =1$, i.e., we consider the FP equation (\ref{tfpeq}) within the framework of Ornstein-Uhlenbeck process \cite{Chen,Risken} with analytic solution:

The source term $g(x,t) = 1/exp(x^2+t)$ and drift
\begin{equation*}
f_1(x,t) = -x
\end{equation*}
and diffusion coefficient $\sigma = 1$ which results in exact solution
\begin{equation*}
  u_e(x,t)= 1/exp(x^2+t).
\end{equation*}
We employ the two-grid Algorithm \ref{tFP_TGalgo} with $3-pre$ and $3-post$ smoothing steps. Results for the discrete $L^1-norm$ and $L^2-norm$ of errors are reported in Table \ref{tab_timesourcefpe1d}), where $$\||u|\|_1 = h^2 \, {\tau}\sum_{m=0}^{N_t} \sum_{i=1}^{N} |u^k_{i}|$$
which we identify with $L^2_{{\tau}}(0,T;L^1)$ and discrete $L^2$-norm $\|\cdot\|_{L^2_{h,{\tau}}(Q)}$ , i.e.,
\begin{equation*}
  \|u\|_2 = {\tau} \, h^2 \sum_{i=1}^{N} |u^2_{i}|.
\end{equation*}

We take the time step size ${\tau} = (\frac{1}{81})(\frac{T}{3^L})$, where $L$ denote the fine level. Second-order accuracy $O(h^2+{\tau})$ is observed in the numerical results for the proposed two-level scheme, see Table \ref{tab_timesourcefpe1d}. In fact, we have a reduction in errors by a factor of nine (i.e., a factor $3^2$) as we refine the mesh by factor of $3$ for the $L^2-norm$ of errors.

In Table \ref{tab_time2sourcefpe1d}, we report results of discrete $L^1-norm$ and $L^2-norm$ of errors with CUP time (second) with second-order time differencing scheme (BDF2) using the proposed two-level scheme to the same numerical example. We use the two-level scheme at $t=t_1$ with first-order time difference scheme $BDF1$. Second-order accuracy $O(h^2+{\tau})$ is achieved , that is, we have a reduction in errors by a factor of $3^2$ as we refine the mesh by factor of $3$ for the $L^2-norm$ of errors.
\begin{table}[h!]
\caption{Convergence of FP equation with BDF1}
\begin{center}
\begin{tabular}{lcccc}
\hline
$N \times N_t$ & $\||u - u_e\||_1$ & $\|u - u_e\|_2$ & $\#TG$ & $CPU$\\
\hline
$81 \times 81$ & $6.2447e-05$ & $1.9392e-06$ & 4 & 0.14 \\
$243 \times 243$& $2.3321e-05$ & $2.4187e-07$ & 2 & 0.22 \\
$729 \times 729$ & $4.1551e-06$ & $1.6076e-08$ & 2 & 0.84 \\
$2187 \times 2187$ & $1.2060e-06$ & $1.4322e-09$ & 2 & 6.13 \\
\hline
\end{tabular}
\end{center}\label{tab_timesourcefpe1d}
\end{table}

\begin{table}[h!]
\caption{Convergence of FP equation with BDF2}
\begin{center}
\begin{tabular}{lcccc}
\hline
$N \times N_t$ & $\||u - u_e\||_1$ & $\|u - u_e\|_2$ & $\#TG$ & $CPU$\\
\hline
$81 \times 81$ & $6.2450e-05$ & $1.9393e-06$ & 4 & 0.08 \\
$243 \times 243$& $2.3330e-05$ & $2.4195e-07$ & 2 & 0.33 \\
$729 \times 729$ & $4.1450e-06$ & $1.6049e-08$ & 2 & 0.89 \\
$2187 \times 2187$ & $1.2077e-06$ & $1.4330e-09$ & 2 & 6.16 \\
\hline
\end{tabular}
\end{center}\label{tab_time2sourcefpe1d}
\end{table}

Next, to have a comparison of the proposed two-level scheme with the Chang-Cooper with first-order time backward difference (CC-BDF1) scheme given by Mohammadi and Borzi \cite{Mohammadi}, we consider the following (Ornstein-Uhlenbeck process) FP equation in $Q:=\Omega \times [0,T]$ and take $B(x,t)=x$, $C(x,t) = \sigma^2$:
\begin{equation}
\frac{\partial u(x,t)}{\partial t} = \partial_x \, \left(B(x,t) \, u(x,t) + C(x,t)\, \partial_x \, u(x,t)\right) + g(x,t), \qquad \mbox{ in } Q
\end{equation}
where the source term is given by
\begin{equation*}
g(x,t)=\frac{(a-x)(2x-a)}{exp((x-a/2)^2+t)}
\end{equation*}
and the drift function $f(x,t) = x$. In partcular, we choose $\Omega = [0,a]$ with $a= 10$ and $\sigma = 1$, $T=1$. Furthermore, the initial condition is given by
\begin{equation*}
u(x,0) = \frac{1}{exp((x-a/2)^2)}
\end{equation*}
with flux zero boundary and the exact solution is given by
\begin{equation*}
u_e(x,t)= \frac{1}{exp((x-a/2)^2+t)}.
\end{equation*}
We employ the two-grid Algorithm \ref{tFP_TGalgo} with $3-pre$ and $3-post$ smoothing steps. Results for the discrete $L^2-norm$ of errors are reported in Table \ref{tab_tfpe1d_expMohBor}. We take the same time step size ${\tau} = 0.01(\frac{1}{3^L})^2$, where $L$ is the fine level in the two-grid Algorithm. Second-order accuracy $O(h^2+{\tau})$ is observed in the numerical results for the proposed two-level scheme, see Table \ref{tab_tfpe1d_expMohBor}. Moreover, we present the numerical results of Chang-Cooper with first-order time difference (BDF1) scheme given by Mohammadi and Borzi \cite{Mohammadi}, in Table \ref{tab_tfpe1d_ComExpMohBor}. From Table \ref{tab_tfpe1d_expMohBor}-\ref{tab_tfpe1d_ComExpMohBor}, we clearly see that our proposed scheme gives better accuracy as compared to the numerical results given by \cite{Mohammadi}. Relative discrete $L^2_h$-norm of errors, on $N=81,N=243,N=729$ are recorded as $1.2253e-1, 8.2669e-2$ and $6.3865e-2$, respectively. In Fig \ref{figs_fpe1dMoh}, we depict the numerical and analytic solution on $N=243$ at $T=1$ with BDF1 and BDF2 to showcase the accuracy of the proposed Two-Level Algorithm.
\begin{table}[h!]
\caption{Convergence of proposed two-level with BDF1 for numerical example c.f. \cite{Mohammadi}}
\begin{center}
\begin{tabular}{lcccc}
\hline
$N \times N_t$  & $\|u - u_e\|_2$ & $\#TG$ & $CPU$\\
\hline
$81 \times 81$ &  $4.6050e-09$ & 2 & 0.16 \\
$243 \times 243$& $1.1507e-10$ & 2 & 0.27 \\
$729 \times 729$& $3.2926e-12$ & 2 & 0.96 \\
\hline
\end{tabular}
\end{center}\label{tab_tfpe1d_expMohBor}
\end{table}

\begin{table}[h!]
\caption{Convergence of proposed two-level with BDF2 for numerical example c.f. \cite{Mohammadi}}
\begin{center}
\begin{tabular}{lcccc}
\hline
$N \times N_t$  & $\|u - u_e\|_2$ & $\#TG$ & $CPU$\\
\hline
$81 \times 81$ &  $4.6051e-09$ & 2 & 0.11 \\
$243 \times 243$& $1.1508e-10$ & 2 & 0.31 \\
$729 \times 729$& $3.2930e-12$ & 2 & 0.91 \\
\hline
\end{tabular}
\end{center}\label{tab_t2fpe1d_expMohBor}
\end{table}

\begin{table}[h!]
\caption{CC-BDF1 scheme for FP equation by Mohammadi and Borzi \cite{Mohammadi}}
\begin{center}
\begin{tabular}{lccc}
\hline
$N \times N_t$   & $\|u - u_e\|_2$\\
\hline
$50 \times 50$   & $1.34e-2$\\
$100 \times 100$ & $3.50e-3$\\
$200 \times 800$ & $8.80e-4$\\
\hline
\end{tabular}
\end{center}\label{tab_tfpe1d_ComExpMohBor}
\end{table}

\begin{figure}[ht]
\centering
\includegraphics[height=2.5in,width=2.5in]{Figs_FPE1d1t_u_ue_N243_expMoh.eps}
\includegraphics[height=2.5in,width=2.5in]{Figs_FPE1d2t_u_ue_N243_expMoh.eps}
\caption{Time-dependent FP equation \cite{Mohammadi}: Numerical solution (dashes) and Analytical solution (solid line) on $N = 243$ with BDF1 (left) and BDF2 (right), respectively.}\label{figs_fpe1dMoh}
\end{figure}
\subsection{FP equation with nonlinear drift \label{exp_tnonlinearfpe}}

In the following, we consider the nonlinear process given by Harrison \cite{Harrison}
\begin{equation}\label{nonlinearStochastic}
dX = (X-X^3)\,dt + \sigma \, dW
\end{equation}
with the corresponding FP equation
\begin{equation}\label{tnonlfpeq}
\frac{\partial u(x,t)}{\partial t} - \frac{\sigma^2}{2} \partial^2_{xx} u(x,t) + \partial_{x} (f(x,t)\,u(x,t)) = 0, \qquad \mbox{ in } Q \end{equation}
with the initial PDF distribution (\ref{inicond_tfpe}) and drift
\begin{equation*}
f(x) = x-x^3.
\end{equation*}
Moreover, we take the diffusion coefficient $\sigma = 0.4$. An analytic solution to (\ref{tnonlfpeq}) is not known. However, the steady state solution is given by
\begin{equation}\label{anasol}
u(x) = C exp((x^2-0.5x^4)/\sigma^2)
\end{equation}
where $C$ is the normalized constant. The numerical solution to this FP equation with the nonlinear drift function is depicted for $T = 0.5, 1.0, 3.0, 5.0, 15.0$ and $T=30$, respectively, on $N=81$ mesh. Moreover, to have a comparison with the results given by \cite{Harrison}, where a numerical solution of the FP equation using {\it{moving finite elements}} is presented, see Fig. \ref{figs_nonlineartfpe1d_N81}. As given by \cite{Harrison}, the deterministic equation $dx/dt=x-x^3$ has two asymptotically stable equilibria at $x = 1$ and $x = -1$ which can been seen in Fig. \ref{figs_nonlineartfpe1d_N81}. In addition, we have a symmetric bimodal distribution as a result of our numerical two-level scheme which is presented in \cite{Harrison}.
\begin{figure}[ht]
\centering
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T0_N81.eps}
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T1_N81.eps} \\
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T3_N81.eps}
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T5_N81.eps} \\
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T15_N81.eps}
\includegraphics[height=2.5in,width=2.5in]{Figs_nonlineartFPE1d_u_T30_N81.eps} \\
\caption{Nonlinear FP equation: Numerical solution when T = 0.5 (top left); T = 1 (top right); T = 3 (middle left); T = 5 (middle right);T = 15 (bottom left); T = 30 (bottom right); on $N = 81 = N_t$ mesh, i.e., with $h=1/81$, and  ${\tau}=(1/81)(T/81)$.}\label{figs_nonlineartfpe1d_N81}
\end{figure}
\section{Conclusions \label{sect_con}}

A two-level scheme with coarsening by a factor-of-three strategy was proposed to solve the Fokker-Planck equation with linear (nonlinear) drift function. The Chang-Cooper scheme was used to discretize the FP equation on staggered grids. Second-order accuracy, that is, $O(h^2+\tau)$ and $O(h^2+\tau^2)$ , was achieved in the numerical results using second-order differences for the spatial variable and first-order (resp. second-order) time differences BDF1 (resp. BDF2). Results of numerical examples outperform the existing numerical works on FP equation (in particular to the Ornstein-Uhlenbeck process given by \cite{Mohammadi} and to the nonlinear FP equation \cite{Harrison}). A natural extension to two-dimensional FP equation is under investigation and more complicated FP equations, that is, FP equations with nonlinear drift functions that depend on space and time both, is also our future work.


text:
\begin{document}

\bibliographystyle{amsplain}

\title{Midpoint sets contained in the unit sphere of a normed space}
\author{Konrad J.\ Swanepoel}
\address{Department of Mathematics,
London School of Economics and Political Science,
Houghton Street,
London WC2A 2AE,
United Kingdom}
\email{k.swanepoel@lse.ac.uk}
\subjclass[2000]{Primary 52A15, Secondary 52A20}
\keywords{Minkowski space, midpoint set, convex position}

\begin{abstract}
The midpoint set $M(S)$ of a set $S$ of points is the set of all midpoints of pairs of points in $S$.
We study the largest cardinality of a midpoint set $M(S)$ in a finite-dimensional normed space, such that $M(S)$ is contained in the unit sphere, and $S$ is outside the closed unit ball.
We show in three dimensions that this maximum (if it exists) is determined by the facial structure of the unit ball.
In higher dimensions no such relationship exists.
We also determine the maximum for euclidean and sup norm spaces.
\end{abstract}

\maketitle

\section{Introduction and main results}

For any set of points $S$ in a vector space, the {\em midpoint set} of $S$ is
\[ M(S) = \{ \half (x+y) : x,y\in S, x\neq y\}.\]
We study sets such that all their midpoints are unit vectors in a finite-dimensional normed space.
This work is contained in the author's PhD thesis \cite{Swanepoel1997}.
Since recently some work appeared \cite{Halman2007, Eisenbrand2008} dealing with related questions, the author decided to write this up as a paper.

First note that if the space is not strictly convex, there exists an infinite set whose midpoint set is a set of unit vectors:
Choose distinct points on a line segment on the boundary of the unit ball.
If however $X$ is strictly convex, it follows from the strict triangle inequality that there is at most one point in $S$ of norm $\leq 1$.

We therefore only look at the cases where all the vectors in $S$ are of norm larger than one.
We say that a set of vectors $S\subseteq X$ is an {\em M-set} if $M(S)$ consists of unit vectors, and each vector in $S$ is of norm $>1$.
Let $m(X)$ be the largest $m$ such that there exists an M-set of cardinality $m$ in $X$, if there is such a largest $m$.
If there is no such largest $m$, we set $m(X)=\infty$.

In two dimensions, $m(X)$ does not distinguish between different norms.

\begin{theorem}\label{helly:th2}
We have $m(X)=3$ for any two-dimensional Minkowski space $X$.
\end{theorem}

In dimensions higher than two, $m(X)=\infty$ is possible.

\begin{theorem}\label{infvectors}
If the unit ball of a Minkowski space has a non-polytopal proper face, then there exists an infinite M-set.
\end{theorem}

However, it is easily seen that an M-set is always countable.

\begin{observation}
Any M-set in a Minkowski space is countable.
\end{observation}
\begin{proof}
It is sufficient to show that a set $S$ of vectors of norm $\geq 1+\epsilon$ satisfying $\norm{x+y}=2$ for distinct $x,y\in S$, is finite (for any $\epsilon>0$).
But such a set is obviously bounded, and for any two vectors $x,y\in S$ we have
$\norm{x-y}+2 = \norm{x-y}+\norm{x+y}\geq \norm{2x}\geq 2+2\epsilon$.
\end{proof}

The following theorem provides a converse for Theorem~\ref{infvectors} in the case of $3$-dimensional spaces.
However, we suspect that there is a $4$-di\-men\-sion\-al strictly convex, smooth space admitting an infinite M-set; see Theorem~\ref{ex4d} for a weaker result.

\begin{theorem}\label{polyth}
Let the unit ball of a $3$-dimensional Minkowski space $X$ have $2$-faces, all of which are polygonal.
Then any M-set is finite.
Furthermore,
\begin{equation*}
m(X) = \sup\{n+1:\text{ The unit ball has an $n$-gonal $2$-face}\}
\end{equation*}
\end{theorem}

The only possibility for three-dimensional spaces not yet covered, is where the unit ball does not have a $2$-face.

\begin{theorem}\label{threelower}
If the unit ball of a $3$-dimensional Minkowski space $X$ does not have a $2$-face, then $m(X)=4$.
\end{theorem}

\begin{conjecture}
 There exists a strictly convex, smooth $4$-dimensional Minkowski space admitting an infinite M-set.
\end{conjecture}

\begin{theorem}\label{ex4d}
For all $n\geq 1$ there exists a $4$-dimensional smooth, strictly convex Minkowski space $X$ such that $m(X)\geq n$.
\end{theorem}

We finally calculate $m(X)$ for euclidean and $\ell_\infty$ spaces of all dimensions.

\begin{theorem}\label{l2m}
For all $d\geq 1$, $m(\ell_2^d)=d+1$.
\end{theorem}

\begin{theorem}\label{inftybound}
For all $d\geq 2$, $m(\ell_\infty^d)=2d-1$.
\end{theorem}

In the next section we esablish some notation and prove some basic properties of M-sets.
In Section~\ref{proofsection} the above theorems are proved.
In the last section we discuss some open questions.

\section{Notation and Basic properties}
We denote the convex hull, affine hull, boundary, interior, and relative interior of a set $S$ in a vector space by $\convex S, \aff S, \bd S, \interior S$ and $\relint S$, respectively.
By {\em polytope} we always mean a convex polytope.
We say that a polytope $P$ is $2$-neighbourly if for any two vertices $x$ and $y$, the segment $[x,y]:=\convex\{x,y\}$ is an edge of $P$, and we say that $P$ is $2$-almost-neighbourly if for any two vertices $x$ and $y$, the segment $[x,y]$ is contained in the boundary of $P$.

We frequently use the following fact:

\begin{lemma}\label{normface}
Let $S$ be a set of unit vectors in a Minkowski space.
If there exists a unit vector $y\in\relint\convex S$, then $\convex S$ is a set of unit vectors.
\end{lemma}
\begin{proof}
By Carath\'eodory's theorem we may assume that $S$ is finite: $S=\{x_1,\dots,x_m\}$.
For any $z\in\convex S$ we have $\norm{z}\leq 1$, since the unit ball is convex.
Alternatively, if $z=\sum_i\lambda_i x_i$ for some $\lambda\geq 0$ such that $\sum_i\lambda_i=1$, then $\norm{z}=\norm{\sum_i\lambda_i x_i}\leq \sum_i\lambda_i\norm{x_i}=1$.

We may assume $z\neq y$.
Then there exist unit vectors $a,b$ on the line $\aff\{y,z\}$ such that $y$ and $z$ are between $a$ and $b$.
Thus there exists $0<\lambda\leq 1$ such that $y=\lambda a + (1-\lambda)z$, and
$1=\norm{y}\leq\lambda\norm{a}+(1-\lambda)\norm{z}=\lambda+(1-\lambda)\norm{z}$, implying $\norm{z}\geq 1$.
\end{proof}

We need the following sharpening of Carath\'eodory's Theorem:

\begin{lemma}\label{caramod}
Let $R\subseteq\R^d$ and $y\in\relint\convex R$.
Then for any $x\in R$ there is an affinely independent set $S\subseteq R$ of at most $d+1$ points such that $x\in S$ and $y\in \relint\convex S$.
\end{lemma}
\begin{proof}
We may assume without loss of generality that $R$ is finite, and that $\convex R$ is $d$-dimensional.
If $x=y$, there is nothing to prove.
Otherwise, let $y'$ be the projection of $y$ along the ray from $x$ through $y$ onto a facet $F$ of $\convex X$.
By Carath\'eodory's theorem, $y'$ is in the relative interior of the convex hull of at most $d$ affinely independent vertices of $F$.
Add the point $x$ to obtain the set $S$.
\end{proof}

\begin{lemma}\label{polytope}
Let $M$ be a finite M-set in a Minkowski space.
Then $M$ is the vertex set of a polytope $P$.

Suppose furthermore that for some two points $x,y\in M$, the open segment $(x,y)$ is contained in the relative interior of some $k$-face $Q$ of $P$, with $k\geq 2$ (thus $[x,y]$ is not an edge).
Then the unit ball has a $k$-face in $\aff Q$.
\end{lemma}

\begin{proof}
We first show that no point is in the convex hull of the remaining points.
Suppose that $x$ is in the convex hull of the remaining vectors.
We may suppose $x = \sum_{i=1}^n \lambda_i x_i$, with $x_i\in M, \lambda_i\geq 0$ and $\sum_{i=1}^n\lambda_i = 1$.
Then
\begin{align*}
 \norm{2x} = & \norm{\sum_{i=1}^n \lambda_i (x+x_i)} \\
		&\leq \sum_{i=1}^n\lambda_i\norm{x+x_i} \\
		&\leq 2,
\end{align*}
and $\norm{x}\leq 1$, a contradiction.

Secondly, let $x$ and $y$ be two vertices of $P$ such that $(x,y)\subseteq\relint Q$.
Note that $x,y\in Q$.
Let $R$ be the set of vertices $z$ of $Q$ such that $[x,z]$ is an edge of $Q$.
Thus $x,y\notin R$.
Let $z\in R$.
By considering the vertex figure of $Q$ at $x$ and Lemma~\ref{caramod}, there is a simplex with vertex set $S$ satisfying $z\in S\subseteq R$, such that $y-x$ is a strictly positive linear combination of $S-x$:
$$y= x+ \sum_{s\in S}\lambda_s (s-x), \qquad \lambda_s>0.$$
Note that $S$ contains at least $2$ and at most $k$ points.
We must have $\lambda:= \sum_s\lambda_s>1$, otherwise $y$ is in the convex hull of $\{x\}\cup S$, contradicting the first part of the proof.
But now $\half (x+y)$ is in the relative interior of the polytope
$$T_z:=\convex\left\{\half(x+s), \half(y+s):s\in S\right\},$$
as can be seen by taking $\mu_s := \lambda_s/\lambda^2$ and $\nu_s := \lambda_s(1-1/\lambda)/\lambda$; then
$$\half (x+y) = \sum_{s\in S}\mu_s \half(x+s) + \sum_{s\in S}\nu_s \half(y+s),\quad\sum_{s\in S}(\lambda_s+\nu_s)=1.$$
It follows that $T_z$ is in the boundary of the unit ball by Lemma~\ref{normface}.
However, it is possible that $T_z$ is not $k$-dimensional, although it is contained in $\aff Q$.

We also have $\half(x+y)\in\relint R'$, where
$$R':=\convex\left\{\half(x+z), \half(y+z):z\in R\right\}$$
is a $k$-polytope contained in $\aff Q$:
if not, then $\half(x+y)$ is contained in some facet $F$ of $R'$.
But choose a $z\in R$ such that $\half(x+z)$ or $\half(y+z)$ is not in $F$.
Then, since $\half(x+y)\in\relint T_z$, we obtain a contradiction.
\end{proof}

\section{Proofs}\label{proofsection}

\begin{proof}[Proof of Theorem~\ref{helly:th2}]
By Lemma~\ref{polytope} we immediately get $m(X)\leq 3$.
For the upper bound, find unit vectors $x,y,z$ such that $x+y+z=0$, by inscribing an affinely regular hexagon in the unit disc.
Then $\{2x,2y,2z\}$ is an M-set.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{infvectors}]
Note that the boundary of a non-polygonal $2$-face $F$ contains either infinitely many straight line segments or a strictly convex arc.
In both cases it is easy to see how to choose infinitely many points in the plane of $F$, but outside $F$, such that the midpoint of any two of the points is in $F$; see Figure~\ref{fig4}.
\begin{figure}
\includegraphics[bb = 0 0 91 59]{midpoint1}
\begin{picture}(0,0)(91,0)
\put(35,25){$F$}
\put(-12,5){$x_1$}
\put(10,56){$x_2$}
\put(61,58){$x_3$}
\put(84,33){$x_4$}
\end{picture}
\caption{Choosing points outside $F$ with midpoint inside $F$}\label{fig4}
\end{figure}
For faces of larger dimension, take an appropriate $2$-dimensional cut to obtain a non-polygonal convex disc $F$.
\end{proof}

\begin{lemma}\label{threelow}
For any Minkowski space of dimension at least three, $m(X)\geq 4$.
\end{lemma}
\begin{proof}
We may assume without loss of generality that the space is three-dimensional.
Consider any plane $H$ through $0$.
Let $x$ be a unit vector such that $x+H$ supports the unit ball at $x$.
For any $\lambda\in\R$, let $H_\lambda = H+\lambda x$.
Now fix $1/3 < \lambda < 1$.
The intersection of the boundary of the unit ball with $H_\lambda$ is a convex curve $C$.
It is possible to inscribe an affinely regular hexagon in $C$ (see e.g.\ \cite[p.\ 242]{Gr3}) with centre $y\in H_\lambda$.
Let the vertices be $\{y\pm a, y\pm b, y\pm c\}$, where $a,b,c\in H$ satisfy $a+b+c=0$.
Then it is easily seen that $\{y+2a, y+2b, y+2c, -3y\}$ is an M-set.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{threelower}]
Let $S\subseteq X$ be an M-set.
By Lemma~\ref{polytope}, $S$ is the vertex set of a polytope of dimension $2$ or $3$, and $2$-neighbourly, implying $m\leq 4$, or of dimension $1$, implying $m\leq 2$.

The lower bound follows from Lemma~\ref{threelow}.
\end{proof}

\begin{lemma}\label{polylower}
If the unit ball of a $3$-dimensional Minkowski space $X$ has an $n$-gonal $2$-face, then $m(X)\geq n+1$.
\end{lemma}
\begin{proof}
Given the polygon $P:= \convex\{a_1,a_2,\dots, a_n\}$ (with vertices in this order) as a face of the unit ball, we construct $n+1$ points outside the unit ball such that the midpoint of any two is in $P$ or $-P$.

As an intermediate step, let $x_i:=\half(a_i+a_{i+1})$ for $i=1,\dots,n$.
Then the midpoint of any two $x_i$'s is obviously in $\relint P$.
We now want to find a point $x_{n+1}$ such that the midpoint of $x_i$ and $x_{n+1}$ is in $\relint -P$ for all $i=1,\dots,n$.
(We need to have the midpoints in the {\em interior} of the reflected polygon in order to later modify $x_i,i=1,\dots,n$ --- at the moment they are still on the boundary of the unit ball.)
See Figure~\ref{proj}.
\begin{figure}
\includegraphics[bb = 0 0 109 150]{helly-type-proj}
\begin{picture}(0,0)(109,0)
\put(23,137){$x_{n+1}$}
\put(40,48){$0$}
\put(-17,17){$P$}
\put(-22,82){$-P$}
\end{picture}
\caption{}\label{proj}
\end{figure}
Note that such an $x_{n+1}$ would be outside the unit ball (in fact $\norm{x_{n+1}}=3$).
For any $x$, the set $\{\half(x+x_i):i=1,\dots,n\}$ is homothetic to $\{x_1,\dots,x_n\}$, with ratio $\half$.
It is therefore sufficient to show that there is a $-\half$-copy of $\{x_1,\dots,x_n\}$ in $\relint P$.
Let $P_i$ be the homothet of $\relint P$ with centre $x_i$ and factor $\tfrac{2}{3}$.
It is easily seen that $c\in P_i$ iff $c$ is a centre of homothety with factor $-\half$ taking $x_i$ into $\relint P$.
We want to find such a $c$ working for all $i$; it is therefore sufficient to show that $\bigcap_{i=1}^n P_i\neq\emptyset$.
By Helly's theorem it is sufficient to show that any three $P_i$'s have non-empty intersection.
This is equivalent to showing that there is a $-\half$-copy of $\{x_1,\dots,x_n\}$ in $\relint P$ in the case $n=3$.
But there is already a $-\half$-copy of $\{x_1,x_2,x_3\}$ in $(\convex\{x_1,x_2,x_3\})\setminus\{x_1,x_2,x_3\}$: for $c$ take the centroid of $\{x_1,x_2,x_3\}$.
(Note that $(\convex\{x_1,x_2,x_3\})\setminus\{x_1,x_2,x_3\}\subseteq \relint P$.)

We now modify $x_1,\dots,x_n$ by slightly shifting each point in the plane $\aff P$ out of $P$ so that $\norm{x_i}>1$.
As the construction was made such that all the midpoints were either in $\relint P$ or $\relint -P$, the midpoints will now still be in $P$ or $-P$.
\end{proof}

The following lemma was first proved by Zamfirescu \cite{Zamfirescu1969}.
For convenience we include a proof.

\begin{lemma}\label{threepolytope}
Let $P$ be a $2$-almost-neighbourly $3$-polytope.
Then $P$ is either combinatorially equivalent to a triangular prism, or $P$ is a polygonal pyramid. \textup{(}See Figure ~\ref{fig1}.\textup{)}
\end{lemma}

\begin{figure}
\includegraphics[bb = 0 0 299 124]{helly-type-prism}
\caption{Combinatorial triangular prism. Polygonal pyramid.}\label{fig1}
\end{figure}

\begin{proof}
Suppose that two faces, an $n$-gon $F$ and an $n'$-gon $G$ , are joined by an edge, with $n,n'\geq 4$.
\begin{figure}
\includegraphics[bb = 0 0 166 106]{helly-type-faces}
\begin{picture}(0,0)(166,7)
\put(25,103){$a$}
\put(130,80){$b$}
\put(163,16){$c$}
\put(23,53){$d$}
\put(112,5){$e$}
\put(-14,28){$f$}
\put(145,50){$G$}
\put(40,11){$F$}
\end{picture}
\caption{}\label{fig2}
\end{figure}
Let the vertices of the joining edge be $a$ and $b$. Let the vertices of $F$ in the neighbourhood of $a$ and $b$ be $f,a,b,e$, and those of $G$ be $d,a,b,c$.
See Figure~\ref{fig2}.
Then $[d,e]$ and $[c,f]$ are on the boundary of $P$ and must be contained in the same face of $P$.
Thus, $\convex\{d,c,e,f\}$ is part of this face.
It follows that $n=n'=4$.
If $\convex\{b,c,e\}$ is not a face of $P$, then $\aff\{b,c,e\}$ strictly separates $a$ from some vertex $x$ of $P$.
Since $x$ cannot be part of $F$ or $G$, it follows that $[a,x]\nsubseteq\bd P$, a contradiction.
Thus, $\convex\{b,c,e\}$ is a face of $P$, and similarly, $\convex\{a,d,f\}$ is a face.
Thus $P$ is a triangular prism.

If $P$ is not combinatorially equivalent to a triangular prism it follows that if two faces meet at an edge, one must be a triangle.
Let $x_1x_2\dots x_n$ be a face with the largest number of vertices $n$.
Then all faces adjacent to $F$ are triangles.
Suppose there are two points $y_1,y_2$ not in the plane $\aff F$.
By Radon's theorem applied to $S:=\{x_1,x_2,x_3,y_1,y_2\}$, there is a partition of $S$ into sets $T,U$ such that $\convex T\cap\convex U\neq\emptyset$.
Since $S$ is the vertex set of a polytope, $T$ and $U$ both contain at least two points.
We may assume that $\#T=2$ points, and $\#U=3$.
It is clear that $T=\{x_i,y_j\}$ for some $i,j$, say $T=\{x_1,y_1\}$.
Thus $[x_1,y_1]$ intersects $\convex U$.
Since $[x_1,y_1]$ is on the boundary of $\convex S$, it cannot intersect $\relint\convex U$.
Therefore, it must intersect the relative interior of some edge of $\convex U$, and we obtain a face of at least $4$ vertices adjacent to $F$, a contradiction.
Thus there is exactly one point not in $\aff F$, and $P$ is a pyramid.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{polyth}]
The lower bound in the equation is proved in Lemma~\ref{polylower}.
By Lemmas~\ref{polytope} and \ref{threepolytope}, any finite M-set is the vertex set of a planar polygon, a combinatorial triangular prism, or a polygonal pyramid.

It is easily seen that a triangular prism is impossible:
Considering the midpoints of a quadrilateral face $F$, we see that the unit ball has a face in the plane of $F$.
Since this holds for all three quadrilateral faces, we contradict the central symmetry of the unit ball.

We may assume that $m\geq 5$, since otherwise
$$m(X)\leq 4\leq \sup\{n+1: \text{there is an $n$-gon}\}.$$
The only remaining possibility is that the convex hull of the M-set is a polygon of $\geq 5$ vertices or a polygonal pyramid with $\geq 4$ vertices.
In both cases, there is an M-set with a polygon $P$ as convex hull.
Thus, by Lemma~\ref{polytope} there is a $2$-face $Q$ on the boundary of the unit ball in the plane of $P$.
By hypothesis, $Q$ is a polygon.
Choose any point in $\relint Q$ and triangulate $\aff Q$ into unbounded regions by joining this point to the vertices of $Q$, as in Figure~\ref{poly}.
\begin{figure}
\includegraphics[bb = 0 0 95 89]{helly-type-poly}
\caption{}\label{poly}
\end{figure}
Obviously, no two points of the M-set can be in the same region.
Thus, $P$ has at most $n$ vertices, and $m\leq n+1$.

It remains to show that an infinite M-set is impossible.
If there is such a set, then, as before, any finite subset must be a polygon or a polygonal pyramid.
It is easily seen that the polygons must all lie in the same plane.
But in this plane the unit ball has some $n$-gon as a $2$-face, and therefore the number of points in this plane is bounded above by $n$, a contradiction.
\end{proof}

\begin{lemma}\label{smoothing}
Let $S$ be the vertex set of a centrally symmetric polytope in $\R^d$.
Then there is a smooth, strictly convex norm $\norm{\cdot}$ on $\R^d$
such that $\norm{x}=1$ for all $x\in S$.\qed
\end{lemma}

\begin{proof}[Proof of Theorem~\ref{ex4d}]
We use Carath\'eodory's moment curve:
For $t\in\R$, define $\phi(t):=(\cos t,\sin t,\cos 2t,\sin 2t)\in\R^4$.
Choose distinct $t_1,\dots,t_n\in[0,\pi/4]$, and let $x_i:=\phi(t_i)$.
The M-set will be $\{x_1,\dots,x_n\}$, which is incidentally the vertex set of a cyclic polytope, which is $2$-neighbourly; see \cite{Gale}.
Let $K:=\convex\{\pm(x_i+x_j):1\leq i<j\leq n\}$.
We now show that for each $i\neq j$, $\pm(x_i+x_j)$ is a vertex of $K$.

Define $p_{ij}(t):=\bigl(1-\cos(t-t_i)\bigr)\bigl(1-\cos(t-t_j)\bigr)$.
Then $p_{ij}(t)\geq 0$ with equality iff $t=t_i$ or $t=t_j$, and also $p_{ij}(t)<1/10$ for all $t\in[0,\pi/4]$.
By expanding $p_{ij}$, it is found that $p_{ij}(t) = \ipr{y_{ij}}{\phi(t)}+c_{ij}$, where
$$y_{ij}:= \left(-\cos t_i-\cos t_j,-\sin t_i-\sin t_j,\half\cos(t_i+t_j),\half\sin(t_i+t_j)\right)$$
and $c_{ij}:= 1+\half\cos(t_i-t_j)$.

We now show that $\{x\in\R^4:\ipr{y_{ij}}{x}+2c_{ij}=0\}$ is a support hyperplane of $K$ containing only $x_i+x_j$.
This follows from the following calculations:
\begin{align*}
\ipr{y_{ij}}{x_i+x_j}+2c_{ij} &= p_{ij}(t_i)+p_{ij}(t_j) = 0,\\
\ipr{y_{ij}}{-x_i-x_j}+2c_{ij} &= -p_{ij}(t_i)-p_{ij}(t_j)+4c_{ij}\\
&= 4c_{ij} > 0,\\
\ipr{y_{ij}}{x_i+x_k}+2c_{ij} &= p_{ij}(t_i)+p_{ij}(t_k)=p_{ij}(t_k)>0,\\
\ipr{y_{ij}}{-x_i-x_k}+2c_{ij} &= -p_{ij}(t_i)-p_{ij}(t_k)+4c_{ij}\\
&= -p_{ij}(t_k)+4+2\cos(t_i+t_j)\\
&\geq -\tfrac{1}{10}+4+\sqrt{2} > 0,\\
\ipr{y_{ij}}{x_k+x_l}+2c_{ij} &= p_{ij}(t_k)+p_{ij}(t_l)>0,\\
\ipr{y_{ij}}{-x_k-x_l}+2c_{ij} &= p_{ij}(t_k)-p_{ij}(t_l)+4c_{ij}\\
&\geq -\tfrac{1}{10}-\tfrac{1}{10}+4+\sqrt{2} > 0,
\end{align*}
for distinct $i,j,k,l$.
By Lemma~\ref{smoothing} there is a strictly convex, smooth norm on $\R^4$ such that $\norm{x_i+x_j}=2$ for all distinct $i,j$.
By strict convexity, at most one of the $x_i$'s can have norm $\leq 1$.
Note that $\dim K=4$.
This follows from the fact that any $5$ of the $x_i$'s are affinely independent, since $\{x_1,\dots,x_n\}$ is the vertex set of a polytope, as noted previously.
\end{proof}

\begin{theorem}\label{l2bound}
Let $x_1,\dots,x_m\in\ell_2^d$ be distinct vectors satisfying
\begin{equation}\label{normeq}
\norm{x_i+x_j}=1 \text{ for all distinct }i,j.
\end{equation}
Then $m\leq d+1$ if $d\neq 2$, and $m\leq 4$ if $d=2$.
Both bounds are sharp.
\end{theorem}
\begin{proof}
Note that by the strict triangle inequality, there is at most one $x_i$ of norm $\leq\half$.
Thus for $d=2$, Theorem~\ref{helly:th2} takes care of the upper bound, and for the lower bound, take $x_1,x_2,x_3$ to be unit vectors at $120^\circ$ angles, and $x_4=0$.
Then \eqref{normeq} is obviously satisfied.

We now assume that $d\neq 2$, and for the sake of contradiction that $m=d+2$.
Then $x_1,\dots,x_m$ are affinely dependent, and there exist $\lambda_1,\dots,\lambda_m\in\R$ such that
$$\sum_{i=1}^m\lambda_i x_i = 0, \sum_{i=1}^m\lambda_i = 0, \quad(\lambda_1,\dots,\lambda_m)\neq 0.$$
From \eqref{normeq} it follows that $2\ipr{x_i}{x_j} = 1-\ipr{x_i}{x_i} - \ipr{x_j}{x_j}$ for all $i\neq j$.
Therefore, for any $j$,
\begin{align*}
0 &= \sum_{i=1}^m\lambda_i\ipr{x_i}{x_j} \\
  &= 2\lambda_j\ipr{x_j}{x_j} + \sum_{i\neq j}\lambda_i\left(1-\ipr{x_i}{x_i} - \ipr{x_j}{x_j}\right) \\
  &= 2\lambda_j\ipr{x_j}{x_j} - \lambda_j\left(1-2\ipr{x_j}{x_j}\right) - \sum_{i=1}^m\lambda_i\ipr{x_i}{x_i} \\
  &= \lambda_j\left(4\ipr{x_j}{x_j} - 1\right) - \sum_{i=1}^m\lambda_i\ipr{x_i}{x_i}.
\end{align*}
Thus
\begin{equation}\label{lambdaj}
\lambda_j(4\ipr{x_j}{x_j} -1) = \sum_{i=1}^m\lambda_i\ipr{x_i}{x_i}.
\end{equation}
Summing over all $j$, we obtain $$(m-4)\sum_{i=1}^m\lambda_i\ipr{x_i}{x_i}=0.$$
Since $d\neq 2$, we have $m\neq 4$, and by \eqref{lambdaj} we obtain that $\lambda_j=0$ for all $j$ such that $\norm{x_j}\neq\half$.
Since there is at most one $x_j$ of norm $\half$, and $\sum_i\lambda_i = 0$, we obtain that $\lambda_j=0$ for all $j$, contradicting affine dependence.

To obtain $d+1$ vectors satisfying \eqref{normeq}, consider the vertices of a regular simplex with centroid at the origin.
By symmetry, the sum of any two vectors will have the same norm.
By a suitable scaling, \eqref{normeq} will be satisfied.
\end{proof}

Theorem~\ref{l2m} now follows immediately from Theorems~\ref{l2bound} and \ref{helly:th2}.

\begin{proof}[Proof of Theorem~\ref{inftybound}]
Let $\{x_1,\dots,x_m\}$ be an M-set in $\ell_\infty^d$.
For each $i$ there is a coordinate of $x_i$ of absolute value larger than $1$.
If $m\geq 2d+1$, then by the pigeon-hole principle, there are three $x_i$'s with the same coordinate of absolute value $>1$.
Some two of these will have the same sign, and therefore the norm of their sum will be $>2$, a contradiction.

Therefore, $m\leq 2d$.
Suppose now that equality holds.
Then, by the same considerations as above, we may reorder the $x_i$'s such that $x_i^{(i)}>1$ and $x_{i+d}^{(i)}<-1$ for $i=1,\dots,d$, and $x_i^{(j)}\leq 1$ in all other cases.
If $x_i^{(i)}+x_{i+d}^{(i)}\geq 2$ for some $i$, then $x_1^{(i)}\geq 2-x_{i+d}^{(i)}>3$.
Thus $x_j^{(i)}<-1$ for all $j\neq i$, since $\norm{x_i+x_j}_\infty =2$.
But then $\norm{x_j+x_{j'}}_\infty >2$ for any two $j,j'\neq i$, a contradiction.
Similarly, $x_i^{(i)}+x_{i+d}^{(i)}\leq -2$ is impossible.
Since $\norm{x_i+x_{i+d}}_\infty=2$, there is a coordinate $j\neq i$ such that $x_i^{(j)} + x_{i+d}^{(j)} = \pm 2$.
Since $\abs{x_i^{(j)}},\abs{x_{i+d}^{(j)}}\leq 1$, we must have $x_i^{(j)}=x_{i+d}^{(j)}=\pm 1$.
But then either $\norm{x_i+x_j}_\infty>2$ or $\norm{x_{i+d}+x_j}_\infty>2$, a contradiction.
Thus $m\leq 2d-1$.

For the lower bound, let $e_1,\dots,e_d\in\R^d$ be the standard unit vectors.
For $i=1,\dots,d-1$, let $x_i:=2e_i+ e_d$ and $x_{i+d-1}:=-2e_1+ e_d$. Finally, let $x_{2d-1}:=-3 e_d$.
\end{proof}

\section{Concluding remarks}

Any centrally symmetric polytope $P$ defines a Minkowski space $X$ for which $m(P)=m(X)$ is finite.
It would be of some interest to study this number for specific polytopes.
In Theorem~\ref{inftybound} we have already determined $M(P)$ for cubes.

We make the following simple observations.

\begin{observation}\label{polytopeupper}
If a Minkowski space $X$ of dimension at least $2$ has a polytope with $f$ facets as unit ball, then $m(X)\leq f-1$.
\end{observation}
\begin{proof}
Imbed the space isometrically into $\ell_\infty^{f/2}$ and apply Theorem~\ref{inftybound}.
Note that since the dimension is at least $2$, $f\geq 4$.
\end{proof}

\begin{observation}\label{polytopelower1}
Let $X$ be a Minkowski space with a polytopal unit ball.
Let $f$ be the largest number of facets that a proper face of the unit ball can have.
Then $m(X)\geq f$.
\end{observation}
\begin{proof}
Let $F$ be a face of the unit ball that has $f$ facets.
Choose points $x_1,\dots,x_f$ just outside the centroid of each facet of $F$ in $\aff F$.
Then $\{x_1,\dots,x_f\}$ is clearly an M-set.
\end{proof}

\begin{observation}\label{polytopelower2}
In any $d$-dimensional Minkowski space with a polytope as unit ball, there exists an M-set of cardinality at least $d+1$.
\end{observation}
\begin{proof}
Choose any facet $F$ of the unit ball.
Let $F$ have $f$ facets.
By the argument of Theorem~\ref{polytopelower1}, $m(X)\geq f$.
If $f\geq d+1$, we are done.
Otherwise $F$ is a simplex with centroid $c$, say.
It is easily seen that we may add $-3c$ to the M-set obtained in Theorem~\ref{polytopelower1}, to obtain an M-set of $d+1$ points.
\end{proof}

These estimates are not sharp.
As a start, consider cross-polytopes.
By Observation~\ref{polytopelower2} we obtain $m(\ell_1^d)\geq d+1$, and Observation~\ref{polytopelower1} doesn't give any larger lower bound.
As upper bound we only have $2^d-1$, by Observation~\ref{polytopeupper}.
Perhaps there is also a linear upper bound?

\begin{conjecture}
 Each $d$-dimensional Minkowski space has an M-set of at least $d+1$ points.
\end{conjecture}
Theorem~\ref{helly:th2} and Lemma~\ref{threelow} show that this conjecture is true for $d=2$ and $3$.
Showing it for $d\geq 4$ would undoubtedly be more difficult than Observation~\ref{polytopelower2}.

\section*{Acknowledgements}
The author thanks Valeriu Soltan for drawing his attention to the work of Zamfirescu, and the anonymous referee for his suggestions on improving the paper.

config: arxiv


TITLE: What do you call geometric patterns like this?
QUESTION [3 upvotes]: What do you call geometric patterns like this ?

REPLY [2 votes]: I'd call this a self-similar logarithmic spiral which is the image of a square under an iterated similarity transformation.
There is a similarity transformation which maps the corners of the outermost square to those of the second-to-largest one, and iterating that transformation will yield an infinite sequence of squares converging to the center. The spiral arms are formed by the corner points of these squares and the outermost bit of edge, while most of the edges contribute to the hatching of the spiral parts. In my own rendering, created with Cinderella, there is even some Moiré effect visible near the center.


TITLE: Existence of a function
QUESTION [0 upvotes]: [also asked here http://math.stackexchange.com/questions/307197]
All arguments are in $\mathbb{R}^3$.
Suppose $n(x)$ is a smooth function where $\mathbf{supp}(n(x)-1)$ is a compact set $\Omega$. i.e. $n(x) = 1$ when $x$ is outside $\Omega$.
Assume there are some points $x_j\in\Omega$, where $j=1,2,\cdots.m$.
Consider Helmholtz equation
$\Delta u + k^2 n(x) u = 0$
And I want to know if there is a function $u$ satisfies the equation and also vanishes at $x_j$. i.e. $u(x_j) = 0$. [Certainly $u$ can have a manifold of zeros, here I just restrict it on $x_j$]
Here you may try to give a method to construct $u$, or prove the existence.
And you may take $m=1$ here at first.
Thanks.

REPLY [1 votes]: I proved the existence here.


TITLE: I need help with this question about convergent series
QUESTION [2 upvotes]: Let $\sum_{n=1}^{∞} a_n$ be a convergent series, consider a growing sequence where $n_1<n_2<...$ , and define
• $b_1=a_1+...+a_{n_1}$;
• $b_2=a_{n_1+1}+...+a_{n_2}$;
• $b_3=a_{n_2+1}+...+a_{n_3}$, and so on.
Prove that $\sum_{n=1}^{∞}b_n=\sum_{n=1}^{∞}a_n$
So, I know that $\sum_{n=1}^{∞} a_n=a_1+a_2+a_3+...$
And I know that
$\sum_{n=1}^{∞}b_n=b_1+b_2+b_3+...=a_1+...+a_{n_1}+a_{n_1+1}+...+a_{n_2}+a_{n_2+1}+...+a_{n_3}+...$
But how can I prove they are equal?

REPLY [1 votes]: Let the series converge to $L$. Therefore
$$
{\exists M\ \ \ ,\ \ \ n>M\implies|\sum_{i=1}^n a_i-L|<\epsilon_1}.
$$
All you need to prove, is that
$$
\exists N\ \ \ , \ \ \ m,n>N\implies |\sum_{i=1}^m a_i-\sum_{i=1}^n b_i|<\epsilon_2,
$$
by expressing $b_n$ in terms of $a_n$ in a bounded series. Can you finish now?


TITLE: What are the Effects of Applying the $\nabla$ Operator to a Scalar Field?
QUESTION [0 upvotes]: For a scalar field $f(\mathbf{r})$, where $\mathbf{r} = (x, y, z)$, what kind of mathematical object results from each of the operations below:
$$ \nabla f(\textbf{r}) $$
$$ \nabla \cdot f(\textbf{r}) $$
$$ \nabla \times f(\textbf{r})? $$

REPLY [2 votes]: The gradient operator $\nabla$ takes a scalar function $f:\mathbf{R}^3 \to \mathbf{R}$ and produces a vector field $(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}).$ The notation $\nabla \cdot F$ denotes the application of $\nabla$ componentwise on the components of a vector field $\mathbf{R}^3\to \mathbf{R}^3$, and describes the divergence. For instance, if the vector field is
$$ F(x,y,z)=P(x,y,z)i+Q(x,y,z)j+R(x,y,z),$$
the divergence operator gives us
$$ \nabla\cdot F(x,y,z)=\frac{\partial P}{\partial x}(x,y,z)+\frac{\partial Q}{\partial y}(x,y,z)+\frac{\partial R}{\partial z}(x,y,z).$$
In particular, this quantity lives in $\mathbf{R}$. However, the divergence operator needs to act on a vector field to be meaningful. In the case where $f:\mathbf{R}^3\to \mathbf{R}$ is a scalar field, there is not a good idea of what that would mean. For the last case, we have the standard curl operator $\nabla\times F$ where $F$ is as above. It is given by
$$ \nabla\times F=\bigg(\frac{\partial R}{\partial y}-\frac{\partial Q}{\partial z}\bigg)i+\bigg(\frac{\partial P}{\partial z}-\frac{\partial R}{\partial x}\bigg)j+\bigg(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}\bigg)k.$$
As we can see, this eats a vector field $F:\mathbf{R}^3\to\mathbf{R}^3$ and produces a vector field in $\mathbf{R}^3$. In your case, it is not obvious how to apply $\nabla\times f$ where $f$ is a scalar field, because the curl depends on the cross product in its definition, which is a construction special to $\mathbf{R}^3$.
In summation, the only one of these expressions that is well-defined is the first.


TITLE: What does negation affect and what does it not?
QUESTION [0 upvotes]: This question came up on Keith Devlin's course on mathematical thinking. I was wondering whether the negation would apply to all terms in this sentence? As in, does the negation apply to 'all' and 'badly' but not the adjective and noun 'foreign cars'?

REPLY [2 votes]: If by 'terms' you mean predicate terms (like 'being well-made' or 'being foreign made') or subject terms (like 'cars'), then the answer is none.
You never negate subject terms or predicate terms.  Likewise, you never negate objects or predicates.  The only thing you ever negate are sentences.
This is because sentences express claims that can be true or false, and that is something we can negate.  But an object is ... well, an object.  What would it mean to 'negate an object'?! Sometimes I see the expression $\neg x = y$ interpreted as $(\neg x)=y$, which supposedly would mean something like 'the object that is the negation of $x$ equals $y$' ... Which is really just unsyntactical and meaningless.
The same goes for predicates. That is, when we write $x \not = y$, we don't mean something like $x (\neg =) y$, but $\neg (x=y)$; again, the only thing you negate are sentences.


TITLE: When integrating by subsitution, once you find the function $u$ and $du$, what function should I integrate?
QUESTION [0 upvotes]: I understand that whenever we have a function to integrate, and we decide to use the method of variable change, what function should we integrate. For example, suppose this problem:
$ \int_{a}^{b}\frac{1}{\sqrt{1-x^{3}}}$
I decided that I will use the new function $u=\sqrt{1-x^{3}}$, and so $du=\frac{-3x^{2}}{2\sqrt{1-x^{3}}}$ or simply $\frac{-3x^{2}}{2u}$
Now, what goal should I keep in my mind when deciding what to put instead of the question mark here:
$\int_{a}^{b}? du$ ? Should I just put $1/u$? Or something else?

REPLY [2 votes]: If you set $u=\sqrt{1-x^{3}}$,
then
$$du=\frac{-3x^{2}}{2\sqrt{1-x^{3}}} \color{red}{dx}.\tag1$$
Don't forget the $dx$ on the right-hand side! Also use $dx$ when writing the integral:
$$\int_{a}^{b}\frac{1}{\sqrt{1-x^{3}}}\color{red}{dx}.
$$
Consistent use of $dx$ will make the formal
manipulations a lot clearer.
From Equation (1) we get
$$dx=\frac{2\sqrt{1-x^{3}}}{-3x^{2}} du,$$
and we can use this equation to simply substitute $\dfrac{2\sqrt{1-x^{3}}}{-3x^{2}} du$ for $dx$ in the integral:
$$
\int_{a}^{b}\frac{1}{\sqrt{1-x^{3}}} dx
= \int_{a}^{b}\left(\frac{1}{\sqrt{1-x^{3}}}\right)
   \frac{2\sqrt{1-x^{3}}}{-3x^{2}}  du
= -\frac23\int_{a}^{b} \frac{1}{x^{2}}  du.
$$
That's correct, but not simple to solve, since we still need to
convert the $x^2$ in the denominator into some function of $u$.
It turns out $x^2 = (1 - u^2)^{2/3}$, so we end up with
$$
-\frac23\int_{a}^{b} \frac{1}{(1 - u^2)^{2/3}}  du,
$$
which is not really an improvement over the original form of the integral.
The obvious substitutions don't always work, and this is one of the
cases where they do not.


TITLE: Summation of gaussians
QUESTION [0 upvotes]: Suppose we have given constants $A_i, x_i (i=1..N)$
Is it possible to approximately calculate the sum of N gaussians in less than N iterations for any x? (may be with some preprocessing)
$$\sum_{i=1}^{N}A_i  e^{-(x-x_i)^2}$$
The same question for 2D case
$$\sum_{i=1}^{N}A_i  e^{-(x-x_i)^2-(y-y_i)^2}$$

REPLY [0 votes]: If the $A_i$, $x_i$ and $y_i$ can take arbitrary values then it is hopeless because you even need $N$ iterations merely to read the input.
If (some of) the Gaussians have the same distribution then their sum can be expressed as a single Gaussian.


TITLE: (Poisson process) probability of least one arrival of $N_2$ and $N_3$, between 2 successive arrivals of $N_1$.
QUESTION [0 upvotes]: Let  $N_1(t),~N_2(t),~N_3(t)$ be independent Poisson processes of rates $\lambda_1,~\lambda_2,~\lambda_3>0$, respectively. Evaluate the probability between two successive arrivals of $N_1$ we get at least one arrival of $N_2$ and at least one arrival of $N_3$.

Attempt. The probability of getting an arrival of $N_i$ (and not of $N_j,~j\neq i$) is $\frac{\lambda_i}{\lambda_1+\lambda_2+\lambda_3}.$ If we already have an arrival of $N_1$, then for
$E_i$= the event we get at least one arrival of $N_i$,
we seek $P(E_2\cap E_3)=P(E_2)P(E_3).$ So:
$$P(E_2)=1-P(E_2^c)=1-\sum_{k=0}^{\infty}\bigg(\frac{\lambda_3}{\lambda_1+\lambda_2+\lambda_3}\bigg)^k\frac{\lambda_1}{\lambda_1+\lambda_2+\lambda_3}$$
and similar we work for $E_3.$ Am I on the right path?
Thanks in advance.

REPLY [1 votes]: After the initial $N_1$ you need to either first get $N_2$ and then $N_3$ before $N_1$, or first $N_3$ and then $N_2$ before $N_1$. Thus the probability is
$$
\frac{\lambda_2}{\lambda_1+\lambda_2+\lambda_3}\frac{\lambda_3}{\lambda_1+\lambda_3}+\frac{\lambda_3}{\lambda_1+\lambda_2+\lambda_3}\frac{\lambda_2}{\lambda_1+\lambda_2}=\frac{\lambda_2\lambda_3}{\lambda_1+\lambda_2+\lambda_3}\left(\frac1{\lambda_1+\lambda_3}+\frac1{\lambda_1+\lambda_2}\right)\;.
$$


TITLE: Proof by Induction Problem
QUESTION [5 upvotes]: I need to prove that for natural numbers $a$, and positive integers $n$, the number $a^{2n+1}-a$ is divisible by $6$.
I have proved the case when $n=1$, that $a^3-a$ is divisible by $6$. I'm having trouble proving that if $a^{2n+1}-a$ is divisible by $6$, then $a^{2(n+1)+1}-a$ is also divisible by $6$.

REPLY [0 votes]: Since you have proved $a^3-a$ is divisive by $6$ that is $a^3-a=6q_1$ where $q_1\in\mathbb{Z}.$ If $a^{2n+1}-a$ is divisible by $6,$ then  $a^{2n+1}-a=6q$ where $q\in \mathbb{Z}.$ This implies $a^{2n+1}=6q+a,$ equivalently $a^{2n+3}=6qa^2+a^3$ or $$a^{2n+3}-a=6qa^2+a^3-a=6qa^2+6q_1=6(qa^2+q_1)$$
That is $a^{2n+3}-a$ is divisible by 6.


TITLE: Component of component of vectors
QUESTION [0 upvotes]: Let's take an example to make you understand my doubt... A force $\vec{F}$ is in +$\hat{y}$ axis. First, I take $\pm 45°$ component of it both the sides against $\hat{y}$ axis... And then again, I repeat the same process, take components of these two, each one their  $\pm 45°$ components ... Resulting in 4 vectors... 2 along y and 2 opposite in directions along $\pm x$ axis... So by this way can we state that taking component of component of vectors is allowed?
Irrespective of this discussion:
Why taking components of a component of a vector is invalid?
And my second doubt was why is it so that we take components inside a quadrant only?... Like 30° or 60°...why not taking components of 130° or greater uncommon?

REPLY [1 votes]: You can write a force vector as a sum of as many vectors as you want. Take the vector $2\hat{j}$, which I will write as $(0,2)$ for clarity. You can write this as two vectors:
$$(0,2) = (1,1) + (-1,1)$$
As three vectors:
$$(0,2) = (1,0.5) + (0,1) + (-1,0.5)$$
As four vectors:
$$(0,2) = (-1,-1/2) + (-1, 1) + (1,3) + (1, -3/2)$$
And so on, in whatever way you like. However, in order to find the resulting motion of a particle being acted upon by the force $(0,2)$, you have to use all of the vectors that sum to $(0,2)$. If you write $(0,2) = (1,1) + (-1,1)$, the particle is not alternately pulled to the left or right. The particle only moves vertically because the total force is only vertical. It is a mathematical convenience to decompose a vector into the sum of other vectors. A vertical force is still a vertical force, no matter how you write it down.
The only reason one vector in the sum would be ignored is if there is another force that cancels it out. For example, I'm dragging a sofa across the floor using a force of $(100, 50)\,\textrm{N}$. The sofa is heavier than 50 newtons, so the sofa stays on the floor. Because the sofa has no vertical motion, I can deduce that gravity and the normal force of the floor act upon the sofa to cancel out the vertical force. This results in a total force of $(100,0) = \vec{F}_{dragging} + \vec{F}_{gravity} + \vec{F}_{normal}$.