A central scheme for coupled hyperbolic systems

Michael Herty^1 Niklas Kolbe^1[Corresponding author] Siegfried Müller^1

A novel numerical scheme to solve coupled systems of conservation laws is introduced. The scheme is derived based on a relaxation approach and does not require information on the Lax curves of the coupled systems, which simplifies the computation of suitable coupling data. The coupling condition for the underlying relaxation system plays a crucial role as it determines the behavior of the scheme in the zero relaxation limit. The role of this condition is discussed, a consistency concept with respect to the original problem is introduced, well-posedness is analyzed and explicit, nodal Riemann solvers are provided. Based on a case study considering the $p$-system of gas dynamics a strategy for the design of the relaxation coupling condition within the new scheme is provided.

35L65, 35R02, 65M08

§ INTRODUCTION

Coupled systems of hyperbolic conservation or balance laws have been discussed in the mathematical literature to address various modeling questions e.g., regarding two-phase dynamics, multi-scale processes and networks, see e.g., <cit.>. The mathematical treatment of interface coupling has been motivated by boundary value problems <cit.> or by applications to traffic flow on road networks <cit.>, or gas dynamics in large-scale pipelines <cit.>. It has been further expanded into a well-posedness result for coupled general gas dynamics, see e.g.  <cit.>. Usually the coupling is imposed at a single spatial interface point as it has been done e.g., to model production systems <cit.> and blood flow through systems of blood vessels <cit.> to mention only a few.
A main challenge has been the derivation of well-posed initial boundary value problems for coupled dynamics out of physically (or problem) induced conditions, see <cit.> for a recent review. A well-known technique by now is the analytical concept of (half-)Riemann problems <cit.> or  solvers <cit.>.
Moreover, the use of wave-front-tracking techniques yields well-posed coupled problems for a variety of mathematical fluid-type models, see e.g. <cit.>.
This analytical tool has given rise to Riemann solver based numerical methods, which have been proposed e.g. in <cit.> and since then extended towards high-order methods <cit.>. Some high-order methods rely on the linearization of the imposed conditions such that the Lax-curves are obtained trivially, see e.g. <cit.>. We focus here on finite-volume schemes —  finite-element based approaches, such as <cit.>, require a different treatment of coupling conditions. Linearization techniques and Roe-type schemes have been used to avoid the explicit computation of eigenvalues in  two-phase problems, see e.g.. <cit.>. A linearization of the nonlinear dynamics at the coupling point is also at the core of an implicit (in space) numerical scheme <cit.>.
Recently using projection of the fluxes in the direction normal to the interface an extension of these techniques to the multi-dimensional case has been introduced in <cit.>. A related problem is given by conservation laws with discontinuous flux at a single point in space. This problem has applications, for instance, in porous media <cit.>, sedimentation processes <cit.>, traffic flow <cit.> or in supply chain models <cit.>. In the scalar, one-dimensional case, if conservation of mass and a modified entropy condition across the interface are imposed, solutions are  constructed <cit.> and numerically approximated  <cit.>. A numerical approximation based on relaxation <cit.> of (scalar) discontinuous fluxes has been proposed in <cit.>. These schemes however, rely on analytical expression of Lax-curves which limits those techniques e.g., when treating phenomena like multi-phase flow. In various cases such curves are not available <cit.> or not explicitly given <cit.>.
In this work a new approach is introduced that does not rely on Lax-curves at the interface. Therefore, building on our prior work on scalar equations on networks <cit.>, the nonlinear systems are rewritten in the relaxation form proposed in <cit.>. Due to the linearity of the relaxed systems the Lax-curves are linear spanned by the constant eigenvectors of the flux matrix. Based on the reformulation we derive an explicit scheme for coupled systems of conservation laws by taking a discretization to the zero relaxation limit as suggested also in <cit.>. A related approach has been followed in <cit.>, where a kinetic relaxation model governs the coupling conditions. Another alternative is the construction of vanishing viscosity solutions, addressed in the schemes in <cit.>, that avoids the use of Lax-curves at the expense of formally treating parabolic systems. The system case that is discussed here requires a detailed analysis of the modified coupling condition imposed to the relaxation formulation, which we will consider.
The outline of this manuscript is as follows: in Section <ref> we introduce the relaxation form for coupled systems of hyperbolic systems as well as a consistency concept for the corresponding coupling condition. In Section <ref> we discuss suitable Riemann solvers of the coupled problem, discuss their well-posedness and provide an explicit form for the case of linear coupling conditions. Section <ref> is concerned with the construction of a finite-volume scheme for the coupled problem. In Section <ref>, we consider the $p$-system of gas dynamics in a case study and conduct numerical experiments for different choices of the relaxation coupling condition and then conclude providing a general strategy for the  construction of this condition.

§ THE COUPLED RELAXATION SYSTEM

We consider a coupled system of hyperbolic conservation laws on the real line of the form

\begin{align}
\ddt U + \ddx F_1(U) &=0,\quad \text{in}(0, \infty) \times (-\infty, 0),\\
\ddt U + \ddx F_2(U) &=0,\quad \text{in} (0, \infty) \times (0,\infty),
\end{align}

with vector valued state variable $U(t,x)\in \R^n$ and flux functions $F_1,F_2:\R^n\rightarrow \R^n$. At the interface located at $x=0$ a transition between the two flux functions takes place and we impose the coupling condition

\begin{equation}\label{eq:originalcoupling}
\Psi_U(U(t, 0^-), U(t, 0^+)) =0 \qquad \text{for a. e. }t>0
\end{equation}

assuming a suitable mapping $\Psi_U:\R^n \times \R^n \rightarrow \R^m$ for $m\in\N$. To obtain a well-posed problem a sufficiently large number of suitable conditions is required and thus $m$ depends on the given flux functions. Moreover, we assume initial data given by the vector valued function $U^0(x)$, which is compatible with the coupling condition (<ref>).

The relaxation system <cit.> was introduced as a dissipative approximation of a hyperbolic system, recovering the original system in the so-called relaxation limit. In this work we use it as a tool to derive coupling information giving rise to a new method.
The relaxation system corresponding to (<ref>) employs the auxiliary state variable $V(t,x)\in \R^n$ and the relaxation rate $\varepsilon>0$ and reads

\begin{align}
\ddt U + \ddx V &= 0, && \text{in }(0, \infty) \times \R\setminus \{0 \}, \label{eq:relsyst11u}\\
\ddt V + A_1 \, \ddx U &= \frac 1 \eps (F_1(U) - V),  && \text{in }(0, \infty) \times (-\infty, 0), \label{eq:relsyst11vl}\\
\ddt V + A_2 \, \ddx U &= \frac 1 \eps (F_2(U) - V),  && \text{in }(0, \infty)\times(0, \infty), \label{eq:relsyst11vr}
\end{align}

where $A_1$ and $A_2$ denote diagonal $n\times n$ matrices with positive diagonal entries $a_j^i$ for $j=1,\dots N$ and $i=1,2$. These matrices satisfy the subcharacteristic condition

\begin{equation}\label{eq:subchar}
A_i - (DF_i(U))^2 \geq 0 \qquad i=1,2
\end{equation}

for all $U$. Here, $DF_i(U)$ refers to the Jacobian of $F_i$ in $U$ and the inequality indicates positive semi-definiteness. We note that both $U$ and $V$ depend on $\varepsilon$. Introducing the vector $Q=(U, V)$ we close system (<ref>) using the relaxation coupling condition

\begin{equation}\label{eq:relaxationcoupling}
\Psiq(Q(t, 0^-), Q(t, 0^+)) =0 \qquad \text{for a. e. }t>0.
\end{equation}

for a mapping $\Psiq:\R^{2n} \times \R^{2n} \rightarrow \R^{l}$, where $l\in\N$. The particular choice of $\Psi_Q$ and $m$ for given $\Psi_U$ is discussed in the following parts of this paper. Initial data for the relaxation system is derived from the original initial data by setting $Q^0=(U^0, V^0)$ such that

\begin{equation}
V^0(x)=F_1(U^0(x))\quad \text{if }x<0, \quad V^0(x)=F_2(U^0(x))\quad \text{if }x>0
\end{equation}

and compatibility of the initial data $Q^0$ with the coupling condition (<ref>) is assumed.

It is well-known that the uncoupled relaxation system recovers the original system of conservation laws as well as $V=F(U)$ as $\varepsilon$ tends to zero. We are interested in this  relaxation limit at the interface as $\varepsilon \rightarrow 0$.

§.§ Continuous relaxation limit

The fact that the uncoupled relaxation system tends to a system of conservation laws in the relaxation limit, see <cit.>, raises the question how the coupled systems and, in particular, the coupling conditions (<ref>) and (<ref>) relate. Here we address this question using asymptotic analysis and suppose that the solution $Q=(U, V)$ of the coupled relaxation system (<ref>) endowed with the coupling condition (<ref>) can be written in terms of a Chapman–Enskog expansion <cit.>, i.e.,

\begin{equation}
Q = \tilde Q^0 + \sum_{k=1}^\infty \varepsilon^k \tilde Q^k.
\end{equation}

As a solution of the relaxation system in the limit $\varepsilon \rightarrow 0$ on both half-axes, the state $\tilde Q^0=(\tilde U^0, \tilde V^0)$ has the property

\begin{equation}
\tilde V^{0}(t, x) = F_1(\tilde U^{0}(t, x)) \text{ if }x<0, \qquad
\tilde V^{0}(t,x) = F_2(\tilde U^{0}(t, x)) \text{ if }x>0
\end{equation}

and $\tilde U^0$ is a solution to the system of conservation laws on both half-axes. If we assume that $\Psi_Q$ is continuously differentiable, we obtain due to Taylor expansion

\begin{equation}
\Psi_Q[Q(t, 0^-), Q(t, 0^+)] =   \Psi_Q\left[\tilde Q^0(t, 0^-), \tilde Q^0(t, 0^+)\right] + \mathcal{O}(\varepsilon)
\end{equation}

and, in particular, in the limit $\varepsilon \rightarrow 0$ we have

\begin{equation}
\Psi_Q\left[\tilde Q^0(t, 0^-), \tilde Q^0(t, 0^+)\right] = 0.
\end{equation}

This observation motivates the following consistency condition.

The coupled relaxation system (<ref>) is consistent with the coupled system of conservation laws (<ref>) iff for a. e. $t>0$ the corresponding coupling conditions (<ref>) and (<ref>) satisfy

\begin{equation}
\Psi_U\left[U(t, 0^-), U(t, 0^+)\right] = 0 \quad \text{iff} \quad \Psiq \left[
\begin{pmatrix} U(t, 0^-) \\
f_1\left(U(t, 0^-) \right)
\end{pmatrix},
\begin{pmatrix} U(t, 0^+) \\
f_2\left(U(t, 0^+) \right)
\end{pmatrix}
\right] = 0.
\end{equation}

Assuming a smaller number of conditions for the system of conservation laws than for the relaxation system, i.e., $m<l$, this definition implies that consistency requires redundancies occurring in the relaxation limit due to (<ref>).

§ RIEMANN SOLVERS FOR THE COUPLED RELAXATION SYSTEM

In this section we consider Riemann solvers for the coupled relaxation system, which identify coupling data that both solve a half-Riemann problem and satisfy the coupling condition. They are an integral component of the finite volume schemes introduced in Section <ref>. We introduce the concept in Section <ref>, discuss their well-posedness with respect to the coupling condition in Section <ref>, derive an explicit form in the linear case in Section <ref> and study their behavior in the relaxation limit in Section <ref>.

§.§ The half-Riemann problem

A finite volume scheme for a coupled problem such as (<ref>) requires suitable coupling/boundary data, which is obtained by solving a half-Riemann problem <cit.>. Later we will introduce the concept of Riemann solvers, which identify coupling data that both solve the half-Riemann problem and satisfy a coupling condition.

To discuss the half-Riemann problem in system (<ref>) we require information about its Lax-curves. The relaxation system can be rewritten on both half-axes ($i=1$ if left and $i=2$ if right) as

\begin{equation}
\ddt Q + S_i \, \ddx Q = \frac{1}{\varepsilon}
\begin{pmatrix}
0 \\
F_i(U) - V
\end{pmatrix}
\end{equation}

where $S_i\in \R^{2n \times 2n}$ refers to the diagonalizable block matrix

\begin{equation}\label{eq:diagonalization}
S_i =
\begin{pmatrix}
0 & I \\
A_i & O
\end{pmatrix}
\begin{pmatrix}
-(\sqrt{A_i})^{-1} & (\sqrt{A_i})^{-1} \\
I & I
\end{pmatrix}
\begin{pmatrix}
-\sqrt{A_i} & 0 \\
0 & \sqrt{A_i}
\end{pmatrix}
\begin{pmatrix}
-\frac 1 2 \sqrt{A_i} & \frac 1 2 I \\
\frac 1 2 \sqrt{A_i} & \frac 1 2 I
\end{pmatrix}
=: R_i \Lambda_i R_i^{-1}
\end{equation}

