Convergence of explicitly coupled Simulation Tools (Cosimulations)

TL;DR

This paper proves the convergence of explicit cosimulation schemes used in engineering, discusses balance errors due to subsystem coupling, and analyzes the stability issues of balance correction methods.

Contribution

It provides a convergence proof for explicit cosimulation schemes and analyzes their stability, including the effects of balance correction techniques.

Findings

Convergence is proven for explicit cosimulation schemes.

Balance errors can accumulate and affect accuracy.

Balance correction methods can mitigate but not eliminate errors.

Abstract

In engineering, it is a common desire to couple existing simulation tools together into one big system by passing information from subsystems as parameters into the subsystems under influence. As executed at fixed time points, this data exchange gives the global method a strong explicit component. Globally, such an explicit cosimulation schemes exchange time step can be seen as a step of an one-step method which is explicit in some solution components. Exploiting this structure, we give a convergence proof for such schemes. As flows of conserved quantities are passed across subsystem boundaries, it is not ensured that systemwide balances are fulfilled: the system is not solved as one single equation system. These balance errors can accumulate and make simulation results inaccurate. Use of higher-order extrapolation in exchanged data can reduce this problem but cannot solve it. The remaining balance error has been handled in past work by recontributing it to the input signal in next coupling time step, a technique labeled balance correction methods. Convergence for that method is proven. Further, a proof for the lack of stability of such methods is given for cosimulation schemes with and without balance correction.

Figures

Tables

Table 1. Table 1: Cosimulation Schemes for the spring-mass system, left constant, right linear extrapolation

Spring	Mass
System States
$x_1:=s=x$	$x_2:=v=\dot{x}$
Outputs
$y_1:=F=-cx$	$y_2:=v=\dot{x}$
Inputs
$u_1:=y_2$	$u_2:=y_1$
Equations
$\begin{array}{c}\dot{x_1}=\mathrm{Ext}\left(u_1\right)=v\end{array}$	$\begin{array}{c}\dot{x_2}=-\frac{1}{m}\mathrm{Ext}\left(u_2\right)\\ =-\frac{F}{m}\end{array}$

Equations

$$S_1:$$

$${\dot{\mathbf{x}}}_1$$

$$0$$

$$S_2:$$

$${\dot{\mathbf{x}}}_2$$

$$0$$

$$S_1:$$

$${\dot{\mathbf{x}}}_1$$

$$0$$

$$S_2:$$

$${\dot{\mathbf{x}}}_2$$

$$0$$

$$\Delta E_i^j={\int }_{T_{j-1}}^{T_j}{\mathbf{u}}_{i}dt-{\int }_{T_{j-1}}^{T_j}{\overline{\mathbf{u}}}_{i}dt$$

$$S_1:{\dot{\mathbf{x}}}_1$$

$$S_2:{\dot{\mathbf{x}}}_2$$

$${\dot{\mathbf{x}}}_{\mathcal{D}}=f(t,\left[{\mathbf{x}}_{\mathcal{D}},{\mathbf{x}}_{\mathcal{I}}\right]).$$

$$ϵ(t,\mathbf{x},h)={\Phi }^{t+h,t}\mathbf{x}-{\Psi }^{t+h,t}\mathbf{x}$$

$$ϵ(t,\mathbf{x},h)=O(h^{p+1})\text{for}h⟶0.$$

$${\mathbf{x}}_{\Delta }(t_0)$$

$${\mathbf{x}}_{\Delta }(t_i)$$

$${ϵ}_{\Delta }(t)=\mathbf{x}(t)-{\mathbf{x}}_{\Delta }(t)$$

$${ϵ}_{\Delta }(t)=O(h^p)\text{for}h⟶0.$$

$$ϵ(t,\mathbf{x},h)=O(h^{p+1})$$

$${ϵ}_{\Delta }(t)=O(h^p).$$

$$\begin{gathered} {ϵ}_{\Delta 1}(t_1)={\mathbf{x}}_{\mathcal{D}}(t_1)-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_1) \\ ={\Phi }^{t_1,t_0}(\left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0}\right])-{\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\Delta ,\mathcal{D},0},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I},0})\right]) \\ =\underset{I}{\underbrace{{\Phi }^{t_1,t_0}(\left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0}\right])-{\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0})\right])}} \\ +\underset{II}{\underbrace{{\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0}\right])-{\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\Delta ,\mathcal{D},0},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I},0})\right])}} \end{gathered}$$

$$\begin{gathered} {\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0}\right])-{\Psi }^{t_1,t_0}(\left[{\mathbf{x}}_{\Delta ,\mathcal{D},0},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I},0})\right]) \\ ={\mathbf{x}}_{\mathcal{D}}(t_0)+h\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},{\mathbf{x}}_{\mathcal{I}}\right](t_0),h\right)-\left({\mathbf{x}}_{\mathcal{D}}(t_0)+h\psi (t_j,\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_0),h)\right) \end{gathered}$$

$$\begin{gathered} \Vert \cdot \Vert \le hL\Vert \left[{\mathbf{x}}_{\mathcal{D},0},{\mathbf{x}}_{\mathcal{I},0}\right]-\left[{\mathbf{x}}_{\Delta ,\mathcal{D},0},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I},0})\right](t_0)\Vert \\ \le hL\Vert {\mathbf{x}}_{\mathcal{I},0}-\mathrm{Ext}({\mathbf{x}}_{\mathcal{I},0})(t_0)\Vert \le hLO(h^{P+1})\le LO(h^{P+2}). \end{gathered}$$

$$\begin{gathered} {ϵ}_{\Delta }(t_{j+1})={\mathbf{x}}_{\mathcal{D}}(t_{j+1})-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_{j+1}) \\ =\underset{=:\mathbf{ϵ}}{\underbrace{{\mathbf{x}}_{\mathcal{D}}(t_{j+1})-{\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j))}} \\ +\underset{=:{\mathbf{ϵ}}_{\text{Prop}}}{\underbrace{{\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j))-{\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\Delta ,\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\Delta ,\mathcal{I}})\right](t_j))}} \end{gathered}$$

