# Numerical Analysis for Real-time Nonlinear Model Predictive Control (NMPC)
of Ethanol Steam Reformers

Robert Joseph
Department of Mathematical and Statistical Sciences

###### Abstract

The utilization of renewable energy technologies, particularly hydrogen, has
seen a boom in interest and has spread throughout the world. Ethanol steam
reformation is one of the primary methods capable of producing hydrogen
efficiently and reliably. This paper provides an in-depth study of the
reformulated system both theoretically and numerically, as well as a plan to
explore the possibility of converting the system into its conservation form.
Lastly, we offer an overview of several numerical approaches for solving the
general first-order quasi-linear hyperbolic equation to the particular model
for ethanol steam reforming (ESR). We conclude by presenting some results that
would enable the usage of these ODE/PDE solvers to be used in non-linear model
predictive control (NMPC) algorithms and discuss the limitations of our
approach and directions for future work.

###### Contents

  1. 1 Acknowledgements
  2. 2 Introduction
  3. 3 Model Description
    1. 3.1 Chemical Reactions
  4. 4 First Order Quasi Linear PDE System
  5. 5 Conservation Law
  6. 6 Solution of the System
    1. 6.1 Method of Characteristics
    2. 6.2 Solution
  7. 7 Uniqueness and Existence Theorem
    1. 7.1 Our Model
    2. 7.2 Uniqueness Theorem for Linear Systems of First Order
  8. 8 Uniqueness Theorem for Quasi Linear Systems of First Order
  9. 9 Existence Theorem for Quasi Linear Systems of First Order
    1. 9.1 Theoretical Analysis of our System
  10. 10 Numerical Analysis
    1. 10.1 Ordinary Differential Equation Solvers
    2. 10.2 Euler Methods
    3. 10.3 Runge-Kutta
    4. 10.4 Partial Differential Equation Solvers
    5. 10.5 Finite Element Method
    6. 10.6 Tools
  11. 11 Discussion
  12. 12 Limitations
  13. 13 Future Work
  14. 14 Conclusion
  15. Appendices
    1. A Eigenvalues and Eigenvectors
    2. B ODE Solver
    3. C Simplified system

##  1 Acknowledgements

I thank professor Xinwei Yu for supervising me this semester; without him, I
wouldn’t have been able to do this. I also want to thank Professor Paul
Buckingham for all of his assistance and criticism during the semester, as
well as for the fantastic course. I’d also like to thank all of the lecturers
that came to class to deliver presentations, as well as my peers and family
for supporting me throughout this semester.

Finally, I just wanted to say thank you to everyone who is a part of my life
and who has stuck by my side and helped me when things became tough this
semester. I thank God for keeping me alive and for providing me with the
strength and wisdom I need to live my life in a fruitful way.
Robert Joseph
08/04/2022
Edmonton, Alberta, Canada

##  2 Introduction

The growth in population has resulted in a considerable increase in global
energy consumption. This is due to the fact that energy is required for
practically all activities like transportation, heating etc. Nowadays,
hydrogen has emerged as an attractive energy source that may be utilized to
store energy from renewable sources. The time has come to capitalize on
hydrogen’s potential to play a significant role in addressing critical energy
challenges. Recent breakthroughs in renewable energy technology have proved
that technical innovation has the potential to build worldwide clean energy
firms that do not rely on hydrocarbon-based fuels, therefore helping to reduce
pollution levels [3].

Because it possesses one of the highest energy densities per unit mass,
hydrogen is appealing as a fuel. According to [5], energy density is defined
as the quantity of energy stored per unit volume in a specific system or
region of space. According to this statistic, hydrogen has an energy density
of 35,000 watts per kilogram (W/kg), but lithium-ion batteries have an energy
density of just 200 watts per kilogram (W/kg) and because hydrogen has a
tenfold higher energy-to-weight ratio than lithium-ion batteries, it can
deliver a greater range while being substantially lighter. Fuel cells when
pure hydrogen is used can be turned into electricity at high efficiency,
durability, and efficiency as needed.

Because of the physical features of hydrogen, it is difficult to transport and
store, prompting us to study safe ways for effectively producing hydrogen
locally. One conceivable method for safe hydrogen transportation that appears
to be the most efficient is to store the fuel as a low-pressure liquid and
then use an on-board ethanol steam reformer (ESR) to produce hydrogen as
needed [1]. This is owing to its mobility, renewable nature, and low toxicity.

The development of real-time efficient and trustworthy control systems to
assure device efficiency while reducing the impact of interruptions such as
significant variations in interior and external temperatures during
transportation, as mentioned in [1] are crucial for the design of ethanol
steam reformers. There have been few mathematically modelled studies on the
best design of control techniques for ethanol steam reformers. Previous
research has examined the steady-state behaviour of ethanol steam reformers in
order to build proportional-integral-derivative (PID)-based decoupled control
loops while disregarding physical and operational restrictions and needs [6]
and [7]. Another study [4] presented the use of a model which manages a mass
flow control of ethanol/water and a temperature regulation of 1 kWe thermal
plasma reformer. Although these work they have certain limitations such as
neglecting physical and operational constraints and do not take maximum use of
the system’s accessible information. Some of the research relies on linear
process models, which might be inaccurate in general as in [4]. Given the
complexity of nonlinear models for ethanol steam reformers, research has been
done to employ nonlinear process models, which have a relatively high
computational cost for computing the control rule.

Model predictive control (MPC) is a sophisticated process control method that
has been utilized extensively in industrial and chemical processes since the
1980s. The MPC approach is a collection of control techniques that employ a
mathematical model of an investigated system to get control actions by
minimizing a cost function connected to specified control objectives while
taking desired system performance into account [2].

The authors of [2] conducted the first research which was a rigorous
examination of the extent of nonlinear model predictive control, in which they
analyze a mechanistic model that is a single distributed parameter system
which in fact makes this physics-based distributed parameter system for
ethanol steam reformer a unique system. In addition, the same paper presented
a method for constructing an approximation reduced-order nonlinear dynamics
model with lower computing cost while preserving the same structure of the
original equations and physical model parameters.

Further extending the study [1] introduced a zero error approached in the
model reformulation. Zero error means that no error is introduced when
reformulating the model and no additional assumptions are made. When putting
the model into a nonlinear model predictive control algorithm, this new
formulation preserves the same advantages of being dependent on physical model
parameters and considerably lowering the real-time computations required.
Figure 1 shows the ESR which is a nonlinear dynamical system comprised of a
catalytic ethanol steam reactor in series with a separation stage that
includes a selective membrane for hydrogen removal as described in [1]. The
mechanistic model described is a function of time and axial direction only
(only one spatial dimension) for the system of partial differential equations.

Our report will provide the first in-depth examination of the aforementioned
reformed system, with the goal of attempting to turn it into some type of
conservation law or similar standard method ( examples include convection
form, diffusion form etc.). We will prove the uniqueness and existence of our
system’s solutions by understanding their characteristic system which is
usually done by the Riemann method. Expanding on this, we look at any
potential singularities, and how to mitigate this during numerical analysis.
This is particularly crucial since any numerical approach we wish to use
requires us to discretize the domain (i.e. discretize both the spatial and
temporal parts of the domain and consider a bounded domain). This is
significant because it allows us to avoid using numerical approaches that
might cause the system to blow up and instead utilize more efficient numerical
methods created expressly for 1D conservation rules (examples include The
finite volume method, Approximate Riemann Solvers etc).

Our main contribution with this paper is to perform a numerical comparison of
the methods applied to our system based on dependability and efficiency, as
defined by computation time and the number of analyses required to achieve a
specific level of accuracy, as well as to guarantee that the algorithm is
stable, converges and is consistent [3]. This would be accomplished by
converting the systems of partial differential equations to a system of
ordinary differential equations, for which there are already various
techniques for solving this system of ODEs (examples include Euler Methods,
Adaptive methods etc). This may also be compared to a direct approach to
solving PDE problems utilizing numerical techniques such as finite difference
methods, finite element methods, finite volume methods, spectral methods etc.
Finally, we take this to the further by determining the best system
formulation and refining the accompanying numerical approach even more to
reduce computing costs.

Aligned with [3] the goals of the study are to enable a mechanistic model to
be employed in real-time control calculations while explicitly accounting for
input, state and output constraints with minimal computation cost. This would
open up a new field of research into nonlinear distributed parameter systems
with additional features such as mass, particle number all of which are common
in other areas of research (examples include thermodynamics, fluid dynamics
etc ). This research may be extended to other reactor systems that are
frequently encountered in chemical process control applications etc. Lastly
this would bring us closer to our goal of manufacturing hydrogen safely which
in turn could be used as green energy.

The remainder of this article is organized as follows: Section 3 describes the
ethanol steam reformer model under consideration, as well as the chemical
processes and assumptions. Section 4 discusses the original system of partial
differential equations as well as the reformed system. Section 5 defines the
conservation law and the convection form of the system. Section 6 summarizes
the efforts made to convert the system into conservation legislation. Section
6 illustrates how to solve a first-order quasi-linear PDE problem. The
Uniqueness and Existence Theorem for First-Order Linear Systems is covered in
Section 7. Sections 8 and 9 explore the Uniqueness and Existence theorem for
First Order Quasi Linear Systems and how it cannot be directly applied to our
model. Section 10 provides a brief explanation of the numerical methods that
are applied to our model. Section 11 summarizes the results as well as the
findings and programming methods employed, whereas Section 12 highlights the
limits of some of the work done in this research. Finally, Sections 13 and 14
discuss future work and the main conclusions, followed by the bibliography and
appendix.

##  3 Model Description

We now proceed to describe the Ethanol Steam Reformer(ESR) as a nonlinear
dynamical system that combines in series two process unit operations:

  1. 1.

a reformer stage comprised of a catalytic ethanol steam reactor

  2. 2.

a separation stage comprised of a selective membrane through which hydrogen
can penetrate.

In the reactor, ethanol is reformed with water to generate a gas mixture from
which hydrogen is extracted. The entire process takes place within a single
integrated module known as a staged membrane reactor which is described in 1.

Figure 1:  Staged membrane reactor scheme

###  3.1 Chemical Reactions

Chemical reactions take place in a tubular packed-bed reactor with a single
intake and output. There are 4 primary chemical reactions over cobalt-based
catalysts that take place in the staged membrane reactor are as follows [1]

$\displaystyle\mathrm{C}_{2}\mathrm{H}_{5}\mathrm{OH}$
$\displaystyle\longrightarrow\mathrm{CH}_{3}\mathrm{CHO}+\mathrm{H}_{2},\hskip
56.9055pt\text{(1a)}$ $\displaystyle\mathrm{C}_{2}\mathrm{H}_{5}\mathrm{OH}$
$\displaystyle\longrightarrow\mathrm{CO}+\mathrm{CH}_{4}+\mathrm{H}_{2},\hskip
52.63759pt\text{(1b)}$ $\displaystyle\mathrm{CO}+\mathrm{H}_{2}\mathrm{O}$
$\displaystyle\rightleftharpoons\mathrm{CO}_{2}+\mathrm{H}_{2},\hskip
85.35826pt\text{(1c)}$
$\displaystyle\mathrm{CH}_{3}\mathrm{CHO}+3\mathrm{H}_{2}\mathrm{O}$
$\displaystyle\longrightarrow 2\mathrm{CO}_{2}+5\mathrm{H}_{2}.\hskip
71.13188pt\text{(1d)}$

These four reactions are occurring in the same location and under the same
circumstances at the same time. To begin, ethanol dehydrogenates to produce
hydrogen and formaldehyde (1a), which is then reformed with water to produce
carbon dioxide (1d). Furthermore, at normal working circumstances, cobalt
catalysts are active for the Water Gas Shift (WGS) process (1c). The
unfavourable process is ethanol breakdown, which produces carbon monoxide and
methane (1b). The Pd-Ag membrane absorbs only hydrogen during the membrane
separation step, leaving waste gases on the retentate side [2].

For the mathematical modelling of the ethanol steam reformer, various
plausible and slightly overlapping assumptions are made which are as follows
[1]

  1. 1.

The concentrations in the ethanol steam reformer are easily characterized as
two plug-flow units connected in sequence.

  2. 2.

The radial dependency of pressure and temperature in the tubular reactor and
membrane separator is insignificant.

  3. 3.

In the tubular reactor and membrane separator, the fluid is entirely radially
mixed at each axial location.

  4. 4.

In comparison to convection in the axial direction, the impact of molecular
diffusion in the axial direction is modest.

  5. 5.

In the operational pressure range, the ideal gas assumption is true.

These assumptions ensure that the mechanistic model for the two process units
is a function of time and axial direction, implying that the partial
differential equations describing the system’s nonlinear dynamics have just
one spatial dimension. This leads to the model being a system of first-order
quasi-linear hyperbolic equations.

##  4 First Order Quasi Linear PDE System

We now formulate the above system as a first order quasi-linear hyperbolic
equation. Recall that the class of nonlinear distributed parameter systems
described by the first-order quasi-linear hyperbolic vector equation is
described as

$A(x,y,u)u_{x}+B(x,y,u)u_{y}=g(x,y,u)$

