TITLE: KPZ relation $\chi = 2 \xi -1$ in a random geometric graph
QUESTION [1 upvotes]: If I have $n$ points uniformly distributed on the surface of a torus, and form a graph by adding an edge between pairs whenever they are within a unit distance (induced by the Euclidean metric), I have a random geometric graph.
The Euclidean geodesic between two vertices is just the shortest path, but where short means the Euclidean lengths of the edges along the path sum to the smallest value available (i.e. over all available paths).
The passage time, as well as the maximum transversal wandering of the path w.r.t. the straight line joining them, are conjectured to follow a power law $L^{\chi}$ and $L^{\xi}$, where $L$ is the Euclidean distance between the endpoints (norm of the straight line segment joining them). The first exponent $\chi$ gives the scaling of the root of the variance of the travel time, while the second $\xi$ gives instead the scaling of the expected wandering. See e.g. "The universal relation between scaling exponents in first-passage percolation", S. Chatterjee, Ann. Maths 177, 2013.
The relation $\chi = 2 \xi -1$ is meant to hold, relating the orders of the expected magnitude of the transversal wandering, and the time fluctuations about their expectation.
Though we don't know the exponents, should such a relation also hold in the spatial setting of the random geometric graph, where the edge lengths are not i.i.d., but induced by a point process?

REPLY [1 votes]: The answer is yes. All that is required for the relation to hold, at least in principle, is that the curvature exponent of the limit shape in direction $u$ to be exactly 2.
See the review of Auffinger and Damron, Section 4.3.


TITLE: Continuity of the restriction of a function 3
QUESTION [0 upvotes]: Let $F$ be a compact of $\mathbb{R}^{n}$, $f$ a real valued function defined on a neighborhood of $F$, and set $$A=\{x\in F: f(x)=0\}.$$ Suppose that for all $y\in A$ we have$$\lim f(y)=0$$ as $y\to x$ form inside $F$.
1) Can we conclude that $A$ is closed in $\mathbb{R}^{n}$?
2) Can we conclude that $A$ is closed in $F$?

REPLY [0 votes]: Let $A=(1,2), F=[0,3]$  and $f(x)=0$ for $x \in A$, $1$ for all other $x$. This shows that $A$ need not be closed either in $F$ or in $\mathbb R$.


text:
\begin{document}

\title{A micro/macro parallel-in-time (parareal) algorithm applied to a climate model with discontinuous non-monotone coefficients and oscillatory forcing}

\author[1]{Giovanni Samaey}
\author[2]{Thomas Slawig}
\affil[1]{
\small
Department of Computer Science, K.U. Leuven,
Celestijnenlaan 200A, 3001 Leuven, Belgium, \tt{giovanni.samaey@kuleuven.be}}
\affil[2]{
\small Department of Computer Science, Kiel University,  24098 Kiel, Germany, \tt{ts@informatik.uni-kiel.de}
}

\maketitle

\begin{abstract}
We present the application of a micro/macro parareal algorithm for a 1-D energy balance climate model with discontinuous and non-monotone coefficients and forcing terms. The micro/macro parareal method uses a coarse propagator, based on a (macroscopic) 0-D approximation of the underlying (microscopic) 1-D model. We compare the performance of the method using different versions of the macro model, as well as different numerical schemes for the micro propagator, namely an explicit Euler method with constant stepsize and an adaptive library routine.
We study convergence of the method and the theoretical gain in computational time in a realization on parallel processors. We show that, in this example and for all settings, the micro/macro parareal method converges in fewer iterations than the number of used parareal subintervals, and that a theoretical gain in performance of up to 10 is possible.
\end{abstract}

\maketitle

\noindent
\textbf{Keywords:} Parallel-in-time algorithm; micro/macro parareal algorithm; energy balance climate model; nonlinear partial differential equation.

\section{Introduction}

Climate simulation is among the most challenging and time-consuming computational tasks, for a number of reasons. First, there is the complexity of the coupled climate system, with interactions between many different components and nonlinearity of many important processes, some of which are not completely understood by now (see e.g. \cite{McgHen14}).
Second, the need for high spatial resolution in global climate models results in a huge dimension of the discretized systems that have to be solved. To reduce the needed computational time, spatial parallelization is a common strategy used in fully coupled high resolution climate models.  Still, for long-time  simulation  runs (e.g., to compute full glacial cycles of  hundreds of thousands of years), spatial parallelization alone is insufficient, since spatial and temporal resolutions are typically coupled via some kind of CFL-type condition. Thus,  such long-time simulation runs are  only feasible with lower complexity models (see, e.g.,  \cite{Ganopolski2017}.  To proceed to long-time \emph{and} high resolution models, an additional parallelization in time becomes very attractive.

In this paper, we present how we generalized and used a micro/macro version of a parallel-in-time algorithm for time integration, called \emph{parareal} \cite{LiMaGa01}, as it was developed in \cite{LeLeSa13}.
 As an exemplary problem for a climate model, we choose an energy balance model (EBM) that describes the evolution of global mean temperature by balancing incoming and outgoing radiation, see, e.g., \cite[Section 3]{McgHen14}). Such models are the simplest way to model the Earth's climate. Incoming radiation is determined by the energy received from the Sun, diminished by a fraction (called albedo) which is reflected by the Earth's atmosphere (e.g., clouds) or surface (e.g., ice).  Outgoing radiation is usually determined by considering the Earth as a perfect black radiating body, for which the Stefan-Boltzmann law gives a relation between the outgoing radiation and the temperature of the body. Since -- due to the greenhouse effect -- the Earth is not perfectly radiating, an emissivity parameter is included in EBMs to take into account that part of the emitted heat is captured in the Earth's atmosphere.
Here, we regard as the ``original'' or ``micro'' model an EBM in one space dimension, using latitude as the spatial coordinate.

The micro/macro parareal method will use an approximate ``macro'' model to obtain a numerically fast predictor, which is iteratively corrected by time-parallel simulations using the original 1-D model on different slices of the time interval. The macro model is obtained by averaging the micro 1-D model in space, which results in a 0-D model that takes the form of an ordinary differential equation (ODE). The resulting macro model is much cheaper to simulation numerically for two reasons. First, due to its low-dimensionality, fewer degrees of freedom need to be accounted for. Second, the macro model only contains the dominant slow time scale, such that larger time steps can be taken.
We mainly use the 0-D model to accelerate simulations. Nevertheless, 0-D EBMs have their own justification,  for educational purposes (see \cite{McgHen14}),
since they include most important features of the Earth's energy balance. They model the  Earth as a point in space, and thus all parameters mentioned above
 (energy received from the Sun as only external forcing of the climate system as well as albedo and emissivity  that enter the differential equations as coefficients)
 represent averaged values over the whole planet.

If only a few parareal iterations are required, the micro/macro parareal algorithm can achieve a significant reduction in the required wall-clock time, compared to a naive simulation using only the micro-simulator by performing the 1-D simulations on different slices of the time interval \emph{in parallel} in each parareal iteration. In that case, one can reach an accuracy that is much higher than that of the macroscopic simulation on the whole time interval of interest, with a wall-clock time that is of the order of the simulation time of a few parareal time slices. In this paper, we study to what extent this potential is realized for climate simulations based on EBMs. In particular, we investigate the ability of the method to deal with oscillatory and abrupt changes in forcing terms, and we study the dependency of the method on the accuracy of the macro model and the number of parareal time slices.