$$\begin{gathered} \mathbf{ϵ}(t_{j+1})={\mathbf{x}}_{\mathcal{D}}(t_{j+1})-{\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j)) \\ ={\mathbf{x}}_{\mathcal{D}}(t_j)+h\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},{\mathbf{x}}_{\mathcal{I}}\right](t_j),h\right)+r(h^{p+1}) \\ -\left({\mathbf{x}}_{\mathcal{D}}(t_j)+h\psi (t_j,\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j),h)\right) \\ =h\left\{\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},{\mathbf{x}}_{\mathcal{I}}\right](t_j),h\right)-\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j),h\right)\right\} \\ +r(h^{p+1}). \end{gathered}$$

$$\begin{array}{r}\Vert ϵ\Vert \le hL\Vert \left[{\mathbf{x}}_{\mathcal{D}},{\mathbf{x}}_{\mathcal{I}}\right](t_j)-\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j)\Vert +\Vert r(h^{p+1})\Vert \\ =hL\Vert ({\mathbf{x}}_{\mathcal{I}}-\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}}))(t_j)\Vert +\Vert r(h^{p+1})\Vert \\ =hLO((jh{)}^{P+1})+\Vert r(h^{p+1})\Vert \\ =j^{P+1}hLO(h^{P+1})+\Vert r(h^{p+1})\Vert \\ \le n^{P+1}O(h^{P+2})+\Vert r(h^{p+1})\Vert \end{array}$$

$$\Vert ϵ\Vert \le O(H^{P+2})$$

$$\Vert ϵ\Vert \le O(H^0)$$

$$\begin{gathered} {\mathbf{ϵ}}_{\text{Prop}}(t_{j+1})= \\ {\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j))-{\Psi }^{t_{j+1},t_j}(\left[{\mathbf{x}}_{\Delta ,\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\Delta ,\mathcal{I}})\right](t_j)) \\ ={\mathbf{x}}_{\mathcal{D}}(t_j)+h\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j),h\right) \\ -\left({\mathbf{x}}_{\Delta ,\mathcal{D}}(t_j)+h\psi (t_j,\left[{\mathbf{x}}_{\Delta ,\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\Delta ,\mathcal{I}})\right](t_j),h)\right) \\ ={\mathbf{x}}_{\mathcal{D}}(t_j)-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_j) \\ +h\left\{\psi \left(t_j,\left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j),h\right)-\psi \left(t_j,\left[{\mathbf{x}}_{\Delta ,\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j),h\right)\right\} \end{gathered}$$

$$\begin{gathered} \Vert {\mathbf{ϵ}}_{\text{Prop}}(t_{j+1})\Vert \le \Vert {\mathbf{x}}_{\mathcal{D}}(t_j)-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_j)\Vert \\ +hL\Vert \left[{\mathbf{x}}_{\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j)-\left[{\mathbf{x}}_{\Delta ,\mathcal{D}},\mathrm{Ext}({\mathbf{x}}_{\mathcal{I}})\right](t_j)\Vert \\ \le \Vert {\mathbf{x}}_{\mathcal{D}}(t_j)-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_j)\Vert +h{L}_{\mathcal{D}}\Vert {\mathbf{x}}_{\mathcal{D}}(t_j)-{\mathbf{x}}_{\Delta ,\mathcal{D}}(t_j)\Vert \\ =(1+h{L}_{\mathcal{D}})\Vert {\mathbf{ϵ}}_{\Delta }(t_j)\Vert , \end{gathered}$$

$$\Vert ϵ\Vert \le O(h^{p+1})$$

Full text

As flows of conserved quantities are passed across subsystem boundaries, it is not ensured that systemwide balances are fulfilled: the system is not solved as one single equation system. These balance errors can accumulate and make simulation results inaccurate. Use of higher-order extrapolation in exchanged data can reduce this problem but cannot solve it. The remaining balance error has been handled in past work by recontributing it to the input signal in next coupling time step, a technique labeled balance correction methods. Convergence for that method is proven.

Further, a proof for the lack of stability of such methods is given for cosimulation schemes with and without balance correction.

Keywords:

Cosimulation, coupled problems, stability, convergence, balance correction, extrapolation of signals,

1 Introduction

Engineers are increasingly relying on numerical simulation techniques. Models and simulation tools for various physical problems have come into existence in the past decades. The desire to simulate a system that consists of well described and treated subsystems by using appropriate solvers for each subsystem and letting them exchange the data that forms the mutual influence is immanent.

The situation usually is described by two coupled differential-algebraic systems $S_{1}$ and $S_{2}$ that together form a system $S$:

[TABLE]

The $(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{1},\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{2})$ are the differential states of $S$, their splitting into $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$ determines the subsystems $S_{i}$ together with the choices of the $\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i}$. In Co-Simulation the immediate mutual influence of subsystems is replaced by exchanging data at fixed time points and subsystems are solved separately and parallely but using the received parameter:

[TABLE]

where $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ are given by coupling conditions that have to be fulfilled at exchange times $T_{k}$

[TABLE]

and are not dependent on subsystem $i$’s states any more, so are mere parameters between exchange time steps.

Full row rank of $d_{\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i}}\mathchoice{\displaystyle\boldsymbol{g}}{\textstyle\boldsymbol{g}}{\scriptstyle\boldsymbol{g}}{\scriptscriptstyle\boldsymbol{g}}_{i}$ can be assumed, such that the differential-algebraic systems are of index 1. This description of the setting is widespread ([3]). With the $\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{ij}$ being solved for $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}$ inside the $S_{j}$ (let solvability be given), for systems with more than two subsystems it is more convenient to write output variable ${\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{j}$ and now redefine $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}$ as the input of $S_{i}$, consisting of some components of the outputs $\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}_{j}$ [2]. This structure is defined as kind of a standard for connecting simulators for cosimulation by the Functional Mockup Interface Standard [1]. It defines clearly what information a subsystems implementation provides.

In Co-Simulation the variables establishing the mutual influence of subsystems are exchanged at fixed time points. This results in continuous variables being approximated by piecewise constant extrapolation, as shown in the following picture:

If one does not want to iterate on those inputs by restarting the simulations using the newly calculated inputs, one just proceeds to the next timestep.

This gives the calculations an explicite component, the mutual influence is now not immediate any more, inducing the typical stability problems, besides the approximation errors.