By taking the root of each diagonal entry in $A$ we obtain $\sqrt{A_i}$, which is invertible and satisfies $\sqrt{A_i} \sqrt{A_i} = A_i$. The decomposition (<ref>) shows that the eigenvalues of the system are given by

\begin{equation}\label{eq:eigenvalues}
\lambda_j = -\sqrt{a^i_j}, \qquad \lambda_{n+j} = \sqrt{a^i_j} \qquad \text{for }j=1,\dots,n.
\end{equation}

The columns of the matrix $R_i$ include the (right) eigenvectors of the system that can be segmented as

\begin{equation}
R_i =
\begin{pmatrix}
-(\sqrt{A_i})^{-1} & (\sqrt{A_i})^{-1} \\
I & I
\end{pmatrix}
\begin{pmatrix}
R_i^- & R_i^+
\end{pmatrix}
\end{equation}

so that $R_i^-\in\R^{2n \times n}$ includes the eigenvectors corresponding to the negative eigenvalues and $R_i^+\in\R^{2n \times n}$ the ones corresponding to the positive eigenvalues. Analogously we define the left eigenvectors as the rows of the matrix

\begin{equation}
L_i \coloneq R_i^{-1}=
\begin{pmatrix}
-\frac 1 2 \sqrt{A_i} & \frac 1 2 I \\
\frac 1 2 \sqrt{A_i} & \frac 1 2 I
\end{pmatrix}
\begin{pmatrix}
L_i^- \\ L_i^+
\end{pmatrix},
\end{equation}

which can again be split with respect to the signs of the corresponding eigenvalues into the matrices $L_i^-$, $L_i^+\in\R^{n \times 2n}$.
The characteristic variables of the system are

\begin{equation}\label{eq:characteristics}
\begin{pmatrix}
W^-_i \\
\end{pmatrix}
:= R_i^{-1}Q = \frac 1 2
\begin{pmatrix}
V - \sqrt{A_i}U \\
V + \sqrt{A_i}U
\end{pmatrix}.
\end{equation}

Starting from a state $Q_0^*=(U_0^*, V_0^*)^T$ the Lax-curves of the systems are lines in the state space given by

\begin{equation}\label{eq:individuallaxcurves}
\mathcal{L}_{\lambda_j}(Q_0^*; \sigma_j^-)= Q_0^* - \sigma_j^-  (a^1_j)^{-1/2} \, e_j + \sigma_j^- e_{n+j}, \quad \mathcal{L}_{\lambda_{n+j}}(Q_0^*; \sigma_j^+)= Q_0^* + \sigma_j^+ (a^2_j)^{-1/2} e_{j} + \sigma_j^+ e_{n+j}
\end{equation}

for $j=1,\dots,n$ with $e_j \in \R^{2n}$ denoting the $j$-th unit vector and parameter $\sigma_j^\pm\in\R$. Moreover, we define the set

\begin{equation}
\mathcal{L}^-(Q_0^*) \coloneq \left\{ Q_0^* + R_1^- \Sigma^-:
\Sigma^- \in \R^{n}
\right\}
\left\{
\begin{pmatrix}
U_0^* -  (\sqrt{A_1})^{-1} \Sigma^- \\
V_0^* + \Sigma^-
\end{pmatrix}:
\Sigma^- \in \R^{n}
\right\}.
\end{equation}

which includes all states that connect to $Q_0^*$ by Lax curves with negative speeds, i.e.,  $\mathcal{L}_{\lambda_1}\dots\mathcal{L}_{\lambda_n}$. Conversely, the set of states connecting to $Q_0^*$ by Lax curves with positive speeds, i.e.,  $\mathcal{L}_{\lambda_{n+1}}\dots \mathcal{L}_{\lambda_{2n}}$, is

\begin{equation}
\mathcal{L}^+(Q_0^*) \coloneq \left\{ Q_0^* + R_2^+ \Sigma^+:
\Sigma^+ \in \R^{n}
\right\}
\left\{
\begin{pmatrix}
U_0^* + (\sqrt{A_2})^{-1} \Sigma^+ \\
V_0^* + \Sigma^+
\end{pmatrix}:
\Sigma^+ \in \R^{n}
\right\}.
\end{equation}

Note that for the following discussion the Lax-curves with negative speeds $\mathcal{L}_{\lambda_1},\dots,\mathcal{L}_{\lambda_n}$ in (<ref>) as well as $\mathcal{L}^-$ in (<ref>) are defined for the relaxation system on the left half-axis, whereas the Lax-curves with positive speeds $\mathcal{L}_{\lambda_{n+1}},\dots,\mathcal{L}_{\lambda_{2n}}$ and $\mathcal{L}^+$ are defined for the relaxation system on the right half-axis.

Wave structure of the coupled relaxation system in the $x$-$t$-plane. The left trace data $Q_0^{-}$ are connected to the coupling data $Q_R$ by the Lax-curves of negative speeds, i.e., $\mathcal{L}_{\lambda_1}\dots\mathcal{L}_{\lambda_n}$ and the outgoing trace data $Q_0^{+}$ are connected to the coupling data $Q_L$ by the Lax-curves of positive speeds, i.e., $\mathcal{L}_{\lambda_{n+1}}\dots \mathcal{L}_{\lambda_{2n}}$. Coupling data on the incoming and outgoing edges are related by the coupling condition $\Psiq$.

We consider a situation, where we are given discrete data next to the coupling interface from a finite volume scheme. We refer to this so-called trace-data by $Q_0^- = (U_0^-, V_0^-)$ left from the interface and $Q_0^+ = (U_0^+, V_0^+)$ right from the interface, respectively. We denote boundary/coupling data for the left half axis $(-\infty, 0)$ by $Q_R=(U_R, V_R)$. Analogously, we denote coupling/boundary data for the right half axis $(0, \infty)$ by $Q_L=(U_L, V_L)$. If $Q_R$ and $Q_L$ constitute a solution of the half-Riemann problem they need to satisfy

\begin{equation}\label{eq:laxcondition}
Q_R \in \mathcal L^-(Q_0^-) \qquad \text{and} \qquad Q_L \in \mathcal L^+(Q_0^+).
\end{equation}

A Riemann solver (RS) for the relaxation system (<ref>) corresponding to the coupling condition (<ref>) is a mapping

\begin{equation}\label{eq:RSdef}
\mathcal{RS}: \R^{2n \times 2n}\rightarrow \R^{2n \times 2n}, \qquad (Q_0^-, Q_0^+) \rightarrow (Q_R, Q_L)
\end{equation}

that assigns coupling data satisfying the coupling condition and solving the half-Riemann problem to given trace data. It is defined by condition (<ref>) and

\begin{equation}\label{eq:psiq}
\Psiq[Q_R, Q_L] = 0
\end{equation}

for the coupling function introduced in (<ref>). We note that the well-posedness of the RS depends on the coupling condition. The relation between traces, coupling data and coupling condition in the relaxation system is illustrated in Figure <ref>.

§.§ Well-posedness

In this section we assume general coupling conditions of the form (<ref>) with $l=2n$, i.e.,

\begin{equation}\label{eq:psigeneral}
\Psiq[Q_R, Q_L] = 0, \qquad \Psi_Q:\R^{2n}\times \R^{2n} \rightarrow \R^{2n}.
\end{equation}

To construct a corresponding RS we aim to identify coupling data that not only satisfies (<ref>) but also solves the half-Riemann problem (<ref>). Condition (<ref>) allows us to parameterize the coupling data as

\begin{equation}\label{eq:couplingdataparameterized}
Q_R(\Sigma; Q_0) = Q_0^- + R^-_1 \Sigma^-, \qquad Q_L(\Sigma; Q_0) = Q_0^+ + R^+_2 \Sigma^+,
\end{equation}

using the notation

\begin{equation}\label{eq:Q0leftright}
Q_0 = \left( Q_0^-, Q_0^+ \right), \quad
\Sigma =
\left( \Sigma^-, \Sigma^+\right)
\end{equation}

and thus to introduce the parametrized form of the coupling function

\begin{equation}
\Psiqp[\Sigma; Q_0] \coloneq \Psi_{Q}[ Q_R(\Sigma; Q_0),  Q_L(\Sigma; Q_0)], \qquad \Psiqp:\R^{2n}\times \R^{4n} \rightarrow \R^{2n}.
\end{equation}

Given this parametrization we can interpret the RS (<ref>) for coupling condition (<ref>) as a routine that first identifies $\Sigma$ by solving the root problem $\Psiqp[\Sigma; Q_0]=0$ for given $Q_0$ and then outputs the data $Q_R(\Sigma; Q_0)$ and $Q_L(\Sigma; Q_0)$.
It is well posed iff for any trace data $Q_0$ the equation $\Psiqp[\Sigma; Q_0]=0$ has a unique solution $\Sigma = \Xi(Q_0)$. We state a result on its local regularity.

Let $\Psiq$ be continuously differentiable. Suppose that for a fixed $Q_0\in\R^{4n}$ there is a unique $\Sigma\in \R^{2n}$ satisfying $\Psiqp[\Sigma; Q_0]=0$ and the Jacobian $D_\Sigma \Psiqp[\Sigma; Q_0]$ is invertible. Then there is an open ball $B_{Q_0}\subset \R^{4n}$ around $Q_0$ such that the RS (<ref>) with coupling condition (<ref>) is well defined and continuously differentiable for all $\tilde Q_0=(\tilde Q_0^-, \tilde Q_0^+)^T\in B_{Q_0}$.

Due to the implicit function theorem there are open balls $B_\Sigma$ and $B_{Q_0}$ around $\Sigma$ and $Q_0$, respectively, and a unique continuously differentiable mapping $\Xi: B_{Q_0} \rightarrow B_{\Sigma}$ such that

\begin{equation}\label{eq:implicitxi}
\Psiqp[\tilde \Sigma; \tilde Q_0] = 0 \quad \text{iff}\quad \tilde \Sigma = \Xi(\tilde Q_0)
\end{equation}

for all $\tilde Q_0 \in B_{Q_0}$ and $\tilde \Sigma \in B_\Sigma$. The statement is obtained taking into account that the output of the RS (<ref>) is uniquely given by $Q_R(\Xi(\tilde Q_0); \tilde Q_0)$ and $Q_L(\Xi(\tilde Q_0); \tilde Q_0)$ for all $\tilde Q_0 \in B_{Q_0}$.

A further result is obtained applying the contraction mapping theorem to a fixed point iteration of the form

\begin{equation}\label{eq:fpiteration}
\Sigma^{n+1} = \Sigma^n - A(Q_0) \Psiqp[\Sigma^n; Q_0]
\end{equation}

with $A(Q_0)\in\R^{2n \times 2n}$ invertible.

Let $\mathcal{S}\subset \R^{4n}$ and $\mathcal{P} \subset \R^{2n}$ denote open sets. Suppose that for any $Q_0 \in \mathcal{S}$ there is a regular matrix $A(Q_0)\in\R^{2n \times 2n}$ such that

\begin{equation}\label{eq:fpselfmapping}
T_{Q_0}(\Sigma) = \Sigma - A(Q_0) \Psiqp[\Sigma; Q_0] \in \mathcal{P}
\end{equation}

for all $\Sigma \in  \mathcal{P}$ and $T_{Q_0}$ is Lipschitz continuous with Lipschitz constant $K < 1$.
Then $\mathcal{RS}:\mathcal{S} \rightarrow \R^{2n}\times \R^{2n}$ defined through (<ref>), (<ref>) and (<ref>) is well defined.

For $Q_0 \in \mathcal{S}$ the function $T_{Q_0}(\Sigma)$ maps from $\mathcal{P}$ to $\mathcal{P}$ due to (<ref>) and is a contraction mapping as it is Lipschitz continuous with sufficiently small Lipschitz constant. Thus, by the contraction mapping theorem, the sequence (<ref>) admits for any arbitrary $\Sigma^0\in \mathcal{P}$ a unique fixed point $\tilde \Sigma$ satisfying $\Psiqp[\tilde \Sigma; Q_0] =0$. This implies well-posedness of $\mathcal{RS}_{\Psiqp}$.

The Lipschitz continuity of $T_{Q_0}$ in Lemma <ref> is implied by the stronger condition

\begin{equation}\label{eq:fpderivativebound}
\sup_{\Sigma \in \mathcal{P}} \|I - A(Q_0)D_\Sigma \Psiqp[Q_0; \Sigma] \| < 1.
\end{equation}

