Hidden convexity in the heat, linear transport, and Euler's rigid body equations: A computational approach

Uditnarayan KouskiyaDepartment of Civil & Environmental Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, email<EMAIL_ADDRESS> Amit AcharyaDepartment of Civil & Environmental Engineering, and Center for Nonlinear Analysis, Carnegie Mellon University, Pittsburgh, PA 15213, email<EMAIL_ADDRESS>
A finite element based computational scheme is developed and employed to assess  a duality based variational approach to the solution of the linear heat and transport PDE in one space dimension and time, and the nonlinear system of ODEs of Euler for the rotation of a rigid body about a fixed point. The formulation turns initial-(boundary) value problems into degenerate elliptic boundary value problems in (space)-time domains representing the Euler-Lagrange equations of suitably designed dual functionals in each of the above problems. We demonstrate reasonable success in approximating solutions of this range of parabolic, hyperbolic, and ODE primal problems, which includes energy dissipation as well as conservation, by a unified dual strategy lending itself to a variational formulation. The scheme naturally associates a family of dual solutions to a unique primal solution; such `gauge invariance' is demonstrated in our computed solutions of the heat and transport equations, including the case of a transient dual solution corresponding to a steady primal solution of the heat equation. Primal evolution problems with causality are shown to be correctly approximated by non-causal dual problems.

§ INTRODUCTION

The goal of this paper is to test a duality based technique proposed in <cit.> for solving differential equations. We implement the proposed scheme on the linear heat and first-order transport equations posed in bounded domains in one space dimension, and the nonlinear system of ODEs of Euler for the motion of a rigid body with a fixed point (which applies no torque to the body). Euler's system comprises three coupled ODEs for the components of the angular velocity vector of the rigid body on the rotating principal axis frame of its moment of inertia tensor. In all cases, the initial value problems are converted to a boundary value problem in time by the scheme. We compute and assess approximate solutions to the dual equations by the finite element method and verify our procedures with exact or approximate solutions of the primal equations, computed directly from the primal descriptions. Obviously, our effort is not meant to be interpreted as the development of competing algorithms for more standard procedures of approximation for the test problems; instead, it is meant to explore the duality scheme in simple settings where the outcomes are clearly understood. The potential utility of the duality scheme is in generating dual variational principles for fairly general classes of nonlinear systems of algebraic, ordinary, and partial differential equations (the primal equations), see <cit.>. Solutions to the Euler-Lagrange equations of the dual variational principle correspond to those of the primal equations in a well-defined, albeit formal, sense. Specifically, the scheme generates an adapted change of variables which constitutes a mapping of a dual to a primal solution. Such a scheme is most useful for primal equations for which standard definitions of a solution and methods for their approximation are not available, e.g. systems of second-order Hamilton-Jacobi equations. In the context of partial differential equations, our scheme also seems to have the unifying property of producing at worst degenerate elliptic dual problems, with the unusual feature of involving oblique natural boundary conditions <cit.>.

The outline of the paper is as follows: The dual formulation of the heat equation is presented in Sec. 2.  In Sec. 3 a dual formulation of the linear transport equation is developed. Sec. 4 describes the dual formulation of Euler's system of ODE for the motion of a rigid body about a fixed point. Sec. 5 describes the finite element implementations of the above problems and presents computational results. Sec. 6 contains concluding remarks.

§ DUAL FORMULATION OF THE HEAT EQUATION

Following the ideas in <cit.>, we write the heat equation in first order form:

\begin{equation}\label{eq:primal_ht}
\p_x \theta - \pi = 0; \qquad \p_t \theta - \p_x (k \pi) = 0; \qquad (x,t)  \in \Omega = (0,L) \times (0,T).
\end{equation}

Using the Lagrange multiplier fields $p,l$ for the two equations, define

\begin{equation*}
\widehat{S}[ \theta, l, \pi, p] = \int_\Omega dtdx \,   - \theta \p_x p - p \pi - \theta \p_t l + k \pi \p_x l + H(\theta, \pi),
\end{equation*}

where the requirement on the function $H$ is that, defining $\dee: = (\p_x p , p, \p_t l, \p_x l)$, the following equations are solvable for $\theta, \pi$ in terms of $\dee$:

\begin{equation}\label{eq:L_trans}
\begin{aligned}
- \p_x p - \p_t l + \p_\theta H & = 0\\
- p + k \p_x l + \p_\pi H & = 0.
\end{aligned}
\end{equation}

\begin{equation}\label{eq:heat_H}
H(\theta, \pi) = \frac{1}{2} (\theta^2 + \pi^2)
\end{equation}

we obtain the following dual-to-primal mapping (DtP)

\begin{equation}\label{eq:DtP}
\theta = \p_x p + \p_t l; \qquad \pi = p - k \p_x l.
\end{equation}

Now eliminate the fields $\theta, \pi$ in $\widehat{S}$ by substituting these functions for $\theta$ and $\pi$ to define the dual functional

\begin{equation}\label{eq:pre_S}
- \frac{1}{2} \int_\Omega  dtdx \, (\p_x p + \p_t l)^2 + (p - k \p_x l)^2.
\end{equation}

The requirement (<ref>) ensures that the dual E-L system of (<ref>) is simply the primal system (<ref>), in which the primal fields $(\theta, \pi)$ are expressed in terms of the dual fields through the DtP mapping (<ref>).

Denoting $(p,l)$ as the dual fields corresponding to the heat equation (<ref>), the `dual heat PDE' system for the quadratic $H$ chosen in (<ref>) is the following constant coefficient second order system of PDE

\begin{equation}\label{eq:dual_heat_pde}
\begin{aligned}
\p_x^2 p  + \p_{xt} l  + k \p_x l  - p & = 0,  \qquad k > 0 \ \mbox{a constant} \\
k^2 \p_x^2 l   + \p_t^2 l     +  \p_{tx} p  - k \p_x p   & = 0.
\end{aligned}
\end{equation}

The quadratic form arising from its principal part is given by

\begin{equation*}
Q_{dual heat} := (\p_x p + \p_t l)^2 + k^2 (\p_x l)^2
\end{equation*}

which shows that the system is degenerate elliptic since $Q_{dualheat}$ is non-negative for all `gradient' matrices $F$ with $F_{11} = \p_x p, F_{12} = \p_t p, F_{21} = \p_x l, F_{22} = \p_t l$, and ellipticity fails on the rank-one direction $(a, 0) \otimes (0, 1)$, for any value of $a$.

§.§ Formal uniqueness of solutions of the dual heat equation

We utilize uniqueness arguments to infer boundary conditions on space-time domains to define well-set dual problems.

The first variation of the dual functional (<ref>) implies that the first dual equation arises from variations $\delta p$ and the second one from variations $\delta l$. Without imposing any boundary conditions, the following statement holds, assuming the dual E-L equations are satisfied:

\begin{equation}\label{eq:weak_dual_heat}
\begin{aligned}
- &  \int_0^T dt \int_0^L dx \, \p_x \delta p \,(\p_x p + \p_t l) - \int_0^T dt \int_0^L dx \, \delta p \,(p - k \p_x l) \\
- & \int_0^T dt \int_0^L dx \, \p_t \delta l\, ( \p_x p + \p_t l) + \int_0^T dt \int_0^L dx \, \p_x \delta l \, (k (p - k \p_xl)) \\
+  & \int_0^T dt \, \delta p (\p_x p + \p_t l)\big|^L_0 + \int^L_0 dx \, \delta l \, (\p_x p + \p_t l) \big|^T_0 - \int_0^T dt \, \delta l (k (p - k \p_x l)) \big|^L_0 = 0.
\end{aligned}
\end{equation}

Let $(p,l)$ be the difference of two (sufficiently smooth, as required) solutions to the system (<ref>) and let each of the solutions $(p^i, l^i)$, $i = 1,2$ satisfy

* either $(\p_x p^i + \p_t l^i)$ or $p^i$ equal specified values on the left and right boundaries of the space-time domain;

* $(\p_x p^i + \p_t l^i)$ equals a specified initial condition on $\theta$ on the bottom boundary, $t = 0$, of the domain;

* $l^i$ equals an arbitrarily specified function of $x$ (w.l.o.g. chosen to  be 0) on the top boundary of the domain, $t = T$;

* either $(p^i - k \p_x l^i)$ or $l^i$ equal specified values on the left and right boundaries of the space-time domain.

* on any point of the left and right boundaries, either $(\p_x p^i + \p_t l^i)$ or $(p^i - k \p_x l^i)$ are specified.

Then, choosing $\delta p = p$ and $\delta l = l$ in (<ref>) results in all the boundary terms vanishing and

\begin{equation*}
- \int_0^T dt \int_0^L dx \, (\p_x p + \p_t l)^2 + (p - k \p_xl)^2 = 0
\end{equation*}

which implies that

\begin{equation}\label{eq:heat_unique}
\p_x p + \p_t l = 0 \ \mbox{and} \  p - k \p_xl = 0  \Longrightarrow \p_x p = k \p_{xx} l \Longrightarrow  \p_tl + k \p_{xx} l = 0   \qquad \mbox{in} \ (0,L) \times (0,T),
\end{equation}

which is the backward heat equation in $l$ for $k > 0$. Since at any left or right boundary point either $l^i$ or $(p^i - k \p_x l^i)$ is specified, the latter implying that at such points $p^i$ is specified (by b.c.s 4.,5.,1. above), (<ref>) (and 4.) imply that either $l = 0$ or $\p_x l = 0$ at such points. Then, defining $\Lambda(t) := \frac{1}{2} \int_0^L dx \,  l^2 (x,t) \geq 0$, (<ref>) gives $\p_t \Lambda(t) \geq 0 \Longrightarrow \Lambda(T) - \Lambda(t) \geq 0$ for all $t \in [0,T]$. But by the `top' boundary condition $\Lambda(T) = 0$, so that $0 \leq \Lambda(t) \leq 0$, and we have $l = 0$ almost everywhere in the space-time domain, and since $p - k \p_xl = 0$ in the same domain, we have uniqueness.