But for good reasons, explicit co-simulation is a widely used method: It allows to put separate submodels, for each of which a solver exists, together into one system and simulate that system by simulating each subsystem with its specialised solver - examination of mutual influence becomes possible without rewriting everything into one system, and simulation speed benefits from the parallel calculation of the submodels. Usually it is highly desirable that a simulation scheme does not require repeating of exchange time intervals or iteration, as for many comercial simulation tools this would already require too deep intrusion into the subsystems method and too much programming in the coupling algorithm.

The following fields of work on explicite co-simulation can be named to be the ones of most interest:

Improvement of the approximation of the exchanged data will most often improve simulation results [4]. This is usually done by higher-order extrapolation of exchanged data, as shown in this plot, where the function plotted with dots is linearly extrapolated:

When the mutual influence between subsystems consists of flow of conserved quantities like mass or energy, it turns out that the improvement of the approximation of this influence by extrapolation of past data is not sufficient to establish the conservation of those quantities with the necessary accuracy. The error that arises from the error in exchange adds up over time and becomes obvious (and lethal to simulation results many times). In a cooling cycle example ([7, Section 6.3]), a gain of 1.25% in coolant mass occurs when simulating a common situation.

It has been tried to meet this challenge by passing the amount of exchanged quantity for the past timestep along with the actual flow on to the receiving system, where then the error that has just been commited is calculated and added to the current flow to compensate the past error. For well damped example problems in fluid circles this method has fulfilled the expectations [10]. It has been labelled balance correction. 3. 3.

There is good reason to prevent jumps in exchanged data by smoothing. Higher order extrapolation polynomials cannot make extrapolated data at the end of the exchange timestep match the newly given value.

There have been some examinations of the convergence of iterative simulator coupling schemes, which can be applied to noniterative coupling as the case where the number of iterations is 0. But to the best of our knowledge, no proof for the convergence of the explicite simulator coupling scheme has been presented so far, despite its importance for industrial applications. Neither has there been an examination of convergence of balance correction methods. This is a considerable gap because the balance correction method is a highly heuristically motivated method with a inherent danger: It means making an error in the exchanged signal $u$ for lowering the accumulated error in the amount of that quantity – it increases the $L1$-error to lower the $L2$-error. This should be considered well.

2 Some preliminaries

Always, a time interval or quantities belonging to one are indexed with the index of its right boundary: $\Delta t_{i}:=[t_{i-1},t_{i})$. Indexing of times begins with 0. Big letters refer to time exchange steps, e.g $\Delta T_{k}:=[T_{k-1},T_{k})$ is an interval between exchange times, whereas above interval might denote a subsystems step. Consistently, $k=1,...,N$, and above $i$ ranges from $1$ to $n$ too. Let $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})_{i}^{k}(t)$ denote the extrapolation of the input variable $i$ in interval $k$, sometimes referred to as hold function.

2.1 Choice of exchanged Variables

The direction in which the variables are exchanged cannot be chosen arbitrarily. Careful consideration is needed to not place contradictorily constraints on the models. If i.e. at one of the model boundaries $S_{1}$ exchanges the value of a differential state, it must only export and not receive this value – otherwise this state is turned into a parameter and extrapolated instead of solved, or one would overwrite an integration result and reinitialize the solver at each exchange. The choice is to receive arguments of the derivative of the state and send the state value itself [10], [7]. Often, the physics of coupling is described by a flow driven by a potential, where the flux depends on the potential by $f\sim\nabla_{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\mathcal{P}$, or in non-spatial context $f=\operatorname{const}\left[\mathcal{P}_{S_{2}}-\mathcal{P}_{S_{1}}\right]$. Then exchange is commonly implemented in the following way: System $S_{1}$ passes the value of the flow-determining potential $\mathcal{P}$ on to $S_{2}$, whereas the latter calculates the flux $f$ and passes it to $S_{1}$.

2.2 Detailed description of Techniques in coupling Simulation Software

After the brief overview in the introduction on the issues that appear when coupling simulation software and means to treat them was given. Those fields are described more precisely here, giving citations.

2.2.1 Determination of consistent Inputs

As $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ fulfils $\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i},\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i})$, where $\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i}$ results from $\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}=\mathchoice{\displaystyle\boldsymbol{g}}{\textstyle\boldsymbol{g}}{\scriptstyle\boldsymbol{g}}{\scriptscriptstyle\boldsymbol{g}}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i},\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{j})$, it in general depends on $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{j}$. So the coupling equations (9) (10) $({\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{i})_{i=1..n}=\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}$ are a coupled system in fact. The $\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{i}$ can be solved with respect to the $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ one after the other only if the directed graph of influence of outputs on each other contains neither bidirectional dependencies nor loops. Otherwise, the coupling conditions remain not fully fulfilled, or one solves them iteratively. Such methods are the Interface-Jacobian based methods, e.g. given in [8], which determine the ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}$ with respect to ${\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{i}=0$ and ${\mathchoice{\displaystyle\boldsymbol{g}}{\textstyle\boldsymbol{g}}{\scriptstyle\boldsymbol{g}}{\scriptscriptstyle\boldsymbol{g}}}_{i}=0$ with the help of a newton solver at exchange times, which means full consistency at those times. But not all tools simulating the subsystems provide the residual of $\mathchoice{\displaystyle\boldsymbol{g}}{\textstyle\boldsymbol{g}}{\scriptstyle\boldsymbol{g}}{\scriptscriptstyle\boldsymbol{g}}_{i}=0$.

2.2.2 Iterative coupled Methods

If restarting the timestep is an option, iterations can be performed and instead of using $\operatorname{Ext}(\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i})$, the time-dependant numerical solution for ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}$ from the old iteration can be used. Schemes are distinguished by the exchange directions of input, as the waveform iteration or the Gauss-Seidel scheme, which can exploit the directions of influences between the subsystems ([9], [3], [4]). Stability and accuracy, also in terms of balance, such can be augmented, of course at computational cost due to the restarts.

In [6], a Newton method is used to solve the coupling equations (9)-(10), but different from [8] not only $({\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}(T_{k}),\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}))_{i=1..n}=\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}$ at the end of the timestep, but $({\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}(\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}),\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}))_{i=1..n}=\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}$. This requires repeated solving of the subsystems alone for Jacobian evaluation.