We note that Lemmata <ref> and <ref> allow to freely choose the vector norm in $\R^n$ determining the induced matrix norm in Remark <ref>.

§.§ Linear coupling conditions

In this section we discuss a special case of (<ref>), namely affine linear coupling functions of the form

\begin{equation}\label{eq:linearcoupling}
\Psiq[Q_R, Q_L] = B_R Q_R - B_L Q_L - P
\end{equation}

for $B_R, B_L \in \R^{2n \times 2n}$ and $P \in \R^{2n}$. Under further conditions on the matrices we can provide explicit formulas of the coupling data  obtained by the corresponding RS (<ref>). We introduce the block matrices

\begin{equation}\label{eq:truncatedR}
\tilde R^- =
\begin{pmatrix}
R_1^- & 0
\end{pmatrix}
\begin{pmatrix}
- (\sqrt{A_1})^{-1} & 0 \\
I & 0
\end{pmatrix}
, \qquad
\tilde R^+ =
\begin{pmatrix}
0 & R_2^+
\end{pmatrix}
\begin{pmatrix}
0 &  (\sqrt{A_2})^{-1} \\
0 & I
\end{pmatrix}
\end{equation}

satisfying $\tilde R^- \Sigma = R_1^- \Sigma^-$ and $\tilde R^+ \Sigma = R_2^+ \Sigma^+$ assuming definition (<ref>).

Suppose the matrix $B_R \tilde R^- - B_L \tilde R^+$ is regular. Then the RS corresponding to the affine linear coupling function (<ref>) is well defined and yields the output

\begin{equation}\label{eq:affinelinearrsoutput}
Q_R = Q_0^- + \tilde R^- \, \Sigma, \qquad Q_L = Q_0^+ + \tilde R^+ \, \Sigma,
\end{equation}

where $\Sigma$ is the unique solution of the linear system

\begin{equation}\label{eq:sigmaexplicit}
\left( B_R \tilde R^- - B_L \tilde R^+\right) \Sigma = P + B_L Q_0^+ - B_R Q_0^-.
\end{equation}

Using the parametrization (<ref>) we rewrite the coupling condition as

\begin{align*}
0 &= \Psiq[Q_R(\Sigma; Q_0), Q_L(\Sigma; Q_0)] = B_R \left( Q_0^- + \tilde R^- \Sigma \right) - B_L \left( Q_0^+ + \tilde R^+ \Sigma \right) - P,
\end{align*}

which implies (<ref>).

In case of block matrices the system (<ref>) can be solved using a block LU decomposion.

Suppose the matrices $B_R$ and $B_L$ are given in block form as

\begin{equation}
B_R =
\begin{pmatrix}
B_R^{11} &  B_R^{12} \\[5pt]
B_R^{21} &  B_R^{22}
\end{pmatrix}, \qquad
B_L =
\begin{pmatrix}
B_L^{11} &  B_L^{12} \\[5pt]
B_L^{21} &  B_L^{22}
\end{pmatrix}
\end{equation}

with each block being a matrix in $\R^{n \times n}$. The system matrix in (<ref>) then takes the form

\begin{equation}
B_R \tilde R^- - B_L \tilde R^+  =
\begin{pmatrix}
B_R^{12} -  B_R^{11}(\sqrt{A_1})^{-1} &   -B_L^{12} - B_L^{11} (\sqrt{A_2})^{-1} \\[5pt]
B_R^{22} - B_R^{21}(\sqrt{A_1})^{-1} &  -B_L^{22} - B_L^{21} (\sqrt{A_2})^{-1}
\end{pmatrix}
\eqcolon
\begin{pmatrix}
B^{11} & B^{12} \\
B^{21} & B^{22}
\end{pmatrix}
= B.
\end{equation}

Suppose now that $B^{11}$ and $S \coloneq B^{22}-B^{21}(B^{11})^{-1} B^{12}$ are regular, then the RS corresponding to the coupling condition (<ref>) is well defined and coupling data can be computed using the inverse system matrix

\begin{equation}\label{eq:luinverse}
B^{-1} =     \begin{pmatrix}
(B^{11})^{-1} \left( I + B^{12} S^{-1} B^{21} (B^{11})^{-1}\right) &  -(B^{11})^{-1} B^{12}S^{-1} \\
-S B^{21} (B^{11})^{-1} & S^{-1}
\end{pmatrix}.
\end{equation}

If $B^{11}$ in Remark <ref> is not regular a block LU decomposition of $B$ and an inverse matrix similar to (<ref>) might still be computable after suitable rotation of the rows and columns.

Note that assuming the coupling condition (<ref>) the relaxation system (<ref>) is consistent (in the sense of Definition <ref>) with system (<ref>) under the coupling condition

\begin{equation}\label{eq:linearcouplinginu}
\Psi_U[U_R, U_L] = B_R
\begin{pmatrix}
U_R \\ f_1(U_R)
\end{pmatrix}
- B_L
\begin{pmatrix}
U_L \\ f_2(U_L)
\end{pmatrix}
- P = 0.
\end{equation}

[Kirchhoff conditions]
A canonical choice of coupling conditions for (<ref>) are the ones due to Kirchhoff, which require equality of the flux left and right from the interface, in all equations of the system, i.e.,

\begin{equation}\label{eq:kirchhoffcouplingcondition}
\Psiq^K[Q_R, Q_L] = S_1 Q_R - S_2 Q_L = 0.
\end{equation}

By choosing the diagonal entries sufficiently large (such that (<ref>) holds) we can achieve $A_1=A_2=:A$, which we assume in the following. Consequently, the quantities in (<ref>), (<ref>) and (<ref>) are independent of the half-axis considered and we therefore neglect the index $i$. Since $S$ is regular, the coupling condition (<ref>) simplifies to

\begin{equation}\label{eq:couplingconditionsimplified}
\Psi_Q^K[Q_R, Q_L] = Q_R - Q_L = 0
\end{equation}

and thus the coupling function is of the form (<ref>) for $B_R=B_L=I$ and $P=0$.

From (<ref>) we obtain

\[
\begin{pmatrix}
-\Sigma^- \\
\end{pmatrix}
\begin{pmatrix}
R^- & R^+
\end{pmatrix}^{-1}
(Q_0^- - Q_0^+)
\begin{pmatrix}
L^- \\
\end{pmatrix}
(Q_0^- - Q_0^+)
\]

and consequently the coupling data is given by

\[
Q_R = Q_0^- - R^- L^-  (Q_0^- - Q_0^+), \quad Q_R = Q_0^+ + R^+ L^+  (Q_0^- - Q_0^+),
\]

or equivalently

\begin{align}
U_R = U_L &= \frac 1 2 (U_0^- + U_0^+) + \frac 1 2 (\sqrt{A})^{-1}(V_0^- - V_0^+), \\
V_R = V_L &= \frac 1 2 (V_0^- + V_0^+) + \frac 1 2 \sqrt{A} (U_0^- - U_0^+).
\end{align}

The relation between trace and coupling data inferred by  (<ref>) and (<ref>) therefore defines a RS with easy-to-compute coupling data. In general the coupling condition (<ref>) does not lead to consistency with the original Kirchhoff coupling problem, in which (<ref>) is endowed with the coupling condition

\[
\Psi_U^K[U(t, 0^-), U(t, 0^+)] = F_1(U(t, 0^-)) -F_2(U(t, 0^+)) = 0
\]

This is due to the lower $n$ components in (<ref>) that in general impose additional conditions in the relaxation limit. Instead concluding from (<ref>) consistency is given in case of the coupling condition

\[
\Psi_U[U(t, 0^-), U(t, 0^+)] =
\begin{pmatrix}
U(t, 0^-) -  U(t, 0^+) \\ F_1(U(t, 0^-)) -F_2(U(t, 0^+)).
\end{pmatrix}
\]

§.§ Discrete relaxation limit

Similarly as for the continuous problem in Section <ref>, we can analyze the relaxation limit of the RS. Therefore, we assume that the trace data $Q_0=(Q_0^-, Q_0^+)$ can be asymptotically expanded around the state $Q_0^0=\left( Q_0^{0,-}, Q_0^{0,+} \right)$ as

\begin{equation}
Q_0^- = Q_0^{0,-} + \sum_{k=1}^\infty \varepsilon^k Q_0^{k,-}, \qquad Q_0^+ = Q_0^{0,+}  + \sum_{k=1}^\infty \varepsilon^k Q_0^{k,+}.
\end{equation}

Moreover, we assume in analogy to Lemma <ref> that $\Psiq$ is continuously differentiable and that for all $Q_0$ there is a unique $(\Sigma^-(Q_0), \Sigma^+(Q_0))=\Sigma\in \R^{2n}$ satisfying $\Psiqp[\Sigma; Q_0]=0$ and the Jacobian $D_\Sigma\Psiqp[\Sigma; Q_0]$ is invertible. From the parametrization (<ref>) we obtain

\begin{equation}\label{eq:limitpara}
Q_R(\Sigma; Q_0) = Q_0^{-} + R_1^- \Sigma^-(Q_0), \qquad Q_L(\Sigma; Q_0) = Q_0^{+} + R_2^+ \Sigma^+(Q_0),
\end{equation}

Due to the implicit function theorem, see the proof of Lemma <ref>, $\Sigma^-$ and $\Sigma^+$ are continuously differentiable functions in $Q_0$ and by Taylor expansion

\begin{equation}
\Sigma^\mp(Q_0) = \Sigma^\mp(Q_0^0) + \mathcal{O}(\varepsilon)
\end{equation}

holds. Hence, we obtain by substitution in (<ref>)

\begin{equation}
Q_R(\Sigma; Q_0) = Q_0^{0,-} + R_1^- \Sigma^-\left(Q_0^0\right) + \mathcal{O}(\varepsilon), \quad Q_L(\Sigma; Q_0) = Q_0^{0,+} + R_2^+ \Sigma^+\left(Q_0^0\right) + \mathcal{O}(\varepsilon)
\end{equation}

or, in different terms

\begin{equation}
\mathcal{RS}(Q_0^-, Q_0^+) = \mathcal{RS}\left(Q_0^{0,-}, Q_0^{0,+}\right) + \mathcal{O}(\varepsilon).
\end{equation}

Using a Taylor expansion of $\Psiq$ we further see that

\begin{equation*}
\Psiq[Q_R(\Sigma(Q_0); Q_0), Q_L(\Sigma(Q_0); Q_0)] = \Psiq[Q_R(\Sigma(Q_0^0); Q_0^0), Q_L(\Sigma(Q_0^0); Q_0^0)] + \mathcal{O}(\varepsilon).
\end{equation*}

and thus, in the limit $\varepsilon \rightarrow 0$ the coupling data satisfies the coupling condition, i.e.,

\begin{equation}
\Psiq[Q_R(\Sigma(Q_0^0); Q_0^0), Q_L(\Sigma(Q_0^0); Q_0^0)] = 0.
\end{equation}

§ THE CENTRAL SCHEME

The relaxation system in the 1-to-1 coupling case on the discretized real line.

We consider a discretization of the real line into uniform mesh cells $I_j=(x_{j-1/2}, x_{j+1/2})$ of width $\Delta x$ with origin located at the cell interface $x_{-1/2}$ as indicated by Figure <ref>. Furthermore, the time line is partitioned into the instances $t^n=\sum_{k=1}^n \Delta t^n$ for some time increments $\Delta t^n>0$, which for brevity we assume uniform and denote by $\Delta t$ throughout this section. Let $Q_j^n$, $U_j^n$, and $V_j^n$ denote a piecewise constant numerical solution of system (<ref>) in terms of cell averages over cell $I_j$ at time $t^n$.

§.§ Relaxation scheme

An asymptotic preserving scheme for the coupled relaxation system (<ref>) and coupling condition (<ref>) is obtained by combining the first order upwind discretization with an implicit-explicit time discretization as in <cit.>. Away from the interface the scheme reads

\begin{align}
U^{n+1}_{j} &= U^{n}_{j} - \frac{\Delta t}{2 \Delta x}\left( V^{n}_{j+1} - V^{n}_{j-1}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_i} \left( U^{n}_{j+1} - 2 U^{n}_{j} + U^{n}_{j-1}\right), \\
V^{n+1}_{j} &= V^{n}_{j} - \frac{\Delta t}{2 \Delta x} A_i \left( U^{n}_{j+1} - U^{n}_{j-1}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_i}\left( V^{ n}_{j+1} - 2 V^{n}_{j} + V^{n}_{j-1}\right) + \frac {\Delta t} \eps \left( F_i(U^{n+1}_{j}) - V^{n+1}_{j} \right),\label{eq:unsplitv0}
\end{align}

for $j\in \Z \setminus \{0, -1\}$, where $i=1$ if $j<0$ and $i=2$ if $j>0$. For the cell that is left next to the coupling interface we have