§.§ Weak formulation for the dual heat equation

For the sake of definiteness, let us assume that the primal heat equation is posed with the following initial and boundary conditions:

\begin{equation}\label{eq:primal_bc}
\theta(x,0) = \theta_0(x); \qquad \theta(0,t) = \theta_l(t); \qquad \pi(L,t) = \pi_r(t).
\end{equation}

Then the following weak form for the problem suffices to compute unique (approximate) solutions to the dual problem, from which the unique primal solution to (<ref>) and (<ref>) can be approximated through the DtP mapping (for arbitrarily specified $T$):

\begin{equation}\label{eq:weak_heat}
\begin{aligned}
-  &  \int_0^T dt \int_0^L dx \, \p_x \delta p \,(\p_x p + \p_t l) - \int_0^T dt \int_0^L dx \, \delta p \,(p - k \p_x l) \\
- & \int_0^T dt \int_0^L dx \, \p_t \delta l\, ( \p_x p + \p_t l) + \int_0^T dt \int_0^L dx \, \p_x \delta l \, (k (p - k \p_xl)) \\
- & \int_0^T dt \, \delta p(0,t) \theta_l(t) - \int^L_0 dx \, \delta l(x,0) \, \theta_0(x) - \int_0^T dt \, \delta l(L,t) (k \pi_r(t)) = 0; \\
& \delta l (x,T) = 0; \qquad \delta p (L,t) = 0; \qquad \delta l (0,t) = 0;\\
& l(x,T) = l_T(x); \qquad p(L,t) = p_r(t); \qquad l(0,t) = l_l(t)
\end{aligned}
\end{equation}

with  $l_T(x), p_r(t), l_l(t)$ being arbitrarily specified continuous functions subject to $l_T(0) = l_l(T)$. This weak statement corresponds to the first variation of the dual functional

\begin{equation*}
\begin{aligned}
S[l,p] & = - \frac{1}{2} \int_\Omega dtdx  \ (\p_x p + \p_t l)^2 + (p - k \p_x l)^2 \\
& \quad - \int_0^T dt \ p(0,t) \theta_l(t)  - \int^L_0 dx \ l(x,0) \, \theta_0(x)  - \int_0^T dt \ l(L,t) (k \pi_r(t)).
\end{aligned}
\end{equation*}

If a Dirichlet boundary condition of the form $\theta(L,t)=\theta_r(t)$ is applied to the right boundary instead of $\pi(L,t)=\pi_r(t)$, the formulation and weak form of the problem remain the same. The only difference is that instead of $p(L,t)$, the dual field $l(L,t)$ is specified arbitrarily as $l_r(t)$ with $\delta l (L,t) = 0$ and the boundary term in (<ref>) for $\pi_r$ on the right boundary is replaced by the boundary term produced based on $\theta_r(t)$, following what is done for the Dirichlet condition for the left boundary.

§ DUAL FORMULATION OF THE LINEAR TRANSPORT EQUATION

An exactly similar procedure as for the heat equation can be followed for the linear, constant coefficient wave equation

\begin{equation}\label{eq:primal_wave}
\p_t u + c \, \p_x u = 0 \  \mbox{in} \  (0,L) \times (0,T), \quad c > 0; \quad u(0,t) = u_l(t); \quad u(x, 0) = u_0(x)
\end{equation}

(with the functions $u_l, u_0$ being specified) to obtain its dual functional with natural boundary conditions

\begin{equation}\label{eq:wave_dual}
S[\lambda] = - \int_0^T dt \int_0^L dx \, \frac{1}{2} (\p_t \lambda + c \p_x \lambda)^2 - \int_0^L dx \, \lambda(x, 0) u_0(x) - \int_0^T dt \, \lambda(0,t) u_l(t) c
\end{equation}

and dual PDE

\begin{equation}\label{eq:dual_wave_EL}
\p_t (\p_t \lambda + c \p_x \lambda) + c \, \p_x (\p_t \lambda + c \p_x \lambda) = 0 = \p_t^2 \lambda + 2c \, \p_{tx} \lambda + c^2 \, \p_x^2 \lambda.
\end{equation}

\[H(u) = \frac{1}{2} u^2
\]

and the DtP mapping emerges as

\begin{equation}\label{eq:DtP_wave}
u = \p_t \lambda + c \p_x \lambda.
\end{equation}

We supplement the dual PDE with the essential boundary conditions

\begin{equation*}
\lambda(x, T) = 0; \qquad \lambda(L, t) = 0,
\end{equation*}

where the specifications of these functions are `arbitrary,' as concerns recovery of the (formally) unique primal solution by the scheme.

As before, the quadratic form in gradients is

\begin{equation*}
Q_{dual1-wave} = (\p_t \lambda + c \p_x \lambda)^2,
\end{equation*}

and this is again positive semi-definite for all gradient vectors $(\p_x \lambda, \p_t \lambda)$ with ellipticity failing for the vector direction $a \, (1, -c)$, for any value of $a$.

\[
grad ({A \, grad \lambda}) =
\begin{bmatrix} \p_x
\left( \begin{bmatrix}
c^2 & c \\
c & 1
\end{bmatrix}
\begin{bmatrix}
\p_x \lambda \\
\p_t \lambda
\end{bmatrix} \right)
& &
\p_t
\left( \begin{bmatrix}
c^2 & c \\
c & 1
\end{bmatrix}
\begin{bmatrix}
\p_x \lambda \\
\p_t \lambda
\end{bmatrix} \right)
\end{bmatrix}
\]

so that $div (A \, grad \lambda)  = 0$ is (<ref>), and ellipticity is equivalent to $n_i A_{ij} n_j \neq 0$ for all $n \neq 0$, which does not hold for the $n$ given.

Weak solutions of this dual system (<ref>) obviously exist, generated explicitly by first integrating the primal equation and boundary and initial conditions by the method of characteristics to obtain, say $u(x,t)$, and then again using the method of characteristics utilizing the DtP mapping, $u(x,t) = \p_t \lambda + c \p_x \lambda$ as the governing equation, and using the `top' and `right' boundary conditions on $\lambda$. Clearly, the smoothness of the solution depends on the smoothness of the specified data $u_l, u_0$ (with the arbitrary data on $\lambda$ always chosen to be smooth, w.l.o.g.).

The degenerate elliptic dual problem allows gradient-discontinuities which translate, through the DtP mapping, to strong discontinuities in the primal $u$. Curiously, the argument in the previous paragraph indicates that  functions $\lambda$ that are discontinuous across characteristics can be constructed as solutions of the method of characteristics, for appropriately chosen data on $u$. Whether such functions can qualify as solutions to the E-L equations of the functional in (<ref>) is a question that remains to be answered.

§.§ Formal uniqueness of solutions of the dual linear transport equation

Let $\lambda$ be the difference of two solutions $\lambda^i$, $i = 1,2$ of (<ref>), each satisfying the following boundary conditions:

* $(\p_t \lambda^i + c \, \p_x \lambda^i)$ equals the specified functions $u_l$ on $x = 0$ and $u_0$ on $t = 0$;

* $\lambda^i(x,T) = \lambda_T(x)$, $\lambda^i(L,t) = \lambda_r(t)$, where $\lambda_T, \lambda_r$ are specified functions.

Then, a direct consequence of the `energy method' (i.e., multiply by $\lambda$, integrate by parts, apply boundary conditions on the space-time domain) is that

\[
\p_t \lambda + c \, \p_x \lambda = 0 \qquad \mbox{in}  \ (0,L) \times (0,T).
\]

Applying the energy method again to the statement above,  and defining $\Lambda(t) := \frac{1}{2} \int_0^L dx \,  \lambda^2 (x,t) \geq 0$, one obtains $\p_t \Lambda (t) + \frac{1}{2} c (\lambda^2(L,t) - \lambda^2(0,t)) = 0$ which, on applying the b.c. at $x = L$, integrating in time from $t$ to $T$, and applying the b.c. at $t = T$ gives $  0 \geq - \Lambda(t) = \frac{1}{2} \int^T_t dt \, c \, \lambda^2(0,t) \geq 0$ since $c > 0$, and we have $\lambda(x,t) = 0$ almost everywhere in $(0,L) \times (0,T)$ and uniqueness.

§.§ Weak formulation for the dual linear transport equation

For the primal first-order wave equation along with the applied boundary conditions stated in (<ref>), we obtain the weak form described as follows. For any field $\delta \lambda$ satisfying the conditions stated below, we intend to find the dual field $\lambda$ which satisfies the following equations:

\begin{equation} \label{eq:weak_dual_wave}
\begin{aligned}
-  &  \int_0^T dt \int_0^L dx \, \p_t \delta \lambda\, (\lambda_t + c \lambda_x) - \int_0^T dt \int_0^L dx \,
c\,\p_x \delta\lambda\, (\lambda_t + c \lambda_x)\\
& - \int_0^T dt \, \lambda(0,t) u_l(t) - \int_0^L dx \, \lambda(x,0) u_0(x) = 0; \\
& \delta \lambda(x,T) = 0; \qquad \delta \lambda(L,t) = 0; \\
& \lambda(x,T) = \lambda_T(T); \qquad
\lambda(L,t) = \lambda_r(t),
\end{aligned}
\end{equation}

where $\lambda_T(\cdot)$ and $\lambda_r(\cdot)$ are arbitrarily specified functions satisfying $\lambda_T(L) = \lambda_r(T)$.
The dual scheme guarantees that the solution to (<ref>) implies the solution to the set of equations (<ref>). We make use of the above weak form to compute an approximate solution for the dual field and utilize the DtP mapping (<ref>) to obtain the corresponding field for the primal problem i.e. the first-order wave equation. The weak form corresponds to the first variation of the dual functional (<ref>).

§ DUAL VARIATIONAL PRINCIPLE FOR EULER'S SYSTEM FOR MOTION OF A RIGID BODY, WITH AND WITHOUT VISCOSITY

Following the derivation in Appendix <ref> of Euler's system for the components of the angular velocity of a rigid body rotating about a fixed point under the possible action of a viscous damping force proportional to its angular momentum, we consider the following primal system of ODEs

\begin{equation}\label{eq:Euler}
I_i \, \dot{\omega}_i + c_i \, \omega_{i+1} \omega_{i+2} + \nu I_i \,\omega_i = 0, \qquad i = 1,2,3, \ \ c_i = I_{i+2} - I_{i+1}, \ \ \nu \geq 0,
\end{equation}

where a superscript dot represents a derivative w.r.t time and all subscript indices in this Section are modulo 3. The system is subject to the initial condition

\[
\omega_i(0) = \omega_i^0.
\]

Using the auxiliary potential

\[
H(\omega, t) = \frac{1}{2} \sum_{i=1}^3 (\omega_i - \tilde{\omega}_i(t))^2,
\]

where $\tilde{\omega_i}$ is an arbitrarily chosen `base state' subject only to facilitating a dual solution <cit.>, denoting the dual functions by $\lambda_i$, and using our standard protocol for generating a dual variational principle, one obtains the DtP mapping for this problem as