2.2.3 Increasing the extrapolation order of inputs

To improve approximation $\operatorname{Ext}{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}^{i}$ has been calculated as an extrapolation polynomial of past ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}^{j}$, $j<i$. In his dissertation [4] Busch examines Lagrangian and Hermite extrapolation polynomials for exchanged quantities using a system concatenated from two coupled spring-mass oscillators, thus a linear second order ODE with four degrees of freedom, as model problem.

2.2.4 Investigations of convergence

Busch examines the effect of the approximation order in exchanged data on global convergence and judges the stability of the method for his problem, assuming exact integration of each subsystem. Convergence to exact solution in this case is limited to $O(h^{2})$ for piecewise constant extrapolation, while it is limited to $O(h^{3})$ for degree one Hermite polynomial extrapolation.

As Busch examines a system of ordinary differential equations without algebraic parts and given explicitely, he feels no need to consider smoothing. Balance errors that occur in his problem remain undiscussed. Besides the aforementioned examination [4] which is limited to linear problems, in [3], similarly in [9], there is a more general examination, treating iterative coupling schemes. Mainly, a fixed point argument is used here.

2.2.5 Smooth switching between old and new input

The smoothing of exchanged data is originally motivated by the need to provide a starting value close to the solution during searching for solutions of the algebraic part of the differential-algebraic equation system.

Another reason for smoothing the inputs can be that it enables the calculation of derivatives resp. difference quotients from them. Although this is unfortunate, such needs sometimes occur in practice.

In [7], two halves of a cooling cycle are co-simulated, the mass flow $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{m}$ being exchanged, so a component of $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}$. The smooth input $\overline{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}^{j}$ is concatenated as a convex sum of $\operatorname{Ext}{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}^{j-1}$ and $\operatorname{Ext}\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}^{j}$ weighted with a sufficiently smooth, here degree 5 polynomial, function switching from 0 at $T_{j-1}$ to 1 at $T_{j}$. Besides from stabilizing the Newton solver as inteded, Kossels work gives an impression on balance errors: During some simulations balance error of about 1.25% in coolant mass is observed.

It is obvious that by this smoothed switching the balance error $\Delta E_{i}^{j}$ possibly increases compared to a unsmoothed switching procedure.

Another work that considers smoothing is the paper [12].

2.2.6 Balance Correction

Classical balance errors of extrapolations is $\Delta E_{i}^{k}=\int_{T_{k-1}}^{T_{k}}\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}-\overline{({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})}_{i}^{k}(t)dt$ is where $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ is a flux of a conserved quantity, but it is defined for arbitrary quantities. Negative and positive contributions from positive and negative intervals partly compensate each other, and as in a typical simulation the system often ends up in a stationary state similar to the one at begin, thus graphs of quantities and their derivatives tend to end at values where they started, the conserved quantities balance error at the end of a simulation may be small. This does not imply it is small during all intervals of the simulation. As a typical such simulation situation think of an automotive driving cycle ([7]): Using piecewise constant extrapolation of $f$, there is a gap in mass after the phase of rising system velocity has passed - this loss remains uncompensated during the (significant) phase of elevated speed.

Scharff, motivated by the loss in balance stated in [7] as a collateral result, proposed in [10] that the errors

[TABLE]

are added in the time step $j+1$. The correction that is applied in the $j$-th interval is then $\Delta E_{i}^{j-1}\phi_{j}(t)$, where the function $\phi_{j}$ is scaled such that its integral is 1 and may be constant but when one wants to preserve smoothness, it is a smooth function smoothly vanishing at the boundaries of the time interval $j$. In spite of the lack of strictness and the increase of errors in derivatives of the exchanged quantities, this method enables coupling of simulations that would be impossible without - Kossels example was recalculated successfully using balance correction. Balance correction methods were applied to nonconserved quantities in [11], as even such quantities have some conservation properties in space and time, and examination will go on here. The method even improved the convergence for a linear such problem.

2.2.7 Reducing exchange induced derivations in input signals

Split systems may have parts that are prone to be excited by quick changes, implying high derivatives, in input signals. An example for such a system is a high-frequency spring-mass system sitting on a slowly moving ground whose displacement is given to the spring as input. In such settings it is essential to avoid the contribution of additional frequencies to the signal that the especially the balance corrected coupling method does. In [11], it was suggested to spread recontribution of balance correction contributions over more than one interval.

2.3 Aim of this work

Most part of our work can be understood as a contribution to Arnold and Guenthers theoretical examination and convergence results. [3], but putting their result about convergence [3, Th. 2.3] into the context of the standard explicit cosimulation scheme from Section 3.1 that is paying respect to practical issues: To the need to run subsystems simultaneously by applying a Jacobi instead of the Gauss-Seidel scheme there, and to the need not to be intrusive into subsystems methods by avoiding iteration. Our convergence results will turn out to be consistent to those presented in [3], but to require simpler derivation as no iterations occur and offer insight into subsystem methods errors effect.

3 Convergence and stability for One-step methods on subsystems

Cosimulation schemes without iteration (repeating timesteps) can be regarded as mixed explicit-implicit methods.

We thus introduce some notations and facts from ODE solving first, mostly the way it is done in [5], before we consider cosimulation methods.

3.1 Scope and notation

The ODE approach that is applied here covers the setting introduced by equations (5) - (8): As $\partial_{{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}}{\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}$ has full rank, equations (9) and (10) can be solved by ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i})={\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{i}^{-1}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i})$. Further, by applying the state space method, which is solving the algebraic equations (6) and (8) at each evaluation of $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}$, the problem is shifted to the solving of the ODE

[TABLE]

The cosimulation scheme for two ODE subsystems whose coupling equations (9) -(10) shall be given explicitely by ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{i}^{-1}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i})$ instead of $\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{i}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i})$, reads:

[TABLE]