\begin{align}
U^{n+1}_{-1} &= U^{n}_{-1} - \frac{\Delta t}{2 \Delta x}\left( V^{n}_{R} - V^{n}_{-2}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_1} \left( U^{n}_{R} - 2 U^{n}_{-1} + U^{n}_{-2}\right), \\
V^{n+1}_{-1} &= V^{n}_{-1} - \frac{\Delta t}{2 \Delta x} A_1 \left( U^{n}_{R} - U^{n}_{-2}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_1}\left( V^{ n}_{R} - 2 V^{n}_{-1} + V^{n}_{-2}\right) + \frac {\Delta t} \eps \left( F_1(U^{n+1}_{-1}) - V^{n+1}_{-1} \right),\label{eq:unsplitvl}
\end{align}

and on the opposite site, for the cell right next to the coupling interface the scheme reads

\begin{align}
U^{n+1}_{0} &= U^{n}_{0} - \frac{\Delta t}{2 \Delta x}\left( V^{n}_{1} - V^{n}_{L}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_2} \left( U^{n}_{1} - 2 U^{n}_{0} + U^{n}_{L}\right), \\
V^{n+1}_{0} &= V^{n}_{0} - \frac{\Delta t}{2 \Delta x} A_2 \left( U^{n}_{1} - U^{n}_{L}\right) + \frac{\Delta t}{2 \Delta x} \sqrt{A_2}\left( V^{ n}_{1} - 2 V^{n}_{0} + V^{n}_{L}\right) + \frac {\Delta t} \eps \left( F_2(U^{n+1}_{0}) - V^{n+1}_{0} \right).\label{eq:unsplitv2}
\end{align}

The discrete coupling data employed in scheme (<ref>) is obtained applying the RS from Section <ref> to the cell averages adjacent to the interface, i.e.,

\begin{equation}
\mathcal{RS}( (U_{-1}^n, V_{-1}^n), (U_0^n, V_0^n)) := ((U_R^n, V_R^n), (U_L^n, V_L^n)).
\end{equation}

Scheme (<ref>) can be rewritten in the form

\begin{align}
U^{n+1}_{j} &= U^{n}_{j} - \frac{\Delta t}{\Delta x}\left( F_{j+1/2}^{n,-} - F_{j-1/2}^{n,+} \right) \label{eq:relaxationschemeu},\\
V^{n+1}_{j} &= V^{n}_{j} - \frac{\Delta t}{\Delta x}\left( G_{j+1/2}^{n,-} - G_{j-1/2}^{n,+} \right) + \frac{\Delta t}{\varepsilon} \left( F_i(U^{n+1}_j) - V^{n+1}_j)\right), \label{eq:relaxationschemev}
\end{align}

where $i=1$ if $j<0$ and $i=2$ if $j\geq 0$. We distinguish between incoming numerical fluxes from the left interface (indicated by a plus sign) and numerical fluxes outgoing through the right interface (indicated by a minus sign). Away from the interface they are equal as we have

\begin{align}
F_{j-1/2}^{n,+} = F_{j-1/2}^{n,-} =  \frac 1 2 \, (V_{j-1}^n + V_{j}^n)  - \frac 1 2  \sqrt{A_i} (U_j^n-U_{j-1}^n),\\
G_{j-1/2}^{n,+} = G_{j-1/2}^{n,-} =  \frac 1 2 A_i  (U_{j-1}^n + U_{j}^n)  - \frac 1 2  \sqrt{A_i} (V_j^n-V_{j-1}^n) \label{eq:vflux}
\end{align}

for $j\in \Z\setminus\{ 0\}$ and $i$ chosen according to the sign of $j$ as in (<ref>). At the interface the scheme is not necessarily conservative since the numerical fluxes are given by

\begin{align}
F_{-1/2}^{n,-} &=  \frac 1 2 \, (V_{-1}^n + V_{R}^n)  - \frac 1 2  \sqrt{A_1} (U_R^n-U_{-1}^n),\\
F_{-1/2}^{n,+} &=  \frac 1 2 \, (V_{L}^n + V_{0}^n)  - \frac 1 2  \sqrt{A_2} (U_0^n-U_{L}^n),\\
G_{-1/2}^{n,-} &=  \frac 1 2 A_1  (U_{-1}^n + U_{R}^n)  - \frac 1 2  \sqrt{A_1} (V_R^n-V_{-1}^n),\\
G_{-1/2}^{n,+} &=  \frac 1 2 A_2  (U_{L}^n + U_{0}^n)  - \frac 1 2  \sqrt{A_2} (V_0^n-V_{L}^n).
\end{align}

§.§ Relaxed scheme

In this section, using asymptotic analysis as done in Section <ref> and <ref> we consider the relaxation limit of the relaxation scheme. The cell averages $U_j^n$ and $V_j^n$ occuring in the schemes above depend on the relaxation rate $\varepsilon$. We assume that for all $j$ and $n$ they can be asymptotically expanded around the zero-relaxation states $U_j^{n,0}$ and $V_j^{n,0}$ for sufficiently small relaxation rates as

\begin{align}
U_j^{n} = U_j^{n,0} + \varepsilon U_j^{n,1} + \mathcal{O}(\varepsilon^2),\\
V_j^{n} = V_j^{n,0} + \varepsilon V_j^{n,1} + \mathcal{O}(\varepsilon^2).
\end{align}

At first we consider the relaxation limit away from the interface. We fix $j\in \Z\setminus \{ -1, 0\}$ and choose $i$ according to the sign of $j$ as in (<ref>) and then substitute the expansions (<ref>) into (<ref>). After Taylor expansion we obtain

\begin{equation}\label{eq:vexpansion}
\begin{split}
V^{n+1, 0}_{j} \left(1 + \frac{\Delta t}{\varepsilon} \right) &= V^{n, 0}_{j} - \frac{\Delta t}{\Delta x}\left( G_{j+1/2}^{n,+, 0} - G_{j-1/2}^{n,-, 0} \right) \\ &\quad + \frac{\Delta t}{\varepsilon} \left( F_i(U^{n+1, 0}_j)  + \varepsilon DF_i(U^{n+1, 0}_j)~U_j^{n+1,1} - \varepsilon V_j^{n,1} \right) + \mathcal{O}(\varepsilon),
\end{split}
\end{equation}

where $DF_i$  denotes the Jacobian of $F_i$ and  $G_{j-1/2}^{n,\mp, 0}$= $G_{j-1/2}^{n,\mp, 0}\left((U^{n,0}_{j-1}, U^{n,0}_{j-1}), (U^{n,0}_{j}, U^{n,0}_{j})  \right)$ the numerical flux equivalents to (<ref>) at the zero-relaxation state. Due to the expansion

\begin{equation}
\left( 1 + \frac{\Delta t}{\varepsilon} \right)^{-1} = \frac{\varepsilon}{\varepsilon + \Delta t} = \frac{\varepsilon}{\Delta t} + \mathcal{O}(\varepsilon^2)
\end{equation}

we obtain by multiplication in (<ref>)

\begin{equation}\label{eq:vlimitpre}
V^{n+1, 0}_{j} = F_i(U^{n+1, 0}_j) + \mathcal{O}(\varepsilon) \qquad \text{for all }n\in \N_0.
\end{equation}

Substituting now the asymptotic expansion (<ref>) into (<ref>) and taking into account (<ref>) we obtain

\begin{equation}\label{eq:asymptoticlimitscheme}
U^{n+1,0}_{j} = U^{n,0}_{j} - \frac{\Delta t}{\Delta x}\left( F_{j+1/2}^{n,+,0} - F_{j-1/2}^{n,-,0} \right)  + \mathcal{O}(\varepsilon),
\end{equation}

where the numerical fluxes take the form

\begin{equation}\label{eq:asymptoticfluxes}
F_{j-1/2}^{n,+,0} = F_{j-1/2}^{n,-,0} =  \frac 1 2 \, (F_i(U_{j-1}^{n,0}) + F_i(U_{j}^{n,0}))  - \frac 1 2  \sqrt{A_i} (U_j^{n,0}-U_{j-1}^{n,0}).
\end{equation}

In the relaxation limit the $\mathcal{O}(\varepsilon)$ terms in (<ref>) and  (<ref>) vanish, which gives rise to the relaxed scheme on both half-axes.

Next, we investigate the relaxation limit at the interface. We assume that in case of the asymptotic expansion (<ref>) the RS (<ref>) satisfies

\begin{equation}\label{eq:RSlimitassumption}
\mathcal{RS}( (U_{-1}^n, V_{-1}^n), (U_0^n, V_0^n)) =  \mathcal{RS}_{0}( (U_{-1}^{n,0}, V_{-1}^{n,0}), (U_0^{n,0}, V_0^{n,0})) + \mathcal{O}(\varepsilon)
\end{equation}

for some mapping

\[
\mathcal{RS}_{0}: \R^{2n} \times \R^{2n} \rightarrow \R^{2n} \times \R^{2n}.
\]

In Section <ref> we have shown that under suitable conditions on $\Psi_Q$ condition (<ref>) holds for $\mathcal{RS}_0=\mathcal{RS}$, i.e., the coupling in the limit is determined by the original RS. We note however, that in general $\mathcal{RS}_{0}$ may not necessarily be a RS in the sense it was introduced in Section <ref>.

Condition (<ref>) allows us to define the limit coupling data

\begin{equation}\label{eq:RSasymptotic}
\mathcal{RS}_{0}( (U_{-1}^{n,0}, V_{-1}^{n,0}), (U_0^{n,0}, V_0^{n,0})) =: ((U_R^{n,0}, V_R^{n,0}), (U_L^{n,0}, V_L^{n,0}))
\end{equation}

and therefore also the interface fluxes

\begin{align}
F_{-1/2}^{n,-,0} &=  \frac 1 2 \, (V_{-1}^{n,0} + V_{R}^{n,0})  - \frac 1 2  \sqrt{A_1} (U_R^{n,0}-U_{-1}^{n,0}),\\
F_{-1/2}^{n,+,0} &=  \frac 1 2 \, (V_{L}^{n,0} + V_{0}^{n,0})  - \frac 1 2  \sqrt{A_2} (U_0^{n,0}-U_{L}^{n,0}),\\
G_{-1/2}^{n,-, 0} &=  \frac 1 2 A_1  (U_{-1}^{n,0} + U_{R}^{n,0})  - \frac 1 2  \sqrt{A_1} (V_R^{n,0}-V_{-1}^{n,0}),\\
G_{-1/2}^{n,+, 0} &=  \frac 1 2 A_2  (U_{L}^{n,0} + U_{0}^{n,0})  - \frac 1 2  \sqrt{A_2} (V_0^{n,0}-V_{L}^{n,0}).
\end{align}

Combining these fluxes with the above arguments and using assumption (<ref>) it is now straightforward to show that (<ref>) and (<ref>) also hold for $j=-1$ and $j=0$.  Eventually we take the relaxation limit $\eps \rightarrow 0$, drop the index $0$ and obtain the coupled relaxed or central scheme for system (<ref>) and coupling condition (<ref>) in the relaxation limit, which reads

\begin{equation}\label{eq:conservativeform}
U_j^{n+1} = U_j^n - \frac{\Delta t }{\Delta x} \left( H_{j+1/2}^{n,-} -  H_{j-1/2}^{n,+}\right) \quad \text{for all }j \in \Z
\end{equation}

and employs the numerical fluxes

\begin{align}
H_{j-1/2}^{n,+} = H_{j-1/2}^{n,-} =  \frac 1 2 \, (F_i(U_{j-1}^n) + F_i(U_{j}^n))  - \frac 1 2  \sqrt{A_i} (U_j^n-U_{j-1}^n)
\end{align}

for $j\in \Z\setminus \{0\}$ and $i$ chosen according to the sign of $j$. At the interface we have

\begin{align}
H_{-1/2}^{n,-} &=  \frac 1 2 \, (F_1(U_{-1}^n) + V_{R}^n)  - \frac 1 2  \sqrt{A_1} (U_R^n-U_{-1}^n)\\
H_{-1/2}^{n,+} &=  \frac 1 2 \, (V_{L}^n + F_2(U_{0}^n))  - \frac 1 2  \sqrt{A_2} (U_0^n-U_{L}^n),
\end{align}

where the variables $V_R^n=V_R^n(U_{-1}^n, U_{0}^n)$, $V_L^n=V_L^n(U_{-1}^n, U_{0}^n)$ depend only on the cell averages of the state variable next to the interface as part of the coupling data