\begin{equation}\label{eq:euler_DtP}
\omega_i - \tilde{\omega}_i =  \sum_j \mathbb{K}^{-1}_{ij} \left(I_j \dot{\lambda}_j - \nu I_j \lambda_j \right),
\end{equation}

\begin{equation*}
\mathbb{K} =
\begin{bmatrix}
a & c_3 \lambda_3 & c_2 \lambda_2 \\
c_3 \lambda_3 & a & c_1 \lambda_1 \\
c_2 \lambda_2 & c_1 \lambda_1 & a
\end{bmatrix},
\end{equation*}

where $a$ is a constant of arbitrary magnitude and we choose $a \gg 1$.

The dual problem is obtained by substituting the change of variables (<ref>) into the primal system (<ref>) appended with the boundary conditions

\begin{equation*}
\label{eq:dual_euler_bc}
\omega^0_i(0) = \tilde{\omega}_i(0) + \sum_j \mathbb{K}^{-1}_{ij}\left (I_j \dot{\lambda}_j (0) - \nu I_j \lambda_j (0) \right); \qquad \lambda_i(T) = 0.
\end{equation*}

The dual problem so obtained is, by design, the E-L equations of the dual functional

\begin{equation}\label{eq:euler_dual_functional}
\begin{aligned}
S[\lambda] &:= \int^T_0  dt \, \sum_i \left( - I_i \, \omega_i \, \dot{\lambda}_i + \lambda_i (c_i \, \omega_{i+1} \omega_{i+2} + \nu I_i \omega_i ) \right) \\
& \qquad - \sum_i I_i \, \omega^0_i \lambda_i (0)  + H(\omega, t),
\end{aligned}
\end{equation}

with all occurrences of $\omega_i$ replaced by its DtP representation (<ref>) in terms of $(\tilde{\omega}_i, \dot{\lambda}, \lambda)$.

It is important to note that the dual variational problem defined by the choice of base state $\tilde{\omega}_i$ given by a solution of the primal ivp guarantees the existence of a solution to the dual problem given by $\lambda_i = 0$. Thus, it may be expected that for base states that are `close' to actual solutions of the primal problem, the dual problem has a solution.

§.§ Weak formulation of the dual Euler rigid body system

For the coupled system of equations (<ref>), we employ the following weak form to obtain the approximate solution for the dual fields and subsequently use the DtP mapping (<ref>) to evaluate the primal field variables:

\begin{equation}\label{eq:weak_euler}
\begin{gathered}
\int_0^T dt \, \sum_i (-I_i \, \omega_i \, \dot{ \delta\lambda_i} + c_i \, \omega_{i+1} \, \omega_{i+2} \, \delta\lambda_i + \nu I_i \, \omega_i \, \delta\lambda_i) - \sum_i I_i\,\omega_i^0\,\delta\lambda_i(0) = 0; \\
\delta\lambda_i(T) = 0; \qquad \lambda_i(T) = \lambda_i^T,
\end{gathered}
\end{equation}

where $\lambda_i^T$ are specified arbitrarily and each of the $\omega_i$ are understood to be expressed in terms of the dual fields and their derivatives using (<ref>). As usual, the solution to the above weak form defines an (approximate) solution to the set of equations (<ref>) through the DtP mapping. This weak form corresponds to the first variation of the functional (<ref>).

§ RESULTS

The following sections describe examples for each of the problems developed in the previous sections. For the heat and transport equations, we solve the dual problems as space-time boundary value problems (bvp), while for Euler's system of ODE (with and without damping), we solve the dual problem as a two-point bvp in time. We use the Finite Element (FE) method to discretize all problems. In the case of heat/transport equations, we define the space-time domain $\Omega=\{(x,t):x \in (0,L), t \in (0,T)\}$, while for Euler's system of ODE, the time domain is defined as $\Omega=\{t:t \in (0,T)\}$. A linear span of globally continuous, piecewise smooth finite element shape functions corresponding to a FE mesh for $\Omega$ is used to achieve this discretization. These shape functions are represented by $N^{(\cdot)}$, where $(\cdot)$ denotes the index of the node under consideration. For the following examples, the error for any field with respect to a reference field within the domain  is given by:

\[
\% \mbox{ error} = \left|\frac{\mbox{dual solution}-{\mbox{reference solution}}}{\mbox{reference solution}}\right| \times 100,
\]

where the `dual solution' above refers to the primal field obtained by the DtP mapping from the corresponding dual approximation. In most cases, this measure is evaluated pointwise. The measure is not implementable at points where the reference field approaches zero. In such cases, alternative measures are used, as described in the relevant sections.

In this work, we use $C^0$ FE shape functions to approximate the dual fields. Hence, the direct evaluation of primal fields through the DtP mapping exhibit discontinuities in the domain in general, even when approximating continuous primal solutions. This is because the DtP mappings involve derivatives of the dual fields. To deal with this feature, our algorithm introduces the following device - essentially, a commonly used $L^2$ projection - to posit the corresponding primal solution, which we then compare with the exact primal solutions to test the accuracy of the overall scheme (including the projection). Let $u$ represent such a discontinuous primal approximation, determined from the dual solution. To achieve a $C^0$ continuous approximation in the domain $\Omega\subset\mathbb{R}^2 \mbox{ or } \mathbb{R}$ (depending on the problem), we employ the following method: Let $u_h$ represent the projection of $u$ onto a space $V_h$ formed by the linear span of globally continuous, piecewise smooth finite element shape functions corresponding to a FE mesh for $\Omega$. We enforce the following condition upon $u_h$:

\begin{equation*}
u_h = \operatorname*{arg\,min}_{v\in V_h}  \int_\Omega \frac{1}{2} |u-v|^2 \, d\Omega.\label{eq:weak_L2}
\end{equation*}

Let $N^A$ represent the basis functions associated with any node with an index $A$. The discrete version of the optimality condition of the above statement is:

\begin{equation}\label{eq:L2_project}
\sum_{A=1}^N\sum_{B=1}^N\delta u^A_h\biggl(\int_\Omega N^A\, N^B \,d\Omega \biggl) u_h^B=  \sum_{A=1}^N\delta u^A_h \int_\Omega N^A\, u\,d\Omega
\end{equation}

where $u_h^B$ denotes the nodal value at node $B$ of the sought continuous projection of $u$, and  $\delta u_h := \delta u^A_h N^A$ is a test function. In solving (<ref>) we impose any known data, e.g. initial and boundary conditions of the primal problem, as known function values $u^A_h$, at corresponding nodes $A$ where the data is known, with $\delta u^A_h = 0$. Using the arbitrariness of the remaining `free' $\delta u^A_h$, we obtain a system of linear equations to solve for the unconstrained nodal values $u_h^B, B \in \{1, \ldots, N\}$.

For several examples presented in the following sections, partitioning the domain into smaller domains (stages) along the time direction simplifies computation and aids in refinement. Each stage, indexed by $s$, uses initial conditions based on the previous stage's results, solving the full problem up to the required final time by concatenating primal fields obtained by DtP mappings for successive stages along increasing time.
We will refer to this process as `time-slicing'. The details of this method are explained in the relevant sections. As also explained later, higher errors in the primal solution occur near the final time of a stage since high gradients in the dual solution have to be resolved, in general, due to the imposition of an arbitrarily fixed final-time boundary condition (this phenomenon also occurs at the outflow boundary of the linear transport equation for the same reason). However,  a portion of the results near the final time of any stage is, and can be, discarded without loss of generality. Thus, the domain is extended beyond the required final time of any specific problem and the obtained results are retained only up to the specified final time.

In the following, we use summation over repeated vector or tensor indices, as well as FE formulation related summations, with summation performed over the appropriate index ranges, unless explicitly mentioned otherwise.

§.§ Heat equation

In order to solve the dual heat equation the weak form (<ref>) is employed. Let $N^A$ represent the shape function corresponding to any node $A$. Using the discrete approximations

\begin{equation*}
\begin{gathered}
p(x,t) = d_1^A N^A(x,t);  \qquad  l(x,t)=d_2^AN^A(x,t), \\
\end{gathered}
\end{equation*}

where $N$ is the total number of nodes in the FE mesh,
we obtain the following set of equations:

\begin{equation*}\label{eq:heat_discrete}
K_{ij}^{AB} d_j^B = R_i^A;
\end{equation*}

\begin{equation*}
\begin{gathered}
\qquad K_{11}^{AB} =
\int_0^T dt \int_0^L dx \Bigl(\,- \p_x N^A(x,t) \,\p_xN^B(x,t) - \,  N^A(x,t)\, N^B(x,t)\Bigl);\\
\qquad K_{12}^{AB} = \int_0^T dt \int_0^L dx \Bigl(\,- \p_x N^A(x,t) \,\p_tN^B(x,t) + \, k N^A(x,t) \,\p_x N^B(x,t)\Bigl);  \\
\qquad K_{21}^{AB} =  \int_0^T dt \int_0^L dx \Bigl(\,- \p_t N^A(x,t) \,\p_xN^B(x,t) + \, k \, \p_x N^A(x,t)\, N^B(x,t)\Bigl);  \\
\qquad K_{22}^{AB} =  \int_0^T dt \int_0^L dx \Bigl(\,- \p_t N^A(x,t) \,\p_tN^B(x,t) - \, k^2 \p_x N^A(x,t) \,\p_x N^B(x,t)\Bigl) ;  \\
R_1^A = \int_0^T dt \, N^A(0,t) \, \theta_l(t); \qquad
R_2^A = \int_0^T dt \, N^A(L,t) \, k\pi_r(t)
+ \int_0^L dx \, N^A(x,0) \, \theta_0(x).
\end{gathered}
\end{equation*}