where $u(x,y)\in\mathbb{R}^{m}$ is an vector of distributed states with each
element being a function of the real independent variables $x$ and
$y,A(x,y,u)\in\mathbb{R}^{m\times m}$ and $B(x,y,u)\in$ $\mathbb{R}^{m\times
m}$ are real matrices, $g(x,y,u)\in\mathbb{R}^{m}$ is a vector with each
element being an algebraic function of its arguments, and
$u_{x}:=\frac{\partial u}{\partial x}$ and $u_{y}:=\frac{\partial u}{\partial
y}$. Each individual equation arises from conservation of some property, such
as mass, moles, particles, energy, or momentum. The model equations in [2] can
be written as [1]

$Mf_{t}+Nf_{z}=g(t,z,f),$

with

$\begin{gathered}f:=\left(\begin{array}[]{c}F\\\ T\end{array}\right),\quad
F=\left(\begin{array}[]{c}F_{1}\\\ \vdots\\\ F_{7}\end{array}\right)\\\
M=\left(\begin{array}[]{cc}I-\frac{1}{\|F\|}Fe^{\top}&-\frac{1}{T}F\\\
0^{\top}&\frac{c_{v}}{RT\|F\|}F\end{array}\right)\\\ \\\
N=\left(\begin{array}[]{cc}\frac{RT}{Ap}I&0\\\
\frac{RT}{Ap}e^{\top}&\frac{c_{p}}{Ap}F\end{array}\right)\\\ \\\
g=\left(\begin{array}[]{cc}RT[F_{1},F_{1},\ldots F_{1}]FV_{1}R^{\prime}\\\
RT[F_{2},F_{2},\ldots F_{2}]FV_{2}R^{\prime}\\\ \vdots\\\
RT[F_{7},F_{7},\ldots F_{7}]FV_{7}R^{\prime}\\\
U\beta(T_{f}-T)-\sum_{i}F_{i}[R^{\prime}H^{\prime}-V_{i}R^{\prime}]^{T}\par\end{array}\right)\end{gathered}$

where

  * •

$I$ is the $7\times 7$ identity matrix,

  * •

$e\in\mathbb{R}^{7}$ is the vector of ones,

  * •

$F_{j}$ is the flow rate of species $j$,

  * •

$\|F\|$ is the Euclidean norm of the vector $F$,

  * •

$T$ is the temperature (in $\mathrm{K}$ ),

  * •

$R$ is the ideal gas constant (in
$\mathrm{Pa}\mathrm{m}^{3}/(\mathrm{mol}\mathrm{K})$ ),

  * •

$\rho$ is the membrane thickness (in $\mathrm{m}$ ),

  * •

$c_{p}$ and $c_{v}$ are the heat capacities (in
$\left.\mathrm{J}/(\mathrm{mol}\mathrm{K})\right)$.

  * •

$p$ is the pressure (in bar),

  * •

$A$ is the cross-sectional area of the tubular reactor (in
$\left.\mathrm{m}^{2}\right)$,

  * •

$U$ is the overall heat transfer coefficient (in
$\mathrm{J}/\left(\mathrm{m}^{2}\mathrm{~{}s}\right.$ $\mathrm{K})$ ).

  * •

$R^{\prime}$ is a vector of the reactions indices ie
$R^{\prime}=[r_{1},r_{2},r_{3},r_{4}]$ where $i\in\\{1,2,3,4\\}$ is the
respective reaction index according to the chemical reactions described
earlier.

  * •

The expressions for the chemical reaction rates are described by elementary
kinetics in [2] as follows

$\displaystyle r_{a}$
$\displaystyle=k_{a}P_{\mathrm{C}_{2}\mathrm{H}_{5}\mathrm{OH}},$
$\displaystyle r_{b}$
$\displaystyle=k_{b}P_{\mathrm{C}_{2}\mathrm{H}_{5}\mathrm{OH}},$
$\displaystyle r_{c}$
$\displaystyle=k_{c}\left(P_{\mathrm{CO}}P_{\mathrm{H}_{2}\mathrm{O}}-\frac{P_{\mathrm{CO}_{2}}P_{\mathrm{H}_{2}}}{k_{WGS}}\right),$
$\displaystyle r_{d}$
$\displaystyle=k_{d}P_{\mathrm{CH}_{3}\mathrm{CHO}}P_{\mathrm{H}_{2}\mathrm{O}}^{3},$
$\displaystyle k_{WGS}$ $\displaystyle=e^{\frac{4577.8}{T}-4.33},$
$\displaystyle k_{i}$ $\displaystyle\left.=k_{\infty
i}e^{-E_{a_{i}}\left(\frac{1}{RT}-\frac{1}{RT_{\text{rel }}}\right.}\right),$

and where

  * •

$P_{j}$ is the partial pressure (in bar) of the $j$ th component,

  * •

$T_{\text{ref }}$ is the reference temperature, equal to $773\mathrm{~{}K}$

  * •

$\beta$ is the area per reactor volume where heat is being transferred (in
$\mathrm{m}^{2}/\mathrm{m}^{3}$ ), which is equivalent to $4/\gamma$, with
$\gamma$ being the catalytic reactor diameter (in $\mathrm{m}$ ),

  * •

$\Delta H_{i}$ is the heat of the $i$ th reaction (in
$\mathrm{J}/\mathrm{mol}$ ).

  * •

$V_{i}$ is a vector consisting of the stoichiometric coefficient
(dimensionless) of all the species in the ith reaction ie
$V_{i}=[\nu_{i,1},\nu_{i,2},\nu_{i,3},\nu_{i,4}]$

  * •

$H^{\prime}$ is the vector consisting of $\Delta H_{i}$ ie $H^{\prime}=[\Delta
H_{1},\Delta H_{2},\Delta H_{3},\Delta H_{4}]$

In the same paper the authors [1] also simplified the above system as by
writing the vector $f$ in terms of an orthonormal basis with
$\|F\|=\sqrt{7}u_{7}$. Define

$U:=\left(\begin{array}[]{c}u_{1}\\\ \vdots\\\ u_{7}\end{array}\right)$

and $C_{v}^{u}$ and $C_{p}^{u}$ to be such that

$C_{v}^{u}\left(\begin{array}[]{c}u_{1}\\\ \vdots\\\
u_{7}\end{array}\right)=C_{v}F\text{ and
}C_{p}^{u}\left(\begin{array}[]{c}u_{1}\\\ \vdots\\\
u_{7}\end{array}\right)=C_{p}F$

Then,

$Me_{i}=\begin{cases}e_{i},&i=1,\ldots,6\\\
-\frac{u_{1}}{u_{7}}e_{1}-\cdots-\frac{u_{6}}{u_{7}}e_{6},&i=7\\\
-\frac{u_{1}}{T}e_{1}-\cdots-\frac{u_{7}}{T}e_{7}+\frac{C_{v}^{u}U}{\sqrt{7}RTu_{7}}e_{8},&i=8\end{cases}$

and

$Ne_{i}=\begin{cases}\frac{RT}{Ap}e_{i},&i=1,\ldots,6\\\
\frac{RT}{Ap}e_{7}+\frac{\sqrt{7}RT}{Ap}e_{8},&i=7\\\
\frac{C_{p}^{u}U}{Ap}e_{8},&i=8\end{cases}$

Hence, the left-hand side of for

$u:=\left(\begin{array}[]{c}u_{1}\\\ \vdots\\\ u_{7}\\\ T\end{array}\right)$

can be written as

$Au_{t}+Bu_{x}=g(t,x,u)$

where

$A=\left(\begin{array}[]{ccc}1&-\frac{u_{1}}{u_{7}}&-\frac{u_{1}}{T}\\\
\ddots&\vdots&\vdots\\\ 1&-\frac{u_{6}}{u_{7}}&-\frac{u_{6}}{T}\\\
&0&-\frac{u_{7}}{T}\\\
&0&\frac{C_{v}^{u}U}{\sqrt{7}RTu_{7}}\end{array}\right)$
$B=\left(\begin{array}[]{ccc}\frac{RT}{Ap}&0&0\\\ \ddots&\vdots&\vdots\\\
0&\frac{RT}{Ap}&0\\\
0&\frac{\sqrt{7}RT}{Ap}&\frac{C_{v}^{u}U}{Ap}\end{array}\right)$

and

$g=\left(\begin{array}[]{cc}RT[u_{1},u_{1},\ldots u_{1}]UV_{1}R^{\prime}\\\
RT[u_{2},u_{2},\ldots u_{2}]UV_{2}R^{\prime}\\\ \vdots\\\
RT[u_{7},u_{7},\ldots u_{7}]UV_{7}R^{\prime}\\\
U\beta(T_{f}-T)-\sum_{i}u_{i}[R^{\prime}H^{\prime}-V_{i}R^{\prime}]^{T}\end{array}\right)$

The latter model is much easier to work on for numerical analysis due to its
simpler matrix structure. This enables us to determine the generalized
eigenvalues and left-generalized eigenvectors due to $A$ and $B$ being almost
diagonal matrices.

##  5 Conservation Law

For a given physical domain, conservation laws are frequently expressed in
integral form [9]. Assume we have a physical domain, $\Omega$, with a domain
boundary, $\partial\Omega$. Then, assuming that the physical domain is fixed,
the canonical conservation equation is of the type

$\frac{d}{dt}\int_{\Omega}Ud\mathbf{x}+\int_{\partial\Omega}\mathbf{F}(U)\cdot\mathbf{n}ds=\int_{\Omega}S(U,t)d\mathbf{x}$

where $U$ is the conserved state, $\mathbf{F}$ is the flux of the conserved
state, $\mathbf{n}$ is the outward pointing unit normal on the boundary of the
domain, and $S$ is a source term. Using Gauss’s theorem, this conservation law
may be represented as a partial differential equation(usually a hyperbolic
system) as follows

$\int_{\partial\Omega}\mathbf{F}\cdot\mathbf{n}ds=\int_{\Omega}\nabla\cdot\mathbf{F}d\mathbf{x}$
$\displaystyle\frac{d}{dt}\int_{\Omega}Ud\mathbf{x}+\int_{\partial\Omega}\mathbf{F}\cdot\mathbf{n}ds$
$\displaystyle=\int_{\Omega}Sd\mathbf{x},$
$\displaystyle\int_{\Omega}\frac{dU}{dt}d\mathbf{x}+\int_{\Omega}\nabla\cdot\mathbf{F}\cdot\mathbf{n}ds$
$\displaystyle=\int_{\Omega}Sd\mathbf{x},$
$\displaystyle\int_{\Omega}\left(\frac{\partial U}{\partial
t}+\nabla\cdot\mathbf{F}-S\right)d\mathbf{x}$ $\displaystyle=0$

Since this last equation must be valid for any arbitrary domain $\Omega$ and
$U_{t}+\nabla\cdot\mathbf{F}$ is assumed to be regular enough this implies
that the integrand must be zero everywhere or equivalently,

$\frac{\partial U}{\partial t}+\nabla\cdot\mathbf{F}=S$

The above equation is the conservation law written as a partial differential
equation. There are many reasons why converting the system into a conservation
law would be useful [9]

  1. 1.

The conservation law’s integral form allows for discontinuous solutions
(leading to Rankine-Hugoniot relation and correct shock speed). In the non-
conservative form, such a character does not exist (it allows only smooth
differentiable solutions).

  2. 2.

Discrete conservation is critical while computing shocks. If the conservation
principle is violated, a shock may travel at an incorrect speed.Successful
non-conservative procedures may fulfill some form of discrete conservation.

  3. 3.

Due to the fact that the integral equation might have non-differentiable
solutions, the equality of both formulations can break down in some
situations, resulting in weak solutions and serious numerical challenges in
simulations of such equations.

#### Example

Burger’s Equation
The inviscid Burgers’ equation is a conservation equation, or more precisely,
a first order quasilinear hyperbolic equation. Let $F(u)=\frac{u^{2}}{2}$ then
we have

$u_{t}+uu_{x}=0.$

Figure 2: Burgers Equation in two space variables

#### Convection

In Convection we have the almost conservation law form following the fact that
if we let $U$ be the ’conserved’ scalar quantity and let the fluxes be given
by [9],

$F=\mathbf{v}U,\quad S=0,$

where $\mathbf{v}(\mathbf{x},t)$ is the known vector of velocity. It should be
noted that a non-zero source term might be added but we will proceed with just
$S$ being zero. The first order partial differential equation may be used to
express this scalar conservation law.

$\frac{\partial U}{\partial t}+\nabla\cdot(\mathbf{v}U)=0$ $\frac{\partial
U}{\partial t}+\mathbf{v}\cdot\nabla U+(\nabla\cdot\mathbf{v})U=0$

Now if the velocity field is divergence free $(\nabla\cdot\mathbf{v}=0)$ then
we arrive at what is commonly referred to as the convection equation,

$\frac{\partial U}{\partial t}+\mathbf{v}\cdot\nabla U=0$

Physically, this equation indicates that the quantity $U$ does not change when
moving along the stream-wise direction (i.e. convecting with the velocity).

#### Transformation