\begin{equation}\label{eq:RSlimit}
\mathcal{RS}_{0}( (U_{-1}^n, F_1(U_{-1}^n), (U_0^n, F_2(U_0^n))) =: ((U_R^n, V_R^n), (U_L^n, V_L^n)).
\end{equation}

Suppose that $\Psi_Q$ is continuously differentiable, for all $Q_0$ there is a unique $\Sigma\in\R^n$ satisfying $\Psiqp[\Sigma; Q_0]=0$ and the Jacobian $D_\Sigma\Psiqp[\Sigma; Q_0]$ is invertible. Then following from Section <ref> the RS for the coupling condition (<ref>) is well defined and $\mathcal{RS}_{0} \coloneqq \mathcal{RS}$ satisfies (<ref>) and defines a limit scheme.

If Kirchhoff conditions are assumed for $\Psi_Q$ together with $A_1=A_2=A$ as discussed in Example <ref> we can substitute the explicit formulas of the coupling data and obtain a conservative scheme of the form (<ref>) with the numerical fluxes

\begin{equation}
H_{j-1/2}^{n,+}=  H_{j-1/2}^{n,-} = \begin{cases}
\frac 1 2 \, (F_1(U_{j-1}^n) + F_1(U_{j}^n))  - \frac 1 2  \sqrt{A} (U_j^n-U_{j-1}^n)  &\text{if }j<0, \\[5pt]
\frac 1 2 \, ( F_1(U_{-1}^n) + F_2(U_{0}^n))  - \frac 1 2 \sqrt{A} (U_0^n-U_{-1}^n) &\text{if }j=0, \\[5pt]
\frac 1 2 \, ( F_2(U_{j-1}^n) + F_2(U_{j}^n))  - \frac 1 2 \sqrt{A} (U_j^n-U_{j-1}^n) &\text{if }j>0.
\end{cases}
\end{equation}

This follows from generalizing the scalar approach in <cit.> to the system case.

§ CASE STUDY: THE $P$-SYSTEM

In this section we apply the developed approach to a problem arising in gas transportation networks considered in <cit.>. Motivated by recent significant increases of gas-fired electrical generation in parts of the USA, see e.g. <cit.>, the modeling of high-pressure gas transportation networks needs to account for the influence of gas-fired power generators. The gas flow through pipelines has been effectively modeled using the compressible Euler equations <cit.>. Neglecting temperature dynamics the $p$-system is used as a simplified model of gas dynamics, see e.g., <cit.>, given by

\begin{align}
\ddt \rho + \ddx(\rho v) &= 0, \\
\ddt (\rho v) + \ddx(\rho v^2 + p(\rho)) &=0
\end{align}

and accounting for the conservation of density $\rho$ and momentum $\rho v$.  As closure for the system we impose the $\gamma$-law for the pressure function

\begin{equation}
p(\rho) = \alpha \rho^{\gamma}
\end{equation}

with parameters $\alpha>0$ and $\gamma\geq0$. To account for the gas-fired power generation we consider (<ref>) on both real half-axes coupled at an interface located at $x=0$ in analogy to (<ref>). In this setting a gas turbine acting as the power generator is assumed at the interface. The impact of the turbine on the gas network is modeled by conditions that leave the pressure of the transmitted gas constant but account for a loss in momentum <cit.> giving rise to the coupling conditions

\begin{align}
\rho_R v_R &= \rho_L v_L + \mathcal{E}, \label{eq:rhovjump}\\
p(\rho_R) &=  p(\rho_L) \label{eq:psystemjumpconditionpressure}
\end{align}

which govern the coupling function $\Psi_U$. In condition (<ref>) the parameter $\mathcal{E}\leq0$ controls the momentum outtake due to the turbine. As the density $\rho$ is assumed nonnegative we note that assuming $\gamma>0$ condition (<ref>) imposing continuity of the pressure at the interface is equivalent to

\begin{equation}\label{eq:psystemcontinuousdensity}
\rho_R = \rho_L.
\end{equation}

§.§ Relaxation form

Introducing the two auxiliary scalar variables $V_1$, $V_2$ and the relaxation rate $\varepsilon$ the relaxation system corresponding to the coupled $p$-system (<ref>) takes the form

\begin{align}
\ddt \rho + \ddx(V_1) &= 0, \\
\ddt (\rho v) + \ddx(V_2) &=0, \\
\ddt V_1 + a \ddx \rho &= \frac 1 \varepsilon (\rho v - V_1), \\
\ddt V_2 + a \ddx (\rho v)  &= \frac 1 \varepsilon (\rho v^2 + p(\rho) - V_2).
\end{align}

As system (<ref>) also its relaxation variant (<ref>) is imposed on both real half-axes coupled at an interface located at $x=0$. Compared to the general relaxation form (<ref>) the parameters $a_i^j$ are chosen independent of both the half-axis and the system component and set to the square of the maximal speed of sound, i.e.,

\[
a = \max_{\rho \in[0, \rho^\text{max}]} p^\prime(\rho).
\]

In vector form system (<ref>) uses the variables $U=(\rho, \rho v)$ and $V=(V_1, V_2)$ and the matrix $A=aI$.

We are interested in identifying a suitable coupling function $\Psi_Q$ for the coupled relaxation system (<ref>). The implied condition should be such that in the limit $\varepsilon \rightarrow 0$ the dynamics of the original coupled $p$-system given by (<ref>) and (<ref>) are attained. At first we look at linear coupling conditions of the type

\begin{align}
\rho_R &= \rho_L, \label{eq:relaxationrhocondition}\\
\rho_R v_R &= \rho_L v_L + \beta_1 \mathcal{E}, \label{eq:rhovjumpagain}\\
V_{1R} &= V_{1L} + \beta_2 \mathcal{E}, \label{eq:v1condition}\\
V_{2R} &= V_{2L} \label{eq:v2continuity}
\end{align}

for $\beta_1$, $\beta_2\in\{0,1\}$. As first approach motivated from (<ref>) and (<ref>) we consider the case $\beta_1=1$ and $\beta_2=0$.
The conditions (<ref>) (with $\beta_2=0$) and (<ref>) for the auxiliary variable $V$ have been chosen assuming continuity of the auxiliary variable.
Noting that in the relaxation limit it holds  $V_{1}= \rho v$, the jump condition (<ref>) can also be implemented in (<ref>). This motivates our second approach, for which we choose $\beta_1=0$ and $\beta_2=1$. Our third approach combines the two former ones by setting $\beta_1=\beta_2=1$.

By our analysis in Section <ref> and the linearity of the coupling approaches 1–3 the auxiliary variables satisfy in the relaxation limit

\begin{equation}
V_{1R} = \rho_R v_R, \quad V_{1L} = \rho_L v_L, \quad V_{2R} = \rho_R (v_R)^2 + p(\rho_R), \quad V_{2L} = \rho_L (v_L)^2 + p(\rho_L).
\end{equation}

This limit condition shows that both our first and our second coupling approach for system (<ref>) lead to a contradiction in (<ref>) and (<ref>) in the relaxation limit due to $\beta_1\neq \beta_2$. Our third approach corrects this, however due to (<ref>) this approach as well as the two former ones impose a condition in the limit that is not implied by the original coupling condition (<ref>). In more details, this additional condition imposes $v_R=v_L$, whenever $\rho_R v_R= \rho_Lv_L \neq 0$. Consequently, our three coupling approaches within (<ref>) cannot lead to consistency of the relaxation system (<ref>) with the p-system (<ref>) and coupling condition (<ref>) in the sense of Definition <ref>.

We aim to design a coupling condition that leads to consistency by taking our third approach as a starting point and deriving a suitable replacement for condition (<ref>). Therefore we assume (<ref>), (<ref>) as well as $\rho_R=\rho_L\neq 0$ and find

\begin{equation}
\rho_R v_R^2 + p(\rho_R) = \frac{\rho_R^2 v_R^2}{\rho_R} + p(\rho_R) = \frac{(\rho_L v_L+ \mathcal{E})^2 }{\rho_R} + p(\rho_L) = \rho_L v_L^2 + p(\rho_L) + \mathcal{E}\, \frac{ 2 \rho_L v_L + \mathcal{E}}{\rho_L}.
\end{equation}

This motivates the new condition

\begin{equation}\label{eq:nonlinearV2condition}
V_{2R} = V_{2L} + \mathcal{E}\, \frac{ 2 \rho_L v_L + \mathcal{E}}{\rho_L},
\end{equation}

which together with (<ref>), (<ref>), (<ref>) and $\beta_1=\beta_2=1$ determines our fourth coupling approach. By construction the relaxation system (<ref>) under this condition is consistent with the coupled $p$-system given by (<ref>) and (<ref>). Unlike the previous approaches this condition is nonlinear.

§.§ Riemann solvers

In this section we analyze the Riemann solvers for the different coupling conditions considered in Section <ref> and provide explicit forms for the implied coupling data,

Coupling approaches 1–3 The first three coupling conditions given by (<ref>) and different choices for $\beta_1$ and $\beta_2$ fit into the linear framework discussed in Section <ref> and correspond to

\begin{equation}
B_R = I, \qquad B_L = I, \qquad P =
\begin{pmatrix}
0 & \beta_1 \, \mathcal{E} & \beta_2 \, \mathcal{E} & 0
\end{pmatrix}^T
\end{equation}

in the coupling function (<ref>). To compute the coupling data from the RS we conclude from Remark <ref> that the system matrix in (<ref>) is of the form

\begin{equation}
B \coloneq B_R \tilde R^- - B_L \tilde R^+ =
- \begin{pmatrix}
( \sqrt{A})^{-1} & ( \sqrt{A})^{-1} \\
-I & I
\end{pmatrix}
\frac 1 2 \sqrt{A} & -\frac 1 2 I \\
\frac 1 2 \sqrt{A} & \frac 1 2 I
\end{pmatrix}^{-1}
\end{equation}

from which we determine the parameter $\Sigma^-$ and $\Sigma^+$ as

\begin{equation}\label{eq:psystemrelaxationcouplingsigma}
\begin{pmatrix}
\Sigma^- \\ \Sigma^+
\end{pmatrix}
= B^{-1} (P + Q_0^+ - Q_0^{-}), \qquad
\Sigma^\mp = - \frac 1 2 \sqrt{A}(P^U + U_0^+ - U_0^-) \pm \frac 1 2 (P^V + V_0^+ - V_0^-).
\end{equation}

Here $P^U=(0, \beta_1 \mathcal{E})$ and $P^V=(\beta_2 \mathcal{E}, 0)$ consist of the upper two and lower two entries in $P$, respectively. The trace data has the form

\begin{equation}\label{eq:psystemrelaxationtracedata}
Q_0^\mp =
\begin{pmatrix}
U_0^\mp \\ V_0^\mp
\end{pmatrix}
\begin{pmatrix}
\rho_0^\mp & \rho_0^\mp v_0^\mp & V_{1,0}^\mp & V_{2,0}^\mp
\end{pmatrix}^T.
\end{equation}

Eventually, the coupling data, for which we assume a notation analogue to (<ref>), is given by (<ref>) implying

\begin{align}
U_R &= \frac 1 2 (U_0^- + U_0^+ + P^U) - \frac 1 2 (\sqrt{A})^{-1}(V_0^+ - V_0^- + P^V), \\
U_L &= \frac 1 2 (U_0^- + U_0^+ - P^U) - \frac 1 2 (\sqrt{A})^{-1}(V_0^+ - V_0^- + P^V), \\
V_R &= \frac 1 2 \sqrt{A} (U_0^- - U_0^+ - P^U) + \frac 1 2 (V_0^+ + V_0^- + P^V),\\
V_L &= \frac 1 2 \sqrt{A} (U_0^- - U_0^+ - P^U) + \frac 1 2 (V_0^+ + V_0^- - P^V).
\end{align}

Due to Remark <ref>, which applies here since $\Psi_Q$ is affine linear, the coupling data in the relaxation limit is determined by the same RS as in case of $\varepsilon>0$. Thus the data (<ref>) also satisfies the coupling condition (<ref>) in the limit $\varepsilon \rightarrow 0$.

Coupling approach 4 We next consider our fourth coupling approach, which introducing the notation $g^V(U_L) = \mathcal{E} (0, (2 \rho_L v_L + \mathcal{E})/\rho_L)^T$ can be written as

\begin{equation}
\begin{pmatrix}
U_R \\ V_R
\end{pmatrix}
\begin{pmatrix}
U_L \\ V_L
\end{pmatrix}
\begin{pmatrix}
P^U \\ P^V
\end{pmatrix}
\begin{pmatrix}
0 \\ g^V(U_L)
\end{pmatrix}.
\end{equation}

To obtain coupling data we substitute the parametrization (<ref>) and obtain

\begin{equation}\label{eq:solve4parametrization}
\begin{pmatrix}
- (\sqrt{A})^{-1} \\ I
\end{pmatrix}
\Sigma^- -
\begin{pmatrix}
(\sqrt{A})^{-1} \\ I
\end{pmatrix}
\Sigma^+ -
\begin{pmatrix}
0 \\ g^V(U_0^+ + (\sqrt{A})^{-1} \Sigma^+)
\end{pmatrix}
\begin{pmatrix}
U_0^+ - U_0^- + P^U \\ V_0^+ - V_0^- + P^V
\end{pmatrix}.
\end{equation}

From the first two components of (<ref>) we get

\begin{equation}\label{eq:solve4Ucomponents}
\Sigma^- =   \sqrt{A}(U_0^- -U_0^+ - P^U) -\Sigma^+,
\end{equation}

which substituted into the components three and four in (<ref>) yields

\begin{equation}\label{eq:solve4Vcomponents}
- 2 \Sigma^+ + \sqrt{A}(U_0^- -U_0^+ - P^U) - g^V(U_0^+ + (\sqrt{A})^{-1} \Sigma^+) =  V_0^+ - V_0^- + P^V.
\end{equation}

Using the notation $\Sigma^+=(\sigma_1^+, \sigma_2^+)^T$ we note that the first component of (<ref>) determines $\sigma_1^+$, which in turn substituted into the second component of (<ref>) determines $\sigma_2^+$. Eventually $\Sigma^-$ is obtained from (<ref>) and therefore a RS for approach four is defined. By the verified well-posedness of the RS and the differentiability of the corresponding coupling function Remark <ref> also applies to our fourth coupling approach and thus the above procedure also determines suitable coupling data in the relaxation limit.

§.§ Numerical experiments

We consider a numerical test case employing constant initial data $\rho_0=\rho_0 v_0=1$ and pressure given by $p(\rho)= 146\,820.4 \rho$, see Section 5.2.1 in <cit.>. As computational domain we choose $[-200,200]$ and impose homogeneous Neumann boundary conditions. The momentum outtake at the interface is assumed time dependent and given by

\begin{equation}
-\mathcal{E}(t) =
\begin{cases}
\min\{ 0.6t, 3t\} & \text{if }0\leq t < 0.
\\ \max \{ 0, -3t + 1.5\} & \text{if } t \geq 0.3.
\end{cases}
\end{equation}

A reference solution for the coupled $p$-system (<ref>) and coupling condition (<ref>) is presented in <cit.>.

We employ the central scheme introduced in Section <ref> to numerically compute solutions for system (<ref>) in the limit $\varepsilon \rightarrow 0$ for the different coupling approaches considered in Section <ref>. In the numerical computations we employ 1000 uniform mesh cells and time increments chosen according to the CFL condition

\begin{equation}
\Delta t = \text{CFL} \frac{\Delta x}{\sqrt{a}}
\end{equation}

and Courant number $\text{CFL}=0.49$. Numerical results are presented in terms of density, momentum and pressure at the time instances $t=0.0716$, $0.2864$ and $0.55$.

group style=group size=3 by 3,
horizontal sep = .8 cm,
vertical sep = .6 cm,
xlabels at=edge bottom,
xticklabels at=edge bottom,
yticklabels at=edge right,
width = .355 ,
height = .22 ,
every tick label/.append style=font=,
title style=font=,
label style=font=,
[title=$t=0.0716$, ylabel=$\rho$, ymin=0.9995, ymax=1.0001]
[] table [x index=0, y index=1] input/psystem11jumpCPL1_2.dat;
[title=$t=0.2864$, ymin=0.9995, ymax=1.0001]
[] table [x index=0, y index=1] input/psystem11jumpCPL1_3.dat;
scaled y ticks=manual:$\cdot 10^{-4} + 1$
ymin=0.9995, ymax=1.0001]
[] table [x index=0, y index=1] input/psystem11jumpCPL1_4.dat;
[ylabel=$\rho v$, ymin=0.8, ymax=1.2]
[] table [x index=0, y index=2] input/psystem11jumpCPL1_2.dat;
[ymin=0.8, ymax=1.2]
[] table [x index=0, y index=2] input/psystem11jumpCPL1_3.dat;
[ymin=0.8, ymax=1.2]
[] table [x index=0, y index=2] input/psystem11jumpCPL1_4.dat;
[ylabel=$p$, ymin=146760, ymax=146830]
[] table [x index=0, y index=3] input/psystem11jumpCPL1_2.dat;
[ymin=146760, ymax=146830]
[] table [x index=0, y index=3] input/psystem11jumpCPL1_3.dat;
[ymin=146760, ymax=146830,
scaled y ticks=manual:$+ 146\,820$
[] table [x index=0, y index=3] input/psystem11jumpCPL1_4.dat;

Numerical solution for system (<ref>) in the relaxation limit using coupling approach 1 given by (<ref>), $\beta_1=1$ and $\beta_2=0$. Peaks appears at the interface ($x=0$). The coupled relaxed scheme on $1000$ mesh cells with $\text{CFL}=0.49$ has been used.

In Figure <ref> we show results obtained by the first coupling approach given by (<ref>), $\beta_1=1$ and $\beta_2=0$. In the numerical results we see a peak in density, momentum and pressure forming and vanishing at the interface as time evolves. Comparing to the reference solution this coupling condition imposed to the coupled relaxation system clearly cannot yield the solution of the original coupled $p$-system (<ref>).

group style=group size=3 by 3,
horizontal sep = .8 cm,
vertical sep = .6 cm,
xlabels at=edge bottom,
xticklabels at=edge bottom,
yticklabels at=edge right,
width = .355 ,
height = .22 ,
every tick label/.append style=font=,
title style=font=,
label style=font=
[title=$t=0.0716$, ylabel=$\rho$, ymin=0.9999, ymax=1.0014]
[] table [x index=0, y index=1] input/psystem11jumpCPL2_2.dat;
[title=$t=0.2864$, ymin=0.9999, ymax=1.0014]
[] table [x index=0, y index=1] input/psystem11jumpCPL2_3.dat;
[title=$t=0.55$, ymin=0.9999, ymax=1.0014,
scaled y ticks=manual:$\cdot 10^{-3} + 1$
, ytick=1, 1.0005, 1.001]
[] table [x index=0, y index=1] input/psystem11jumpCPL2_4.dat;
[ylabel=$\rho v$, ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL2_2.dat;
[ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL2_3.dat;
[ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL2_4.dat;
[ylabel=$p$, ymin=146810, ymax=147000]
[] table [x index=0, y index=3] input/psystem11jumpCPL2_2.dat;
[ymin=146810, ymax=147000]
[] table [x index=0, y index=3] input/psystem11jumpCPL2_3.dat;
[ymin=146810, ymax=147000,
scaled y ticks=manual:$+ 146\,820$
, ytick=146820, 146895, 146970]
[] table [x index=0, y index=3] input/psystem11jumpCPL2_4.dat;

Numerical solution for system (<ref>) in the relaxation limit using coupling approach 2 given by (<ref>), $\beta_1=0$ and $\beta_2=1$. Artifacts at the interface ($x=0$) are visible. The coupled relaxed scheme on $1000$ mesh cells with $\text{CFL}=0.49$ has been used.

Figure <ref> shows the numerical solution in case of the second coupling approach given by (<ref>), $\beta_1=0$ and $\beta_2=1$.
While this solution qualitatively coincides with the reference solution for $x\neq0$ it exhibits layer artifacts at the interface at the time instances $0.0716$ and $0.2864$.

[spy using outlines=rectangle, red, magnification=5,
size=0.7cm, connect spies]
group style=group size=3 by 3,
horizontal sep = .8 cm,
vertical sep = .6 cm,
xlabels at=edge bottom,
xticklabels at=edge bottom,
yticklabels at=edge right,
width = .355 ,
height = .22 ,
every tick label/.append style=font=,
title style=font=,
label style=font=
[title=$t=0.0716$, ylabel=$\rho$, ymin=0.9999, ymax=1.0014]
[] table [x index=0, y index=1] input/psystem11jumpCPL3_2.dat;
[title=$t=0.2864$, ymin=0.9999, ymax=1.0014]
[] table [x index=0, y index=1] input/psystem11jumpCPL3_3.dat;
[gray, dotted] table [x index=0, y index=1] input/psystem11jumpCPL4_3.dat;
(spypoint) at (axis cs:0,1.00077);
(magpoint) at (axis cs:100,1.001);
[red] on (spypoint) in node at (magpoint);
[title=$t=0.55$, ymin=0.9999, ymax=1.0014,
scaled y ticks=manual:$\cdot 10^{-3} + 1$
, ytick=1, 1.0005, 1.001]
[] table [x index=0, y index=1] input/psystem11jumpCPL3_4.dat;
[ylabel=$\rho v$, ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL3_2.dat;
[ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL3_3.dat;
[gray, dotted] table [x index=0, y index=2] input/psystem11jumpCPL4_3.dat;
[ymin=0.65, ymax=1.35]
[] table [x index=0, y index=2] input/psystem11jumpCPL3_4.dat;
[ylabel=$p$, ymin=146810, ymax=147000]
[] table [x index=0, y index=3] input/psystem11jumpCPL3_2.dat;
[ymin=146810, ymax=147000]
[] table [x index=0, y index=3] input/psystem11jumpCPL3_3.dat;
[gray, dotted] table [x index=0, y index=3] input/psystem11jumpCPL4_3.dat;
[ymin=146810, ymax=147000,
scaled y ticks=manual:$+ 146\,820$
, ytick=146820, 146895, 146970]
[] table [x index=0, y index=3] input/psystem11jumpCPL3_4.dat;

Numerical solution for system (<ref>) in the relaxation limit using coupling approaches 3 (solid line) given by (<ref>) and $\beta_1=\beta_2=1$ and 4 (dotted line) given by (<ref>), (<ref>), (<ref>), (<ref>) and $\beta_1=\beta_2=1$. The coupled relaxed scheme on $1000$ mesh cells with $\text{CFL}=0.49$ has been used. Except for a jump of small magnitude occuring in approach 3 at the interface both solutions are indistinguishable.

Eventually, numerical results for the third coupling approach given by (<ref>) and $\beta_1=\beta_2=1$ are shown in Figure <ref>. Using this condition the solution of the coupled relaxation system qualitatively matches the one of the original system (<ref>) with coupling condition (<ref>). This indicates that the intrinsic contradiction in the relaxation limit between (<ref>) and (<ref>) in both coupling approaches one and two due to $\beta_1\neq \beta_2$ deteriorates the quality of the relaxation approach as an approximation of the original problem. We note hence, that consistency of the relaxation approach (<ref>) to system (<ref>) as we have verified for our fourth coupling approach has not been necessary to achieve an accurate numerical solution. This is presumably due to the small magnitude of the additional term $\mathcal{E} (2 \rho_L v_L + \mathcal{E})/\rho_L$ in (<ref>) that has been added in the fourth coupling approach, which compared to the magnitude of $V_{2R}$ and $V_{2L}$ is smaller by a factor of $10^{-6}$. A magnification at the interface of the solution corresponding to approach 3 reveals a jump of low magnitude not present in the reference solution. The solution for the fourth coupling approach given by (<ref>), (<ref>), (<ref>), (<ref>) and $\beta_1=\beta_2=1$ also shown in Figure <ref> and being qualitatively almost similar to the one of coupling approach 3 does not exhibit this artifact and therefore achieves the best fit to the reference solution.

group style=group size=2 by 1,
horizontal sep = 1.2 cm,
vertical sep = .6 cm,
xticklabels at=edge bottom,
xlabels at=edge bottom,
width =.5 ,
height = .22 ,
xmin=0, xmax=0.55,
every tick label/.append style=font=,
label style=font=,
title style=font=,
legend columns=-1,
legend style=font=, fill=none, /tikz/every even column/.append style=column sep=.25cm
[title=coupling error concerning (<ref>), ylabel=$E^1$]
[color=cyan, thick] table [x index=0, y index=2] input/psystem11jumpCPL_cplerror_1.dat;
[color=olive, thick, dashed] table [x index=0, y index=2] input/psystem11jumpCPL_cplerror_2.dat;
[color=red, thick] table [x index=0, y index=2] input/psystem11jumpCPL_cplerror_3.dat;
[color=black, thick, dashed] table [x index=0, y index=2] input/psystem11jumpCPL_cplerror_4.dat;
[title=coupling error concerning (<ref>), ylabel=$E^2$, legend to name=leg:cplapproaches]
[color=cyan, thick] table [x index=0, y index=1] input/psystem11jumpCPL_cplerror_1.dat;
[color=olive, thick, dashed] table [x index=0, y index=1] input/psystem11jumpCPL_cplerror_2.dat;
[color=red, thick] table [x index=0, y index=1] input/psystem11jumpCPL_cplerror_3.dat;
[color=black, thick, dashed] table [x index=0, y index=1] input/psystem11jumpCPL_cplerror_4.dat;
coupling 1, coupling 2, coupling 3, coupling 4

Coupling errors $E^1$ (left) and $E^2$ (right) defined in (<ref>) over the time interval $[0,0.55]$ for the four considered coupling approaches. Errors for numerical solutions computed by the central scheme on 1000 mesh cells using $\text{CFL}=0.49$ are considered. Both coupling approaches 3 and 4 achieve significantly lower errors than approaches 1 and 2 concerning (<ref>). Approach 4 reduces the error concerning (<ref>) over approaches 1–3.

For an insight into the approximation quality of the central scheme at the interface we consider the coupling errors with respect to the original coupling conditions (<ref>) and (<ref>) defined as

\begin{equation}\label{eq:couplingerrors}
E^1(t^n) = \left| (U_{-1}^n - U_{0}^n)_2 - \mathcal{E} \right|, \qquad E^2(t^n) = \left|(U_{-1}^n - U_{0}^n)_1 \right|
\end{equation}

Using the numerical solutions by the central scheme over 1000 mesh cells and employing Courant number $\text{CFL}=0.49$ we present these errors over time in Figure <ref>. In case of the error concerning the density condition (<ref>) the temporal errors follow the monotony of the momentum outtake $\mathcal{E}$. Approaches 4 and 1 here yield lower errors than approaches 2 and 3. Concerning the momentum condition (<ref>) both approaches 1 and 2 yield coupling errors in the magnitude of the momentum outtake and are clearly outperformed by approaches 3 and 4 both achieving similarly small errors.

Mesh convergence of the coupling error $L^1(E^1)$ over time for the considered coupling approaches. First order error convergence for approaches 3 and 4 is indicated.
cells     coupling 1     EOC     coupling 2     EOC     coupling 3     EOC     coupling 4     EOC
100     \(8.999 \times 10^{-2}\)         \(9.155 \times 10^{-2}\)
\(1.278 \times 10^{-2}\)         \(1.278 \times 10^{-2}\)

200     \(8.999 \times 10^{-2}\)     0.00     \(9.040 \times 10^{-2}\)     0.02
\(6.325 \times 10^{-3}\)     1.01     \(6.324 \times 10^{-3}\)     1.01

400     \(8.999 \times 10^{-2}\)     0.00     \(9.011 \times 10^{-2}\)     0.00
\(3.147 \times 10^{-3}\)     1.01     \(3.146 \times 10^{-3}\)     1.01

800     \(8.999 \times 10^{-2}\)     0.00     \(9.004 \times 10^{-2}\)     0.00
\(1.569 \times 10^{-3}\)     1.00     \(1.569 \times 10^{-3}\)     1.00

1600     \(8.999 \times 10^{-2}\)     0.00     \(9.003 \times 10^{-2}\)     0.00
\(7.838 \times 10^{-4}\)     1.00     \(7.837 \times 10^{-4}\)     1.00

To investigate the scaling of the coupling errors with the mesh resolution we compute numerical solutions for increasing numbers of mesh cells and consider the errors (<ref>) in the discrete $L^1$ norm over time, i.e.,

\[
L^1(E^i)= \sum_{k=0}^{k_T} \Delta t E^i(t^k),\qquad k_T = \lceil 0.55/\Delta t\rceil, \qquad i=1,2.
\]

Table <ref> shows the error $L^1(E^1)$ for the considered coupling approaches after mesh refinement along with the corresponding experimental order of convergence (EOC)[The EOC is computed by the formula $\text{EOC}=\log_2(E^1/E^2)$ with $E^1$ and $E^2$ denoting the error in two consecutive lines of the table.]. In case of coupling approaches 1 and 2 the errors are not significantly decreasing when refining the mesh. Conversely mesh convergence of order one can be clearly observed for the coupling approaches 3 and 4. Moreover, a minor advantage of coupling approach 4 over coupling approach 3 can be seen from the errors in Table <ref>.

Mesh convergence of the coupling error $L^1(E^2)$ over time for the considered coupling approaches. First order error convergence is only indicated for approach 4.
cells     coupling 1     EOC     coupling 2     EOC     coupling 3     EOC     coupling 4     EOC
100     \(3.131 \times 10^{-7}\)         \(9.173 \times 10^{-7}\)
\(1.224 \times 10^{-6}\)         \(8.685 \times 10^{-8}\)

200     \(3.083 \times 10^{-7}\)     0.02     \(9.179 \times 10^{-7}\)     0.00
\(1.224 \times 10^{-6}\)     0.00     \(4.307 \times 10^{-8}\)     1.01

400     \(3.070 \times 10^{-7}\)     0.01     \(9.179 \times 10^{-7}\)     0.00
\(1.224 \times 10^{-6}\)     0.00     \(2.143 \times 10^{-8}\)     1.01

800     \(3.067 \times 10^{-7}\)     0.00     \(9.179 \times 10^{-7}\)     0.00
\(1.224 \times 10^{-6}\)     0.00     \(1.069 \times 10^{-8}\)     1.00

1600     \(3.066 \times 10^{-7}\)     0.00     \(9.179 \times 10^{-7}\)     0.00
\(1.224 \times 10^{-6}\)     0.00     \(5.336 \times 10^{-9}\)     1.00

The error $L^1(E^2)$ after mesh refinement is presented in Table <ref>. In case of the linear coupling approaches no significant variation with the mesh resolution can be observed while the errors are of small magnitude. Convergence with respect to the mesh resolution is only observed in case of coupling approach 4 and consistency with the original problem.

§ CONCLUSION

We have presented a novel framework to handle systems of conservation laws coupled at an interface using a relaxation approach. When systems of conservation laws are coupled usually information on the Lax-curves of the underlying systems as well as complex and error-prone computations are required to identify suitable coupling data. The relaxation approach avoids this complication by replacing the nonlinear systems of conservation laws by systems of linear balance laws. We have introduced Riemann solvers for the new approach, discussed their well-posedness and proposed a finite volume scheme that can be generally applied to coupled hyperbolic problems.

The approach requires a suitable adjustment of the original coupling condition to account for the modified systems. Considering the $p$-system in a case study we have illustrated how the choice of this condition affects the numerical solutions. While accurate results could be obtained using simple linear ad-hoc conditions, best results that also exhibit convergence at the interface have been achieved constructing a coupling condition both leading to a well-defined Riemann solver and being in the relaxation limit consistent with the original coupling condition. Motivated by these results we suggest the same strategy for the design of the relaxation coupling condition also for different coupled  hyperbolic systems.

§ STATEMENTS AND DECLARATIONS

Compliance with Ethical Standards The authors state that they comply with the ethical standards.

Funding The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for the financial support through 320021702/GRK2326,  333849990/IRTG-2379, B04, B05 and B06 of 442047500/SFB1481, HE5386/18-1,19-2,22-1,23-1,25-1, ERS SFDdM035 and under Germany’s Excellence Strategy EXC-2023 Internet of Production 390621612 and under the Excellence Strategy of the Federal Government and the Länder. Support through the EU DATAHYKING is also acknowledged.

Conflict of Interest On behalf of all authors, the corresponding author states that there is no conflict of interest.

Ethical Approval Not applicable.


# Knowledge-Infused LLM-Powered Conversational Health Agent: A Case Study for
Diabetes Patients

Mahyar Abbasian1, Zhongqi Yang1, Elahe Khatibi1, Pengfei Zhang1, Nitish
Nagesh1,
Iman Azimi1, Ramesh Jain1, and Amir M. Rahmani1,2
1Department of Computer Science, University of California, Irvine
2School of Nursing, University of California, Irvine
{abbasiam, zhongqy4, ekhatibi, pengfz5, nnagesh1, azimii, rcjain,
<EMAIL_ADDRESS>
###### Abstract

Effective diabetes management is crucial for maintaining health in diabetic
patients. Large Language Models (LLMs) have opened new avenues for diabetes
management, facilitating their efficacy. However, current LLM-based approaches
are limited by their dependence on general sources and lack of integration
with domain-specific knowledge, leading to inaccurate responses. In this
paper, we propose a knowledge-infused LLM-powered conversational health agent
(CHA) for diabetic patients. We customize and leverage the open-source openCHA
framework, enhancing our CHA with external knowledge and analytical
capabilities. This integration involves two key components: 1) incorporating
the American Diabetes Association dietary guidelines and the Nutritionix
information and 2) deploying analytical tools that enable nutritional intake
calculation and comparison with the guidelines. We compare the proposed CHA
with GPT4. Our evaluation includes 100 diabetes-related questions on daily
meal choices and assessing the potential risks associated with the suggested
diet. Our findings show that the proposed agent demonstrates superior
performance in generating responses to manage essential nutrients.

###### Index Terms:

LLMs, Knowledge Graph, Diabetes, Nutrition Therapy, Health Agents.

##  I Introduction

Effective diabetes management plays a pivotal role in maintaining optimal
health in individuals diagnosed with diabetes. Given the widespread prevalence
of diabetes and its significant impact on both global healthcare
infrastructures and personal health outcomes, the availability and utilization
of robust management strategies are essential [1, 2, 3]. Such services require
holistic approaches integrating a healthy lifestyle, nutritious diet, and
physical activity assessment [4]. Particularly, a critical aspect is dietary
regulation, which plays a direct role in controlling blood glucose levels and
influencing the progression of the disease [5, 6]. The recent surge in
advanced technologies, particularly in Large Language Models (LLMs), has
substantially enhanced the accessibility and efficacy of diabetes management
as transformative educational and assistive tools, enabling patients to
acquire comprehensive knowledge and bridge the gap in diabetes self-management
[7, 8, 9].

Recent studies have investigated the implementation of LLMs on diabetes
management [7, 8, 9]. For example, [7] evaluated the use of ChatGPT for
Diabetes Self-Management and Education. They leveraged ChatGPT to answer a
range of diabetes self-management questions, covering four areas:
diet/exercise, glucose level education, insulin storage, and administration.
ChatGPT responded to the inquiries, showcasing a structured response style and
offering guidance that is understandable to laymen. However, in some specific
contexts, such as devising dietary plans, it still requires extra prompts to
comprehensively produce guidelines for tasks like insulin administration.
Additionally, Yang et al. [8] introduced ChatGLM to provide diabetes treatment
strategies. The proposed model was fine-tuned through P-tuning [10] and LoRA
[11] techniques with Electronic Health Record (EHR) of patients with diabetes.
ChatGLM showed proficiency in generating treatment recommendations for most
cases. However, its dependency on the specific data used for fine-tuning could
result in misleading outputs, potentially harmful in the absence of certain
essential external diabetes knowledge.

The lack of integration with verified domain-specific knowledge underscores a
significant gap in existing approaches, compromising the accuracy and
reliability of their outputs. Unlike systems anchored in specific knowledge
bases, existing LLM-based approaches fall into limitations in their reliance
on general information sources rather than specialized for diabetes. For
instance, ChatGPT’s foundation on a general database – rather than medically
specialized and validated ones – might account for its shortcomings in
grasping medical nuances. Furthermore, due to the model’s inability to
validate the reliability of its information, they may cause “hallucination,”
generating inaccurate responses presented persuasively and fluently. This
problem could easily mislead individuals lacking prior knowledge on the topic
[7].

Figure 1: An overview of openCHA [12] framework.

We believe that the integration of external knowledge bases with LLMs presents
a robust solution to current challenges by enhancing the models’ capacity to
access reliable information. To achieve this, Conversational Health Agents
(CHAs) can play an important role [13]. CHAs are conversational systems that
provide healthcare services, such as assistance and diagnosis, leveraging
agents as the decision-making core. openCHA [12], an open-source CHA
framework, provides a flexible framework for implementing LLM-based diabetes
management through its adept integration of diverse external sources. openCHA
possesses the capability to integrate and infuse knowledge into chatbots,
while also integrating external analytical methods for data analysis. Figure 1
shows an overview of this integration process, wherein an orchestrator –
powered by an LLM-based planner and an LLM-based response generator – engages
with various sources to gather the required information. This collected
information is then used to craft a response to the user’s query.

This paper presents a knowledge-based diabetes management system that is
implemented through the utilization of LLM-powered CHAs. We customize the
open-source openCHA framework for diet assessment by integrating external
knowledge – including the diet guidelines from the American Diabetes
Association reports and the nutritional data from the Nutritionix knowledge
base – into LLMs. We also incorporate an analysis tool to precisely calculate
the total daily nutritional intake and compare it against the guidelines. We
evaluate the proposed approach with real-world diabetes-related questions
focusing on daily meals in comparison to GPT4. The evaluation is centered on
the assessment of the potential risk associated with the recommended dietary
options.

##  II Method

We develop an LLM-based CHA, designed to assess the risk associated with daily
food intake according to recommended nutritional thresholds. We leverage
openCHA [12], as a foundational framework, for our development. The proposed
CHA aims to interact with users with diabetes, get their daily food intake
information in the form of a conversation, and use a food knowledge base as
well as guidelines to answer users in a more reliable way.

Our proposed CHA contains three main parts as Interface, Orchestrator, and
External Sources (see Figure 2). The Interface serves as a connection point
between users and our framework, facilitating interactions through text using
a web chat interface and relaying user queries to the Orchestrator.

Figure 2: LLM-based CHA for diabetes management enabled by the openCHA
framework.

The Orchestrator serves as the core component of the CHA and is responsible
for problem-solving, planning, executing actions, and generating appropriate
responses based on user queries. It operates on the principles of the
Perceptual Cycle Model [14] and involves perceiving, transforming, and
analyzing the world (i.e., input query and metadata). It interacts with
external sources, acquires information, performs data integration, and
extracts insights. The Orchestrator comprises four major components: 1) Task
Planner: to perform decision-making and planning, 2) Task Executor: to enable
actuation and data conversion, 3) Data Pipe: to act as a repository for
acquired metadata and intermediate data, 4) and Response Generator: to refine
information and delivers conclusive responses. For the Orchestrator setup, we
use the OpenAI’s [15] GPT-3.5-turbo model as our base LLM, and the Tree of
Thought [16] prompting technique for Planning.

External Sources are crucial for acquiring valuable information from diverse
origins. Through the openCHA framework [12], we integrate two primary external
sources to the CHA as: 1) Knowledge Base and 2) AI and Analysis Models (see
Figure 2).