Additionally, we apply the Dirichlet Boundary conditions as

\begin{equation*}
\begin{gathered}
l(0,t) = l_l(t);\qquad p(L,t) = p_r(t); \qquad l(x,T) = l_T(x).
\end{gathered}
\end{equation*}

If the boundary condition at $x=1$ is altered from $\pi(L,t)=\pi_r(t)$ to $\theta(L,t)=\theta_r(t)$ in (<ref>), $R_i^A$ gets modified in the obvious manner following the left boundary condition:

\begin{equation*}
R_1^A = \int_0^T dt \, \bigl(N^A(0,t) \, \theta_l(t)
-  \, N^A(L,t) \, \theta_r(t)\bigl); \qquad
R_2^A =
\int_0^L dx \, N^A(x,0) \, \theta_0(x),
\end{equation*}

and the Dirichlet boundary condition $p(L,t)=p_r(t)$ changes to  $l(L,t)=l_r(t)$ on the right boundary. Each of the integrals in the above expressions are evaluated using a two-point Gauss quadrature scheme in each direction. Once the discretized dual fields $p(x,t)$ and $l(x,t)$ are obtained, the primal field $\theta(x,t)$ is evaluated at the Gauss points and the scheme explained in the preamble of Sec. <ref> is used to generate its continuous projection.

§.§.§ Steady solution of the Heat equation

A steady solution of the heat equation (<ref>) for $k = 1$ given by

\begin{equation*}\label{eq:heat_p1_exact_sol}
\theta(x,t) = 3x+1 %%\qquad \forall t\in[0,\infty)
\end{equation*}

is considered on a domain $\Omega$ with $L=1$, satisfying the following initial and boundary conditions:

\begin{equation}\label{eq:heat_p1_primal_bcs}
\theta_0(x) = 3x+1; \qquad \theta_l(t) = 1; \qquad \theta_r(t) = 4.
\end{equation}

In order to evaluate our dual scheme, we approximate the solution to the current problem by applying the boundary conditions

\begin{equation}\label{eq:dual_heat_hom_l_p}
l_l(t)=l_r(t)=0; \qquad l_T(x)=0
\end{equation}

in the weak form (<ref>). The Dirichlet boundary conditions at $x=0$ and $x=L$ in the primal problem are transformed into oblique natural boundary conditions in the dual formulation. For a final time of $T=1.1$ and a mesh of $100\times110$ elements in the space-time domain, the field $\theta(x,t)$ and the error obtained with respect to the exact solution are shown in Fig. <ref> and Fig. <ref>,  respectively.

$\%$ error in $\theta(x,t)$

$\%$ error in $\theta(x,t)$
up to $T=1$

DtP mapping generated primal field $\theta$ for the steady state heat equation. The mesh is $100\times110$, $T=1.1$ and $k=1$.

As evident from <ref>, while the errors are minimal, they are higher at at $T=1.1$. Given the arbitrariness in the choice of imposition of b.c.s for the dual problem at the final time and the latter's value, given a final time of interest for the primal problem, we can always choose $T$ appropriately and discard a reasonable portion of data from the dual solution near the final time and focus on results within the remaining domain. For illustration, Fig. <ref> displays the results until $T=1$, which has a significantly lower maximum error compared to until $T=1.1$.

Transient Dual field $l(x,t)$

Transient Dual field $p(x,t)$

Steady Dual field $l(x,t)$

Steady Dual field $p(x,t)$

Dual fields for steady state heat equation. (a) and (b) represent dual fields when $l_l(t)=l_r(t)=l_T(x)=0$ are imposed. (c) and (d) represent the dual fields when the boundary condition set (<ref>) is imposed.

The dual fields for this problem are presented in Fig. <ref>-<ref>. Although these fields are transient in nature, the corresponding primal solution is steady in time.
As explained in <cit.>, when the primal problem has a unique solution, any solution to the dual problem, regardless of the choice of $H$ and the space-time b.cs imposed on the dual problem beyond the imposition of the primal set, must recover the unique primal solution through the DtP mapping.  As an illustration of this invariance of the DtP mapping generated primal solution from two different dual solutions, we set the following Dirichlet Boundary condition on the dual fields:

\begin{equation}\label{eq:heat_p1_steady_bc}
l_l(t)=0;\qquad l_r(t)=-\frac{13}{6};\qquad l_T(x)= \frac{x^3}{3} + \frac{x^2}{2}-3x.
\end{equation}

The computed dual solution with these b.c.s are shown in Fig. <ref>-<ref>. It can be easily verified that

\begin{equation}
\label{eq:heat_p1_steady_duals}
p(x,t) = \frac{3x^2}{2}+2; \qquad l(x,t) =\frac{1}{k}\left(\frac{x^3}{3} + \frac{x^2}{2}-3x\right)
\end{equation}

satisfies the boundary conditions (<ref>) and the dual heat equation (<ref>). Furthermore, on application of the DtP mapping (<ref>), it results in the expression $\theta(x,t) = 3x+1$ which is the unique solution to the primal problem, as must be. We note that dual solution (<ref>)  corresponding to b.cs (<ref>) is steady, unlike the dual solution corresponding to the b.cs (<ref>). This fact is also evident in the simulation results Fig. <ref>-<ref> (except close to $T = 1$, without any essential loss of generality, as explained earlier).

§.§.§ Transient solution of the Heat equation by duality

The purpose of the following example is to demonstrate an approximation of a time-dependent solution of the heat equation for the following combination of Dirichlet and Neumann boundary conditions:

\begin{equation*}
\theta_l(t) = 1; \qquad \pi_r(t) = 0,
\end{equation*}

and the initial condition

$$\theta_0(x) = \sin{\frac{\pi x}{2}} + 1.$$

For $L=1$, the exact solution in this case is given by:

$$\theta(x,t) = \sin\left(\frac{\pi x}{2}\right) e^{-\frac{\pi^2 kt}{4}}+1. $$

Applying $l_T(x)=p_r(t)=l_l(t) =0$ as the boundary conditions along with $k=0.2$, we utilize (<ref>) and simulate the problem up to $T=1.1$. The results for a mesh of $100\times110$ are presented in Fig. <ref> and Fig. <ref>. The error plot up to a time $T=1$ is also shown in Fig. <ref> which illustrates that the error in the interior of the domain is smaller relative to that near the top boundary.

$\%$ error in $\theta(x,t)$

$\%$ error in $\theta(x,t)$
upto $T=1$

DtP mapping generated primal field $\theta$ for the transient heat equation. Mesh is $100\times110$, $T=1.1$ and $k=0.2$.

§.§.§ Discontinuous initial data for the Heat equation

As a final example involving the heat equation, we investigate how the dual methodology deals with any problem which involves a discontinuity in initial conditions, especially since primal Dirichlet b.cs and initial conditions get imposed as natural boundary conditions in the dual problem. Consider the following set of initial and boundary conditions for $L=1$:

\begin{equation}\label{eq:heat_p3_initial_condition}
\theta_0(x) = \begin{cases}
10+2x &  \mbox{for } 0 \leq x< 0.5 \\
8+2x &  \mbox{for }  0.5 < x \leq 1
\end{cases};
\qquad \theta_l(t) = 10;
\qquad \theta_r(t) = 10.
\end{equation}

In this case, separation of variables along with the application of Fourier series yields the solution

\begin{equation}\label{heat_p3_exact}
\theta(x,t) = 10 + 2\sum_{m=1}^\infty \frac{(-1)^{m+1}}{\pi m} \sin{(2\pi mx)}e^{-(2m)^2\pi^2kt}\qquad \text{ for }0\leq x\leq1.
\end{equation}

It has been verified that the series converges in the $L^2$ sense to a value within machine precision on summing the first $10^5$ terms and beyond, and we retain $10^5$ terms for generating the exact solution. Additionally, at the point of discontinuity, the function value approximates to  the average of the values obtained from the left and the right of the jump. Consequently, for the initial condition in the simulation, we provide a value of $10$, the average jump, at $x=0.5$. The results are presented in Fig. <ref>.

$\%$ Error for $\theta(x,t)$

Line plots for $\theta(x,t)$

DtP mapping generated primal field $\theta$ for a discontinuous initial condition. Mesh is $200\times50$, $T=0.25$ and $k=0.1$. The line plots for $\theta(x,t)$ close to the initial time are illustrated in (c).

Even though Fig. <ref> provides a satisfactory result in attempting to capture a discontinuous problem, achieving convergence to (<ref>) on a finite mesh through refinement is practically not possible.
This is mainly due to the fact that the linear span of continuous FE basis functions can represent a discontinuity only in the limit of infinite mesh refinement. To test our code when convergence on a finite mesh is expected for a sharp transition, we modify the initial condition <ref> to the following:

\begin{equation*}
\theta_0(x) = \begin{cases}
\beta+2x &  \mbox{for } 0 \leq x< 0.5-\veps  \\
kx+c    &  \mbox{for } 0.5 - \veps \leq x \leq 0.5+ \veps \\
\beta-2+2x &  \mbox{for }  0.5+\veps < x \leq 1
\end{cases};    \qquad k = \frac{2\veps-1}{\veps}; \qquad c = \beta-\frac{2\veps-1}{2\veps}.
\end{equation*}

