## 5 Numerical Examples  
### The Viscous Burgers' Equation  
We start with an illustrative example to show the quadratic and bilinear modeling and reduction concepts proposed in the previous sections: _the viscous Burgers' equation_. The equation was first introduced by [PERSON] in 1915 and later studied by [PERSON] in 1948 [47]. In the sequel, we follow the description from [42]. As was stated there, the dynamics are characterized by the following PDE:  
\[\frac{\partial v(x,t)}{\partial t}+v(x,t)\frac{\partial v(x,t)}{\partial x}= \frac{\partial}{\partial x}\Big{(}v\frac{v(x,t)}{\partial x}\Big{)},\quad(x,t )\in(0,1)\times(0,T)\, \tag{58}\]  
subject to the initial and boundary conditions given (in our application) by:  
\[v(x,0)=f(x),\;x\in[0,1],\;v(0,t)=u(t),\;v(1,t)=0,\;t\geqslant 0\.\]  
This PDE commonly occurs in the area of fluid dynamics where it is used for modeling gas dynamics and traffic flow. The solution \(v(x,t)\) can be interpreted as a function describing the velocity at \((x,t)\). In general, the viscosity coefficient (here denoted by \(v(x,t)\)) might depend on space and time as well. In what follows, it is assumed that the viscosity coefficient \(v\) is constant w.r.t. the time variable. Furthermore, consider a zero initial condition on the system, i.e., \(f(x)=0\). Finally, we assume that the left boundary is subject to a control. A spatial discretization of Equation (58), is performed, by means of a classical first-order Euler approximation scheme with equidistant step size \(h=\frac{1}{n+1}\). Here, n denotes the number of interior points of the interval \((0,1)\). The following spacial derivative approximations are used:  
\[\frac{\partial v}{\partial x}=\frac{v(x+h)-v(x)}{h},\;\;\frac{\partial^{2}v}{ \partial x^{2}}=\frac{v(x+h)-2v(x)+v(x-h)}{h^{2}}.\]  
By substituting these formulas into (58), we manage to obtain the following QBTI characterizing the original dynamics:\[\dot{v}_{k}(t)= \begin{cases}-\frac{v_{1}(t)v_{2}(t)}{2h}+\frac{v}{h^{2}}(v_{2}(t)-2 v_{1}(t))+\Big{(}\frac{v_{1}(t)}{2h}+\frac{v}{h^{2}}\Big{)}u(t),\ \ k=1,\\ -\frac{v_{k}(t)}{2h}(v_{k+1}(t)-v_{k-1}(t))+\frac{v}{h^{2}}(v_{k+1}(t)-2v_{k}(t )+v_{k-1}(t)),\ \ 2\leqslant k\leqslant n-1,\\ -\frac{v_{n}(t)v_{k-1}(t)}{2h}+\frac{v}{h^{2}}(-2v_{n}(t)+2v_{n-1}(t)),\ \ k=n. \end{cases}  \tag{59}\]  
The observed output is chosen to be the average o fall states, i.e., \(y(t)=\frac{1}{n}\sum_{k=1}^{n}v_{k}(t)\). More precisely, for \(n=3\), we can write that  
\[\dot{v}_{1}(t) =-\frac{v_{1}(t)v_{2}(t)}{2h}+\frac{v}{h^{2}}(v_{2}(t)-2v_{1}(t) )+\bigg{(}\frac{v_{1}(t)}{2h}+\frac{v}{h^{2}}\bigg{)}u(t), \tag{60}\] \[\dot{v}_{2}(t) =-\frac{v_{2}(t)}{2h}(v_{3}(t)-v_{1}(t))+\frac{v}{h^{2}}(v_{3}(t )-2v_{2}(t)+v_{1}(t)),\] \[\dot{v}_{3}(t) =-\frac{v_{3}(t)v_{2}(t)}{2h}+\frac{v}{h^{2}}(-2v_{3}(t)+2v_{2}(t )),\ \ \text{and}\ \ y(t)=\frac{v_{1}(t)+v_{2}(t)+v_{3}(t)}{3}.\]  
In what follows, for the numerical experiments, we use \(n=90\), and we put together a QBTI system in 90 state variables described by equations in (59). The time-domain full-order solution of this system is depicted in Figure 1 on the temporal-spatial grid \([0,1]\times[0,0.2]\). Next, by making use of the Carleman linearization technique, we approximate the above \(n\)th QBTI system \(\Sigma_{\text{QB}}\) described by Equation (59), with a bilinear system of order \(\mathcal{N}=n^{2}+n=8190\). We follow the procedure discussed in the previous sections. We denote with \(\Sigma_{\text{B}}^{\otimes}\) the 8190 th order initial bilinear system obtained by means of the Carleman linearization. The model \(\Sigma_{\text{B}}^{\otimes}\) is reduced using the bilinear Loewner framework to obtain \(\hat{\Sigma}_{\text{B}}^{(1)}\) of dimension \(r=34\) by following the procedure in Section 4. As discussed there, one needs to start with a grid of interpolation points; in this case, we choose 300 logarithmically spaced points in the interval \([10^{-2},10^{2}]\imath\) with \(\imath=\sqrt{-1}\).