Knowledge Base: serves to fetch current healthcare information from credible
sources. We developed a task focused on retrieving nutritional information
from Nutritionix for various foods, enhancing the up-to-dateness of the
responses. This task takes a food query as input and uses the NutritioniX API
[17] to retrieve a list of foods along with their corresponding nutrient
information, based on standard measurements. This is achieved by leveraging
the API’s Natural Language for Nutrients feature, which allows for parsing of
natural language requests and returning detailed nutrition information for
common foods.

AI and Analysis Models: provide data analysis, enabling the extraction of
insights and associations. This is particularly valuable in intricate
healthcare contexts where LLMs may not efficiently handle intensive
computations. We develop a task focused on daily food risk assessment into the
openCHA framework [12]. This task follows established guidelines for daily
nutrient intake to perform food risk assessments. The task also leverages the
knowledge extracted from the NutritioniX using Knowledge Base developed task.
This algorithm evaluates potential risks associated with food intake by
comparing extracted nutritional data to recommended thresholds for diabetics,
accounting for factors like carbohydrate content, sugar levels, protein
amount, and dietary fiber.

Figure 3: A sample question and responses from the proposed CHA and GPT4. The
green text indicates matching the identified risk with the ground truth, the
red text indicates mismatching.

The threshold levels for the risk assessment algorithm are selected based on
commonly recognized and accepted standards by medical associations for
nutrient intake. According to an American Diabetes Association report [18],
individuals suffering from Type 2 Diabetes should eat 20 to 35 grams (g) of
fiber from raw vegetables and unprocessed grains per day and limit sodium
intake to 2,300 milligrams (mg) per day, both recommendations aligning with
those for the general population. Although there are no overarching
recommended ranges for carbohydrate intake, based on [19], patients with Type
1 and type 2 diabetes in the US consumed about 45% of their total energy
intake through carbohydrates. While there is no concrete evidence about a
certain range of protein intake aiding in glycemic control, we set protein
intake to be 15-20% of total calorie intake based on ADA standards of care
[18]. The recommended intake for fat is 20-35% of total calories according to
[20, 21] while saturated fat is best restricted to less than 10% of total
calories [22]. Based on the American Heart Association (AHA) guidelines [23],
we set stringent restrictions on sugar intake due to its negative impact on
diabetes with a maximum of 6 teaspoons or 25g of sugar.