With $\beta=10$, the above-stated initial condition tends to (<ref>) for $\veps \to 0$.
To obtain the dual results, we set $\veps = 0.01$, $l_l(t)=l_r(t)=0, l_T(x)=0$, $k =0.1$, and perform the simulation on a $200\times25$ and a $400\times50$ mesh up to $T=0.125$. For this case, no data is discarded near the final time. The results for these simulations are shown in Fig. <ref>. The solution is assessed by comparing it to the analytical result, which is obtained by considering the Fourier series of the modified initial condition, similar to the discontinuous initial condition. This analytical result has been stated in Appendix <ref>. Comparing Fig. <ref> with Fig. <ref>, it is evident that refinement improves the results.

$\%$ Error for $\theta(x,t)$

Line plot for $\theta(x,t)$

$\%$ Error for $\theta(x,t)$

Line plots for $\theta(x,t)$

[short 1]
Errors and line plots for DtP mapping generated primal field $\theta$ with a smoothed jump in the initial condition.  $\beta=10$. (a) and (b) are obtained on a mesh of $200\times25$. (c) and (d)  are obtained on a mesh of $400\times50$.

The parameter $\beta$ determines the constant factor by which the function profile is shifted with respect to $x$-axis. By increasing $\beta$, one can artificially create a smaller error field. Furthermore, when $\beta$ is set close to zero, there are regions in the domain where the function value approaches zero, making it difficult to evaluate errors in these regions. To address these issues, we compute the following unbiased error measures ($err_1$ and $err_2$) as well,

\begin{equation*}
\% \,err_1(x,t) = \frac{|u(x,t) - u^e(x,t)|}{\mbox{rms}(u^e,t)}\times 100;
\qquad \% \,err_2(t) = \frac{\mbox{rms}(u-u^e,t)}{\mbox{rms}(u^e,t)}\times 100,
\end{equation*}

where we define rms$(u,t)$ to be the spatial root mean square value of any function given by

\begin{equation*}
\mbox{rms}(u,t) = \frac{1}{L}\sqrt{\int_0^L\bigl(u(x,t)\bigl)^2\,dx}
\end{equation*}

and $u^e(x,t)$ represents the exact solution which can be obtained using the expression presented in Appendix <ref>. $err_1(x,t)$ represents a local measure of error while $err_2(t)$ represents a global measure of error. For $\beta=10$, we also find that the maximum value of $err_2(t)$  decreases from $0.93\%$ to $0.82\%$.
when the mesh is refined from $200\times25$ to $400\times50$. Furthermore, we compute the results for a value of $\beta=0$ for which only $err_1(x,t)$ and $err_2(t)$ are evaluated. The results for this problem along with the error profiles are shown in Fig. <ref>.

Line plots for $\theta(x,t)$

$\% \,err_1(x,t)$

Errors and line plots for DtP mapping generated primal field $\theta$ with a smoothed jump in the initial condition. Mesh is  $200\times50$, $\beta=10$, $k=0.1$ and $T=0.25$.

§.§ Linear transport equation

In order to solve the dual wave equation (<ref>), the weak form (<ref>) is employed. Using the discrete approximation

\begin{equation*}
\lambda(x,t) = \lambda^A N^A(x,t),
\end{equation*}

where $N$ is the total number of nodes in the FE mesh,
we obtain the following set of equations:

\begin{equation*}\label{eq:weak_dual_wave_discrete}
K^{AB} \lambda^B = R^A,
\end{equation*}

\begin{multline*}
K^{AB} = - \int_0^T dt \int_0^L dx \, \Bigl( \p_t N^A(x,t)\, (\p_tN^B(x,t) + c \p_xN^B(x,t))  \\ -
c\,\p_x N^A(x,t)\, (\p_tN^B(x,t) + c \p_xN^B(x,t)) \Bigl);
\end{multline*}

\begin{equation*}
\begin{gathered}
R^A =
\int_0^T dt \, N^A(0,t) u_l(t) +
\int_0^L dx \, N^A(x,0) u_0(x).
\end{gathered}
\end{equation*}

Additionally the Dirichlet boundary conditions are given by
$$\lambda(L,t) = \lambda_r(t); \qquad \lambda(x,T) = \lambda_T(x).$$

To approximate the integrals in each of the above expressions, we use a two-point Gauss quadrature rule in each direction. Upon obtaining the discretized field $\lambda$, we evaluate the primal field, $u$, at the Gauss points of each element and compute its continuous projection as described at the beginning of Sec. <ref>.

§.§.§ A Discontinuous Wave

Consider the wave equation (<ref>) with the following boundary conditions  on a domain $\Omega$ with $L=2$:

\begin{equation*}
u_0(x) = \begin{cases}
2 &  \mbox{for } x< 0.2 \\
4 &  \mbox{for } x>0.2
\end{cases};  \qquad u_l(t) = 2; \qquad c = 0.25.
\end{equation*}

The exact solution to this problem is given by

\begin{equation} \label{eq:wave_p1_exact}
u(x,t) = \begin{cases}
2 &  \mbox{for } x< 0.2 +t \\
4 &  \mbox{for } x>0.2+t.
\end{cases}
\end{equation}

We utilize the weak form (<ref>) along with the Dirichlet boundary conditions $\lambda_T(x) = \lambda_r(t) = 0$ to compute the dual field in $\Omega$.

$\%$ error in $u(x,t)$

DtP mapping generated primal field $u$ for a discontinuous i.c. in wave equation. Mesh is $200\times55$ per stage, $T_e^{(s)}=0.5$, $T=5$ and $c=0.25$.

In order to simulate this problem for a longer period of time, we utilize the `time-slicing' strategy as alluded to in the preamble of Sec. <ref>. We set a time per stage, $T^{(s)}$, and discard the results beyond  $T_e^{(s)}, 0 < T_e^{(s)} < T^{(s)}$. For any stage, $s$, we refer to the dual solution at $(x,t)$ as $\lambda^{(s)}(x,t), t \in [0,T^{(s)}]$, and the corresponding solution for $u$, generated through the DtP mapping, as $u^{(s)}(x,t)$. The global-in-time $u$ solution corresponding to $t \in [0,T^{(s)}]$ is then given by $u\left(x, t + \sum^{s-1}_{i = 1} T^{(i)}_e\right) = u^{(s)}(x,t)$ (with the sum taking the value $0$ for $s = 1$). For any stage, $s$, we set the initial condition, $u^{(s)}(x,0) = u^{(s-1)}\left(x,T_e^{(s-1)}\right) $, with $u^{(1)}(x,0) = u_0(x)$. We repeat this procedure to obtain the results until the simulation reaches a required time, $T$, such that

$$ T = \sum_s T_e^{(s)}.$$

We assess our scheme on this problem by setting $T^{(s)}=0.55$,  $T_e^{(s)}=0.5$ and $T=5$.
Fig. <ref> and <ref> illustrate the field $u(x,t), t \in [0,T]$ and its associated error, respectively.
These figures have been produced by `stitching' together the results for consecutive stages in time, as defined in the previous paragraph.

Line plot for $u(x,t)$

Overshoot values

Undershoot values

(a) illustrates DtP mapping generated primal field $u$ for a discontinuous i.c. in wave equation (Line plots). The error blips near $x=2$ are a consequence of the Dirichlet boundary conditions on $\lambda$. (b) and (c) illustrate the overshooting, $h_t$, and the undershooting, $h_b$, values at the top and the bottom of the jump, respectively.

$\lambda^{(s)}(x,t), t \in [0,T^{(s)}_e], s \in \{1,\ldots, N^{*}\}$, where $N^{*} = 10$ is the total number of computed stages for this problem.

Line plots for $\lambda^{(1)}(x,t)$

Dual field $\lambda$ for a discontinuous i.c. in the wave equation. Mesh is $200\times55$ per stage, $T_e^{(s)}=0.5$, $T=5$ and $c=0.25$. For every stage, $\lambda_T^{(s)}(x)=\lambda_r^{(s)}(t)=0$. (b) Illustrates $\lambda^{(1)}$ at different times.

Fig. <ref> shows localized errors near the right boundary throughout time. As explained in the preamble to this section, this is a mesh resolution related error arising from the imposition of the b.c. $\lambda_T(x)=\lambda_r(t)=0$ at each stage, resulting in high gradients of the dual $\lambda$ field, as shown in Fig. <ref>.  The primal field $u$ is a function of the gradients in the $\lambda$ field, and errors in the latter transfer to the primal solution. In this specific problem an exact solution for the dual problem can be generated by the method of characteristics and b.cs for the dual field at the top and right boundaries can be determined to minimize such gradients, given the primal initial and inflow boundary conditions. However, this is not a procedure that generalizes to more involved, e.g. nonlinear, problems and hence we have not adopted such an approach. We have also checked that such errors decrease with mesh refinement.

The exact solution (<ref>) involves a spatial discontinuity persisting in time. While the use of strictly piecewise smooth, globally $C^0$ finite element interpolation for the dual problem allows gradient discontinuities, such discontinuities are not necessarily aligned with the characteristic lines of the primal problem, which results in approximation errors.
These errors manifest as overshooting and undershooting of the value of $u$ near the jump, as evident in Fig. <ref>.
Fig. <ref> and Fig. <ref> show that these approximations remain stable as time progresses and, as shown in Fig. <ref>, these errors are localized only in the vicinity of the discontinuity.

In any case, this approximated discontinuity at different times strengthens the claim that our dual scheme can deal with discontinuous functions given that appropriate computational methods are employed.

As a final remark on this problem, we note that while $u^{(s)}(x,0) = u^{(s-1)}\left(x,T_e^{(s-1)}\right)$ by design, $\lambda^{(s)}(x,0) \neq \lambda^{(s-1)}\left(x,T_e^{(s-1)}\right)$, in general. Nevertheless, in any stage $s$, the $u^{(s)}$ and $\lambda^{(s)}$ functions are consistently related by the DtP mapping, exemplifying the many-to-one nature of the DtP change of variables. The latter is a general feature of the duality scheme not particular to this problem (cf. the distinct dual solutions to the steady heat equation in Sec. <ref>) which we exploit in our algorithm.

§.§ Euler's Rigid Body system, with and without damping