Since its introduction in~\cite{LiMaGa01}, the parareal strategy has been applied to a wide range of problems, including fluid-structure interaction~\cite{farhat2003time}, Navier--Stokes equation simulation~\cite{fischer2005parareal}, and reservoir simulation~\cite{garrido2005}. We refer to ~\cite{maday2002parareal,maday2005parareal} for further analysis, and to \cite{bal2005convergence,staff2005stability} for stability results. In~\cite{gander2007analysis}, the method is reformulated in a more general setting that relates the parareal strategy to earlier time-parallel algorithms, such as multiple shooting (see e.g.~\cite{keller1968numerical,nievergelt1964parallel}) or multigrid waveform relaxation (see e.g.~\cite{lubich1987multi,vandewalle1992efficient}).  The micro-macro parareal method in this paper is a direct generalization of the method in \cite{LeLeSa13}, in which the micro model was a high-dimensional stiff ODE and the macro model was an approximate, low-dimensional ODE for a limited set of slow degrees of freedom. Similar micro-macro parareal methods have been considered in the literature. The authors of~\cite{BBK,maday41parareal} consider a singularly perturbed system of ordinary differential equations (ODEs) at the microscopic level and the limiting differential-algebraic equation at the macroscopic level. In~\cite{engblom2009parallel}, a parareal algorithm for multiscale stochastic chemical kinetics is presented, in which the macroscopic level uses the mean-field limiting ODE. In~\cite{mitran2010time}, the parareal algorithm is used with kinetic Monte Carlo at the macroscopic level and molecular dynamics at the microscopic level.

This text is organized as follows. In Section~\ref{sec:model}, we describe the model problem in its 0-D and 1-D versions. We also describe  how the 1-D coefficients are captured in the 0-D version. In Section \ref{sec:parareal}, we describe the micro/macro parareal algorithm that is the focus of the present paper. In Section \ref{sec:numexp}, we describe in detail the numerical experiments that we performed. We present the numerical results in Section~\ref{sec:results}, discussing both the effect of the choice of coefficients in the 0-D model and the choice of micro and macro time integration methods. We end the paper with a summary and conclusions in Section \ref{sec:discussion}.

\section{Model problem\label{sec:model}}

In this Section, we discuss the energy balance models (EBMs) that will be used throughout the paper.  In the 1-D model, we have a single spatial coordinate $\phi$, which varies from $\phi=0$ at the north pole to $\phi= \pi/2$ at the equator, thus named colatitude. Due to a symmetry assumption (following \cite{Bay91,Ghi75}), only half of the sphere is modeled.   The state variable, temperature, then is a function $T=T(\phi,t)$ of colatitude and time. In this model, spatial redistribution of heat is included via diffusion, with a potentially temperature-dependent diffusion coefficient. Incoming radiation, albedo and thermal capacity depend on the colatitude $\phi$. Albedo and emissivity (which models the greenhouse effect) also depend on temperature, to take into account ice melting or increase of water vapor in the atmosphere, respectively, leading to additional nonlinearities in the model.
By averaging over the colatitude $\phi$, a 0-D model can be obtained, in which the Earth is seen as a point in space, and thus all parameters (energy received from the Sun, albedo, emissivity) are  averaged over the whole planet. The state variable  is again temperature $T$, but now as a global mean value, i.e., a scalar function $T$ of time $t$.

For clarity of exposition, we first discuss the 0-D model in Section~\ref{sec:model-0-D}, after which we elaborate on the 1-D model in Section~\ref{sec:model-1-D}.  The choice of the parameters in the 1-D model is discussed in Section~\ref{sec:model-params}, along with the connection between the 0-D and 1-D models.  We conclude this Section with some comments on existence and uniqueness of solutions (Section~\ref{sec:model-exist}).

\subsection{0-D model\label{sec:model-0-D}}
In a 0-D EBM \cite[Section 3.2]{McgHen14},
we consider the instantaneous change of the time-dependent temperature $T(t)$ due to the difference between incoming and outgoing radiation.
Any difference between ingoing and~outgoing, energy -- $R_{in}(t)$, resp.,~$R_{out}(t)$ -- induces a temporal change of thermal energy,
\begin{eqnarray*}
\left(4\pi r^{2}h\;c\;\rho\right) \; T'(t) = R_{in}(t)-R_{out}(t)
\end{eqnarray*}
in which  $4\pi r^{2} h$ is the volume of the considered spherical shell, $c$ is the specific heat of the fluid, and $\rho$ the respective density.

The total amount of incoming radiant energy per  time unit for the whole Earth  is given as
\begin{eqnarray}
\label{eq:Rin}
R_{in}&=&(1-\alpha)\pi r^{2} S
\end{eqnarray}
 with  $\alpha\in[0,1]$ the albedo, $\pi r^{2}$ the area of the Earth's two-dimensional projection  ($r$ being the Earth's radius), and $S\approx 1367\,{\rm Wm}^{-2}$ the amount of energy per second and area (often denoted as solar ``constant'', but it fact not constant due to temporal variations of solar activity). In general, the albedo $\alpha$ may be temperature-dependent. We will specify its value in Section~\ref{sec:model-params}.

The outgoing radiant energy per unit time is given by the Stefan-Boltzmann law,
\begin{eqnarray}
\label{eq:Rout}
R_{out}&=&4\pi r^{2}\;\epsilon\;\sigma T^{4}
\end{eqnarray}
 with $4\pi r^{2}$ the Earth's surface and $\sigma= 5.67\times 10^{-8}{\rm Wm}^{-2}\,{\rm K}^{-4}$ the Stefan-Boltzmann constant, and including the emissivity $\epsilon$, i.e., the fraction of outgoing radiation that is not captured by the atmosphere.  Like the albedo $\alpha$, the emissivity $\epsilon$ usually depends on the temperature, see~Section~\ref{sec:model-params}.

Since $\pi r^{2}$ cancels out in the balance equation, we obtain
as resulting ODE
\begin{eqnarray}
\label{eq:0-D}\label{eq:macro}
T'(t) &=&\frac1C\left((1-\alpha)Q- \epsilon\sigma T(t)^4\right),
\end{eqnarray}
where we introduced the symbol $Q:=S/4$ for notational convenience.

In the easiest case of constant solar insolation $Q:=S/4$, heat capacity $C=hc\rho$, emissivity $\epsilon$, and albedo $\alpha$, a stationary solution can be  computed from \eqref{eq:macro}
as
\begin{eqnarray}
\label{eq:stat0-D}
 T_{stat}&=&\sqrt[4]{\frac{(1-\alpha)Q}{ \epsilon\sigma}}.
\end{eqnarray}

\subsection{1-D model\label{sec:model-1-D}}

The following 1-D model, based on \cite{Bay91,Ghi75} and \cite[Section 10]{GhiChi87}, includes variation of the temperature in the latitudinal direction
 $\phi \in [0,\pi/2]$, i.e., we model a half sphere from north pole to equator. To this end, we add diffusive term of the form $\nabla\cdot (k\nabla T(\phi,t))$. Using spherical coordinates and taking into account that colatitude is the only spatial coordinate, this term reduces to
\begin{eqnarray}
\label{eq:diffterm}
\nabla\cdot (k\nabla T(\phi,t))&=&\frac{1}{\sin\phi}\frac{\partial}{\partial \phi}\left(k\sin \phi\frac{\partial T(\phi,t)}{\partial \phi}\right).
\end{eqnarray}

The coefficient $k$ may, in principle, depend on space and temperature. We also include a dependence on the spatial coordinate for the solar radiation $Q$, the albedo $\alpha$, and the heat capacity $C$ a on the spatial coordinate $\phi$.

In the most general setting, the complete 1-D equation of energy balance gives:
\begin{eqnarray}
\label{eq:1}
C(\phi,T)\frac{\partial T}{\partial t}=
\frac{1}{\sin\phi}\frac{\partial}{\partial \phi}\left(k(\phi,T)\sin \phi\frac{\partial T}{\partial \phi}\right)+\left[1-\alpha(\phi,T)\right]Q(\phi)
 - \epsilon(T)\sigma T^4.