By shifting the evaluation of ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{i}^{-1}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i})$ to the receiving subsystem $S_{j}$, ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{i}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{i}$ can be written in the $S_{1}-S_{2}$ system (this straightforwardly generalizes to more subsystems), and a notation is achieved that is better suitable for applying techniques from ODE. Moreover, for readability the subsystems differential states set of indices shall now be $\mathcal{D}$, and that of the input variables, which shall be differential states of other subsystems, shall be $\mathcal{I}$. Thus e.g. ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}}:={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2}$ in first subsystem, and ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}=\left[{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}_{\mathcal{D}},{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}}\right]$ and the ODE governing the subsystem can be written as

[TABLE]

In cosimulation methods, input variables ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}}$ are extrapolated from past data into the current timestep. This extrapolation is denoted $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}})$.

3.2 Consistency

Let $\left\{t_{j}\right\}$, $j=1,...n$ be the time steps taken by the subsystems method in one step. It is $h_{j}=t_{j+1}-t_{j}$.

As in [5], the evolution $\Phi^{t+\tau,t}x$ of some ODE through $x$ at $t$ until $t+\tau$ is the solution of the IVP $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=f(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})$, ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}(t)={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$. Similarly, the discrete evolution $\Psi^{t+\tau,t}x$ is the numerical solution of that IVP produced by a given method, $\Psi^{t+h,t}$ thus the evolution on one time step.

The following common definitions are used:

Definition 3.1** (Consistency).**

[TABLE]

is called the consistency error of a one-step method. A discrete evolution $\Psi^{t+h,t}x$ has consistency order $p$ if

[TABLE]

The following is given in brief, indeed it is a consequence of the two consistency criteria $\Psi^{t,t}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}=x$ and $\frac{d}{d\tau}\Psi^{t+\tau,t}x\mid_{\tau=0}=f(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})$.

Lemma 3.2**.**

Let $\Psi^{t+h,t}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$ be differentiable with respect to $h$. Then the following is equivalent:

• •

The discrete evolution $\Psi^{t+h,t}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$ is consistent.

• •

The discrete evolution can be written as $\Psi^{t+h,t}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}+h\psi(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}},h)$ where $\psi$ is called the increment function and is continuous w.r.t $h$.

• •

$\epsilon(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}},h)=o(h)$* for $h\longrightarrow 0$, i.e. the consistency error vanishes near $t$.*

Proof.

See [5, lemma 4.4] ∎

3.3 Convergence

We call the set of time points between which one step of the method is executed the grid: $\Delta:=\{t_{i}\}$. Each one-step method defines a grid function or numerical solution ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta}$ by solving the IVP $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=f(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})$, ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}(t_{i})={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta}(t_{i})$ on $t_{i},t_{i+1}$, in other words, by the recursion

[TABLE]

Definition 3.3** (Convergence).**

The convergence error of a discrete solution ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta}$ is defined as

[TABLE]

A discrete solution ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta}$ has convergence order $p$ if

[TABLE]

The following result for ODEs is well-known (see, e.g., [5], theorem 4.10:)

Theorem 3.4**.**

Let the consistency error of an evolution $\Psi^{t+h,t}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}+h\psi(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}},h)$ given by an one-step method with $\psi$ Lipschitz-continuous in ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$ satisfy

[TABLE]

Then for all grids $\Delta$ with sufficiently small width $h$ the evolution defines a numerical solution ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta}$ for the initial value ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{0}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}(t_{0})$. The numerical solution converges of order $p$ towards ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$, equivalently,

[TABLE]

In short, the theorem says that if a method has consistency order $p$, it follows that it is also convergent of order $p$.

3.4 Convergence of Cosimulation schemes

Now times of data exchange are denounced as $\left\{T_{i}\right\}$, $i=0,...m$, the exchange step width as $H_{i+1}=T_{i+1}-T_{i}$, whereas the $\left\{t_{j}\right\}$, $j=0,...n$ are the time steps taken by the subsystems method, counting beginning again with 0 in each timestep (indexing by macro and local timestep is not needed). It is $h_{j+1}=t_{j+1}-t_{i}$. For convienience, it is assumed that $H_{i}=\sum_{j=1}^{n}h_{j}$.

As in [5], the evolution $\Phi^{t+\tau,t}x$ of some ODE through $x$ at $t$ until $t+\tau$ is the solution of the IVP $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=f(t,{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})$, ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}(t_{0})={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$.

Proposition 3.5**.**

There exists a consistent and thus convergent cosimulation scheme, given that derivatives for all subsystems are lipschitz-continuous for all subsystems.

This holds for the following: If Euler-forward scheme with one step between exchange times $T_{i}$ is executed, this is equivalent to executing Euler-forward on the whole system $S$ and thus, for Euler-forward being consistent of order 1, this cosimulation method is convergent of order 1 due to theorem 3.4.

Motivated by this simple result, we try to transfer the methods used to prove Theorem 3.4 to realistic subsystem methods as convergence should improve with use of higher-order subsystem methods and higher-order extrapolation of exchanged data. The proof also sketches which steps are to be taken:

Consistency on the subsystems for step $[T_{0},T_{1})$ is nearly trivially given only for $[t_{0},t_{1})\subset[T_{0},T_{1})$ by the consistency of the subsystems methods. For the other subintervals $[t_{i-1},t_{i})\subset[T_{0},T_{1})$ a consistency-like estimate has to be proven – due to the contributions from extrapolated exchanged data the convergence on $[T_{0},T_{1})$ can not be concluded from consistency of subsystem methods using theorem 3.4 in an obvious way. 2. 2.

If an estimate on $[t_{i-1},t_{i})$ is given, convergence of subsystems methods on $[T_{0},T_{1})$ can be proven in a similar way as theorem 3.4. This then means consistency of a step of the cosimulation method on the whole system $S$ is provided. 3. 3.

If by (2) consistency of a step $[T_{0},T_{1})$ of the cosimulation method is provided, again by applying theorem 3.4 to cosimulation steps global convergence of the cosimulation follows.

We have to show items (1) and (2).

3.4.1 Consistency of Subsystems steps inside $\left[T_{0},T_{1}\right)$

For readability, we calculate for the first timestep

[TABLE]

The difference $I$ is of the order consistency of the subsystems method, thus $O(h^{p+1})$. The difference $II$ is

[TABLE]

and is estimated using the Lipschitz continuity of $\psi$ and ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{D},0}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{D},0}$:

[TABLE]

Remind that capital letters denote exchange parameters: $P$ is the degree of the extrapolation of exchanged data, $H$ the exchange time stepwith.

Now consider a subsystem step $j$ inside $[T_{0},T_{1})$. The grid error is

[TABLE]

Expression ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}$ resembles the usual local cutoff error made by $\Psi$ in timestep $j$, but has additional contributions from extrapolations, whereas ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{\text{Prop}}$ describes the propagation of previous errors. It is

[TABLE]

Again using the Lipschitz continuity of $\psi$, we get

[TABLE]

where $n$ is the number of subsystem steps per exchange step. Using exchange stepsize $H=nh$ and provided that the subsystems step size $h$ is bounded by $H$, $h<H$, and that $p>P$,

[TABLE]

formally holds. Unfortunately, in situations of interest $n$ is bigger than one, thus $h$ in general remains unchanged when $H$ is reduced, so strictly seen

[TABLE]

for $H=nh$, $n>1$ in typical application. Fortunately, this still states that $\epsilon$ approaches the error of the subsystems methods, and further, this most common situation will turn out to be covered.

Consistency as required for theorem 3.4 is $\Phi^{t,t+h}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}-\Psi^{t,t+h}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}=O(h^{p+1})$, so there is no difference between exact, numerical and extrapolated solution at $t$, so $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{I}})=\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}})$. Using this, expression ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{Prop}$ of equation (24) is calculated as

[TABLE]

Thus, using the Lipschitz continuity of $\psi$ again

[TABLE]

writing ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{D}}(t_{j})-{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{D}}(t_{j})=\epsilon_{\Delta}(t_{j})$ again.

3.4.2 Consistency of a cosimulation step

First, the situation that exchange and subsystem stepsize are reduced together, $H_{i}=ch_{i}$, is treated, a situation that is occuring when $H_{i}=nh_{i}$ is established, e.g. if $h$ is already reduced by $H>h$ when $H\longrightarrow 0$ or by reducing subsystems solvers tolerances together with $h$.

This situation again splits into the situation when $p=P$ and when $p>P$. In the first case,

[TABLE]

holds, in the second

[TABLE]

Equations (30) and, depending on situation, (28), (31) or (32) now give different recurrence schemes for the propagated error.

First, let $H=ch$ be established and $p=P$. In this case (31) yields the scheme

[TABLE]

where $\tilde{C}$ is a constant.

The solution of this recurrence scheme is

[TABLE]

which is shown in the proof of 3.4, see [5, Th. 4.10]. The reader is also referred to the solving of the next recurrence scheme, where the arguments are repeated. As

[TABLE]

and $H=ch$

[TABLE]

is proven for $p\leq P$.

Now, let $H=ch$ be established and $p>P$. The recurrence scheme, using (32), now reads

[TABLE]

It is shown that this error is of order $H^{P+1}$ by transferring the arguments from the proof of 3.4 as given in [5, Th. 4.10] to this situation. Inductively,

[TABLE]

where $C$ is a constant and $h_{\Delta}$ the maximum stepwidth in the interval of interest is shown:

First, (46) holds for $t_{0}$. Now assume it holds for $t_{j}$. Then

[TABLE]

where $(1+h_{j}L)\leq e^{h_{j}L}$ has been applied. Again using $\left(e^{L(t_{n}-t_{0})}-1\right)=(LH+(LH)^{2}+O(H^{3}))$,

[TABLE]

is proven for $H=ch$ and $p>P$.

The preceding two results provide a solution for the recurrence scheme without restriction by relating $h$ and $H$ to each other, using full equation (26):

[TABLE]

Now

[TABLE]

is proven inductively. Assume it holds for $t_{j}$. Then

[TABLE]

Expression (49) equals (33) (ii) and so is estimated using (34), and expression (62) equals (37) and has been estimated using (46). This finally is, as before,

[TABLE]

In the case where the refinement of the exchange time grid does not (yet) influence $h$, this means, after the last refinement $H$ still is a multiple of $h$, the preceding result applies with $P=0$, giving

[TABLE]

3.4.3 Global convergence

Consistency of a cosimulation step for all subsystems in the cosimulation scheme is given by (53), which is in fact estimating the error of the $i$th subsystems differential states, by reintroducing the subsystems index that could be ommitted in the previous section ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{\Delta,i}(t_{j+1})={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{D},i}(t_{j}+1)-{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{D},i}(t_{j}+1)$,

[TABLE]

So for the consistency error $\left\|{\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{S}(T)\right\|$ of all states of the system,

[TABLE]

so the cosimulation method applied to overall system $S$ is consistent according to definition (3.1) with respect to methods stepwith $H$.

Using this, one proceeds in the straightforward way: The recurrence inequality now reads

[TABLE]

where the index $S$ indicates that the error of whole system $S$ is considered. Calculation is performed in the appendix and is completely analogous to equations (24)-(30).

The recurrence scheme is solved again analogously to the proof of (46):

[TABLE]

is proven inductively. Assume it holds for $t_{j}$. Then

[TABLE]

By this result, refining the subsystems step size alone does not ensure convergence, and refinement of exchange step sizes alone while it does not bound $h$ has an offset in error:

[TABLE]

Thus unique rates might be less clearly observable as usual in practice due to the presence of contributions of different order and of the two variables $H$ and $h$. See section 4.2 for examples.

The above clear cases are consistent with the conclusions one can draw from theorem 3.4 using consistency orders

[TABLE]

Moreover, this is consistent with estimates given in [2, Example 3]. Finally, the result is given as a

Theorem 3.6**.**

Let S be a set of ODE which is split into disjoint subsystems $S_{k}$ of the shape

[TABLE]

the $\mathcal{D}_{k}$ and $\mathcal{I}_{k}$ denoting index sets. Let ${T_{i}}$ be a time grid with width $H$, and let the inputs of all subsystems ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}_{k}}$ be extrapolated at $T_{i}$ with polynomial order $P$ and then be solved with an one-step method of order $p$ and maximal stepwidth $h\leq H$. Then for the error of the numerical solution the estimate

[TABLE]

holds.

Remark