To solve the system of equations (<ref>), we employ the Newton Raphson scheme. Correspondingly, we generate a residual $R[\lambda,\delta \lambda] $ given by

\begin{equation}\label{eq:euler_resi}
R\,[\lambda;\delta \lambda] =
\int_0^T dt \, \sum_i (-I_i \, \hat{\omega}_i \, \dot{ \delta\lambda_i} + c_i \, \hat{\omega}_{i+1} \, \hat{\omega}_{i+2} \, \delta\lambda_i + \nu I_i \, \hat{\omega}_i \, \delta\lambda_i) - \sum_i I_i\,\omega_i^0\,\delta\lambda_i(0),
\end{equation}

where $\delta \lambda_i$ represents the test functions corresponding to each of the equations in (<ref>) and for each $i$ and $j$, $\hat{\omega}_i(\lambda_j,\dot{\lambda}_j) = \omega_i(t)$. Consequently,
the Jacobian matrix is generated by considering the variation of the residual (<ref>) in the direction $d\lambda$ and is given by

\begin{multline*}
J\,[d\lambda;\delta \lambda] =
\int_0^T dt \, \sum_i  \left(\bigl(-I_i (\dot{ \delta\lambda_i}\bigl) + \nu \, I_i (\delta \lambda_i)\bigl ) \left(\parderiv{\hat{\omega}_i}{\lambda_j}d\lambda_j+\parderiv{\hat{\omega}_i}{\dot{\lambda_j}}d\dot{\lambda_j}\right)\right)  \\+ \int_0^T dt \, \sum_i  c_i (\delta \lambda_i)\, \hat{\omega}_{i+2}  \left(\parderiv{\hat{\omega}_{i+1}}{\lambda_j}d\lambda_j+\parderiv{\hat{\omega}_{i+1}}{\dot{\lambda_j}}d\dot{\lambda_j}\right)
\\
+\int_0^T dt \, \sum_i c_i (\delta \lambda_i) \, \hat{\omega}_{i+1}  \left(\parderiv{\hat{\omega}_{i+2}}{\lambda_j}d\lambda_j+\parderiv{\hat{\omega}_{i+2}}{\dot{\lambda_j}}d\dot{\lambda_j}\right),
\end{multline*}

where $d\lambda_i$ represent the variations in each of the $\lambda_i$ respectively.

Next, we discretize the domain and approximate the dual fields as

$$\lambda_i(t) = \lambda_i^A N^A(t).$$
The subscript index values are interpreted as modulo 3 below, and the index $A$ runs from $1$ to $N$, where $N$ represents the total number of nodes generated upon discretization of time. Using these approximated fields, the value of $\hat{\omega_i}(t)$ can be evaluated using (<ref>), which is being reproduced below:

\begin{equation*}
\omega_i(t) - \tilde{\omega}_i(t) =  \sum_j \mathbb{K}^{-1}_{ij} \left(I_j \dot{\lambda}_j(t) - \nu I_j \lambda_j(t) \right); \qquad \hat{\omega}_i(\lambda_j,\dot{\lambda}_j) = \omega_i(t)
\end{equation*}

\begin{equation*}
\mathbb{K} =
\begin{bmatrix}
a & c_3 \lambda_3 & c_2 \lambda_2 \\
c_3 \lambda_3 & a & c_1 \lambda_1 \\
c_2 \lambda_2 & c_1 \lambda_1 & a
\end{bmatrix}.
\end{equation*}

We utilize a base state $\tilde{\omega_i}$ (see Sec. <ref>) in this nonlinear problem. The specific choice of $\tilde{\omega}_i$ is subsequently explained in the discussion of Algorithm <ref>.

We utilize the residual (<ref>), along with the approximated dual fields to generate a discrete residual $R_i^A$
at each node $A$ corresponding to the degree of freedom pair $(A,i)$ given by

\begin{equation*}\label{eq:euler_resi_discrete}
R_i^A = \int_0^T dt \, (-I_i \, \hat{\omega}_i \, \dot{N^A}(t) + c_i \, \hat{\omega}_{i+1} \, \hat{\omega}_{i+2} \, N^A(t) + \nu I_i \, \hat{\omega}_i \, N^A(t)) - I_i\,\omega_i^0\,\delta_{i0}),
\end{equation*}

no sum on $i$.

Initialization: Set $s=1$, $t_f^{(0)} = 0$ and $\bfomega_f^{(0)}=\boldsymbol{\omega_0}$. Choose a value for $N_e^{(s)}$ and $tol$.

$\boldsymbol{s_{th}}$ stage:

*  Set $\bflambda^{(s)(-1)}= \bf0$, $t_i^{(s)} = t_f^{(s-1)}$ and  $\bfomega_0^{(s)} = \bfomega_f^{(s-1)}$. Over the domain $\bfOmega^{(s)}$, set $\boldsymbol{\tilde{\omega}}^{(s)} = \bfomega_0^{(s)}$

* Newton's Method is used to evaluate the current-stage results. Set $i=0$

For $i \geq 0$:

*  Set $\bflambda^{(s)(i)}=\bflambda^{(s)(i-1)}$.

* Evaluate $\bfR^{(s)(i)}$ and $\bfJ^{(s)(i)}$.

* Solve for $\boldsymbol{d\lambda}^{(i)(s)}$ using (<ref>) and set $\bflambda^{(s)(i+1)}=\bflambda^{(s)(i)} + \boldsymbol{d\lambda}^{(s)(i)}$.

* Evaluate $d^{(s)(i)} = \max_{(A,j)}\,|(\lambda^A_j)^{(s)(i+1)}-(\lambda^A_j)^{(s)(i)}|$.

if $d^{(s)(i)}<tol$ then go to step <ref>

else do $i=i+1$ and go to step <ref>

*  Evaluate the values of $\bfomega$ at the nodes using the projection (<ref>). Discard the last $N_c$ elements along with the results obtained on them. The retained nodes and their values establish the solution for the  current stage  $\bfomega^{(s)}$.

* Set $\bfomega_f^{(s)}$ and $t_f^{(s)}$ based on the retained values of $\bfomega$.

* Set $s=s+1$ and repeat steps 1-4 until $t_f^{(s)} \geq T$

Algorithm to solve Euler's Rigid Body system

Additionally, $R_i^N=0$ is a consequence of the Dirichlet boundary condition applied as in (<ref>). The discrete version of the Jacobian corresponding to the degree of freedom pair $(A,i),(B,j)$ is given as follows:

\begin{multline*}
J_{ij}^{AB} = \int_0^T dt \,  \bigl(-I_i \, \dot{ N^A}(t) + \nu \,I_i \,N^A(t) \bigl) \left(\parderiv{\hat{\omega}_i}{\lambda_j}N^B(t)+\parderiv{\hat{\omega}_i}{\dot{\lambda_j}}\dot{N^B}(t)\right)  \\
+ \int_0^T dt \,  c_i \,N^A(t) \,\hat{\omega}_{i+2} \left(\parderiv{\hat{\omega}_{i+1}}{\lambda_j}N^B(t)+\parderiv{\hat{\omega}_{i+1}}{\dot{\lambda_j}}\dot{N^B}(t)\right) \\
+ \int_0^T dt \, c_i \, N^A(t) \,\hat{\omega}_{i+1} \left(\parderiv{\hat{\omega}_{i+2}}{\lambda_j}N^B(t)+\parderiv{\hat{\omega}_{i+2}}{\dot{\lambda_j}}\dot{N^B}(t)\right),
\end{multline*}

no sum on $i$, and the derivatives of $\hat{\omega_i}$ can be evaluated by differentiating the DtP mapping (<ref>) with respect to $\lambda$ and $\dot{\lambda}$, respectively. The final expressions for these derivatives are shown below:

\begin{gather*}
\parderiv{\hat{\omega_i}}{\lambda_k}=  \mathbb{K}^{-1}_{ij} [\bff_k]_j;
\qquad \parderiv{\hat{\omega_i}}{\dot{\lambda_k}} =  \mathbb{K}^{-1}_{ij} [\bfg_k]_j,
\end{gather*}

\begin{gather*}
\bff_1 = \begin{bmatrix}
-\nu I_1 \\
-c_1 \omega_3 \\
-c_1 \omega_2
\end{bmatrix}; \qquad
\bff_2 = \begin{bmatrix}
-c_2 \omega_3 \\
-\nu I_2 \\
-c_2 \omega_1
\end{bmatrix} ;\qquad
\bff_3 = \begin{bmatrix}
-c_3 \omega_2 \\
-c_3 \omega_1 \\
-\nu I_3
\end{bmatrix}; \\
\bfg_1 = \begin{bmatrix}
I_1 \\
0 \\
\end{bmatrix}; \qquad
\bfg_2 = \begin{bmatrix}
0 \\
I_2 \\
\end{bmatrix} ;\qquad
\bfg_3 = \begin{bmatrix}
0 \\
0 \\
\end{bmatrix}. \\
% \mathbb{K} =
% \begin{bmatrix}
% a & c_3 \lambda_3 & c_2 \lambda_2 \\
% c_3 \lambda_3 & a & c_1 \lambda_1 \\
% c_2 \lambda_2 & c_1 \lambda_1 & a
% \end{bmatrix}
\end{gather*}

Let $(\cdot)^{(n)}$ denote the value of $(\cdot)$ at the $n^{th}$ iterate of a NR solve. Using this notation, at each NR iterate $n$, the corrections can be evaluated as:

\begin{equation*}
J^{AB\,(n)}_{ij} d\lambda^{B\,(n)}_j = -
\end{equation*}

The corrected $\lambda_i$ for the next iteration is obtained via:

\begin{equation}\label{eq:increment_vector}
\bflambda^{(n+1)}=\bflambda^{(n)} + \boldsymbol{d\lambda}^{(n)},
\end{equation}

where $\boldsymbol{\lambda}$ and $\boldsymbol{d\lambda}$ represent collective terms for $\lambda_i^A$ and $d\lambda_i^A$, respectively.