\end{eqnarray}

We use homogeneous Neumann boundary conditions
\begin{eqnarray*}
\frac{\partial T}{\partial \phi}(\phi,t)&=&0\quad\text{for }\phi\in\left\{0,\frac\pi2\right\},t\ge0,
\end{eqnarray*}
where the condition at $\phi=\pi/2$ (equator)  comes from the  assumed symmetry.  The precise choice of the coefficients is the subject of the next Section~\ref{sec:model-params}.

In the numerical experiments, we discretize equation~\eqref{eq:1} in space. We introduce a spatial grid for $\phi\in[0,\pi/2]$ with stepsize $\Delta\phi=\pi/(2I),I\in\mathbb{N}$. The gridpoints then are $\phi_{i}=i\Delta\phi,i=0,\ldots,I$. On these points, we denote the approximate  solution as $\bT_{i}(t)\approx T(\phi_{i},t)$. We then use standard finite differences, see Appendix~\ref{sec:discr}, to obtain a semi-discretization,
\begin{eqnarray}
\bT'(t)&=&f(t,\bT(t)),\quad t\ge0\label{eq:micro},
\end{eqnarray}
with $\bT(t),f(t,\bT(t))\in \mathbb{R}^{I-1}$ and some initial value $\bT(0)\in \mathbb{R}^{I-1}$. Equation~\eqref{eq:micro} is subsequently discretized in time, either by a forward Euler method with fixed time step or by a variable-step variable order linear multistep method, see Section~\ref{sec:numexp}.

\subsection{Spatial and temperature-dependent modeling\label{sec:model-params}}

 In the 1-D model, the parameters solar radiation $Q$ and albedo $\alpha$ can depend on colatitude $\phi$. The heat capacity $C$, albedo $\alpha$ and emissivity $\epsilon$ can depend on temperature, leading to nonlinearity in the models.  We first discuss the heat capacity, diffusion and emissivity (Section~\ref{sec:hde}), after which we go into more depth into the albedo (Section~\ref{sec:alb}) and the solar radiation (Section~\ref{sec:rad}).

\subsubsection{Heat capacity, diffusion and emissivity\label{sec:hde}}

\paragraph{Heat capacity.} For the heat capacity $C(T)$, we use
the temperature-dependent function  from \cite[Section 4.5]{Bay91}:
\begin{eqnarray*}
C(T)&=&\left(C_{1}+C_{2}\tanh\left(C_{3}(T-T_{s})\right)\right)\times 10^{8}\,{\rm Jm}^{-2}{\rm K}^{-1}\end{eqnarray*}
with $C_{1}=3.14,C_{2}=1.15,C_{3}=0.08$, and $T_{s}=263.15$. It is shown in Figure \ref{fig:emmi}.

\paragraph{Diffusion coefficient.} For the diffusion coefficient $k(\phi,T)$, Ghil \cite[Eqn. (2e)]{Ghi75} used
the following nonlinear function of colatitude and temperature
\begin{eqnarray*}
k(\phi,T)&=&k_{1}(\phi)+k_{2}(\phi)g(T),\quad g(T)=\frac{c_{4}}{T^{2}}\exp\left(-\frac{c_{5}}{T}\right),
\end{eqnarray*}
with
 given  data for the coefficient functions  $k_{1},k_{2}$ at the grid points (see \cite{Ghi75}) and  parameters
\begin{eqnarray*}
c_{4}

&=&6.105\times0.75\times\exp(19.6)\times10^{2}{\rm N  K  m}^{-2}
 \approx1.4891\times 10^{11} {\rm NKm}^{-2},\\
c_{5}&=&5.35\times10^{3}\,{\rm K}.
\end{eqnarray*}

Because this choice (especially the coefficients suggested in \cite{Ghi75}) does not guarantee positivity of $k$, we use a constant value, namely $k=0.591\,{\rm W m}^{-2}{\rm K}^{-1}$, taken from \cite[Section 4.5]{Bay91}.

Note that, also in the literature cited above, mostly only one source of nonlinearity (either $C$ or $k$) has been chosen.

\paragraph{Emissivity.} For the emissivity $\epsilon$, we use
\begin{eqnarray*}
\epsilon(T)&=&1-\epsilon_{1}\tanh\left(\left(\frac{T}{T_{\epsilon}}\right)^{6}\right)
\end{eqnarray*}
with $\epsilon_{1}  = 0.5, T_{\epsilon}=284.15\,{\rm K}$, which  is  Ghil's   \cite{Ghi75} suggestion.

 \begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/emissivity}
	\includegraphics[scale=0.475]{figuren/heat_capacity}
	\caption{Left: Emissivity  $\epsilon=\epsilon(T)$. Right: Heat capacity $C=C(T)$. \label{fig:emmi}}
\end{figure}

\subsubsection{Albedo\label{sec:alb}}
The albedo $\alpha(\phi,T)$ describes the fraction of solar radiation that is absorbed by the earth. Since this fraction depends on features such as clouds or snow and ice on the surface, the temperature-dependence is obvious.
Piecewise constant and linear functions as well as smoothed counterparts are used in the literature. For the 1-D model,  spatial dependency reflects the bigger area covered by ice or snow (with higher albedo) near the poles. We use
 the model of North \cite{Nor75}, see also \cite[Section 4.5]{Bay91}:
\begin{eqnarray}
\label{eq:albedo}
\alpha(\phi,T)&=&
\left\{
\begin{array}{ll}
\alpha_{\rm max},& T\le T_{s},\\
\alpha_{1}+\alpha_{2}(\alpha_{3}\cos^{2}(\phi)-\alpha_{4}),&T> T_{s},
\end{array}
\right.
\end{eqnarray}
with $\alpha_{\rm max}=0.62,\alpha_{1}=0.303,\alpha_{2}=0.0779,\alpha_{3}=1.5,\alpha_{4}=0.5,$ and $T_{s}=263.15.$
Figure \ref{fig:albedo} shows the resulting non-discontinuous function.
 \begin{figure}
	\centering
	\includegraphics[scale=0.525]{figuren/albedo}
	\includegraphics[scale=0.45]{figuren/insolation}
	\caption{Left: Albedo  $\alpha=\alpha(\phi,T)$. Right: Solar radiation $Q=Q(\phi)$. \label{fig:albedo}}
\end{figure}

In the 0-D version, we restrict to a simple step function
\begin{eqnarray}
\label{eq:albedo-step}
\alpha(T)&=&\left\{
\begin{array}{ll}
\alpha_{\rm max},&T\le T_{s} ,\\[0.2em]
\alpha_{\rm min},\quad&T> T_{s},
\end{array}\right.
\end{eqnarray}
which is also mentioned in \cite{Bud69}, \cite[Section 10.2]{GhiChi87},
and \cite[Section 3.2.2]{McgHen14}. The parameter $\alpha_{\rm min}$ in the 0-D model is not present in the 1-D model. We chose $\alpha_{\rm min}=0.275$, which gives a good approximation of the long-time steady state obtained by the micro model.

\subsubsection{Solar radiation\label{sec:rad}}

\paragraph{Space dependence.} Clearly, the incoming solar radiation $Q(\phi)$ depends on the colatitude of the considered position on the Earth's surface.
We use a second order polynomial in the variable $x=\sin\phi$, see \cite[Section 4.3]{Sto14} or \cite[Section 4.2]{Bay91}),
 \begin{eqnarray}
 \label{eq:Q}
Q(\phi) &=& \frac{S}{4}\left(Q_{1} + Q_{2} \sin^2\phi\right)
\end{eqnarray}
with $Q_{1}=0.5294$ and $Q_{2}=0.706$, see Figure \ref{fig:albedo}.
Ghil \cite{Ghi75} took data of which the spatial distribution is similar to the one of the polynomial above.

 For the 0-D model, we use the average of~\eqref{eq:Q} over the half-sphere.
The spatial mean of a quantity depending on colatitude $\phi$ is obtained by integrating it over $\phi\in[0,\pi/2]$ and dividing by the measure of the area. In spherical coordinates the arc length of the circle at constant colatitude $\phi$ is given by $\pi r\sin\phi$ ($r$ being again the Earth's radius). This gives
\begin{eqnarray}\label{eq:spatial_mean}
\overline{Q}&=&\left(\int_{0}^{\frac\pi2}Q(\phi)\sin\phi \,d\phi\right)\left(\int_{0}^{\frac\pi2}\sin\phi \,d\phi\right)^{-1}
=\int_{0}^{\frac\pi2}Q(\phi)\sin\phi \,d\phi,
\end{eqnarray}
since the integral in the denominator equals 1. The integral for $Q$ given by the formula \eqref{eq:Q} can be evaluated exactly.
Since we also made experiments with the date provided by Ghil, we use a numerical approximation by the trapezoidal quadrature rule on the numerical grid on which equation~\eqref{eq:1} is solved.

\paragraph{Temporal variation.}
Temperature data reconstructions over one or more glacial cycles typically show relatively sharp gradients which are superimposed by small-scale fluctuations \cite{Brook2008}. Note that, in this context,  ``sharp'' has to be interpreted in relation to  the considered time ranges of  several hundred thousands of years. Furthermore,  a global warming trend
can be seen in temperature observations in the last  decades.
Our motivation for the use of the parareal method is to capture such multi-scale behavior in climate  models. Thus we
additionally varied $Q$ in the 1-D model by adding some jumps, a linear trend and small-scale random fluctuations, described by a time-dependent function $\Delta Q$. We write
\begin{align}
\label{Qdist0}
Q_{\Delta}(\phi,t)&=Q(\phi)(1 + \Delta Q(t)),
\end{align}
in which the perturbation $\Delta Q(t)$ has the following form:
\begin{align}
\label{Qdist1}
\Delta Q(t)&=\sum_{i=1}^{2}q_{i}\chi_{[a_{i},b_{i}]}(t)+q_{3}\chi_{[a_{3},b_{3}]}(t)(t-a_{3})+q_{4}R(t).
\end{align}
Here $\chi_{[a,b]}$ is the characteristic function on the interval $[a,b]$ and $\{R(t):t\in\{1,\ldots,t_{\rm end}\}\}$ is a set of uniformly distributed random numbers in $[-1,1]$. These were computed once and  then fixed for all experiments. For our tests, we use  the values,
\begin{align}
\label{Qdist2}
\begin{array}{llllll}
q_{1}=0.05,&a_{1}=283,&b_{1}=335,&
q_{2}=- 0.03,&a_{2}=487,&b_{2}=564,\\
q_{3}=0.0001,&a_{3}=700,&b_{3}=1000,&q_{4}=0.05.
\end{array}
\end{align}
In the time integration schemes, the values of $R(t)$ were interpolated linearly  whenever needed at non-integer time instants. Figure \ref{fig:fluct-1} shows the relative deviation from the constant incoming solar radiation $Q$ as function of time.
 \begin{figure}
	\centering
	\includegraphics[scale=0.5]{figuren/solar_forcing}
	\caption{Introduced variation of solar forcing with four jumps, a linear trend after year $700$, and a random variation in the whole time interval $[0,1000]$, compare  \eqref{Qdist0}, \eqref{Qdist1} and \eqref{Qdist2}. \label{fig:fluct-1}}
\end{figure}

For the 0-D model, we tested two choices to handle the time-dependence of $Q_{\Delta}(t)$:
\begin{enumerate}[(1)]
\item a spatially averaged value that takes into account the time-dependence of $Q_{\Delta}(\phi,t)$ (ignoring the random fluctuations);
\[
\overline{Q}_{\Delta}(t)= \int_{0}^{\frac\pi2}Q_{\Delta}(\phi,t)\sin\phi \,d\phi
\]
\item a constant value $\overline{Q}_{\Delta}(t)=\overline{Q}$ for all $t$.
\end{enumerate}

\subsection{Existence and uniqueness of solutions\label{sec:model-exist}}

In this section, we briefly summarize existence and uniqueness results of the two models.
With the temperature-dependent modeling introduced above, the 0-D model takes the form
\begin{eqnarray}
\label{eq:0-D2}
T'=
\frac{(1-\alpha(T))\overline{Q}
 - \epsilon(T)\sigma T^4}{C(T)}.
\end{eqnarray}
With the albedo being the step-function \eqref{eq:albedo-step}, the right-hand side does not fulfill the classical assumptions of the theorems of Peano or Picard-Lindel\"of for existence and uniqueness. On the other hand, we have $\tanh(x)\in[-1,1], x\in\mathbb{R}$, and thus
$C(T)>0$ for $ T\in\mathbb{R}_{\ge0}.$ Hence, the sign of $T'$ can be deduced from the  numerator of the right-hand side of \eqref{eq:0-D2}, compare Figure \ref{fig:ana}: For $T<T_{s}$ we have $T'<0$ and for $T>T_{s}$ we have $T'>0$. Thus any solution of the 0-D initial value problem for \eqref{eq:0-D2} with $T(0)=T_{0}$ will always remain in either $(0,T_{s})$ or $(T_{s},\infty)$, depending in which interval $T_{0}$ lies. In each of these intervals,  the right-hand side of the equation is continuously differentiable  w.r.t. $T$ and thus locally Lipschitz continuous. Thus the Picard-Lindel\"of  Theorem is applicable and gives existence and uniqueness of the solution.
 \begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/rhd-0d}
	\includegraphics[scale=0.475]{figuren/monotone_withC}
	\caption{Left: Numerator of the right-hand side of the 0-D model. Right: Zero-order term of the 1-D model, here for spatially averaged values of $\alpha,Q$. \label{fig:ana}}
\end{figure}

The 1-D model is a semilinear parabolic PDE of the general form
$$
C(T)\,T_{t}+\left(kT_{\phi}\right)_{\phi}+g(\phi,T)=f,
$$
where subscripts denote partial derivatives.
Since $C$ is positive it can be re-written as
$$
T_{t}+\frac{1}{C(T)}\left(kT_{\phi}\right)_{\phi}+\tilde g(\phi,T)=0.$$
 It is not the aim of this paper to perform an analysis of the equation. We just  note that the zero-order term
$$
\tilde g(\phi,T)=\frac{(\alpha(\phi,T)-1)Q(\phi)
 +\epsilon(T)\sigma T^4}{C(T)}
 $$
 shows -- besides the discontinuity at $T=T_{s}$ coming from the albedo function -- a non-monotonicity w.r.t. $T$ for $T\ge T_{s}$, compare Figure \ref{fig:ana}. Classical existence theorems (see e.g. \cite[Section 9.2, Theorem 2]{Eva98}, \cite[Section 7.3, Lemma 5.3]{Tro10}) require either Lipschitz continuity or monotonicity of $\tilde g$. This lack of theoretical results shows the typical analytical properties of nonlinear climate models, even if they are low-dimensional as in this case.

\section{Micro/macro parareal algorithm\label{sec:parareal}}

In this section, we propose the micro/macro parareal algorithm that will be studied further on.  It is a straightforward generalization of the parareal algorithm proposed in \cite{LiMaGa01}; see also \cite{LeLeSa13}, where a similar algorithm was proposed in the context of singularly perturbed ordinary differential equations.  We first introduce the necessary notation (Section~\ref{sec:parareal-notation}), after which we outline the algorithm (Section~\ref{sec:parareal-alg}).

\subsection{Notations\label{sec:parareal-notation}}

We introduce a time discretization $(t_n)_{n=0}^{N}$, with $t_n=n\Delta t$, as well as the numerical approximations $\bT^n\approx \bT(\phi,t_n)$ of the 1-D model \eqref{eq:1} (or, more precisely, its spatial discretization~\eqref{eq:micro}), and $T^n\approx T(t_n)$ of the macroscopic model \eqref{eq:macro}, respectively.

\paragraph{Fine-scale and macro propagators.}
The micro/macro parareal algorithm makes use of a \emph{micro-scale propagator}, that advances the microscopic model \eqref{eq:micro} over a time-step $\Delta t$,
\begin{equation}\label{eq:micro-propagator}
\bT^{n+1}=F_{\Delta t}(\bT^n),
\end{equation}
and, similarly, a \emph{macro propagator} for the macroscopic model
\eqref{eq:macro},
\begin{equation}\label{eq:macro-propagator}
T^{n+1}=G_{\Delta t}(T^n).
\end{equation}
For now, we consider a forward Euler time discretization with time step $\Delta t$ for the macro propagator and with time step $\delta t \ll \Delta t$ for the micro-scale propagator.

\paragraph{Lifting, restriction and projection.} Furthermore, we introduce operators that connect the microscopic and macroscopic levels of description. The \emph{restriction} operator
\begin{equation}
\R:  T \mapsto T = \R(\bT),
\end{equation}
maps a microscopic state to the corresponding macroscopic state.
For the model problem \eqref{eq:micro}-\eqref{eq:macro}, this restriction operator simply becomes the averaging formula~\eqref{eq:spatial_mean}.
For notational convenience, we also introduce the complement of the restriction operator,
 \[
 \R^\perp(\bT):= \dfrac{\bT}{\R(\bT)},
 \]
such that we can write $T(\phi,t)=\R(T(\phi,t))\cdot\R^\perp(T(\phi,t))$.

Conversely, to reconstruct a microscopic state from a given macroscopic state, we distinguish between a \emph{lifting} operator, and a \emph{projection} operator.  A lifting operator $\L$ initializes a microscopic temperature profile $T(\phi,t)$ that is uniquely determined by the given spatially averaged temperature $T(t)$,
\begin{equation}
\L: T \mapsto \bT = \L(T).
\end{equation}
For the model problem \eqref{eq:micro}-\eqref{eq:macro}, for instance, we could choose
\begin{equation}\label{eq:lift}
	T(\phi,t)=\L(T(t)):=T(t)\Psi(\phi),
\end{equation}
with $\Psi(\phi)$ an arbitrary function such that $\R(\Psi(\phi))=1$.
Clearly, one requires $\R \circ \L=\Id$.

In contrast, one may also \emph{match} a ``nearby'' (\emph{prior}) temperature profile $\bT^*(\phi)$ with a desired macroscopic spatially averaged temperature $T$. Then, the result is not uniquely determined by the macroscopic state $T$, but depends also on the prior $\bT^*(\phi)$. We call the resulting reconstruction operator a matching operator,
\begin{equation}
\P: T,\bT^*(\phi) \mapsto \bT = \mathcal{P}(T,\bT^*(\phi)).
\end{equation}
This operator projects a microscopic temperature profile $\bT^*(\phi)$ onto the manifold of microscopic temperature profiles consistent with the macroscopic spatially averaged temperature $T$.
Here, we require $ T = \left(\R \circ \P\right)(T,\bT^*(\phi))$, for any $\bT^*(\phi)$.
Additionally, a so-called \emph{self-consistency} property is of particular importance.
\begin{defn}[Self-consistency] \label{def:self-consistent} A projection operator $\P : T,\bT^*(\phi) \mapsto \bT(\phi) = \mathcal{P}(T,\bT^*(\phi))$, is called self-consistent if, and only if,
\begin{equation}
\forall \bT(\phi): \R(\bT(\phi))=T \;\; \Rightarrow \;\; \bT(\phi)=\P(T,\bT(\phi)).
\end{equation}
\end{defn}
 When this property holds, a microscopic temperature profile is not altered if it is projected onto a macroscopic spatially averaged temperature with which it is already consistent.
As a guideline, the matching should be such that $\mathcal{R}(\bT(\phi))=T$, while requiring $\bT(\phi)$ to be as close to $\bT^*(\phi)$ as possible, in a sense to be made precise for the problem at hand.
For the model problem \eqref{eq:micro}-\eqref{eq:macro}, we choose the matching operator as
\begin{equation}\label{eq:proj}
	\bT(\phi)=\P(T,\bT^*(\phi)):=T\dfrac{\bT^*(\phi)}{\R(\bT^*(\phi))}=T\cdot \R^\perp(\bT^*(\phi)).
\end{equation}

\subsection{Algorithm\label{sec:parareal-alg}}

The parareal algorithm iteratively constructs approximations on the whole time domain.  We denote by $T_k^n$ the approximate solution at time $t_n$, obtained during the $k$-th parareal iteration.
We start from an initial condition $\bT(\phi,0)=T_0$, and create an initial approximation on the whole time interval by using the macro propagator, i.e. we restrict $T_0 = \R(\bT_0)$, and compute
\begin{equation}
T_0^{n+1}=G_{\Delta t}(T_0^n), \qquad T^0_0 = T_0.
\end{equation}
We then lift this macro approximation to the micro scale,
\begin{equation}
\label{eq:lift1}
\bT_0^n = \L(T_0^n).
\end{equation}
We now have an initial approximation of the microscopic solution at each of the time instants $t_n$, $1\le n \le N$.
The parareal iterations then proceed as follows.
\begin{enumerate}[a)]
\item Compute (in parallel) the time propagation at each time instance, using the macro and micro-scale propagators,
\begin{align}
\bar{T}_k^{n+1}&=G_{\Delta t}(T_k^n), \\
\bar{\bT}_k^{n+1}&=F_{\Delta t}(\bT_k^n).\label{eq:micro_prop}
\end{align}
\item Compute the jumps (the difference between the two propagated values) \emph{at the macro level},
\begin{align}\label{eq:parareal-jumps}
J_k^{n+1}= \R(\bar{\bT}^{n+1}_k)-\bar{T}^{n+1}_k.
\end{align}
\item Propagate (serially) the macro jumps over the whole time domain using the macro propagator,
\begin{align}\label{eq:parareal-propagate}
T_{k+1}^{n+1} &= G_{\Delta t}(T^n_{k+1})+J^{n+1}_k,
\end{align}
and reconstruct the micro-scale state from the corrected macro state,
\begin{align}\label{eq:parareal-projection}
\bT_{k+1}^{n+1}=\P(T_{k+1}^{n+1},\bar{\bT}_k^{n+1}).
\end{align}
\end{enumerate}
Remark that the micro-scale state is reconstructed by projecting the intermediate value $\bar{\bT}_k^{n+1}$ onto the corrected macro value $T_{k+1}^{n+1}$, using the matching operator~\eqref{eq:proj}.