For extrapolation of input data to a certain order, exact values and derivatives of it are required. Those values and derivatives may depend on up-to-date input. If dependencies on input are circular, it thus may be impossible to get exact derivatives without some iterating procedure for solving the coupling equations (9) and (10), as e.g. presented in [8]. When estimating the consistency error (eq. (24)), thus order loss may occur.

3.5 Stability

Aim of this section is to examine the stability of the overall system.

According to the established notions of A- and B-stability in standard ODE theory, a method is called A- (resp. B-) stable if it conserves the stability properties of the ODE.

First, A-stability is a property of performance of methods on the scalar problem $\dot{x}=\lambda x$ and as such of no use for examination of coupling schemes. We instead examine linear stability, which shall be the conservation of linear vector valued problems

[TABLE]

stability by the method in question.

Consider a system $S$ governed by $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=\tensor[]{B}{}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$, which is split into subsystems $S_{i}$ with variables $x_{i}$, $i=1,2$ governed by

[TABLE]

The last summands in this setting are inputs and in previously used notation one would write $({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{D}},{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}})=(x_{1},x_{2})$. Cosimulation method means

[TABLE]

The method is not stable if it does not preserve the stability of the ODE, i.e. if there is a matrix $\tensor[]{B}{}$ with $\rho(\tensor[]{B}{})<1$ but

[TABLE]

Choose $B_{11}=0$ and $B_{12}\neq 0$ but subject to $\rho(\tensor[]{B}{})<1$. Then the solution of the method induced ODE (73) at $t_{k}\in[T_{j},T_{j+1}]$ is

[TABLE]

For order $n$ extrapolation the integral is a polynomial of degree $n+1$, for constant extrapolation that is $x_{2}t$ - the solution of the method induced ODE such is not stable for any extrapolation, so the numerical solution cannot be stable for any subsystems integration scheme.

We state that from B-stability linear stability follows: Consider above linear ODE, which shall be stable, and thus $\operatorname{Re}(\rho(\tensor[]{B}{}))\leq 0$. From these properties $\left<\tensor[]{B}{}x,x\right>\leq 0$, which is $\left<\tensor[]{B}{}(y_{1}-y_{2}),(y_{1}-y_{2})\right>\leq 0$ and so the dissipativity. So the discrete solution of a B-stable method applied to a stable linear ODE will be stable. As linear stability is not given for explicite cosimulation schemes, explicite cosimulation is neither B-stable.

These results could be expected, if one sees cosimulation methods as mixed explicit-implicit methods.

4 Convergence of balance correction methods

Balance correction methods seek to establish conservation of an input quantity ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}}$ by adding its extrapolation error to the input at a later timestep, when it can be calculated. Denoting

[TABLE]

where

[TABLE]

and $g(t,t_{i},h)$ is a smooth function with a support beginning at $t_{i}$, extending to 1-4 data exchange intervals and fulfilling $\int g(t,t_{i},h)\,dt=1$, the general one-step balance correction method reads

[TABLE]

Our examination now follows exactly the arguments outlined in the introduction of 3.4 and treats the extrapolation error contribution in the method as an input error and shows that this input error does not spoil validity of any of the arguments used.

4.1 Consistency of subsystem step inside $\left[T_{j},T_{j+1}\right)$

One may split $\epsilon_{\Delta}(t_{j})$ for a balance correction method according to (24) and further, such that

[TABLE]

• •

The first difference $\epsilon$ is the same as in equation (24), thus named like it and thus yields

[TABLE]

which then is estimated analogously to (26) as

[TABLE]

• •

Expression ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{\text{Prop}}$ is the propagated error again and is treated as in equation (29) using the Lipschitz continuity. The fact that consistency as required for theorem 3.4 is $\Phi^{t,t+h}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}-\Psi^{t,t+h}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}=O(h^{p+1})$ justifies the use of $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}})=\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{I}})$ here again. This results in equation (30) again:

[TABLE]

• •

Finally, ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{\text{bal}}$ becomes

[TABLE]

Mind that $\operatorname{bal}_{j}$ is part of the method and has to be considered during consistency analysis even if otherwise $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}})=\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\Delta,\mathcal{I}})$ is assumed on $[T_{0},T_{1}]$. By the Cauchy-Schwarz inequality

[TABLE]

Together,

[TABLE]

If $H$ while $h$ remains constant, this is $H^{P+1}$, if $h$ is reduced while $n$ constant, this is $h^{P+3}$, so higher order compared to ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}$ in any case. The recurrence scheme can be set up the same way as in section 3.4,

[TABLE]

with the same result (46)

[TABLE]

and same implications (53) for consistency and (66) for global convergence:

[TABLE]

Finally, the result is formulated as

Theorem 4.1**.**

Let $S$, $S_{k}$, ${{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}_{\mathcal{D}_{k}}$, ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}_{k}}$ and ${T_{i}}$ be as in theorem 3.6, but balance correction contributions be added to the extrapolated variables. Then the estimate from theorem 3.6 still holds.

This results can be interpreted in the following way: If we regard the balance correction as an errorous disturbation, it still is too small to disturb convergence.

It should be possible to show that for a sufficiently small macro time step size, the error of a cosimulation method with balance correction is smaller than one without balance correction.

4.2 Numerical results

4.2.1 Linear Problem for convergence examination

For examination of convergence, a two-component linear system split into two subsystems of one variable each is solved, which written in FMI style reads:

General scheme

$S_{1}$

$S_{2}$

System States

${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1}$

${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2}$

Outputs

${\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{1}$

${\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{2}$

Inputs

${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{1}={\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{2}$

${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{2}={\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{1}$

Equations

$\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1}}={\mathchoice{\displaystyle\boldsymbol{f}}{\textstyle\boldsymbol{f}}{\scriptstyle\boldsymbol{f}}{\scriptscriptstyle\boldsymbol{f}}}_{1}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1},\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{1}))$

$\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2}}={\mathchoice{\displaystyle\boldsymbol{f}}{\textstyle\boldsymbol{f}}{\scriptstyle\boldsymbol{f}}{\scriptscriptstyle\boldsymbol{f}}}_{2}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2},\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{2}))$

$\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=\tensor[]{A}{}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$

$S_{1}$