We take our original model and notice that it breaks up into two sets of
equations

  1. 1.

The Molar Conservation Equations

  2. 2.

The Energy Conservation Equation

Now the Molar Conservation Equations can easily be converted into a
conservation law. For the reforming stage, the molar conservation equations
for the seven $j\in(1,\ldots,7)$ species
$\left(\mathrm{C}_{2}\mathrm{H}_{5}\mathrm{OH},\mathrm{H}_{2}\mathrm{O},\mathrm{CH}_{4},\mathrm{H}_{2}\right.$,
$\left.\mathrm{CO},\mathrm{CO}_{2},\mathrm{CH}_{3}\mathrm{CHO}\right)$ are
described by a set of seven nonlinear differential equations, initial
conditions, and boundary conditions [1]:

$\displaystyle p\frac{\partial F_{j}}{\partial t}-p$
$\displaystyle\frac{F_{j}}{\sum_{i}F_{i}}\sum_{i}\frac{\partial
F_{i}}{\partial t}-\frac{p}{T}F_{j}\frac{\partial T}{\partial
t}=RT\sum_{i}F_{i}\left(\sum_{i}\nu_{j,i}r_{i}\left(\left\\{F_{j}\right\\}\right)-\frac{1}{A}\frac{\partial
F_{j}}{\partial z}\right)$ $\displaystyle F_{j}(z,0)$
$\displaystyle=F_{j,0}(z),\quad\forall z\in\left[0,\ell_{1}\right]$
$\displaystyle F_{j}(0,t)$ $\displaystyle=F_{j,\text{ in }}(t),\quad\forall
t\geq 0$

where most of the constants were defined as before in section 4 and the
following boundary and initial conditions variables are as follows

  * •

\- $F_{j}$ is the molar flow (in $\mathrm{mol}/\mathrm{min}$ ) of the $j$ th
species,

  * •

$F_{j,0}$ is the molar flow of the $j$ th species at time $t=0$,

  * •

$F_{j,\text{ in }}$ is the molar flow of the $j$ th species at the reactor
inlet $(z=0)$,

  * •

$\ell_{1}$ is the axial length (in m) of the reactor,

At any given point, the concentration of the $j$ th species, $C_{j}$, and the
flow’s velocity, $v$, can be computed as [1]

$\displaystyle C_{j}$ $\displaystyle=\frac{F_{j}}{\sum_{i}F_{i}}\frac{p}{RT}$
$\displaystyle v$ $\displaystyle=\frac{1}{A}\sum_{i}F_{i}\frac{RT}{p}$

Now notice we can transform the Molar Conservation Equation into a
conservation equation by using $C_{j}$ and $v$ as defined above.

$\displaystyle\frac{\partial\mathrm{C}_{\mathrm{j}}}{\partial
t}+\frac{\partial\left(vC_{j}\right)}{\partial z}=\sum_{i}v_{j,i}r_{i},$
$\displaystyle\mathrm{C}_{\mathrm{j}}(0,z)=\mathrm{C}_{\mathrm{j},0}(z),\quad\forall
z\in[0,l_{1}]$
$\displaystyle\mathrm{C}_{\mathrm{j}}(t,0)=\mathrm{C}_{\mathrm{j},\text{ in
}}(t),\quad\forall t\geq 0$

which is in a conservation form. The spatiotemporal dynamics of the ESR
temperature are described by the energy conservation equation [2]

$\begin{array}[]{r}\frac{p}{RT}\frac{\sum_{j}c_{v_{j}}F_{j}}{\sum_{i}F_{i}}\frac{\partial
T}{\partial t}=U\beta\left(T_{f}-T\right)-\frac{RT}{A}\sum_{j}\frac{\partial
F_{j}}{\partial
z}-\frac{1}{A}\left(\sum_{j}c_{p_{j}}F_{j}\right)\frac{\partial T}{\partial
z}-\sum_{j}\sum_{i}r_{i}\left(\left\\{F_{j}\right\\}\right)\left(\Delta
H_{i}-\nu_{j,i}RT\right),\end{array}$

with initial and boundary conditions

$\begin{array}[]{ll}T(z,0)=T_{0}(z),&z\in\left[0,\ell_{1}\right]\\\
T(0,t)=T_{\text{in }}(t),&\forall t\geq 0\end{array}$

Now the above set of equation was particularly derived from these equations by
using $C_{j}$ and $v$ ie

$\displaystyle\rho_{g}c_{p_{g}}\left(v\frac{\partial T}{\partial
z}+\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)\frac{\partial
T}{\partial t}\right)=Ua\left(T_{f}-T\right)+H_{r}$
$\displaystyle\begin{aligned} &H_{r}=\sum_{j=1}^{N}-\Delta H_{j}r_{j}\quad
j=1,2,\ldots,4\\\ &T(0,z)=T_{0}(z)\quad z\in[0,L]\\\
&T(t,0)=T_{in}(t)\quad\forall t>0\end{aligned}$

Now notice that this is not in a conservation form. An attempt was made to
convert the above equation into a conservation form.

#### Convection Form

Let us now take the above equation and try to convert it into a conservation
form as a function of $F^{-1}$ where $F=\sum_{i}F_{i}$ and $v=\frac{FRT}{Ap}$.
Let $k=\frac{R}{Ap}$ where $k$ is just a constant then we get
$T=\frac{v}{Fk}$. Now substituting this in our equation we get

$\rho_{g}c_{p_{g}}\left(v\frac{\partial\frac{v}{Fk}}{\partial
z}+\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)\frac{\partial\frac{v}{Fk}}{\partial
t}\right)=Ua\left(T_{f}-T\right)+H_{r}$