\section{Setup of numerical experiments}
\label{sec:numexp}
In this section, we describe the numerical experiments that we performed.  In Section~\ref{sec:time-choices}, we discuss the two options for the time integration of the micro and macro models.  In Section~\ref{sec:macro-choices}, we discuss the two versions of the macro model that will be considered.
The source code is available at \url{https://doi.org/10.5281/zenodo.1287561}.

\subsection{Time integration\label{sec:time-choices}}

\paragraph{Two time discretization methods.} We simulate up to time $t_{\rm end}=1000$, using an initial value of $285\,{\rm K}$ for all computations, in 0-D as well as uniform in space for the 1-D model.
For the micro-scale propagator $F_{\Delta t}$, we use the spatial discretization~\eqref{eq:micro} for the 1-D model, using $I=45$. For the time discretization in the micro 1-D model, we use the built-in \verb#odeint# routine in Python's \verb$scipy$ library, which is the routine \verb$lsoda$ \cite{Pet83} from the \verb$odepack$ library \cite{Hin83}. This method automatically switches between an Adams-Bashforth method (for non-stiff problems) or a backward Differentiation Formula (BDF, for stiff ones). The 1-D model is stiff, thus nearly in all cases and at all time instants the implicit version was chosen by the algorithm.
The \verb#odeint# routine also uses a variable time step, based on a user-prescribed error tolerance. We use $ \varepsilon_{\rm rel}= \varepsilon_{\rm abs}=10^{-6}$  for both relative and absolute tolerances. This means that a  time-step is accepted if the predicted error $e_{n}$ in time instant $t_{n}$ satisfies $$\|e_{n}\|_{\infty}\le \varepsilon_{\rm rel}\|\bT_{n}\|_{\infty}+\varepsilon_{\rm abs}.$$

Note that, because of the time-step adaptation strategy, we have no control over the time steps taken during time integration. In particular, it is very unlikely that the same time step will be taken in subsequent parareal iterations (or in a serial micro simulation over the full time interval). As a consequence, the exactness property of the parareal method becomes very hard to check numerically.  We therefore also perform numerical experiments in which the micro model is integrated using the explicit Euler method with fixed stepsize, chosen small enough (here $0.005$) to give reasonable results for all experiments. Specifically the trajectory of the solution after the second jump ($t> a_{2}$, compare \eqref{Qdist2}), is difficult to obtain with larger stepsizes.

\paragraph{Reference solutions.} We compute wo 1-D micro solutions sequentially using each of the above-described time discretization methods, and use them as reference solutions to compute the differences to the respective parareal solutions. These reference solutions are are denoted by ${\bf T}^{*}(t)$. and shown in figure \ref{fig:refsol}.
\begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/solution_ref_Euler5e-3}
        \includegraphics[scale=0.475]{figuren/solution_ref_odeint1e-4}
	\caption{1-D macro reference solutions obtained by the Euler method with constant stepsize 0.005 (left) and the adaptive  {\tt lsoda} library routine. \label{fig:refsol}}