$S_{2}$

System States

$x_{1}$

$x_{2}$

Outputs

$y_{1}:=x_{1}$

$y_{2}:=x_{2}$

Inputs

$u_{1}:=y_{2}$

$u_{2}:=y_{1}$

Equations

$\dot{x_{1}}=a_{1,1}x_{1}+a_{1,2}\operatorname{Ext}(u_{1})$

$\dot{x_{2}}=a_{2,2}x_{2}+a_{2,1}\operatorname{Ext}(u_{2})$

Efficiently, the coupled equations are:

[TABLE]

The matrix entries are chosen such that

• •

an unidirectional dependency on input is given: $a_{11},a_{21},a_{22}\neq 0$, $a_{12}=0$,

• •

a mutual dependency is given: $a_{12},a_{21}\neq 0$, $a_{11}=a_{22}=0$.

In the unidirectional flow of data case, one expects the convergence rates from (69) to become observable. In the example where the data flows in both directions, the linear extrapolation of input data case is an example for a problem where the coupling equations (9)-(10) are not exactly fulfilled when calculating the $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ and where one would have to apply a method to make outputs consistent as discussed in section 2.2.1; it might illuminate the effects of output dependency on past-timestep input data: In this case, the output row in above scheme reads

[TABLE]

and as

[TABLE]

and either $\operatorname{Ext}(x_{i})(T_{N})$ or $\operatorname{Ext}(x_{j})(T_{N})$ in the corresponding equation of the subsystem $j$ has to be the extrapolant from $(T_{N-1},T_{N}]$, one expects a bigger error for this problem. On the other hand, with $\operatorname{Ext}^{(T_{N-1},T_{N}]}(x_{i})(T_{N})$ denoting the extrapolant for the interval given in superscript,

[TABLE]

and so an order loss is not to be expected.

Figure 6 shows the convergence result for the four situations. To show the predictions made by (69) in terms of $H$, it is necessary that the subsystems methods contribution $\tilde{\tilde{C}}h_{\Delta}^{p}$ is of higher order than the extrapolation and that the method used is a one-step method, which will be discussed in Section 4.2.3. Thus dopri5, an explicite Runge-Kutta method, was chosen, using the built-in stepsize control with default absolute tolerance $10^{-12}$.

The figure shows that convergence is of order 1 for constant extrapolation and of order 2 for linear extrapolation, as predicted by (69). As discussed, there is no order loss for linear extrapolation and circular dependency of inputs, but not even an higher error, in spite of the negative effects that should occur.

The error of the first component of the lower triangular, thus unidirectionally coupled system is very low as it has no extrapolation contribution, thus indicating that the error made by $dopri$ is low enough to allow for judgement of the effect of extrapolation error.

4.2.2 Spring-Mass system for convergence and stability

To numerically examine the stability of the method, we chose the linear spring-mass oscillator as discussed in [11] and given by

[TABLE]

with mass $m$, spring constant $c$, and in which the damping constant $d$ shall vanish, for being simple while still showing stability problems. It is indeed a physical interpretation of the above system with offdiagonal matrix.

This system was treated in the cosimulation scheme 1. Output of the spring is the force $F=-cx$, that of the mass is the velocity $v=\dot{x}$. As ODE solver on subsystems, any solver that does not dominate the convergence and stability behavior of the cosimulation scheme could be used. The plots 8 and 9 show simulations done with vode and zvode from the numpy Python numerics library, which both implement implicit Adams method if problem is nonstiff and BDF if it is and behave stable due to their step size adjustment. The convergence plot 7 was made with dopri5 (see Section 4.2.3).

The numerical convergence examination backs up the results from section 3.4, estimate (69), and equation (88) – even more, balance corrected scheme for this problem converges of one order higher than proven there. A sharper theoretical result should be achievable.

But the method is unstable for its explicite contributions, as proven in section 3.5. This means $\left\|{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\right\|\longrightarrow\infty$ for $t\longrightarrow\infty$. The energy of our system is $E=\frac{1}{2}mv^{2}+\frac{1}{2}cs^{2}=\left<{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}},{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\right>_{\frac{1}{2}\operatorname{diag}(m,c)}$, which is an equivalent norm, so lack of stability is equivalent to energy augmentation.

This lack of stability can be interpreted in physics as a consequence of extrapolation errors in factors of power acting on subsystems boundaries: In [11] the problem arose that errors in the force $y_{1}$ made during data exchange lead to errors in the power that acts on the mass. The system picks up energy and behaves unstable (See figure 8). This led to a reclassification to balance errors: Cumulated extrapolation errors in input variables that in fact are conserved quantities, and those in input variables that are factors of conserved quantities and such disturb a balance [11, Section 3.2].

4.2.3 Pitfall Subsystems methods

Multistep methods need additional initial values: an $n$-step methods needs $n$ of them. Those initial values are usually calculated using one-step methods, possibly of lower order. Its error usually is not visible in the result, as it occurs only once in the calculation. But as subsystem solver restarts after data exchange, the error occurs $N$ times in cosimulation. Assume that the one-step method is of order $p^{\prime}$. Then even if $H=ch$ holds, sum of its error is $NCh^{p^{\prime}}=NHCh^{p^{\prime}-1}/c=(T_{N}-T_{0})Ch^{p^{\prime}-1}=O(h^{p^{\prime}-1})$. In figure 10, this is $0(h)$, as visible for the last two refinements. The local minimum in error is due to different signs of the two contributions $O(h^{p^{\prime}-1})$ and $O(H^{P+1})$.

5 Discussion and Conclusion

By the estimates given in section 3.4, theorem 3.6 cosimulation methods are proven to be efficient methods, if the overall system is not too vulnerable to stability problems (see Section 3.5). The stability issue will be tackled in a separate publication.

The convergence result is consistent with those in [2] and in [3], but include the error of the subsystems methods. Moreover, with the provided methodology it can be proven that balance correction method converges at least with the same convergence order as the cosimulation scheme without balance correction, see theorem 4.1 from section 4. For most problems it should be the case that convergence order rises by one if balance correction is applied, as the test problems and superficial considerations indicate.

Finding a sharper theoretical result for balance correction techniques seems a realistic future task.