$=\rho_{g}c_{p_{g}}\left(vk\frac{\partial\frac{v}{F}}{\partial
z}+\left(k+\frac{k\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)\frac{\partial\frac{v}{F}}{\partial
t}\right)=Ua\left(T_{f}-T\right)+H_{r}$

$=\rho_{g}c_{p_{g}}kv\frac{\partial\frac{v}{F}}{\partial
z}+\rho_{g}c_{p_{g}}k\frac{\partial\frac{v}{F}}{\partial
t}+k\rho_{s}c_{p_{s}}\frac{\partial\frac{v}{F}}{\partial
t}=Ua\left(T_{f}-T\right)+H_{r}$

$=v\frac{\partial\frac{v}{F}}{\partial z}+\frac{\partial\frac{v}{F}}{\partial
t}+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\frac{\partial\frac{v}{F}}{\partial
t}=\left(Ua\left(T_{f}-T\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{1}{k}$

$=v\frac{\partial\frac{v}{F}}{\partial
z}+\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\frac{\partial\frac{v}{F}}{\partial
t}\right)=\left(Ua\left(T_{f}-T\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{1}{k}$

Let $\beta=\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)$

Then we get the final equation

$=v\frac{\partial\frac{v}{F}}{\partial
z}+\beta\frac{\partial\frac{v}{F}}{\partial
t}=\left(Ua\left(T_{f}-T\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{1}{k}$

$=v\frac{\partial v{F^{-1}}}{\partial z}+\beta\frac{\partial
v{F^{-1}}}{\partial
t}=\left(Ua\left(T_{f}-T\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{1}{k}$

Let $g=v{F^{-1}}$ Then we get

$=v\beta^{-1}\frac{\partial g}{\partial z}+\frac{\partial g}{\partial
t}=\left(Ua\left(T_{f}-T\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{\beta^{-1}}{k}$

Which is not exactly in the conservation form that we needed but more of in a
convection form, which is in terms of a function $F^{-1}$.

#### Final System

Finally collecting both the attempts from above, we have two sets of systems
of equations one of them being

$\displaystyle\frac{\partial\mathrm{C}_{\mathrm{j}}}{\partial
t}+\frac{\partial\left(vC_{j}\right)}{\partial z}=\sum_{i}v_{j,i}r_{i},$
$\displaystyle\mathrm{C}_{\mathrm{j}}(0,z)=\mathrm{C}_{\mathrm{j},0}(z),\quad\forall
z\in[0,l_{1}]$
$\displaystyle\mathrm{C}_{\mathrm{j}}(t,0)=\mathrm{C}_{\mathrm{j},\text{ in
}}(t),\quad\forall t\geq 0$

and the other system where $T_{f}=vF^{-1}_{f}=g_{f}$ is just a constant

$\displaystyle
v\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)^{-1}\frac{\partial\left(v{F^{-1}}\right)}{\partial
z}+\frac{\partial\left(v{F^{-1}}\right)}{\partial
t}=\left(Ua\left(F^{-1}_{f}-\frac{v}{Fk}\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{Ap}{R}\left(1+\frac{\rho_{s}c_{p_{s}}}{\rho_{g}c_{p_{g}}}\right)^{-1}$
$\displaystyle H_{r}=\sum_{j=1}^{N}-\Delta H_{j}r_{j}\quad j=1,2,\ldots,4$
$\displaystyle vF^{-1}(0,z)=vF^{-1}_{0}(z)\quad z\in[0,L]$ $\displaystyle
vF^{-1}(t,0)=vF^{-1}_{in}(t)\quad\forall t>0$

which is equivalent to

$\displaystyle\frac{\partial g}{\partial
t}+v\beta^{-1}\frac{\partial\left(g\right)}{\partial
z}=\left(Ua\left(g_{f}-\frac{g}{k}\right)+H_{r}\right)\frac{1}{\rho_{g}c_{p_{g}}}\frac{1}{k}\beta^{-1}$
$\displaystyle H_{r}=\sum_{j=1}^{N}-\Delta H_{j}r_{j}\quad j=1,2,\ldots,4$
$\displaystyle g(0,z)=g_{0}(z)\quad z\in[0,L]$ $\displaystyle
g(t,0)=g_{in}(t)\quad\forall t>0$

As a result, our PDE system would be a tuple
$\left(C_{1},C_{2},C_{3},C_{4},C_{5},C_{6},C_{7},g\right)$, which is not
entirely in the conservation form but rather a hybrid of the conservation and
convection forms. As a consequence, the entire system is not in conservation
form. We could further decompose it by decomposing the 8 by 8 system into the
7 by 7 system of Conservation form and the last one on its own, but because
this is a coupled first order partial differential equation, we would have to
treat the system as a whole rather than independently.

##  6 Solution of the System

We now proceed to derive the solution of the above system. We use the method
of characteristics to help aid us in solving this which is introduced as
follows

###  6.1 Method of Characteristics

The primary concept of the method of characteristics is to convert a partial
differential equation to an ordinary differential equations using
characteristic curves which enables us to then solve the PDE.
A parametric curve is a the image of a mapping $\gamma(t)\in\mathbb{R}^{2}$
where $t$ ranges in some open interval, $\gamma(t):=(x(t),y(t))$. Assume
$u=u(x,y)$ is an unknown function in the $(x,y)$-plane [13]. The restriction
of $u$ to $\gamma$ is just $u(\gamma(t)):=u(x(t),y(t))$. An ordinary
differential equation of $u$ along $\gamma$ if of the form

$\frac{du}{dt}(\gamma(t))=f(t,u(\gamma(t)))$

where $\frac{du}{dt}(\gamma(t))$ is the derivative of $u$ along $\gamma(t)$.
The characteristic technique is based on locating a family of curves known as
the characteristic curves associated with a particular partial differential
equation, such that the partial differential equation reduces to an ordinary
differential equation along the curves which then enables us to solve the PDE.
For example, consider the family of curves $\gamma(t;c):=(t+c,t)$ for a
parameter $c$. If $u=u(x,y)$ is a continuously differentiable function, then
$\frac{d}{dt}u(\gamma(t))$ by the chain rule has the following form

$\frac{du}{dt}(\gamma(t))=\partial_{x}u\frac{dx}{dt}+\partial_{y}u\frac{dy}{dt}=\partial_{x}u+\partial_{y}u$

Therefore, the partial differential equation $\partial_{x}u+\partial_{y}u=0$
reduces to the ordinary differential equation $\frac{du}{dt}=0$ along
$\gamma(t)$. The solution is

$u(\gamma(t))=u(\gamma(0))=u(c,0),$

for the parameter $c$. If $u(c,0)$ is provided, for example, $u(c,0)=h(c)$,
then by the relation $c=x-y$, we obtain $u(x,y)=h(x-y)$ that is the general
solution of the equation. This can easily be extended to systems of PDE by
considering the characteristic lines in a vectorized format which will be
presented in the solution of our above system.

###  6.2 Solution

This goes as follows [1]. The first step is to multiply both sides of by a row
vector $\boldsymbol{l}^{\top}$ with $\boldsymbol{l}\in\mathbb{R}^{m}$ :

$\boldsymbol{l}^{\top}u_{t}+\boldsymbol{l}^{\top}Bu_{x}=\boldsymbol{l}^{\top}g.$

The method of characteristics analyzes the dynamics in terms of a
characteristic curve, which is the ordered pair $(x(\omega),t(\omega))$
parameterized by the real scalar $\omega$. Multiplying each term on the left-
hand side by $x(\omega)$ and $t(\omega)$ gives

$\left(\frac{1}{t_{\omega}}\boldsymbol{l}^{\top}\right)u_{t}t_{\omega}+\left(\frac{1}{x_{\omega}}\boldsymbol{l}^{\top}B\right)u_{x}x_{\omega}=\boldsymbol{l}^{\top}g.$

Now if a vector $l$ and ordered pair $(x(\omega),t(\omega))$ can be found such
that

$\frac{1}{t_{\omega}}\boldsymbol{l}^{\top}=\frac{1}{x_{\omega}}\boldsymbol{l}^{\top}B,$

then this can be written as
$\left(\frac{1}{t_{\omega}}\boldsymbol{l}^{\top}\right)\left(u_{x}x_{\omega}+u_{t}t_{\omega}\right)=\boldsymbol{l}^{\top}g.$
Application of the chain rule and multiplication by a scalar simplifies
$u_{\omega}=u_{x}x_{\omega}+u_{t}t_{\omega}$ to
$\left(\boldsymbol{l}^{\top}u_{\omega}\right)=\boldsymbol{l}^{\top}gt_{\omega}.$
Then this gives the equivalent expression

$\left(\boldsymbol{l}^{\top}Bu_{\omega}\right)=\boldsymbol{l}^{\top}gx_{\omega}$

The above expressions can be written using the chain rule as

$\quad l^{\top}\frac{du}{dt}=l^{\top}g$

and

$\quad l^{\top}B\frac{du}{dx}=l^{\top}g$

The next step is to characterize the set of vectors $l$ and ordered pairs
$(x(\omega),t(\omega))$ that satisfy the equation. From 6.2 we define the left
generalized eigenvalue decomposition, that is, multiplication by $t_{\omega}$
gives the equivalent expression
$\boldsymbol{l}^{\top}B\frac{du}{dt}=\boldsymbol{l}^{\top}g$ characterize the
$t(\omega))$ that satisfy generalized eigenvalues by $x_{\omega}$ gives the
$\boldsymbol{l}^{\top}(B-\sigma I)=0$ where

$\sigma=\frac{x_{\omega}}{t_{\omega}}=\frac{dx}{dt}$

The following equation has $m$ solutions (possibly repeated) defined by the
generalized eigenvalues $\sigma_{k}$ and left generalized eigenvectors
$\boldsymbol{l}_{k}$

$l_{k}^{\top}\left(B-\sigma_{k}I\right)=0,\quad k=1,2,\ldots,m$

The hyperbolicity of $\boldsymbol{l}$ is defined by the requirement that all
the roots $\sigma_{k},\;\forall k\in 1,2\ldots$ of $(B-\sigma I)=0$ are real
and that in addition there exists $m$ linearly independent eigenvectors. The
last step is to collect the above set of equations into a system of ODEs,

$\displaystyle\frac{\mathrm{d}x}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{k}(x,t,u)}$
$\displaystyle\boldsymbol{l}_{k}(x,t,u)^{\top}A(x,y,u)\frac{\mathrm{d}x}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{k}(x,t,u)^{\top}g(x,t,u)$ $\displaystyle
k=1,2,\ldots,m.$

which can easily be solved by standard techniques either by hand or
numerically.

##  7 Uniqueness and Existence Theorem

###  7.1 Our Model

We have our system of Partial Differential Equations which is given by

$Au_{t}+Bu_{x}=g(x,t,u)$

Now a problem arises here as we have the fact that $\operatorname{det}A=0$ and
so we cannot transform it into the equation

$u_{t}+Bu_{x}=g(x,t,u)$

We will thereby first prove the uniqueness and existence theorem to the above
system and then show how it can be applied to our system.

###  7.2 Uniqueness Theorem for Linear Systems of First Order

The method of ”energy integrals” permits a transparent uniqueness proof for
hyperbolic systems of first order [14]. We write the in- homogeneous system in
the vector and matrix form

$u_{t}+Au_{x}+Bu+c=0$

where $u$ is the unknown vector function, $A(x,t)$ and $B(x,t)$ are given
matrices, and $c(x,t)$ is a given vector. We assume, moreover, that $A$ has
continuous derivatives ie ($A\in C^{1}$ ) and that $B$ and $c$ are continuous.
Without restriction of generality we assume for the proof that
$A=\left(a_{ik}(x,t)\right)$ is a symmetric matrix and we can convert a linear
hyperbolic system into a symmetric form (e.g., into characteristic normal
form). The following is done for a quasi linear hyperbolic system( Similarly
follows for a linear system)

#### Symmetric Form

We have that

$u_{t}+B(x,t,u)u_{x}+g(x,t,u)=0$

is a first-order quasi linear system. We now set $U=Hu$, where the rows of the
matrix $H$ are the eigenvectors $l_{k}$, and substitute

$u=H^{-1}U,\quad u_{x}=H^{-1}U_{x}+H_{x}^{-1}U,\quad
u_{t}=H^{-1}U_{t}+H_{t}^{-1}U$

into (6.2). Then we obtain

$BH^{-1}U_{x}+H^{-1}U_{t}=-g-\left(BH_{x}^{-1}+H_{t}^{-1}\right)U.$

Now, multiplying by $H$ from the left, we obtain the diagonal form

$U_{t}+HBH^{-1}U_{x}=G,$

where $G=H\left(-g-BH_{x}^{-1}U-H_{t}^{-1}U\right)$ is a known vector,
depending linearly on $U$ but not depending on the derivatives of $U$. Thus
setting $B^{\prime}=HBH^{-1}$ yields the final equation

$U_{t}+B^{\prime}U_{x}=G$

which is in the symmetric form.

#### Characteristic Lines

Through a point $P$ with coordinates $\left(\xi,\tau\right)$, we draw the
characteristics $C_{1},C_{2},\cdots,C_{k}$ backward, obtaining the points
$P_{1},P_{2},\cdots,P_{k}$ of intersection with the initial line $t=0$. This
is drawn in 3

Figure 3: Characteristic Lines

#### Uniqueness Theorem Restated

###### Theorem 1

If $c=0$ and $u(x,0)=0$ in a closed interval of the line $t=0$ including all
the points $P_{x}$, then $u(\xi,\tau)=0$.

The smallest such interval, i.e., the interval cut out by the outer
characteristics through $P$, contains the domain of dependence of $P$. As a
matter of fact, if $u(x,0)=0$ on $[P_{1},P_{k}]$, then $u(x,t)=0$ in the
triangular region $\Gamma_{P}$ 3 cut out by the outer characteristics through
$P$ and having the interval $[P_{1},P_{k}]$ of the line $t=0$ as base.

#### Proof

We use Green’s identity which is equivalent to

$(u,Au)_{x}=\left(u_{x},Au\right)+\left(u,A_{x}u\right)+\left(u,Au_{x}\right)$

since $A$ is symmetric, $\left(u,Au_{x}\right)=\left(Au,u_{x}\right)$, this
identity reduces to

$2\left(u,Au_{x}\right)=(u,Au)_{x}-\left(u,A_{x}u\right).$

Taking the inner product of the differential equation (7.2) with the vector
$u$ and using the preceding relation we have (for $c=0$ )

$\frac{1}{2}(u,u)_{t}+\frac{1}{2}(u,Au)_{x}-\frac{1}{2}\left(u,A_{x}u\right)+(u,Bu)=0.$

We now introduce with a constant $\mu$ instead of $u$ the unknown vector

$v=e^{-\mu t}u$

which leads to the differential equation

$\quad v_{t}+Av_{x}+\left(B+\mu I\right)v=0$

and the initial condition $v(x,0)=0$. Let us denote $B*$ with the unit matrix
$I$ to be the following,

$B^{*}=B+\mu I.$

The preceding quadratic identity becomes, if we write $u$ again instead of
$v$,

$\frac{1}{2}(u,u)_{t}+\frac{1}{2}(u,Au)_{x}=(u,\hat{B}u).$

The quadratic form on the right-hand side is formed with the matrix

$\hat{B}=-B^{*}-\frac{1}{2}A_{x},$

which is the aim of the construction of our constant $\mu$. By choosing $\mu$
sufficiently large, we can insure that $(u,\hat{B}u)\leq 0$.

Figure 4: Trapezoid region

Now we integrate over the trapezoid-like domain $\Gamma_{h}$ in 4 which has
the boundary $\beta_{h}=P_{1}P_{k}A_{k}A_{1}$ denoting the components of the
outward unit normal vector by $x_{\nu},t_{\nu}$, and obtain Green’s formula

$\displaystyle\frac{1}{2}\iint_{\Gamma_{h}}\left[(u,u)_{t}+(u,Au)_{x}\right]dxdt\leq
0$
$\displaystyle=\frac{1}{2}\int_{\Gamma_{h}}\left[(u,u)t_{\nu}+(u,Au)x_{\nu}\right]ds\leq
0$
$\displaystyle=\frac{1}{2}\int_{A_{1}A_{k}+P_{k}P_{1}}(u,u)dx+\frac{1}{2}\int_{C_{1}+C_{k}}x_{\nu}\left(u,\left[A+\frac{t_{\nu}}{x_{\nu}}I\right]u\right)ds$

With the notation

$E(h)=\frac{1}{2}\int_{A_{1}}^{A_{k}}(u,u)dx$

this means

$E(h)-E(0)\leq-\frac{1}{2}\int_{C_{1}+C_{k}}x_{\nu}\left(u,\left[A+\frac{t_{\nu}}{x_{\nu}}I\right]u\right)ds$

We now show that the right-hand side is non positive. For this purpose we
recall that $C_{1}$ and $C_{k}$ are characteristic curves
$\left(\phi^{1}(x,t),\ldots\phi^{k}(x,t)\right)=\vec{0}^{t}$ and that,
therefore, they satisfy the differential equations

$-[\frac{\phi_{t}^{1}}{\phi_{x}^{1}},\dots\frac{\phi_{t}^{k}}{\phi_{x}^{k}}]=[l_{1},\ldots
l_{k}]$

where $\left(\l_{1},\dots\l_{k}\right)$ are eigenvalues of the matrix $A$,
i.e., values for which $A-\l I$ is singular.

Now along $C_{1}=\phi^{1}$, the outward normal has components proportional to
$\phi_{x}^{1},\phi_{t}^{1}$; hence $-t_{\nu}/x_{\nu}=\tau^{1}$. Since $C_{1}$
is the outer characteristic at the left, $\l_{1}$ is the largest eigenvalue of
$A$. Similarly, along $C_{k}=\phi^{k}$, the components of the normal are
proportional to $\phi_{x}^{k},\phi_{t}^{k}$ and $-t_{v}/x_{v}=\l_{k}$ is the
smallest eigenvalue of $A$. We recall that, for a symmetric matrix $A$, the
extreme eigenvalues can be characterized by

$\l_{1}=\max\frac{(u,Au)}{(u,u)},\quad\l_{k}=\min\frac{(u,Au)}{(u,u)}$

Consequently

$(u,u)\l_{1}\geq(u,Au),\quad(u,u)\l_{k}\leq(u,Au)$

or

$\left(u,\left[A-\l_{1}I\right]u\right)\leq
0,\quad\left(u,\left[A-\l_{k}I\right]u\right)\geqq 0.$

Since $x_{\nu}<0$ on $C_{1}$ and $x_{\nu}>0$ on $C_{k}$, we have

$-\frac{1}{2}\int_{C_{i}}x,\left(u,\left[A+\frac{t_{\nu}}{x_{\nu}}I\right]u\right)ds\leq
0\quad\forall i\in\\{1,\ldots k\\}$

This shows that the right hand side of 7.2 is non-positive and thus we can
write it as

$E(h)\leq E(0).$

Since by assumption $E(0)=0$, we have $0\leq E(h)\leq 0$; hence $E(h)=0$ for
$h>0$; therefore $u=0$ for $t=h$. Hence the uniqueness theorem is proved for
First Order Linear Systems.

##  8 Uniqueness Theorem for Quasi Linear Systems of First Order

The following uniqueness theorem remains valid for quasi-linear systems of
first order [14]

$u_{t}+B(x,t,u)u_{x}+g(x,t,u)=0$

even though the characteristic lines $C_{k}$ of the above equation depend on
the solution $u$.

We assume here that the matrices $B$ and $g$ possess continuous derivatives
with respect to $x,t$, and $u$ in the region under consideration. Let $u$ and
$v$ be two solutions of the above system be defined in a domain $D$ and having
the same initial values on the initial interval. Then we consider the
following equations

$\displaystyle u_{t}+B(x,t,u)u_{x}+g(x,t,u)=0$ (1) $\displaystyle
u(x,0)=\psi(x)$ (2)

and

$\displaystyle v_{t}+B(x,t,v)v_{x}+g(x,t,v)=0$ (3) $\displaystyle
v(x,0)=\psi(x)$ (4)

Subtracting the differential equation for $v$ from that for $u$ and denoting
this by $z$ ie ($z(x,t)=u(x,t)-v(x,t)$) we get

$\displaystyle z_{t}+B(v)z_{x}+[B(u)-B(v)]u_{x}+g(u)-g(v)=0$ (5)
$\displaystyle z(x,t)=0$ (6)

Because both $B$ and $g$ are differentiable and continuous we may apply the
mean value theorem

$B(u)-B(v)=H(u,v)z;\quad g(u)-g(v)=K(u,v)z$

where $H,K$ are continuous functions. We now consider $u,v,u_{x}$ as known
expressions in $x,t$ and substitute these expressions in $H$ and $K$ as well
as in $B(v)$; thereby becomes a linear homogeneous differential equation for
$z$ of the form

$z_{t}+bz_{x}+gz=0$

with initial values zero. The theorem proved in the above section asserts now
that $z$ is identically zero, and thus the uniqueness theorem is proved also
for quasi-linear equations.

##  9 Existence Theorem for Quasi Linear Systems of First Order

We shall now prove the existence theorem for first order Quasi linear systems

$u_{t}+B(x,t,u)u_{x}+g(x,t,u)=0$

Note that the characteristics $C_{k}$, eigenvalues $\sigma_{k}$ and
eigenvectors $l_{k}$ $\forall k\in\\{1,2,\dots m\\}$ now depend on the
specific function $u$ [14].

#### Assumptions

  1. 1.

The matrix $B(x,t;v)$ has $k$ real eigenvalues $\sigma_{k}(x,t;v)$, briefly
$\sigma(v)$, and a matrix $\Lambda(v)$ of independent eigenvectors $l_{k}$
with Lipschitz-continuous derivatives with respect to all arguments in a fixed
domain $G_{h}$ to be described presently provided we substitute for $v$ any
”admissible” function restricted by the above conditions

  2. 2.

The matrix $g$, and the initial values $\psi(x)$ possess Lipschitz-continuous
first derivatives with respect to $x,t,u$, then in a suitable neighbourhood
$0\leq t\leq h$ of a portion $g$ of the $x$-axis a uniquely determined
solution with Lipschitz-continuous first derivatives exists, provided that the
system is hyperbolic for the given initial values $u(x,0)=\psi(x)$.

  3. 3.

We now consider functions $v$, not necessarily solutions, with the fixed
initial values $v(x,0)=$ $\psi(x)$, with Lipschitz-continuous first
derivatives and satisfying the inequalities

$\|v\|\leq M,\quad\max\\{|v_{x}|,|v_{t}|\\}\leq M_{1}$

with fixed $M$ and $M_{1}$.

#### Proof

The solutions of the ordinary differential equations $dx/dt=\sigma_{k}$ are
called the characteristics $C_{k}$ of the $v$-field. For admissible functions
their slopes are uniformly bounded: $\left|\sigma_{k}\right|<\mu$, and we now
specify that a closed domain $G$ and in it the strip $G_{h}$ consists of
points such that the $v$-characteristics $C^{*}$. followed backwards towards
$t=0$ remain in $G$ and intersect the portion $g$ of the $x$-axis. The set or
”space” of all such functions $v$ is again called $S_{h}$. The eigenvectors
$l$ and eigenvalues as well as their derivatives depend on the variables
$x,t$-and $v$-and are Lipschitz continuous.

We define the characteristic differential operators in the $v$-field by
$D=\partial/\partial t+l(\partial/\partial x)$. Let us now proceed by a
natural iteration scheme.

#### Iteration

To allow enough leeway for $v$, we may set $M=N+1$ assuming a bound
$\|\psi(x)\|<N$. Then, after substituting an admissible function $v(x,t)$ in
$B$ and $g$, the above equation becomes linear. The solution $u$ of Cauchy’s
problem with the given initial values $u(x,0)=\psi(x)$ is called $u=Tv$, and
we obtain the solution $u$ of the above equation as a fixed element of the
transformation $T$, specifically as the uniform limit for $n\rightarrow\infty$
of iterations $Tu_{n}=u_{n+1}$ with $u_{0}=\psi(x)$, where $T=T_{n}$ depends
on $n$.

First we obtain $u=TV$ according to the procedure, based on the introduction
of a new function $U=\Lambda u,V=\Lambda v,\Psi=\Lambda\psi$, etc. Then in the
basic formulas in [1] we have merely to observe that the differentiation with
respect to $x$ must account for the dependence of $l,\cdots$, on $v$; that is,
we must interpret $(d/dx)l=l_{x}+l_{v}v_{x},\cdots$, thus introducing
derivatives of $v$. This leads quite directly to the lemma:

#### Lemma

For $M=$ $N+1$ we can choose a sufficiently large bound $M_{1}$ and
sufficiently small $h$ such that any function $v$ of $S_{h}$ is transformed
into another function $w=Tv$ of $S_{h}$.

Furthermore, the sequence $u_{n}$ converge uniformly in $G_{h}$ to a limit
function $u$. This function is uniquely determined and obviously satisfies the
differential equation in the characteristic form $lDu+$ $lg=0$ where

$D=\partial/\partial t+l(\partial/\partial x)$

and $l$ can be chosen $l=l_{k}$ along any characteristic line $C_{k}$. The
contracting property and uniqueness follow directly by considering for the
difference $z=u-u^{*}$ of two admissible functions $u=Tv$ and $u^{*}=Tv^{*}$
with the difference $\zeta=v-v^{*}$, the differential equation

$z_{t}+B(v)z_{x}+\left(B(v)-B\left(v^{*}\right)\right)z_{x}^{*}+g(v)-g\left(v^{*}\right)=0.$

Using the mean value theorem we have

$z_{t}+B(v)z_{x}+\zeta K=0,$

where $K$ is a bounded function of $x,t$, and obtain an inequality of the form

$\|z\|<M_{3}hk\|K\|\text{. }$

With sufficiently small $h$ we assure the inequality
$\|z\|\leq\frac{1}{2}\|\zeta\|$, and hence the contracting character of the
transformation $T$. Therefore $u_{n}$ converges uniformly to a function $u$ in
$G_{h}$. Hence the uniqueness theorem is proved.

###  9.1 Theoretical Analysis of our System

We have the reformulated system as

$Au_{t}+Bu_{x}=g(x,t,u)$

Now as stated previously a problem arises here as we have the fact that
$\operatorname{det}A=0$ and so we cannot transform it into the equation

$u_{t}+Bu_{x}=g(x,t,u)$

An attempt was made to prove the uniqueness theorem for the system by first
trying to first collect the set of equations
$u=(u_{1},u_{2},u_{3},u_{4},u_{5},u_{6},T)$ into a 7 by 7 system. We then try
to rewrite $u_{7}$ in terms of the other $u_{i},\forall i\in{1,\ldots 6}$ and
$T$. This enables $u_{7}$ being written in terms of an integral form and then
an analysis could be made, but due to lack of time we could not proceed. Hence
we take a different approach and try to use the original system formulation ie

$\displaystyle Mf_{t}+Nf_{z}=g(t,z,f)$ $\displaystyle
f:=\left(\begin{array}[]{c}F\\\ T\end{array}\right),\quad
F=\left(\begin{array}[]{c}F_{1}\\\ \vdots\\\ F_{7}\end{array}\right)$
$\displaystyle
M=\left(\begin{array}[]{cc}I-\frac{1}{\|F\|}Fe^{\top}&-\frac{1}{T}F\\\
0^{\top}&\frac{C_{v}}{RT\|F\|}F\end{array}\right)\text{, }$ $\displaystyle
N=\left(\begin{array}[]{cc}\frac{RT}{Ap}I&0\\\
\frac{RT}{Ap}e^{\top}&\frac{C_{p}}{Ap}F\end{array}\right)$

with initial and boundary conditions as follows

$\displaystyle F_{j}(z,0)=F_{j,0}(z),\quad\forall z\in\left[0,\ell_{1}\right]$
$\displaystyle F_{j}(0,t)=F_{j,\text{ in }}(t),\quad\forall t\geq 0$
$\displaystyle T(z,0)=T_{0}(z),\quad z\in\left[0,\ell_{1}\right],$
$\displaystyle T(0,t)=T_{\text{in }}(t),\quad\forall t\geq 0$

Notice we have that $\operatorname{det}(M)\neq 0$ which was also confirmed by
using Mathematica. Therefore we may solve the above system with respect to
$u_{t}$ and write it in the equivalent form

$f_{t}+M^{-1}Nf_{z}=M^{-1}g(z,t,f)$ $=f_{t}+N^{\prime}f_{z}=g^{\prime}(z,t,f)$

with $N^{\prime}=M^{-1}N$ and $g^{\prime}=M^{-1}g(z,t,f)$, which implies the
line $t=\text{constant}$ are now non characteristic or free. Hence by using
this system we just apply the above uniqueness and existence theorem that were
proved earlier to the above system.

##  10 Numerical Analysis

We take the reformulated system $Au_{t}+Bu_{x}$ to obtain the generalized
eigenvalues and left eigenvectors as mentioned in Section 6. Notice over here
that we need to find the generalized eigenvalues and left eigenvectors of
$\operatorname{det}(A-\sigma B)$. This is solved in [1] and the solution is as
follows. In particular, the analytical solution to
$\operatorname{det}(A-\sigma B)=0$ is

$\begin{gathered}\sigma_{k}=\frac{Ap}{RT},\quad k=1,\cdots,6,\\\
l_{1}=\left(\begin{array}[]{c}1\\\ 0\\\ \vdots\\\
0\end{array}\right),\ldots,\boldsymbol{l}_{6}=\left(\begin{array}[]{c}0\\\
\vdots\\\ 0\\\ 1\\\ 0\\\ 0\end{array}\right)\end{gathered}$

and $\sigma_{7}$ and $\sigma_{8}$ solve

$\operatorname{det}\left(\begin{array}[]{cc}-\sigma\frac{RT}{Ap}&-\frac{u_{7}}{T}\\\
-\sigma\frac{\sqrt{7}RT}{Ap}&\frac{C_{v}^{u}U}{\sqrt{7}RTu_{7}}-\sigma\frac{C_{p}^{u}U}{Ap}\end{array}\right)=0$

which can be written as

$\frac{RTC_{p}^{u}U}{(Ap)^{2}}\sigma^{2}-\left(\frac{C_{v}^{u}U}{\sqrt{7}Apu_{7}}-\frac{\sqrt{7}Ru_{7}}{Ap}\right)\sigma=0$

whose solutions are

$\sigma_{7}=0,\quad\sigma_{8}=\frac{Ap}{RTC_{p}^{u}U}\left(\frac{C_{v}^{u}U}{\sqrt{7}u_{7}}-\sqrt{7}Ru_{7}\right)$

Since the eigenvalues are known, computing the left generalized eigenvectors
$l_{7}$ and $l_{8}$ is straightforward. We obtain the following set of
eigenvectors for $l_{7}$ and $l_{8}$ when computed, this can be done
numerically or can be found out as follows

$l_{7}^{T}A=0$
$(l_{7}^{1},l_{7}^{2},l_{7}^{3},l_{7}^{4},l_{7}^{5},l_{7}^{6},l_{7}^{7},l_{7}^{8})\left(\begin{array}[]{ccc}1&-\frac{u_{1}}{u_{7}}&-\frac{u_{1}}{T}\\\
\ddots&\vdots&\vdots\\\ 1&-\frac{u_{6}}{u_{7}}&-\frac{u_{6}}{T}\\\
&0&-\frac{u_{7}}{T}\\\
&0&\frac{C_{v}^{u}U}{\sqrt{7}RTu_{7}}\end{array}\right)=0$
$(l_{7}^{1},l_{7}^{2},l_{7}^{3},l_{7}^{4},l_{7}^{5},l_{7}^{6},l_{7}^{7},l_{7}^{8})=\left\\{\frac{{u}_{1}}{{u}_{7}},\frac{{u}_{2}}{{u}_{7}},\frac{{u}3}{{u}_{7}},\frac{{u}_{4}}{{u}_{7}},\frac{{u}_{5}}{{u}_{7}},\frac{{u}_{6}}{{u}_{7}},1,0\right\\}$

which implies

$\boldsymbol{l}_{7}=\left(\begin{array}[]{c}\frac{{u}_{1}}{{u}_{7}}\\\
\frac{{u}_{2}}{{u}_{7}}\\\ \frac{{u}3}{{u}_{7}}\\\ \frac{{u}_{4}}{{u}_{7}}\\\
\frac{{u}_{5}}{{u}_{7}}\\\ \frac{{u}_{6}}{{u}_{7}}\\\ 1\\\
0\end{array}\right)$

Similarly we can find out the left eigenvector $l_{8}$ as follows

$l_{8}^{T}(A-\sigma_{8}B)=0$
$(l_{8}^{1},l_{8}^{2},l_{8}^{3},l_{8}^{4},l_{8}^{5},l_{8}^{6},l_{8}^{7},l_{8}^{8})\left(\left(\begin{array}[]{ccc}1&-\frac{u_{1}}{u_{7}}&-\frac{u_{1}}{T}\\\
\ddots&\vdots&\vdots\\\ 1&-\frac{u_{6}}{u_{7}}&-\frac{u_{6}}{T}\\\
&0&-\frac{u_{7}}{T}\\\
&0&\frac{C_{v}^{u}U}{\sqrt{7}RTu_{7}}\end{array}\right)-\frac{Ap}{RTC_{p}^{u}U}\left(\frac{C_{v}^{u}U}{\sqrt{7}u_{7}}-\sqrt{7}Ru_{7}\right)\left(\begin{array}[]{ccc}\frac{RT}{Ap}&0&0\\\
\ddots&\vdots&\vdots\\\ 0&\frac{RT}{Ap}&0\\\
0&\frac{\sqrt{7}RT}{Ap}&\frac{C_{u}^{p}U}{Ap}\end{array}\right)\right)=0$

This was extremely hard to solve by hand, but after some simplifications we
finally arrive with the eigenvector $l_{8}$ which is given by

$\boldsymbol{l}_{8}=\left(\begin{array}[]{c}\frac{7u_{1}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{7u_{2}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{7u_{3}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{7u_{4}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{7u_{5}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{7u_{6}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{7}+\sqrt{7}{C_{u}^{p}}^{2}UU^{\top}u_{7}^{2}-7C_{u}^{p}RUu_{7}^{3}\right)}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)\left(\sqrt{7}C_{u}^{v}U-7C_{u}^{p}Uu_{7}+7\sqrt{7}Ru_{7}^{2}\right)}\\\
{}\hfil\\\
\frac{\sqrt{7}C_{u}^{p}Uu_{7}^{2}}{T\left(C_{u}^{v}U+7Ru_{7}^{2}\right)}\\\
{}\hfil\\\ 1\end{array}\right)$

We can simply determine all of the system’s eigenvalues and use the
eigenvalues to get the system’s left eigenvectors. We then gather it into the
ODE system specified in Section 6, i.e.

$\displaystyle\frac{\mathrm{d}x}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{k}(x,t,u)}$
$\displaystyle\boldsymbol{l}_{k}(x,t,u)^{\top}A(x,t,u)\frac{\mathrm{d}u}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{k}(x,t,u)^{\top}g(x,t,u)$ $\displaystyle
k=1,2,\ldots,m.$

Note that we have a point of singularity due to the fact that $\sigma_{7}=0$.
We discuss how we deal with this in Section 11 and how the solution changes as
we change the value of $\epsilon$. The above reformulated system is easier to
deal with numerically just for the fact that both $A$ and $B$ are almost in
diagonal forms and that the eigenvalues and left-generalized eigenvectors are.

We could also proceed with the original system ie $Mf_{t}+Nf_{z}=g(t,z,f)$.
First we try to find the generalized eigenvalues and left eigenvectors of this
system $\boldsymbol{l}_{k}^{\top}\left(M-\sigma_{k}N\right)=0,\;\forall
k=1,2,\ldots,m$. In particular, the analytical solution to
$\operatorname{det}(M-\sigma N)=0$ is done as follows

$\begin{gathered}\operatorname{det}\left(\left(\begin{array}[]{cc}I-\frac{1}{\|F\|}Fe^{\top}&-\frac{1}{T}F\\\
0^{\top}&\frac{c_{v}}{RT\|F\|}F\end{array}\right)-\sigma\left(\begin{array}[]{cc}\frac{RT}{Ap}I&0\\\
\frac{RT}{Ap}e^{\top}&\frac{c_{p}}{Ap}F\end{array}\right)\right)=0\end{gathered}$

$\begin{gathered}=\operatorname{det}\left(\left(\begin{array}[]{cc}I-\frac{1}{\|F\|}Fe^{\top}-\sigma\frac{RT}{Ap}I&-\frac{1}{T}F\\\
-\sigma\frac{RT}{Ap}e^{\top}&\frac{c_{v}}{RT\|F\|}F-\sigma\frac{c_{p}}{Ap}F\end{array}\right)\right)=0\end{gathered}$

Now using properties of block matrix we get that

$\begin{gathered}=\operatorname{det}\left(\left(\begin{array}[]{cc}\left(I-\frac{1}{\|F\|}Fe^{\top}-\sigma\frac{RT}{Ap}I\right)\left(\frac{c_{v}}{RT\|F\|}F-\sigma\frac{c_{p}}{Ap}F\right)-\left(-\frac{1}{T}F\right)\left(-\sigma\frac{RT}{Ap}e^{\top}\right)\end{array}\right)\right)=0\end{gathered}$

From this easily we can figure out all the eigenvalues of the system and use
the eigenvalues to figure out the left eigenvectors for this system. We then
collect it into a system of ODE’s as mentioned in Section 6 ie

$\displaystyle\frac{\mathrm{d}z}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{k}(z,t,f)}$
$\displaystyle\boldsymbol{l}_{k}(z,t,f)^{\top}M(z,t,f)\frac{\mathrm{d}f}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{k}(z,t,f)^{\top}g(z,t,f)$ $\displaystyle
k=1,2,\ldots,m.$

Both the systems can therefore then easily be solved by any numerical method
that exists out there for ODE’s or even PDE’s. An attempt is also made to try
to figure out a closed form expression for a much simpler subset of the ODE
system.

###  10.1 Ordinary Differential Equation Solvers

In the following section we discuss some of the common techniques to solve the
above system of ODE. We will proceed with trying to numerically solve 11.

###  10.2 Euler Methods

We start of with one of the most common approaches to solve ODE’s. The Euler
method is a first-order numerical algorithm for solving ordinary differential
equations (ODEs) using a known starting point. As described in [11], suppose
we have the initial value problem to be solved on an interval $[t_{0},t_{f}]$

$y^{\prime}(t)=f(t,y(t)),\quad y\left(t_{0}\right)=y_{0}$

Now divide this interval by the mesh points such that the value for
$h=\frac{(t_{f}-t_{0})}{N}$ is the size of every step and set
$t_{n}=t_{0}+nh,\quad\forall n\in\\{0,1,2,\ldots N\\}$. Now, one step of the
Euler method from $t_{n}$ to $t_{n+1}=t_{n}+h$ is

$y_{n+1}=y_{n}+hf\left(t_{n},y_{n}\right),\quad\forall n\in\\{0,\dots N-1\\}$

The value of $y_{n}$ is an approximation of the solution to the ODE at time
$t_{n}:y_{n}\approx y\left(t_{n}\right)$. The Euler method is explicit, i.e.
the solution $y_{n+1}$ is an explicit function of $y_{i}$ for $i\leq n$ and
this is known as the Euler’s method. We may now get an approximation of a
nearby point on a curve by going a short distance down a line tangent to the
curve from any point on the curve. We begin by replacing the derivative in the
differential equation above $y^{\prime}$ by the finite difference
approximation and then integrating the differential equation between two
consecutive mesh points $t_{n}$ and $t_{n+1}$

$y\left(t_{n+1}\right)=y\left(t_{n}\right)+\int_{t_{n}}^{t_{n+1}}f(t,y(t))\mathrm{d}t,\quad
n=0,\ldots,N-1,$

and then applying the numerical integration rule

$\int_{t_{n}}^{t_{n+1}}g(t)\mathrm{d}t\approx hg\left(t_{n}\right),$

with $g(t)=f(t,y(t))$, to get

$y\left(t_{n+1}\right)\approx
y\left(t_{n}\right)+hf\left(t_{n},y\left(t_{n}\right)\right),\quad n=0,\ldots
N-1,\quad y\left(t_{0}\right)=y_{0}.$

To solve the mathematical model presented in the previous section ie section
6, a discretization of the space in the axial direction has been performed, as
represented in the following equations [6]:

$\displaystyle\frac{\partial u}{\partial x}$
$\displaystyle=\frac{U(x)-U(x-1)}{\Delta x}$ $\displaystyle\frac{\partial
T}{\partial x}$ $\displaystyle=\frac{T(x)-T(x-1)}{\Delta x}$

where $u,T$ are the molar concentration and the temperature, respectively.
Specifically, the discretization process results in a 320 order system (seven
chemical components and the temperature per 40 differential volumes) which has
been implemented Mathematica (®) and MATLAB (®) and
$\operatorname{SIMULINK}^{(B)}$ and integrated with the ode15s (stiff/NDF)
solver and NDSolve feature. An example on how to use Mathematica’s NDSolve
method or even an implementation of the above method is as follows

⬇

(* steps = N *)

(* initial = initialvalue *)

myEuler[t0_, t1_, initial_, steps_] :=

With[{h = (t1 - t0)/steps}, FoldList[#1 + h f[#2, #1] &, initial, Range[t0 +
h, t1, h]]]

(* or we can use the NDSolve method*)

NDSolve[\\{y’[t] == y[t], x[0] == 1\\}, x, \\{t, 0, 10\\}, StartingStepSize ->
1, Method -> ‘‘ExplicitEuler’’]

###  10.3 Runge-Kutta

We take a look at a generalization of the Euler technique. One of the most
popular approaches for solving ODEs is the Runge-Kutta method. This method is
a higher-order approximation of the halfway approach. Instead of shooting to
the midpoint, computing the derivative, and then firing across the full
interval, the Runge-Kutta technique takes four steps: shooting across one-
quarter of the interval, estimating the derivative, shooting to the midway,
and so on. To summarize, the goal of Runge–Kutta techniques is to approximate
a Taylor series by executing successive (weighted) Euler steps. In this case,
function evaluations (rather than derivatives) are used. This is simple to
implement in Mathematica, as seen below.

⬇

NDSolve[{y’[t] == -y[t], y[0] == 1}, y[t], {t, 0, 1},

Method -> "ExplicitMidpoint", "StartingStepSize" -> 1/10]

We can also extend this to incorporate a higher order method by repeated
function evaluation.

#### Remarks

As a result, these approaches may be utilized to solve this system of $2m$
equations. The following equations are solved with an initialization of the
state at the non-time coordinate’s inlet boundary condition, and the results
are gathered to generate ordered triplets $(x,t,u)$, which are vector $u$
values for each ordered pair $(x,t)$. The beginning values may be determined
by selecting suitable values while keeping the molar and energy conservation
equations in mind. In our case we randomize the initial values to solve the
above equation numerically. A more controlled approach would require the
presence of a domain expert to let us know what are the appropriate initial
values for the system.

###  10.4 Partial Differential Equation Solvers

In the following section we discuss the most common technique to solve PDEs

###  10.5 Finite Element Method

Explicit closed-form solutions for partial differential equations (PDEs) are
rarely available. The finite element method (FEM) is a technique to solve
partial differential equations numerically.

It is important for at least two reasons. First, the FEM is able to solve PDEs
on almost any arbitrarily shaped region. Second, the method is well suited for
use on a large class of PDEs. While it is almost always possible to conceive
better methods for a specific PDE on a specific region, the finite element
method performs quite well for a large class of PDEs.

It should be noted that the FEM can only provide an approximate answer. As a
result, it is not the preferred method of resolving a physical condition. A
closed form analytical solution is the best technique to address a physical
issue controlled by a differential equation. Unfortunately, in many real
circumstances, the analytical answer is impossible to find or does not exist.
This method can also be easily implemented in Mathematica as follows

⬇

Needs["NDSolve‘FEM‘"]

(*Load the package FEM and try to solve the laplacian with boundary
conditions*)

pde = {Laplacian[u[x, y, z], {x, y, z}] == 1,

DirichletCondition[u[x, y, z] == 0, z == -1],

DirichletCondition[u[x, y, z] == 1, z == 1]};

(*Solve the PDE using FEM*)

AbsoluteTiming[

MaxMemoryUsed[

ILUOSolution =

NDSolveValue[pde, u, Element[{x, y, z}, mesh],

Method -> {"PDEDiscretization" -> {"FiniteElement",

"PDESolveOptions" -> {"LinearSolver" -> {Automatic,

Method -> {"Krylov", Method -> "BiCGSTAB",

Tolerance -> 10^-3, "Preconditioner" -> "ILU0"}}}}}]]/1024.^2]

###  10.6 Tools

#### Mathematica

The primary programming used for all of the analysis here is Mathematica.
Mathematica is a very powerful language that can solve even handle symbolic
calculations. More details of the various methods can be found in [15]

#### Matlab

Matlab has a variety of methods that can solve not only systems of ODE but
also systems of PDE. The primary two are the built-in functions ode23 and
ode45, which implement Runge–Kutta 2nd/3rd-order and Runge–Kutta 4th/5th-order
variants, respectively. The methods described are found in Matlab’s
documentation. The Matlab results were consistent with the results using
Mathematica.

##  11 Discussion

We now proceed to numerically solve our reformulated system. We first consider
the simplified version of the system thereby considering only 2 molar
equations ie Molar Flow Rate of $\mathrm{C_{2}H_{5}OH}$ and $\mathrm{CH_{4}}$
and the same temperature equation. Therefore we have our simplified system as
$U=(u_{1},u_{2})$ and $u=(u_{1},u_{2},T)$

$\begin{gathered}A=\left(\begin{array}[]{ccc}1&-\frac{u_{1}}{u_{2}}&-\frac{u_{1}}{T}\\\
0&0&-\frac{u_{2}}{T}\\\
0&0&\frac{C_{v}^{u}U}{\sqrt{2}RTu_{2}}\end{array}\right)B=\left(\begin{array}[]{ccc}\frac{RT}{Ap}&0&0\\\
0&\frac{RT}{Ap}&0\\\
0&\frac{\sqrt{2}RT}{Ap}&\frac{C_{v}^{u}U}{Ap}\end{array}\right)\\\
g=\left(\begin{array}[]{cc}RTu_{1}UV_{1}R^{\prime}\\\
RTu_{2}UV_{2}R^{\prime}\\\
U\beta(T_{f}-T)-\sum_{i}u_{i}[R^{\prime}H^{\prime}-V_{i}R^{\prime}]^{T}\end{array}\right)\end{gathered}$

Again before attempting to solve the system we must first try to find the
generalized eigenvalues and left eigenvectors so that we could transform our
system of PDE into a system of ODEs. The eigenvalues of the above system are

$\sigma_{1}=0,\quad\sigma{2}=\frac{Ap}{RT},\quad\sigma_{3}=\frac{Ap}{RTC_{p}^{u}U}\left(\frac{C_{v}^{u}U}{\sqrt{2}u_{2}}-\sqrt{2}Ru_{2}\right)$

and

$\begin{gathered}l_{1}=\left(\begin{array}[]{c}\frac{u_{1}}{u_{2}}\\\ 1\\\
0\end{array}\right),l_{2}=\left(\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right),l_{1}=\left(\begin{array}[]{c}\frac{2u_{2}\left(-C_{u}^{v}C_{u}^{p}UU^{\top}u_{2}+\sqrt{2}{C_{u}^{p}}^{2}UU^{\top}u_{2}^{2}-2C_{u}^{p}RUu_{2}^{3}\right)}{T\left(C_{u}^{v}U+2Ru_{2}^{2}\right)\left(\sqrt{2}C_{u}^{v}U-2C_{u}^{p}Uu_{2}+2\sqrt{2}Ru_{2}^{2}\right)}\\\
\frac{\sqrt{2}C_{u}^{p}Uu_{2}^{2}}{T\left(C_{u}^{v}U+2Ru_{2}^{2}\right)}\\\
1\end{array}\right)\end{gathered}$

Figure 5: Eigenvalues and Eigenvectors of our system

This is also confirmed in 5 where the constants $c$ and $c1$ are $C_{u}^{v}=c$
and $C_{u}^{p}=c1$. Now that we have our eigenvalues and eigenvectors we have
the following ODE system

$\displaystyle\frac{\mathrm{d}x}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{k}(x,t,u)}$
$\displaystyle\boldsymbol{l}_{k}(x,t,u)^{\top}A(x,t,u)\frac{\mathrm{d}u}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{k}(x,t,u)^{\top}g(x,t,u)$ $\displaystyle
k=1,2,3.$

which when explicitly written out can be written as

$\displaystyle\frac{\mathrm{d}x_{3}}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{1}(x,t,u)}$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{3}(x,t,u)}$
$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{~{}d}t}$
$\displaystyle=\frac{1}{\sigma_{2}(x,t,u)}$
$\displaystyle\boldsymbol{l}_{1}(x,t,u)^{\top}A(x,t,u)\frac{\mathrm{d}u_{1}}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{1}(x,t,u)^{\top}g(x,t,u)$
$\displaystyle\boldsymbol{l}_{2}(x,t,u)^{\top}A(x,t,u)\frac{\mathrm{d}u_{2}}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{2}(x,t,u)^{\top}g(x,t,u)$
$\displaystyle\boldsymbol{l}_{3}(x,t,u)^{\top}A(x,t,u)\frac{\mathrm{d}T}{\mathrm{~{}d}t}$
$\displaystyle=\boldsymbol{l}_{3}(x,t,u)^{\top}g(x,t,u)$

Notice that we have $\sigma_{1}=0$ which implies that
$\frac{\mathrm{d}x_{1}}{\mathrm{~{}d}t}=\frac{1}{\sigma_{1}(x,t,u)}$ would
blow up and is a cause of singularity in our system. We replace the $0$ by
$\epsilon>0$ just to notice how the solution changes as $\epsilon$ increases.
Finally we make the following set of assumptions

  1. 1.

We now proceed to put the values for the constants ($R$, $T_{f}$, $\beta$,
stoichiometric coefficients and etc) where some of the constants have the true
values and others have been randomized.

  2. 2.

We should also note that one of the eigenvalues in our system is $0$ which is
a point of singularity. This was replaced by $\epsilon=0.00001$ just to
analyze the system better.

  3. 3.

Finally we set all the boundary conditions and initial conditions to small
non-zero values but randomized.

  4. 4.

We set the independent variable $x$ to be solved in the range $[\epsilon,10]$

We now use Mathematica’s NDSolve feature to numerically solve the system by
using the following code

⬇

(*Equations are defined in the appendix*)

(*Initial conditions are all set to epsilon or 1*)

(*Explictly without mentioning the Runge - Kutta Method was used*)

eqns = {eqns1,eqns2,eqns3,eqns4,eqns5,eqns6,y1[0] == eps,y2[0] == 1,y3[0] ==
eps,u1[0] == eps,u2[0] == eps,T[0] == eps};

sol = NDSolve[eqns,{y1,y2,y3,u1,u2,T},{x,eps,10}];

The above numerical method took exactly 113 steps to evaluate and converge to
a solution. This was calculated as follows

⬇

(*Counter to keep track of the number of steps to converge*)

Module[{c = 0},

NDSolve[qns,{y1,y2,y3,u1,u2,T},{x,eps,10},

EvaluationMonitor :> c++]; c]

We also notice that during the evaluation of the numerical method none of
evaluations failed at each step as noted in 6.

Figure 6: Successive differences vs steps

Finally we plot all the solutions $x_{1},x_{2},x_{3},u_{1},u_{2},T$ both for
$x>0$ and $x<0$ as shown in figures 7 8 9 10 11 12 13 14 15 16. Both the ODE
Solvers Runge-Kutta and Euler method result in the same numeric solutions and
almost the same figures. We also consider changing the value of $\epsilon$ to
see how the solution changes (all solutions plotted ie
$x_{1},x_{2},x_{3},u_{1},u_{2},T$ represented by the blue, yellow, green, red,
purple, orange lines respectively. We now note the following observations from
the figures

  1. 1.

As $\epsilon$ approaches $0$, the solution line $x_{3}$ approaches the $x$
axis.

  2. 2.

All the other solutions relatively remain the same and steady.

  3. 3.

When extrapolated into the negative $x$ axis we see that $x_{3}$ is highly
random and discontinuous compared to the other solutions

  4. 4.

The solution has no closed form solution. This is also confirmed in 25 26.

  5. 5.

Increasing the number of evaluation steps and working precision gives a more
accurate solution.

  6. 6.

We cannot gain a grasp of the numerical error of our system as there is no
closed form solution to compare it to, but we do know the truncation error
(error made in a single evaluation) for Euler’s Method is $O(h^{2})$ and for
Runge-Kutta Method is $O(h^{p})$ where $p$ is the order of the method chosen.

  7. 7.

Overall the solution of the system is stable and this can be fed now into a
Non Linear Model Predictive control algorithm.

Figure 7: $x_{1}$ solution for $x>0$ Figure 8: $x_{1}$ solution for $x<0$
Figure 9: $x_{2}$ solution for $x>0$ Figure 10: $x_{2}$ solution for $x<0$
Figure 11: $x_{3}$ solution for $x<0$ Figure 12: $x_{3}$ solution for $x<0$
Figure 13: $u_{1}$ solution for $x>0$ Figure 14: $u_{1}$ solution for $x<0$
Figure 15: $u_{2}$ solution for $x>0$ Figure 16: $u_{2}$ solution for $x<0$
Figure 17: $T$ solution for $x>0$ Figure 18: $T$ solution for $x<0$ Figure 19:
All solutions plotted $x>0$ and $\epsilon=0.00001$ Figure 20: All solutions
plotted $x<0$ and $\epsilon=0.00001$ Figure 21: All solutions plotted $x>0$
and $\epsilon=1$ Figure 22: All solutions plotted $x<0$ and $\epsilon=1$
Figure 23: All solutions plotted $x<0$ and $\epsilon=100$ Figure 24: All
solutions plotted $x>0$ and $\epsilon=100$ Figure 25: DSolve method to find
the closed form solution Figure 26: Error Message

##  12 Limitations

There were certain limitations and assumptions that were made due to lack of
time, some of them were

  * •

The whole 16 equations were programmed, but we need a more detailed analysis
and possibly a simplification of the equations for Mathematica to even find a
numerical solution.

  * •

Mathematica is similarly unable of handling an 8 by 8 system of PDEs,
therefore we utilize the characteristics approach to convert it back to ODEs.

  * •

Our reformulated system has a singularity point due to the fact that
$\operatorname{det}(A)=0$, and so we need to do more theoretical analysis on
how it behaves.

  * •

All of the constants were randomized; a more regulated setup is required to
set out the constants as precisely as feasible in relation to the reactants.

##  13 Future Work

Although we did get some solid results in this project over the past semester,
some of the future work for this research would include

  * •

An in depth study of the system (more theoretical Analysis). This would
involve proving the uniqueness theorem for the 7 by 7 system and how the $T$
equation behaves as $t\rightarrow\infty$. Hypothetical shock developments are
also crucial to explore as well as more attempts to convert the whole system
into a conservation law.

  * •

More research on the various numerical methods which can help us to solve this
system faster as well as more numerical analysis.

  * •

This technique may be generalized by combining these ODE solvers into
Nonlinear Model Predictive Control algorithms that can be utilized for various
Nonlinear systems.

  * •

In this paper we provided a static analysis of the our model and in the future
hope to provide a dynamic analysis of the model ie( when time changes how does
our solution evolve).

  * •

Lastly to completely numerically solve the 16 equations (Although I predict
that it would just be the same as our simplified system with $u_{7}$ being our
source of singularity and all the other solutions staying the same).

##  14 Conclusion

This paper provides a detailed analysis of the reduced-order model for an
ethanol steam reformer that employs no approximations when converting a system
of singular nonlinear distributed parameter systems to a system of nonlinear
ordinary differential equations for computing the states required in nonlinear
model predictive control formulations. We also go over the theoretical
analysis of both the original and reformed systems, which is accomplished by
establishing the uniqueness and existence theorem for quasi-linear and linear
first order systems. Despite the fact that we only have one point of
singularity, we explain how to deal with it numerically.

We also examined the potential of attempting to turn the system into a
conservation law, which appeared to be impossible, but an attempt to convert
the system into a conservation/convection law was successful. Finally, we
offer some numerical analysis of the system and some techniques to solving it.
Our results are consistent with those of prior works, allowing us to feed
these ODE solvers into a Nonlinear Model Predictive Control method.

Because of the cheap online compute cost of the model equations (about $2m$
for ordinary differential equations), a mechanistic model may be employed in
real-time control computations while explicitly accounting for input, state,
and output limits.

In conclusion this research might be expanded to generalize the approach so
that we could include additional reactor systems that are often used in
chemical process control applications, for example. Finally, this would move
us closer to our aim of producing hydrogen safely, which might then be
utilized as green energy.

## Bibliography

  * [1] Pedro Reyero , Xinwei Yu, Carlos Ocampo-Martinez, Richard D. Braatz, _A Reduced-order Model for Real-time NMPC of Ethanol Steam Reformers_ , 2020
  * [2] Pedro Reyero-Santiago, Carlos Ocampo-Martinez, Richard D. Braatz, _Nonlinear Dynamical Analysis for an Ethanol Steam Reformer: A Singular Distributed Parameter System_ , 2020
  * [3] Robert Joseph, _Numerical Analysis for Real-time NMPC of Ethanol Steam Reformers_ , Proposal, 2022
  * [4] M. Torchio, C. Ocampo-Martinez, L. Magni, M. Serra, R. D. Braatz and D. M. Raimondo, _Fast Model Predictive Control for hydrogen outflow regulation in Ethanol Steam Reformers_ , 2016
  * [5] Kasper T. Møller, Torben R. Jensen, Etsuo Akiba, Hai-wen Li, _Hydrogen - A sustainable energy carrier, Progress in Natural Science: Materials International_ , ISSN 1002-0071
  * [6] S. Biset, L. N. Deglioumini, M. Basualdo, V. M. Garc’ia, and M. Serra, _Analysis of the control structures for an integrated ethanol processor for proton exchange membrane fuel cell systems_ , Journal of Power Sources, vol. 192, pp. 107–113, 2009.
  * [7] V. M. Garc´ıa, E. Lopez, M. Serra, J. Llorca, and J. Riera, _Dynamic modeling and controllability analysis of an ethanol reformer for fuel cell application_ , International Journal of Hydrogen Energy, vol. 35, pp. 9768–9775, 2010
  * [8] Wen Shen, _A Very Brief Introduction to Conservation Laws_ , Department of Mathematics, Penn State University, 2013.
  * [9] Yu Xinwei, Intermediate PDE (Fall 2012), _Method of Characteristics_ , 2012
  * [10] MIT, _Conservation Laws_ , Partial Differential Equations, 2020
  * [11] Michael Zeltkevic, _Forward and Backward Euler Methods_ , Differential Equation Notes, 1998
  * [12] Qingkai Kong, Timmy Siauw, Alexandre M. Bayen, _Chapter 22 - Ordinary Differential Equations (ODEs) Initial-Value Problems_ , Python Programming and Numerical Methods, Academic Press, 2021, Pages 371-398, ISBN 9780128195499,
  * [13] Mohammad Niksirat, _Partial Differential Equations_ , Chapter 2 - First Order PDEs, 2022
  * [14] Richard Courant, Hilbert David, _Methods of mathematical physics. / Volume II, Partial differential equations_ , pp. 551-798, 1989
  * [15] Wolfram Mathematica Tutorial Collection, _Advanced Numerical Differential Equation_ , 2008

## Appendices

###  A Eigenvalues and Eigenvectors

⬇

(*Define the matrices A,B,g and the constants values have been randomized*)

A = 100;

V1 = V2 = V3 = v4 = V5 = V6 = V7 =1;

c = c1 = 100;

R’ = R = 10;

H = 10;

A =
{{1,0,0,0,0,0,-u1[x]/u7[x],-u1[x]/T[x]},{0,1,0,0,0,0,-u2[x]/u7[x],-u2[x]/T[x]},{0,0,1,0,

0,0,-u3[x]/u7[x],-u3[x]/T[x]},{0,0,0,1,0,0,-u4[x]/u7[x],-u4[x]/T[x]},{0,0,0,0,1,0,

-u5[x]/u7[x],-u5[x]/T[x]},{0,0,0, 0,0,1,-u6[x]/u7[x],-u6[x]/T[x]},{0,0,0,0,0

,0,0,-u7[x]/T[x]},{0,0,0,0,0,0,0,(c U)/(Sqrt[7] R T[x] u7[x])}}

B = {{(R T[x])/(A p), 0,0,0,0,0,0,0},{0,(R T[x])/(A p), 0,0,0,0,0,0}, {0,0, (R
T[x])/(A p), 0,0,0,0,0}, {0,0,0,(R T[x])/(A p), 0,0,0,0}, {0,0,0,0,(R T[x])/(A
p), 0,0,0}, {0,0,0,0,0,(R T[x])/(A p), 0,0}, {0,0,0,0,0,0,(R T[x])/(A p),
0},{0,0,0,0,0,0,(Sqrt[7] R T[x])/(A p),(c1 U)/(A p)}}

g = {{R T[x] (u1[x] u1[x] + u1[x] u2[x] +u1[x] u3[x] + u1[x] u4[x] + u1[x]
u5[x] + u1[x] u6[x] +u1[x] u7[x]) V1 R’},{R T[x](u2[x] u1[x] + u2[x] u2[x]
+u2[x] u3[x] + u2[x] u4[x] + u2[x] u5[x] + u2[x] u6[x] +u2[x] u7[x]) V2 R’},{R
T[x] (u3[x] u1[x] + u3[x] u2[x] +u3[x] u3[x] + u3[x] u4[x] + u3[x] u5[x] +
u3[x] u6[x] +u3[x] u7[x]) V3 R’},{R T[x] (u4[x] u1[x] + u4[x] u2[x] +u4[x]
u3[x] + u4[x] u4[x] + u4[x] u5[x] + u4[x] u6[x] +u4[x] u7[x]) V4 R’},{R T[x]
(u5[x] u1[x] + u5[x] u2[x] +u5[x] u3[x] + u5[x] u4[x] + u5[x] u5[x] + u5[x]
u6[x] +u5[x] u7[x]) V5 R’},{R T[x] (u6[x] u1[x] + u6[x] u2[x] +u6[x] u3[x] +
u6[x] u4[x] + u6[x] u5[x] + u6[x] u6[x] +u6[x] u7[x]) V6 R’},{R T[x] (u7[x]
u1[x] + u7[x] u2[x] +u7[x] u3[x] + u7[x] u4[x] + u7[x] u5[x] + u7[x] u6[x]
+u7[x] u7[x]) V7 R’},{U1 B (Tf - T[x]) -u1[x] (R’ H - V1 R T[x]) - u2[x](R’ H
- V2 R T[x]) - u3[x](R’ H - V3 R T[x]) - u4[x](R’ H - V4 R T[x]) - u5[x](R’ H
- V5 R T[x]) - u6[x](R’ H - V6 R T[x]) - u7[x](R’ H - V7 R T[x])}}

(*Find the generalized Eigenvalues and Eigenvectors of our system*)

l = Eigenvectors[{a,b}]

s = Eigenvalues[{a,b}]

###  B ODE Solver

⬇

(*Define the 16 equations of our system *)

eqns1 = y1’[x] == 1/(s[[2]]);

eqns2 = y2’[x] == 1/(s[[3]]);

eqns3 = y3’[x] == 1/(s[[4]]);

eqns4 = y4’[x] == 1/(s[[5]]);

eqns5 = y5’[x] == 1/(s[[6]]);

eqns6 = y6’[x] == 1/(s[[7]]);

eqns7 = y7’[x] == 1/(s[[1]] + 1);

eqns8 = y8’[x] == 1/(s[[8]]);

eqns9 = u1’[x] == (R T[x] (u1[x] u6[x] + u2[x] u6[x] +u3[x] u6[x] + u4[x]
u6[x] + u5[x] u6[x] + u6[x] u6[x] +u7[x] u6[x]) V6 R’);

eqns10 = u2’[x] == (R T[x] (u1[x] u5[x] + u2[x] u5[x] +u5[x] u5[x] + u4[x]
u5[x] + u5[x] u5[x] + u6[x] u5[x] +u7[x] u5[x]) V5 R’);

eqns11 = u3’[x] == (R T[x] (u1[x] u4[x] + u2[x] u4[x] +u5[x] u4[x] + u4[x]
u4[x] + u5[x] u4[x] + u6[x] u4[x] +u7[x] u4[x]) V4 R’);

eqns12 = u4’[x] == (R T[x] (u1[x] u3[x] + u2[x] u3[x] +u5[x] u3[x] + u4[x]
u3[x] + u5[x] u3[x] + u6[x] u3[x] +u7[x] u3[x]) V3 R’);

eqns13 = u5’[x] == (R T[x] (u1[x] u2[x] + u2[x] u2[x] +u5[x] u2[x] + u4[x]
u2[x] + u5[x] u2[x] + u6[x] u2[x] +u7[x] u2[x]) V2 R’);

eqns14 = u6’[x] == (R T[x] (u1[x] u1[x] + u2[x] u1[x] +u5[x] u1[x] + u4[x]
u1[x] + u5[x] u1[x] + u6[x] u1[x] +u7[x] u1[x]) V1 R’);

array1 = l[[8]]*g;

ufinal = {u1’[x],u2’[x],u3’[x],u4’[x],u5’[x],u6’[x],u7’[x],T’[x]};

new1 = ArrayReshape[array1, Dimensions[array1] ~DeleteCases~ 1];

eqns15 = Total[l[[8]].a ufinal] == Total[new1];

eqns16 = -((u1[x]u1[x]/(u7[x] u7[x])) +(u2[x]u2[x]/(u7[x] u7[x])) +
(u3[x]u3[x]/(u7[x] u7[x])) + (u4[x]u4[x]/(u7[x] u7[x])) + (u5[x]u5[x]/(u7[x]
u7[x])) + (u6[x]u6[x]/(u7[x] u7[x]))) u7’[x] == (((R T[x] u6[x] (u1[x] u6[x] +
u2[x] u6[x] +u3[x] u6[x] + u4[x] u6[x] + u5[x] u6[x] + u6[x] u6[x] +u7[x]
u6[x]) V6 R’)/(u7[x])) + ((R T[x] u5[x] (u1[x] u5[x] + u2[x] u5[x] +u3[x]
u5[x] + u5[x] u6[x] + u5[x] u5[x] + u6[x] u5[x] +u7[x] u5[x]) V5 R’)/(u7[x]))
+ ((R T[x] u4[x] (u1[x] u4[x] + u2[x] u4[x] +u3[x] u4[x] + u4[x] u4[x] + u5[x]
u4[x] + u6[x] u4[x] +u7[x] u4[x]) V4 R’)/(u7[x])) + ((R T[x] u3[x] (u1[x]
u3[x] + u2[x] u3[x] +u3[x] u3[x] + u4[x] u3[x] + u5[x] u3[x] + u6[x] u3[x]
+u7[x] u3[x]) V3 R’)/(u7[x])) + ((R T[x] u2[x] (u1[x] u2[x] + u2[x] u2[x]
+u3[x] u2[x] + u4[x] u2[x] + u5[x] u2[x] + u6[x] u2[x] +u7[x] u2[x]) V2
R’)/(u7[x])) + ((R T[x] u1[x] (u1[x] u1[x] + u2[x] u1[x] +u3[x] u1[x] + u1[x]
u1[x] + u1[x] u6[x] + u1[x] u1[x] +u1[x] u6[x]) V1 R’)/(u7[x])) + R T[x]
(u1[x] u7[x]+u2[x] u7[x]+u3[x] u7[x]+u4[x] u7[x]+u5[x] u7[x]+u6[x] u7[x]+u7[x]
u7[x]) V7 R’);

eps = 10^-5;

(*Solve the system numerically with initial conditions*)

qns = {eqns1,eqns2,eqns3,eqns4,eqns5,eqns6,eqns7, eqns8, eqns9, eqns10,
eqns11, eqns12, eqns13, eqns14, eqns15, eqns16,y1[0] == 0,y2[0] == 0,y3[eps]
== 0,y4[eps] == eps,y5[eps] == eps,y6[eps] == eps,y7[0] == 0,y8[0] == 0,
u1[eps] == 1,u2[eps] == eps, u3[eps] == eps,u4[eps] == eps, u5[0] == eps,u6[0]
== eps,u7[eps] == eps,T[0] == eps};

sol
=NDSolve[qns,{y1,y2,y3,y4,y5,y6,y7,y8,u1,u2,u3,u4,u5,u6,u7,T},{x,eps,1000}];

###  C Simplified system

⬇

(*Definining our system*)

a1= {{1,-u1[x]/u2[x],-u1[x]/T[x]},{0,0 ,-u2[x]/T[x]},{0,0,(c U)/(Sqrt[2] R
T[x] u2[x])}}

b1= {{(R T[x])/(A p), 0,0},{0,(R T[x])/(A p), 0},{0,(Sqrt[2] R T[x])/(A p),(c1
U)/(A p)}}

g1 = {{R T[x] (u1[x] u1[x] + u1[x] u2[x] ) V1 R’},{R T[x] (u2[x] u1[x] + u2[x]
u2[x] ) V2 R’},{U1 B (Tf - T[x]) -u1[x] (R’ H - V1 R T[x]) - u2[x](R’ H - V2 R
T[x]) }}

(*Eigenvalues and Eigenvectors of our system*)

s1 = Eigenvalues[{a1,b1}]

l1 = Eigenvectors[{a1,b1}]

(*Define constants and our 6 equations*)

array1 = l1[[3]]*g1;

new1 = ArrayReshape[array1, Dimensions[array1] ~DeleteCases~ 1];

new1;

epsilon =1000000;

R = 6.022*10;

R’ = 1;

V1 = 2;

eps = 10^-5;

inf = 5;

V2 = 3;

c1 = 10;

U =1 ;

H = 1;

Tf = 232;

B = 10;

U1 = 10;

c = 10;

A = 343;

p = 10;

u10 = {u1’[x],u2’[x],T’[x]};

eqns1 = y1’[x] == 1/(s1[[2]]);

eqns2 = y2’[x] == 1/(s1[[3]]);

eqns3 = y3’[x] == 1/(s1[[1]]+epsilon);

eqns4 = u1’[x] + u1[x] (u2’[x]/u2[x]) + u1[x] (T’[x]/T[x]) == (R T[x]
(u1[x]u1[x] + u1[x]u2[x]) V1 R’);

eqns5 =(-(u1[x]^2/(T[x] u2[x])) - u2[x]/T[x]) Derivative[1][T][x] + (u1[x]
Derivative[1][u1][x])/u2[x] - (u1[x]^2 Derivative[1][u2][x])/u2[x]^2 == (R
T[x] u1[x] (u1[x]^2 + u1[x] u2[x]) V1 Derivative[1][R])/u2[x] + R T[x] (u1[x]
u2[x] + u2[x]^2) V2 Derivative[1][R];

eqns6 = Total[l1[[3]].a1 u10] == Total[new1];

(*Solution of our system*)

qns = {eqns1,eqns2,eqns3,eqns4,eqns5,eqns6,y1[0] == eps,y2[0] == 1,y3[0] ==
eps,u1[0] == eps,u2[0] == eps,T[0] == eps};

sol = NDSolve[qns,{y1,y2,y3,u1,u2,T},{x,eps,10},Method ->
"ExplicitRungeKutta"];

(*2D Plot of all our solutions*)

DynamicModule[{xc,dx},Manipulate[Plot[Evaluate[{y1[x],y2[x],y3[x],u1[x],u2[x],T[x]}/.
sol[[1]]],{x,xc-dx,xc+dx}],

{{xc,1000001/200000,"center"},-(4999993/200000),1399999/40000},{{dx,5.,"zoom"},2999997/100000,0.3}]

,DynamicModuleValues:>{}]