\end{figure}

\paragraph{Macro propagator.} As the macro-scale propagator $G_{\Delta t}$ for the 0-D model, we again used the explicit Euler method, but now with the constant stepsize  $10.0$.

 \subsection{Two versions of macro model\label{sec:macro-choices}}
Besides varying the microscopic time discretization, we also use two versions of the  0-D macro model. One has a similar temporally varying forcing as the micro model, but without the small-scale random fluctuations. This means we use  $\Delta Q$ as in \eqref{Qdist1} and \eqref{Qdist2}, but with $q_{4}=0$. The other version of the macro model has a constant solar forcing, i.e., $\Delta Q=0$. The  1-D micro reference solution (with varying forcing) obtained with the Euler method (the same as in Figure \ref{fig:refsol} on the left, now averaged in space) is compared in Figure \ref{fig:fluct-2} with these two macro solutions.  As can be seen and is obvious by construction, the macro model with constant solar forcing is an even more coarse and inaccurate approximation of the 1-D micro model. It does not follow neither the jumps nor the linear trend of the forcing (which is catched by the micro model), whereas the macro model with variable forcing does, at least to some extend. Our aim was to see how much this fact influences  the parareal convergence.
 \begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/solution_mean_Euler10_var_forcing_Euler5e-3}
		\includegraphics[scale=0.475]{figuren/solution_mean_Euler10_const_forcing_Euler5e-3}
	\caption{Spatially averaged 1-D micro solution $\overline{\bT}$ with introduced perturbation compared to macro solution $T$ with variable (left) and constant forcing (right). \label{fig:fluct-2}}
\end{figure}

\subsection{Summary of numerical setup}

We thus end up with four configurations, namely using either Euler or an adaptive library routine as micro solver and using either the macro model with or without added perturbation $Q_{\Delta}$ in solar forcing.
For all these four settings, we varied the parareal time-step $\Delta t\in\{10,20,25,40,50,100\}$ (corresponding to $N\in\{100,50,40,25,20,10\}$) and compared the results w.r.t. difference of the parareal solution $\bT_{k}$ in iteration $k$ and the serially computed micro reference solution  $\bT^{*}$. The differences were always evaluated  at all integer time instants $t\in \{1,\ldots,t_{\rm end}\}$.

With the tolerances set as in Section~\ref{sec:time-choices}, the library routine \verb#odeint# and the forward Euler method require approximately the same time for a computation. In this setting, the  0-D macro model with constant forcing requires (approximately and on average (over all performed runs) only a fraction $1/700$ of the micro computation time. The  macro model with temporally varying forcing takes about 2-3 times as much time as macro model with constant forcing. Recall that the 1-D micro model always uses the varying forcing.

\section{Numerical results\label{sec:results}}

\subsection{Results using Euler method as micro propagator}
Figures \ref{fig:conv-1} and \ref{fig:conv-2} show the convergence of the parareal method to the micro reference solution using the Euler method  for the micro model.
\begin{figure}
	\centering
	\includegraphics[scale=0.3]{figuren/conv_Euler10_var_forcing_Euler5e-3}
	\caption{Convergence of parareal method using Euler method for the micro model and  the macro model with  temporally varying solar forcing (depicted in left plot of Figure \ref{fig:fluct-2}) for different numbers $N$ of parareal subintervals. Plotted are maximum norms over space and time.  \label{fig:conv-1}}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[scale=0.3]{figuren/conv_Euler10_const_forcing_Euler5e-3}
	\caption{Same as Figure \ref{fig:conv-1}, but   using  the macro model with constant solar forcing (depicted in right plot of Figure \ref{fig:fluct-2}).\label{fig:conv-2}}
\end{figure}
It can be seen that the differences reach machine precision, see  also Table \ref{table:1}. A reasonable difference level of $10^{-2}$ (talking about temperature which is the range of 200-300 K)  is reached even faster.

To investigate the theoretical gain in computational effort, we assume that  $N$ processors are used and discard the effort for communication and macro propagator. It can be seen from the values in the table that  a  reduction by  a factor of approximately  10    is reached. This refers to a computation up to a difference of the given reasonable tolerance of $10^{-2}$. For computations up to machine precision the factor is approximately 3 to 5.

\begin{table}
\centering
\begin{tabular}{|c|c|r|r|r|r|r|r|}
\hline
Micro&$N$ (number of parareal subintervals)&100&50&40&25&20&10\\
\cline{2-8}
propagator&$\Delta t$ (length of parareal subintervals)&10&20&25&40&50&100\\
\hline
explicit&\multicolumn{7}{|c|}{macro model: temporally varying forcing}\\
\cline{2-8}
Euler&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<eps\}$&17&12&9&7&6&3\\\cline{2-8}
&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<10^{-2}\}$&9&5&4&3&2&1\\\cline{2-8}
&\multicolumn{7}{|c|}{macro model: constant forcing}\\\cline{2-8}
&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<eps\}$&21&11&10&7&6&3\\\cline{2-8}
&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<10^{-2}\}$&9&5&4&3&2&1\\\hline
{\tt lsoda}&\multicolumn{7}{|c|}{macro model: temporally varying forcing}\\\cline{2-8}
library&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<10^{-2}\}$&12&5&4&3&3&1\\\cline{2-8}
routine&\multicolumn{7}{|c|}{macro model: constant forcing}\\\cline{2-8}
&$\min \{k: \|\bT_{k}-\bT^*\|_\infty<10^{-2}\}$&10&5&4&3&3&1\\\hline
\end{tabular}
\caption{Convergence of the micro/macro parareal method to a reasonable tolerance of $10^{-2}$, using the Euler method and the adaptive library method in the micro model. For the Euler method also the convergence  to machine precision (in double precision IEEE arithmetic, $eps\approx 2.2\times 10^{-16}$) is shown.  $\bT_{k}$ denotes the parareal solution in the $k$-th iteration, $\bT^{*}$ the reference solution.
\label{table:1}}
\end{table}