##  III Results

TABLE I: Sample evaluation for one question.

| Carbohydrate | Fat | Saturated Fat | Protein | Sodium | Sugars | Dietary Fiber
---|---|---|---|---|---|---|---
Grount Truth | R | R | NR | R | NR | R | NR
Proposed CHA | R | R | NR | R | NR | R | NR
GPT4 | R | NR | NR | NR | NR | R | NR
R=Risky, NR=Not Risky

In this section, we assess the effectiveness of our proposed CHA in evaluating
food-related risks for individuals with diabetes. Our objective is to
determine the CHA’s ability to gauge the risk associated with daily food
intake when the nutritional content deviates from recommended ranges. We
compare the performance of our CHA, equipped with comprehensive food knowledge
and analytical capabilities, with the OpenAI’s GPT4 model [15].

We collected 100 sample questions related to individuals inquiring about their
daily meal intake. For the ground truth, we manually extract the nutritional
data for each question and compute the overall nutritional values for
parameters, such as fat, saturated fat, protein, carbohydrates, sodium,
sugars, and dietary fiber. Subsequently, we evaluate the alignment of the
nutrients with established guidelines. If any of these nutrients fall outside
the recommended range, we categorize them as ”Risky,” while those within the
range are labeled as ”Not Risky.” Ultimately, we generate a table of 100 rows
and 7 columns (as our ground truth), representing Risky or Not Risky for each
nutrient.