To obtain a solution to the set of equations (<ref>) at a final time $T$, we proceed incrementally by solving the equations on smaller subdomains. Each of these subdomains is referred to as a stage. We select a time $T^{(s)} \leq T$ as the length of each stage. After each stage is completed, we accumulate its time and continue running the simulation until the cumulative time reaches the final time $T$.

We denote the domain for the $s^{th}$ stage as $\Omega^{(s)}$. We further discretize each of the stages and refer to the number of corresponding elements in stage $s$ as $N_e^{(s)}$. Within each stage, the Newton Raphson (NR) algorithm is employed to solve the set of equations (<ref>). The algorithm employed to achieve this, along with the motivation for a choice of $\tilde\omega$, has been explained in the following paragraph.

An initial guess for the dual field $\bflambda$ is required by the NR algorithm within each stage. It is practical, and more feasible, to make an initial guess for a primal field $\bfomega$ with physical meaning, and infer a dual field $\bflambda$ consistent with it as a solution of the DtP mapping. By (<ref>), a primal initial guess of $\bfomega^{(s)}(t) := \tilde{\bfomega}^{(s)}(t)$ is consistent with the dual initial guess of $\bflambda^{(s)}(t) = \bfzero$, for $t$ belonging to the $s^{th}$ stage. We choose this initial guess for $\bfomega$ within each stage $s$ as $\bfomega^{(s)}(t)=\tilde{\bfomega}^{(s)}(t) = \bfomega_0^{(s)}$, where $\bfomega_0^{(s)}$ denotes the initial condition of $\bfomega$ in that stage.
For the first stage, we set the specified initial condition for the primal problem, $\bfomega_0$, as this constant and for all subsequent stages, we set the constant equal to the value of $\bfomega$ obtained at the end of the previous stage.
The algorithm capturing these ideas has been outlined in Table <ref>.

Thus, our overall algorithm in effect designs different dual functionals for each stage, the set of functionals parametrized by the piecewise constant choice of a base state defining the potential, say $H^{(s)}$, for that stage, as well as the initial condition for the stage.

The following notation is used for the description of the algorithm:
$s$ denotes the active stage. The collective vectors, $\bfomega$ and $\boldsymbol{\tilde{\omega}}$, consist of $\omega_i^A$ and the base state values $\tilde{\omega_i}^A$, respectively, where the index $i$ takes values from $1$ to $3$ and  the index $A$ ranges over all the nodes in the current stage. $\bfomega_0^{(s)}$ and $\bfomega_f^{(s)}$ denote the values of $\omega_i$ at the first and the last node under the current stage, respectively. For Newton's method, let $(\cdot)^{(s)(i)}$ denote the value of $(\cdot)$ at a stage $s$ and an iteration number $i$. It is important to note that, typically, each stage generates an error towards the end of the stage. To address this problem, we discard the results obtained for a specific number of elements at the end of each stage, and then start the next stage from the final node of the last retained element. We denote the number of discarded elements as $N_c$. Also, $t_i^{(s)}$ and $t_f^{(s)}$ denote the time at the start and end of any stage.

§.§.§ Example 1: Free Rotation

To demonstrate the algorithm, consider the set of equations (<ref>) with the following conditions:

\begin{equation*}
I = [1,2,3];  \qquad
\omega_0 =
\qquad \nu = 0.
\end{equation*}

We run the problem with $T^{(s)}=0.5$ for every stage until $T=3$ is reached. The number of elements considered per stage is 20. The results obtained for $\omega_i$ are plotted in Fig. <ref>.

$E(t)$ and $L(t)$

DtP mapping generated results for a free rotation. Mesh is 20 elements per stage in (a) and (b) and $N_c=5$ for all stages. (c) illustrates error for $\omega_i$ obtained on additional meshes. The maximum error drops quadratically when compared against the mesh size.

Since this example represents a free rotation, the magnitude of angular momentum and kinetic energy are conserved. The obtained results approximate these features well, as shown in Fig. <ref>.
We evaluate the error obtained in $\bfomega(t)$  by comparing it against a reference solution $\bfomega^\bfe(t)$.
The percentage error for $\bfomega$ (represented by $err(\omega)$) is given by:

err(\omega) = 100 \times\sqrt{\frac{\mathstrut (\omega_1 - \omega^e_1)\strut^2+(\omega_2 - \omega^e_2)\strut^2+(\omega_3 - \omega^e_3)\strut^2}{ (\omega^e_1)\strut^2+(\omega^e_2)\strut^2+(\omega^e_3)\strut^2}}.

For the current example, we set the reference solution to be the exact solution obtained in Appendix <ref>.
Similar to the cases of the heat and wave equation examples, error concentrations appear at the end of every stage.
To overcome this shortcoming, we discard the results obtained for a few elements (denoted by $N_c$) near the end of the stage as explained at the start of the section <ref>.

Additionally, small errors that are relatively larger than the domain errors are found near the beginning of each stage. Predictably, the refinement of the mesh per stage reduces these errors. The error obtained for the current case of 20 elements along with a few different types of meshes ($N_c=5$ for all the stages) are shown in Fig.<ref>.

This example demonstrates that the dual methodology has the potential to solve nonlinear problems, preserving conservation properties when expected.

§.§.§ Example 2: Dissipative System

The following example demonstrates the algorithm with an active viscous damping term. Consider the following:

\begin{equation*}
I =
[1,2,5];  \qquad
\omega_0 = [5,3,0];
\qquad \nu = 0.4.
\end{equation*}

These parameter values have been chosen solely for the purpose of demonstrating an example. In fact, the algorithm remains applicable to any value other than the selected values. To evaluate the dual scheme on the current problem, we set the solution obtained from MATLAB's ode45 tool with the default settings <cit.> as the benchmark solution. The plots for $\omega_i$ along with the error obtained for $\bfomega$ against the benchmark are shown in Fig. <ref> and <ref>, respectively.
Additionally, plots for the magnitude of angular momentum $L$ along with its respective error $err(L)$ when compared against the exact expression obtained in Appendix <ref> are produced.
These have been illustrated in Fig. <ref> and Fig. <ref>, respectively.

$\%$ Error in $L(t)$

DtP mapping generated results for a dissipative rigid body rotation. The mesh is 20 elements per stage. $N_c=5$, $T^{(s)}=0.3$ per stage and $T=5$.

Surprisingly, these results are based on the Euler-Lagrange equations of the variational principle (<ref>), even though it is believed that, in general, this is not possible for dissipative (P)ODE (in this regard, it is perhaps worth noting that a postulated gradient flow is not an E-L equation of a functional defined over (space)-time fields).

§ CONCLUDING REMARKS

We have provided proof-of-principle demonstrations of a duality-based scheme for generating solutions to differential equations. Non-standard dual formulations for linear parabolic and hyperbolic PDE, and a nonlinear ODE system have been approximated by the finite element method and shown to reproduce correct primal responses. Interestingly, our approach provides a somewhat unifying view of converting initial value problems to well-set degenerate elliptic boundary value problems in (space)-time. This works because the order of the equation in time derivatives (as well as space derivatives when present) increases, allowing for the imposition of an extra, final-time boundary condition on the dual variables which, however, does not get in the way of correct, `free' evolution of the primal field at the considered final time (i.e., constrained only by the primal equation and its initial and boundary conditions), as conjectured in <cit.>. A curious observation is that in computationally solving problems with wave-propagation, no explicit account of `upwinding' or domains of dependence of the solution has been involved in the example considered, but much more study to precisely understand the details of such features is required.

Our approach falls within the umbrella of ideas of `Hidden Convexity' in PDE advanced by Yann Brenier <cit.>, with some key differences of interpretation. In particular, our understanding and allowance for

* the use of a large class of potentials $H$ in formulating the dual problem, which includes the judicious use of `base states' in designing the dual variational principle for nonlinear problems;

* final-time dual boundary conditions as a necessary, routine matter for solving strict primal initial value problems without invoking the degenerate ellipticity of the dual E-L equations for associating a well-set dual variational principle with a primal initial value problem; and

* recovery of the primal equations as a critical point or Euler-Lagrange equations of an appropriate dual functional as the principal guiding requirement of the scheme

serve as the main conceptual points of departure from Brenier's approach.

§ ACKNOWLEDGMENT

We thank Vladimir Sverak for helpful discussions. This work was supported by the grant NSF OIA-DMR #2021019 and was also supported by a Simons Pivot Fellowship to AA. It was completed while he was on sabbatical leave at a) the Max Planck Institute for Mathematics in the Sciences in Leipzig, and b) the Hausdorff Institute for Mathematics at the University of Bonn funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC-2047/1 – 390685813, as part of the Trimester Program on Mathematics for Complex Materials. The support and hospitality of both institutions is acknowledged.

§ APPENDIX: EXACT SOLUTION OF THE HEAT EQUATION WITH AN ALMOST DISCONTINUOUS INITIAL CONDITION

For given parameters $\veps$ and $\beta$, the solution of the heat equation (<ref>), along with the boundary conditions $\theta_l(t)=\theta_r(t)=\beta$ and the initial condition

\begin{equation*}
\theta_0(x) = \begin{cases}
\beta+2x &  \mbox{for } 0 \leq x< 0.5-\veps  \\
kx+c    &  \mbox{for } 0.5 - \veps \leq x \leq 0.5+ \veps \\
\beta-2+2x &  \mbox{for }  0.5+\veps < x \leq 1,
\end{cases};    \qquad k = \frac{2\veps-1}{\veps}; \qquad c = \beta-\frac{2\veps-1}{2\veps},
\end{equation*}

is given by:

\begin{equation*}
\theta(x,t) = \beta + \sum_{m=1}^\infty a_m \sin(2\pi mx)e^{-(2m)^2\pi^2 kt},
\end{equation*}

$$a_m = a_{m_1} + a_{m_2} + a_{m_3};$$

\begin{equation*}
a_{m_1} = \frac{1}{2\pi^2m^2}\sin(2\pi ml) - \frac{l}{\pi m}\cos(2 \pi ml) + \frac{\beta}{2\pi m}\bigl(1-\cos(2\pi ml)\bigl);
\end{equation*}