Figure \ref{fig:plots} shows as example the solution obtained by the parareal method  for the setting with $N=100$ subintervals (Euler method in the micro, constant forcing in the macro model). This choice of subintervals showed the worst convergence behavior w.r.t. the needed number of parareal iterations, compare Figure \ref{fig:conv-2} and Table \ref{table:1}. It can be seen that even in this case already after 4 iterations there is no  visible difference between parareal and reference solution in the spatial mean, and also only  a small one in the spatial state.
\begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/solution_mean1_Euler}
	\includegraphics[scale=0.475]{figuren/solution_mean4_Euler}
	\includegraphics[scale=0.475]{figuren/diff_pr_ref1_Euler}
	\includegraphics[scale=0.475]{figuren/diff_pr_ref4_Euler}
	\caption{Top: Spatial mean of parareal and reference and macro solution after first (left) and fourth (right) iteration for $N=100$ subintervals, using Euler method in the micro model and constant forcing in the macro model. Bottom: pointwise relative difference between parareal and reference solution.\label{fig:plots}}
\end{figure}
 The main  difficulty for the method is to capture the steep gradient at the time instant close to $t=b_{2}$ (compare \eqref{Qdist2}), which is the second upward jump in solar  forcing.  This can be seen in Figure \ref{fig:maxerror}. It shows the difference between reference and parareal solution over time. Between  6th and 8th iteration, a significant reduction (by a factor $\approx$ 100) can be seen.
 \begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/diff_max6}
	\includegraphics[scale=0.475]{figuren/diff_max8}
	\caption{Temporal variation of maximum difference  (in space) between parareal and reference  solution after 6th (left) and 8th (right) iteration for $N=100$ subintervals, using Euler method in the micro model and constant forcing in the macro model. \label{fig:maxerror}}
\end{figure}

No big difference can be seen  in the convergence behavior between the use of both macro model versions. Both versions of the   macro model  give the same precision after nearly the same number of parareal ierations.

\subsection{Results using adaptive method as micro propagator}
 Figures \ref{fig:conv-3} and \ref{fig:conv-4} show the convergence of the parareal method using the \verb$lsoda$ library routine for the micro model. Here, the choice of the macro model leads to no relevant differences in any case.
Due to the reasons already mentioned at the beginning of this section,  a reduction of the difference between reference and parareal solution up to machine precision could not be reached in this setting.  We have chosen the tolerances of the method such that a difference of less than $10^{-2}$ was reached. This allows us to compare  the results, shown in Table \ref{table:1}, to those of the Euler method. No big difference in the number of iterations needed to give a maximal error less then $10^{-2}$ in both methods can be seen.

 \begin{figure}
	\centering
	\includegraphics[scale=0.3]{figuren/conv_Euler10_var_forcing_odeint1e-6}
	\caption{Same as Figure \ref{fig:conv-1}, but  using {\tt lsoda} library routine for the micro model, and the macro model with  varying forcing. \label{fig:conv-3}}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[scale=0.3]{figuren/conv_Euler10_const_forcing_odeint1e-6}
	\caption{Same as Figure \ref{fig:conv-3}, but   using  the macro model with constant forcing. \label{fig:conv-4}}
\end{figure}

 Figure \ref{fig:plots2} shows again as example the solution obtained  for the setting with $N=100$ subintervals with constant forcing in the macro model. Also here, this choice of subintervals showed the worst convergence  w.r.t. the needed number of parareal iterations, compare Figure \ref{fig:conv-4}. The behavior is comparable to the results obtained using the Euler method: After four iterations no big difference is visible.

\begin{figure}
	\centering
	\includegraphics[scale=0.475]{figuren/solution_mean1_odeint1e-6}
	\includegraphics[scale=0.475]{figuren/solution_mean4_odeint1e-6}
	\includegraphics[scale=0.475]{figuren/diff_pr_ref1_odeint1e-6}
	\includegraphics[scale=0.475]{figuren/diff_pr_ref4_odeint1e-6}
	\caption{Same as Figure \ref{fig:plots}, but now for the {\tt lsoda} library routine as micro model propagator after the first (left) and the fourth  parareal iteration.\label{fig:plots2}}
\end{figure}

\section{Summary and conclusions \label{sec:discussion}}

 We applied the micro/macro parareal method to a 1-D  climate model with  temporally multi-scale forcing.
As macro model, we use two 0-D versions of the model with spatially averaged coefficients.

We investigated four configurations (micro model with explicit Euler method with constant stepsize or  adaptive library routine on one hand and macro model with varying or constant forcing on the other) and  a variety of number of parareal subintervals.

In all configurations, a reasonable difference to the respective reference solution  (obtained serially with the same time integrator) was obtained after quite a low number of parareal iterations.
The computational gain (discarding macro solver and communication), was about a factor of 10 when using any of the two  micro propagators.
Using the Euler method for the micro propagator, machine precision was reached in less iterations than numbers of used subintervals. For the adaptive library routine, this was not the case due to the different internal time-grids used  on the whole time interval and when restarted on each subinterval.

There was no relevant difference between the two versions of the macro model: The more inaccurate version with constant forcing was able to predict the solution in all cases but one as well as the one with time-varying forcing.

We conclude that the applied micro/macro parareal algorithm is appropriate for this kind of problem, even if the macro model does not include the temporal multi-scale  features. The computational gain
could be even increased when lower (but still reasonable) accuracy requirements are used.
Our results motivate the application of the micro/macro parareal method to more realistic climate models with higher spatial dimension and resolution.

\section*{Acknowledgements}
GS acknowledges the support of the Research Council of the University of Leuven through grant 'PDEOPT', and of the Research Foundation – Flanders (FWO –- Vlaanderen) under grant G.A003.13.
TS acknowledges the support of the German Federal Ministry of Education and Research (BMBF) as part of the Research for Sustainability Initiative (FONA) in the Palmod project,  under grant 01LP1514A.

\bibliographystyle{plain}
\bibliography{literature}

\appendix

\section{Discretization of equation~\eqref{eq:1}\label{sec:discr}}

\subsection{Spatial discretization}
In this appendix, we treat the spatial discretization of the right-hand side of \eqref{eq:1}, which basically means the discretization of the diffusion term  \eqref{eq:diffterm}. Here, we omit the temporal dependency of $T$ in the notation. All indices referring to spatial discretization are denoted as subscripts, whereas indices coming from the temporal discretization will (later on) be denoted as superscripts.

We introduce a spatial grid for $\phi\in[0,\pi/2]$ with stepsize $\Delta\phi=\pi/(2I),I\in\mathbb{N}$.
The gridpoints then are $\phi_{i}=i\Delta\phi,i=0,\ldots,I$. On these points we compute the approximate  solution denoted by
$\bT_{i}\approx T(\phi_{i})$ for every discrete time-step. We also need intermediate points
$\phi_{i+\frac12}=(i+\frac12)\Delta\phi,i=0,\ldots,I-1$.

We discretize the diffusion term \eqref{eq:diffterm}
by applying central finite differences with stepsize $(\Delta\phi)/2$ twice:
We approximate the inner derivative on the intermediate grid points by
\begin{eqnarray}
\label{eq:theta}
\frac{dT}{d \phi}(\phi_{i+\frac12})&\approx&\frac{\bT_{i+1}-\bT_{i}}{\Delta\phi}
=:\Theta_{i+\frac12},\quad i=0,\ldots,I-1.
\end{eqnarray}
Then, the  derivative  of
\begin{eqnarray*}
 F(\phi)