These questions are posed to both our CHA and GPT4 [15] in order to compare
their responses with the ground truth. Figure 3 provides an illustrative
instance of a question and responses from the proposed CHA and GPT4. The green
text highlights the model’s accurate identification of risk for specific
nutrients based on the ground truth, while the red text indicates instances
where the chatbot incorrectly assesses the risk for the corresponding
nutrient. Table I shows the ground truth values and summarizes the assessments
of the two chatbots for each nutrient.

Similarly, we replicate this procedure for each set of 100 questions and then
compare the quantity of correctly identified risks against the ground truth.
For GPT4, we formulate a customized prompt to instruct GPT4 to estimate the
nutritional information of the mentioned foods and assess their alignment with
recommended guidelines and rules.

TABLE II: The risk assessment accuracy of the 100 sample questions.

| Carbohydrate | Fat | Saturated Fat | Protein | Sodium | Sugars | Dietary Fiber
---|---|---|---|---|---|---|---
GPT4 | 49% | 52% | 68% | 17% | 73% | 58% | 46%
Proposed CHA | 84% | 94% | 99% | 90% | 95% | 91% | 92%

Table II indicates the responses to the 100 questions of both GPT4 and the
proposed CHA. Our CHA outperforms GPT4 across all seven nutrients. This shows
the significance of incorporating guidelines, knowledge bases (e.g.,
NutritioniX), and analytical tools into LLMs for health management tasks. It
is worth noting that, in certain instances, the nutrient calculations for the
proposed CHA bordered closely to the established guidelines, resulting in an
increased error rate in some cases. For instance, when manually calculated,
Carbohydrate content was found to be 44%, which falls within the recommended
range of less than 45%, whereas our CHA calculated it to be 45%. This might be
due to minor variations in portion sizes or rounding of fractional numbers.

##  IV Discussion

The proposed CHA offers a significant degree of flexibility for the
integration of LLMs with external health data sources, knowledge bases, and
analytical tools. This flexibility provides opportunities for mitigating
hallucination problems while enhancing the personalization and reliability of
the CHA. For instance, it allows for the incorporation of diverse elements
such as food-related knowledge graphs, personal health biomarkers, individual
demographics, and food preferences into the existing CHA framework, thereby
enabling more refined and tailored recommendations.

Additionally, this development emphasizes a strong focus on explainability. It
enables users to inquire about the sequence of tasks and actions involved in
generating a response. For example, when users seek information about specific
food nutrients, data sources, or risk calculation methods, our CHA can provide
insights into the underlying tasks and display results (in contrast to the
state-of-the-art chatbots). This heightened level of transparency enhances the
trustworthiness of the CHAs, fostering user confidence in the responses.

##  V Conclusion

In this paper, we proposed an LLM-based CHA for diabetes management enabled by
knowledge-infused LLMs. To achieve this, We leveraged openCHA for our
development. We developed a nutrition information retrieval task to integrate
the Nutritionix knowledge base into the CHA. Moreover, we developed a food
risk assessment tool based on the American Diabetes Association dietary
guidelines. We evaluated the effectiveness of our proposed CHA, in comparison
to GPT4, in food-related risks for individuals with diabetes. Our findings
showed that the proposed CHA outperforms the general-purpose GPT4 in
responding to 100 diabetes-related queries regarding daily meal choices. This
advancement highlights the potential of LLM-powered conversational agents in
improving the accessibility and effectiveness of diabetes management,
addressing a crucial aspect of healthcare for diabetic patients.

## References

  * [1] J. U. Poulsen et al., “A diabetes management system empowering patients to reach optimised glucose control: from monitor to advisor,” in 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, pp. 5270–5271, IEEE, 2010.
  * [2] L. Mertz, “Automated insulin delivery: taking the guesswork out of diabetes management,” IEEE pulse, vol. 9, no. 1, pp. 8–9, 2018.
  * [3] H. A. Klein et al., “Self management of medication and diabetes: Cognitive control,” IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, vol. 34, no. 6, pp. 718–725, 2004.
  * [4] G. Fico et al., “Integration of personalized healthcare pathways in an ict platform for diabetes managements: a small-scale exploratory study,” IEEE Journal of Biomedical and Health Informatics, vol. 20, no. 1, pp. 29–38, 2014.
  * [5] I. Shalahuddin et al., “Blood sugar levels regulation in diabetes mellitus type 2 patients through diet management,” Jurnal Aisyah : Jurnal Ilmu Kesehatan, 2022.
  * [6] W. Russell et al., “Nutritional management of blood glucose levels,” Annals of Nutrition and Metabolism, vol. 63, pp. 1339–1340, 2013.
  * [7] G. G. R. Sng et al., “Potential and pitfalls of chatgpt and natural-language artificial intelligence models for diabetes education,” Diabetes Care, vol. 46, no. 5, pp. e103–e105, 2023.
  * [8] H. Yang et al., “Exploring the potential of large language models in personalized diabetes treatment strategies,” medRxiv, pp. 2023–06, 2023.
  * [9] H. Sun et al., “An ai dietitian for type 2 diabetes mellitus management based on large language and image recognition models: Preclinical concept validation study,” Journal of Medical Internet Research, vol. 25, p. e51300, 2023.
  * [10] X. Liu et al., “P-tuning: Prompt tuning can be comparable to fine-tuning across scales and tasks,” in Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers), pp. 61–68, 2022.
  * [11] E. J. Hu et al., “Lora: Low-rank adaptation of large language models,” arXiv preprint arXiv:2106.09685, 2021.
  * [12] M. Abbasian et al., “Conversational health agents: A personalized llm-powered agent framework,” arXiv preprint arXiv:2310.02374, 2023.
  * [13] Y. Li et al., “Personal llm agents: Insights and survey about the capability, efficiency and security,” arXiv preprint arXiv:2401.05459, 2024.
  * [14] U. Neisser, “Perceiving, anticipating, and imagining,” in Perception and Cognition: Issues in the Foundations of Psychology, University of Minnesota Press, Minneapolis, 1978.
  * [15] OpenAI, “Chatgpt: Openai’s conversational ai model.” https://openai.com/chatgpt, Accessed: January 2024.
  * [16] S. Yao et al., “Tree of thoughts: Deliberate problem solving with large language models,” arXiv preprint arXiv:2305.10601, 2023.
  * [17] Syndigo Company, “Nutritionx: The world’s largest verified nutrition database.” https://www.nutritionix.com/.
  * [18] “5. Lifestyle Management: Standards of Medical Care in Diabetes-2019,” Diabetes care, vol. 42, pp. S46–S60, 1 2019.
  * [19] A. Gray and R. J. Threlkeld, “Nutritional Recommendations for Individuals with Diabetes,” Diabetologia, vol. 54, 10 2019.
  * [20] B. E. Millen et al., “2013 American Heart Association/American College of Cardiology Guideline on Lifestyle Management to Reduce Cardiovascular Risk: practice opportunities for registered dietitian nutritionists,” Journal of the Academy of Nutrition and Dietetics, vol. 114, pp. 1723–1729, 11 2014.
  * [21] J. K. Snell-Bergeon et al., “Adults with type 1 diabetes eat a high fat, atherogenic diet which is associated with coronary artery calcium,” Diabetologia, vol. 52, p. 801, 5 2009.
  * [22] L. Van Horn et al., “The evidence for dietary prevention and treatment of cardiovascular disease,” Journal of the American Dietetic Association, vol. 108, pp. 287–331, 2 2008.
  * [23] R. K. Johnson et al., “American heart association nutrition committee of the council on nutrition, physical activity, and metabolism and the council on epidemiology and prevention. dietary sugars intake and cardiovascular health: a scientific statement from the american heart association,” 2009.