\begin{multline*}
a_{m_2} = \frac{1}{4 \pi^2 m^2}\Bigl(2\pi c m\bigl(\cos(2\pi ml) - \cos(2\pi rm)\bigl) + k\bigl(-2\pi rm\cos(2\pi rm) \\ + \sin(2\pi rm)
+ 2\pi lm \cos(2\pi lm) - \sin(2\pi lm)\bigl)\Bigl);
\end{multline*}

\begin{equation*}
a_{m_3} = \frac{1}{2\pi^2m^2}\bigl(-1\beta\pi m + (\beta-2+2r)\pi m \cos(2\pi mr) - \sin(2\pi mr) \bigl);
\end{equation*}

\begin{equation*}
l = 0.5 - \veps; \qquad r =0.5+\veps.
\end{equation*}

§ APPENDIX: DERIVATION OF EULER'S SYSTEM FOR RIGID BODY MOTION WITH VISCOSITY

Consider the motion of a rigid body with one of its points fixed. In the following part, $\dot{()}$ represents the time derivative $\deriv{()}{t}$. Let $\bfI$ represent the moment of Inertia tensor for the body with $\bfa_i$ and $I_i$ representing its principal directions and principal values, respectively. The indices presented in this section are modulo 3. The $\bfa_i$ vectors form an orthonormal triad. Let

\begin{equation*}
\bfI = \sum_i I_i \, \bfa_i \otimes \bfa_i.
\end{equation*}

Let $\bfomega$ denote the angular velocity of the body and

\begin{equation*}
\bfomega =  \omega_i \, \bfa_i.
\end{equation*}

While $\bfa_i$ and $\omega_i$ vary with time, $I_i$ remain constant due to the rotational invariance of eigenvalues. The angular momentum of the body can be represented by $\bfL = \bfI \bfomega$. Denoting $\bfN$ as the external torque acting on the body, the balance of angular momentum leads to

\begin{equation}\label{eq:angular_mom}
\frac{d\bfL}{dt} = \frac{d(\bfI\bfomega)}{dt} = \bfN.
\end{equation}

Let $\bfa_i(t) = \bfR(t) \bfa_i(0)$ for some reference basis vectors being fixed at time zero and $\bfR(t)$ represent the rotation tensor of the triad at time $t$. The time derivative of $\bfL$ can be represented ain terms of a convected derivativ as follows:

\begin{equation}   \label{eq:angular_simplified}
\frac{d\bfL}{dt}  =  \frac{d(L_i\bfa_i)}{dt}
= \dot{L_i}\bfa_i + L_i\dot{\bfa_i}
= I_i \dot{\omega_i} + L_i \dot{\bfR} \bfa_i(0)
= I_i \dot{\omega_i} + L_i \dot{\bfR} \bfR^T \bfa_i(t).
\end{equation}

Let $\bfW$ be the unique skew tensor corresponding to the angular velocity $\bfomega$ satisfying $\bfomega \times \bfv = \bfW\bfv$ $\forall \bfv \in \mathbb{R}^3$ and $\bfW(t) = \dot{\bfR(t)} \bfR(t)^T$. Using (<ref>) and (<ref>),

\begin{equation*}
\frac{d\bfL}{dt} = I_i \dot{\omega_i} + \bfomega \times (L_i \bfa_i) =
I_i \dot{\omega_i} + (\omega_i \bfa_i) \times (L_i \bfa_i)
= \bfN.
\end{equation*}

Writing out in components on the $\bfa_i(t)$ basis, Euler’s equations of motion for rigid-body rotation are obtained as

\begin{equation}\label{eq:euler_original}
\begin{aligned}
I_1 \frac{d\omega_1}{dt} + (I_3 - I_2)\,\omega_2 \,\omega_3 = N_1, \\
I_2 \frac{d\omega_2}{dt} + (I_1 - I_3)\,\omega_3 \,\omega_1 = N_2, \\
I_3 \frac{d\omega_3}{dt} + (I_2 - I_1)\,\omega_1 \,\omega_2 = N_3, \\
\end{aligned}
\end{equation}

where $N_i$ are the components of torque on the $\bfa_i(t)$ basis.

Free rotation of a rigid body about a fixed point is governed by the system (<ref>) for $N_i = 0$. (Rayleigh) damped rotations of the body about a fixed point can be represented by considering viscous damping torques given by $\bfN= - \nu\bfL, \ \nu > 0$ in (<ref>). In the rotating frame defined by the principal axes of inertia, the governing system reduces to

\begin{equation}\label{eq:euler_free_equations}
I_i \dot{\omega_i} + (I_{i+2} - I_{i+1})\,\omega_{i+2} \,\omega_{i+1} + \nu I_i \omega_i = 0,
\end{equation}

no sum on i,
with the system for free-rotation corresponding to $\nu = 0$.

An analytical solution for the free rotation problem can be found in <cit.>, and is reproduced in Appendix <ref> for the sake of this paper being self-contained.
The dissipative system given by $\frac{d\bfL}{dt} = - \nu \bfL$ implies that the magnitude of the angular momentum $L = |\bfL|$ evolves as $L(t) = L(0)e^{-\nu t}$, where $L(0)$ is its initial value.

For $\nu = 0$, angular momentum is conserved and $L(t) = L(0)$.

§ APPENDIX: ANALYTICAL SOLUTIONS TO FREE ROTATIONS OF A RIGID BODY <CIT.>

For no external torque $\bfN$, angular momentum is conserved (Appendix <ref>) and

\[
|\bfL(t)|^2 = L^2(0).
\]

\begin{equation*}
\deriv{\bfL}{t} = 0 \Rightarrow \bfomega\cdot\deriv{(\bfI\bfomega)}{t} = 0\Rightarrow \frac{1}{2}\deriv{}{t}\left(\bfomega\cdot\bfI\bfomega\right) = 0\Rightarrow \frac{1}{2} \bfomega(t)\cdot\bfI(t)\bfomega(t) =: E(t) = E(0),
\end{equation*}

which represents conservation of kinetic energy.
In the following, we refer to the constants $L(0), E(0)$ as simply $L, E$ for notational convenience.

In the rotating frame defined by the principal axes of inertia, the above conditions reduce to:

\begin{equation*}
\begin{aligned}\label{eq:euler_conservation}
L_1^2 + L_2^2 &+ L_3^2 = L^2,\\
I_1\omega_1^2 +I_2\omega_2^2&+I_3\omega_3^2=2E.
\end{aligned}
\end{equation*}

Using (<ref>) with $\nu=0$, $\omega_1$ and $\omega_3$ can be expressed  in terms of $\omega_2$ as

\begin{equation}\label{eq:euler_substitution}
\begin{gathered}
\omega_1^2 = \frac{(2EI_3-L^2)-I_2(I_3-I_2)\omega_2^2}{I_1(I_3-I_1)}; \qquad
\omega_3^2 = \frac{(L^2-2EI_1)-I_2(I_2-I_1)\omega_2^2}{I_3(I_3-I_1)}.
\end{gathered}
\end{equation}

Also, when $i=2$ and $\nu=0$, the equation obtained from (<ref>) is:

\begin{equation*}
\left(\deriv{\omega_2}{t}\right)^2 = \left(\frac{(I_3-I_1)\omega_1\omega_3}{I_2} \right)^2.\nonumber
\end{equation*}

Using the substitution (<ref>) in the last equation results in:

\begin{equation*}
\left(\deriv{\omega_2}{t}\right)^2
= \frac{[(2EI_3-L^2)-I_2(I_3-I_2)\omega_2^2]\times[(L^2-2EI_1)-I_2(I_2-I_1)\omega_2^2]}{I_2^2 I_1I_3}.
\end{equation*}

\begin{equation}\label{eq:euler_new_variables}
\begin{gathered}
\tau = t\sqrt{\frac{(I_3-I_2)(L^2-2EI_1)}{I_1I_2I_3}}; \qquad
s = \omega_2\sqrt{\frac{I_2(I_3-I_2)}{2EI_3-L^2}}; \qquad
k^2 = \frac{(I_2-I_1)(2EI_3-L^2)}{(I_3-I_2)(L^2-2EI_1)}.
\end{gathered}
\end{equation}

where $k^2<1$ forms a positive parameter for the elliptic function. This leads to the equation

\begin{equation}\label{eq:app2_elliptic_sn}
\left(\deriv{s}{\tau}\right)^2 = (1-s^2)(1-k^2s^2).
\end{equation}

The value of $\varphi$ satisfying:

$$ \tau = \int_0^\varphi \frac{d\theta}{\sqrt{1-k^2\sin^2(\theta)}}$$
is known as the Jacobi amplitude  and denoted by am$(u,k^2)$. Correspondingly, the Jacobi Elliptic functions are defined as

\begin{equation*}
\mbox{sn}(\tau,k^2) = \sin(\text{am}(\tau,k^2));
\qquad
\mbox{cn}(\tau,k^2) = \cos(\text{am}(\tau,k^2));
\qquad
\mbox{dn}(\tau,k^2) =
\frac{d}{du}(\text{am}(\tau,k^2)).
\end{equation*}

The Jacobi Elliptic function $\mbox{sn}(\tau,k^2)$ satisfies the equation (<ref>) exactly. Using (<ref>), one can obtain $\omega_2$ explicitly in terms of the elliptic Jacobi functions. Also, $\omega_1$ and $\omega_3$ can algebraically obtained through equations (<ref>).  The final expressions for each of the $\omega_i$ are shown below:

\begin{equation*}
\begin{aligned}
\omega_1 &= \sqrt{\frac{2EI_3-L^2}{I_1(I_3-I_1)}}\mbox{cn}\,(\tau,k^2); \quad
\omega_2 &= \sqrt{\frac{2EI_3-L^2}{I_2(I_3-I_2)}}\mbox{sn}\,(\tau,k^2); \quad
\omega_3 &= \sqrt{\frac{2EI_3-L^2}{I_3(I_3-I_1)}}\mbox{dn}\,(\tau,k^2).
\end{aligned}
\end{equation*}