&:=&k(\phi,T(\phi))\sin (\phi)\frac{\partial T}{\partial \phi}(\phi)
\end{eqnarray*}
at the points $\phi_{i}$ is  computed,  again by central finite differences:
\begin{eqnarray}
\label{eq:4}
\left.\nabla\cdot (k\nabla T(\phi))\right|_{\phi=\phi_{i}}=\frac{dF}{d \phi}(\phi_{i})
\approx
\frac{F(\phi_{i+\frac12})-F(\phi_{i-\frac12})}{\Delta\phi},
\quad i=1,\ldots,I-1.
\end{eqnarray}
We need the values of the function $k=k(\phi,T(\phi))$ on the intermediate points:
\begin{eqnarray*}
\kappa_{i+\frac12}&:=&k(\phi_{i+\frac12},T(\phi_{i+\frac12}))
\;\approx \;k\left(\phi_{i+\frac12},\frac{\bT_{i}+\bT_{i+1}}{2}\right),\quad i=0,\ldots,I-1,
\end{eqnarray*}
where the unknown temperature values  at the intermediate points are interpolated linearly.
For the outmost points $\phi_{\frac12},\phi_{I-\frac12}$ this gives, using the homogenous Neumann boundary conditions:
\begin{eqnarray*}
\kappa_{\frac12}
\approx k(\phi_{\frac12},\bT_{1}),&&
\kappa_{I-\frac12}
\approx k(\phi_{I-\frac12},\bT_{I-1}).
\end{eqnarray*}
Eventually, needed values of the coefficient $k$ (if not constant) at the intermediate points $\phi_{i+\frac12}$  have to be interpolated accordingly.
We then obtain for the  terms on the right-hand side of  \eqref{eq:4}:
\begin{eqnarray}
\label{eq:5}
F(\phi_{i+\frac12})
\approx
\kappa_{i+\frac12}\sin (\phi_{i+\frac12})\Theta_{i+\frac12},&&
F(\phi_{i-\frac12})\approx\kappa_{i-\frac12}\sin (\phi_{i-\frac12})\Theta_{i-\frac12}.
\end{eqnarray}

\subsection{Time discretization}
We now obtain a nonlinear nonautonomous system of ODEs
\begin{eqnarray}
\bT'(t)&=&f(t,\bT(t)),\quad t\ge0\label{eq:micro_app}
\end{eqnarray}
with $\bT(t),f(t,\bT(t))\in \mathbb{R}^{I-1}$ and some initial value $\bT(0)\in \mathbb{R}^{I-1}$. The $i$-th component of the function $f$ is given by \eqref{eq:1} using \eqref{eq:theta}, \eqref{eq:4}, and \eqref{eq:5}.
For both models, we re-scaled the time from seconds  to years, i.e. we set
\begin{eqnarray*}
\tilde t&:=& \frac{t}{s_{year}}\quad \text{with } s_{year}=60\times60\times24\times 365 = 3.1536\times10^{7}.
\end{eqnarray*}
The time derivative on the left-hand side of both models can be transformed using the definition $\tilde T(\tilde t):=T(s_{year}\tilde t)$ and the formula
\begin{align*}
\frac{dT}{dt}(t)=\frac{dT}{d(s_{year}\tilde t)}(s_{year}\tilde t)=
\frac{1}{s_{year}}\frac{dT}{d\tilde t}(s_{year}\tilde t)
=\frac{1}{s_{year}}\frac{d\tilde T}{d\tilde t}(\tilde t).
\end{align*}
The dependency on $\phi$ in the 1-D model was suppressed here for simplicity.
Then  both models, \eqref{eq:0-D} and \eqref{eq:1}, retain their formulation (omitting the tildes) when
the respective right-hand sides are multiplied by $s_{year}$.

\end{document}

config: arxiv


TITLE: Does sampling a distribution linearly at random preserve that distribution?
QUESTION [0 upvotes]: I have a set of values, S, that have a gaussian distribution.
I order the set in some way, assigning an integer "position" to each value in S.
I then sample the set S, by choosing values based on a position linearly at random (so the distribution of my position sampling is linear). My sampled values form a new set, Q.
Do the values in the new set, Q, also have a gaussian distribution? Do the relative sizes of the sets Q and S affect this? Does the order I choose to assign integer "positions" to values in S prior to sampling affect this?
I am not sure what concepts from mathematics this even deals with, beyond "stats". I am working with generating synthetic data, and I want to manipulate sets of values while preserving underlying statistical properties (in this case distributions). What should I search for to learn more about this?

REPLY [1 votes]: I would search for the keyword "bootstrap" which is any kind of sampling with replacement to enlarge a dataset. In general I think it preserves the distribution if S and Q are large enough. See a picture of a mixture Gaussian distribution of 9,900 data points on the left and a sample in the middle of size 5000 without replacement, and a bootstrapped sample of size 10,000 on the right with replacement.


TITLE: Combinatorial inequality involving alternating signs
QUESTION [15 upvotes]: I would like to prove the following inequality. It arises from my study of random matrices.
I have verified the inequality for $q\in \{0.01,0.02, \ldots, 0.99\}$ and $1\le n\le 100$.
Let $n$ be any positive integer and $0\le q\le 1$. Then the following inequality is true.
$$\sum_{k=0}^n(-1)^k\binom{n}{k}(q^k-q^n)^n\le (1-q^n)^n.$$
I suspect that the Wilf-Zeilberger method applies here. But I can't get it to work.

REPLY [21 votes]: Mark each box of an $n\times n$ table with probability $q$. By inclusion-exclusion the difference RHS-LHS equals to the probability that there exists a full row (with all boxes marked) but there does not exist a full column: that's because for given $k$ rows the probability that (they are full but no column is full) equals $(q^k-q^n)^n$.


TITLE: Geometry problem involving a cyclic quadrilateral and power of a point theorem?
QUESTION [6 upvotes]: Convex cyclic quadrilateral $ABCD$ are inscribed in circle $O$. $AB,CD$ intersect at $E$, $AD,BC$ intersect at $F$. Diagonals $AC, BD$ intersect at $X$. $M$ is midpoint of $EF$. $Y$ is midpoint of $XM$. circle $Y$ with diameter $XM$ intersects circle $O$ at $P,Q$. Prove that $PY$, $QY$ are tangent to circle $O$
It looks like a pretty interesting problem that could be solved by power of a point theorem, because they are a lot of line segments we can use to compute. But I didn't go very far.

REPLY [0 votes]: As can be seen in the figure:
$YP=YQ$
$OP=OQ$
Therefore triangles YPQ and OPQ are isosceles and consequently YO is perpendicular on PQ. PQ is the common chord of circles centered at Y and O, but YO is not perpendicular on PQ. Now we use the fact that if two Chords of two intersecting circles are perpendicular the line connecting their centers is perpendicular on their common chord. This means PQ can be the common chord of a circle centered at O and another circle passing P, Q and O. Now if Y is on this circle then YO or PQ must be the diameter of this circle. PQ can not be the diameter because  circles centered at O and Y are not equal. This is a fact that if two equal circles intersect such that each one passes the center of the other ,then  the measure of their common chord is equal to the radius of circles and also equal to the measure of line connecting their centers, in this case both common chord and center connector can be the diameter of third circle passing intersections and centers.This is possible if the ABCD is regular.We can conclude that YO is the diameter of a circle passing P, Q, Y and O, hence $\angle YPO=\angle YQO=90^o $ and YP and YQ are tangents on circle centered at O.