**Numerical and Symbolic Computation Developments and Applications**

Maria Amélia Ramos Loja and Joaquim Infante Barbosa Edited by

> Printed Edition of the Special Issue Published in *Mathematical and Computational Applications*

www.mdpi.com/journal/mca

## **Numerical and Symbolic Computation**

## **Numerical and Symbolic Computation Developments and Applications**

Special Issue Editors **Maria Am´elia Ramos Loja Joaquim Infante Barbosa**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Maria Am´elia Ramos Loja CIMOSM, ISEL, Instituto Superior de Engenharia de Lisboa Portugal

Joaquim Infante Barbosa IDMEC, IST, Instituto Superior Técnico, Universidade de Lisboa Portugal

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Mathematical and Computational Applications* (ISSN 2297-8747) (available at: https://www.mdpi.com/ journal/mca/special issues/SYMCOMP2019).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-302-5 (Pbk) ISBN 978-3-03936-303-2 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Special Issue Editors**

**Maria Am ´elia Ramos Loja** is presently an Adjunct Professor in the Mechanical Engineering Department of the Engineering Institute of Lisbon (ISEL, IPL), an invited Associate Professor of the Physics Department of the University of Evora (UEvora) and a Senior Researcher of the Mechanical ´ Engineering Institute (IDMEC, IST). Her academic background integrates a BSc with honors in Marine Engineering from the Portuguese Nautical School and a BSc in Computer Science. Her MSc and Ph.D. degrees in Mechanical Engineering were conferred by the IST, University of Lisbon and the Habilitation in Mechatronic Engineering by the University of Evora. Her principle areas ´ of interest include the scientific areas of Computational Solids Mechanics, Composite Materials, Smart Materials, Optimization, and Reverse Engineering. She is Chairperson of the ECCOMAS thematic series of conferences SYMCOMP (International Conference on Numerical and Symbolic Computation: Developments and Applications) and she coordinates the Research Centre on Modelling and Optimization of Multifunctional Systems (CIMOSM, ISEL). Since 2017 she has been invited by European Commission Research Agencies to evaluate project proposals in different subjects related to her competences.

**Joaquim Infante Barbosa** is presently a Jubilated Full Professor for the Mechanical Engineering Department of the Engineering Institute of Lisbon (ISEL, IPL, Polytechnic Institute of Lisbon) and Senior Researcher of the Mechanical Engineering Institute (IDMEC, IST, University of Lisbon). His academic background integrates a degree in Mechanical Engineering by IST, University of Lisbon. His MSc and Ph.D. degrees in Mechanical Engineering were conferred by IST, University of Lisbon and the Habilitation in Structural Mechanics by the University of Evora. His major areas of ´ interest include the scientific areas of Computational Mechanics, Structural Optimization, Composite Materials and Vibration Suppression, among others. He is a researcher on several Portuguese and European projects and a reviewer of various engineering journals. He is a member of the organizing committee of the ECCOMAS thematic series of conferences SYMCOMP (International Conference on Numerical and Symbolic Computation: Developments and Applications) and a senior member of APMTAC—Portuguese Association of Theoretical, Applied and Computational Mechanics.

## *Editorial* **Preface to** *Numerical and Symbolic Computation: Developments and Applications—2019*

**Maria Amélia R. Loja 1,2,3,\* and Joaquim I. Barbosa 1,3**


Received: 11 May 2020; Accepted: 11 May 2020; Published: 11 May 2020

This book constitutes the printed edition of the Special Issue *Numerical and Symbolic Computation: Developments and Applications—2019*, published by *Mathematical and Computational Applications* (MCA) and comprises a collection of articles related to works presented at the 4th International Conference in Numerical and Symbolic Computation—SYMCOMP 2019—that took place in Porto, Portugal, from April 11th to April 12th 2019.

This conference series has a multidisciplinary character and brings together researchers from very different scientific areas, aiming at sharing different experiences, in a cross-fertilization perspective. Therefore, the articles contained in this book, although sharing a common characteristic related to the use of numerical and/or symbolic methods and computational approaches, also present an overview of their use in a transversal way to science and engineering fields.

In the first contribution *Bridging Symbolic Computation and Economics: A Dynamic and Interactive Tool to Analyze the Price Elasticity of Supply,* from Andraz et al. [1], the authors propose a new dynamic and interactive tool, created in the computer algebra system Mathematica and available in the Computable Document Format. This tool can be used as an active learning tool to promote better student activity and engagement in the learning process, among students enrolled in socio-economic programs.

The second article of the book is authored by Escobar et al. [2] and has the title *The Invariant Two-Parameter Function of Algebras* ψ. In this article, it is proven that the five-dimensional classical-mechanical model built upon certain types of five-dimensional Lie algebras cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra. Other applications to physical problems are also addressed.

Gavina et al. [3], in their article *Solving Nonholonomic Systems with the Tau Method*, propose a numerical procedure based on the spectral tau method to solve nonholonomic systems. The Lanczos' spectral tau method is used to obtain an approximate solution to these nonholonomic problems exploiting the tau toolbox software library, adding to the ease of use characteristics and providing accurate results.

The contribution of Matos and Rodrigues [4], *Almost Exact Computation of Eigenvalues in Approximate Di*ff*erential Problems*, addresses differential eigenvalue problems that arise in many fields of Mathematics and Physics. These authors present a method for eigenvalues computation following the Tau method philosophy and using Tau Toolbox tools, wherein the eigenvalue differential problem is translated into an algebraic approximated eigenvalues problem, after which by making use of symbolic computations, they arrive at the exact polynomial expression of the determinant of the algebraic problem matrix, allowing us to get high accuracy approximations of differential eigenvalues.

In a different area, Monteiro et al. [5], through their article *Factors for Marketing Innovation in Portuguese Firms CIS 2014*, aim at understanding which factors influence marketing innovation and also aim to establish a business profile of firms that innovate or do not in marketing. These authors used multivariate statistical techniques, such as, multiple linear regression (with the Marketing Innovation Index as dependent variable) and discriminant analysis where the dependent variable is a dummy variable, indicating if the firm innovates or not in marketing.

The sixth article *Numerical Optimal Control of HIV Transmission in Octave*/*MATLAB,* from to Campos et al. [6], provides a GNU Octave/MATLAB code for the simulation of mathematical models described by ordinary differential equations and for the solution of optimal control problems through Pontryagin's maximum principle. A control function is introduced into the normalized HIV model and an optimal control problem is formulated, where the goal is to find the optimal HIV prevention strategy that maximizes the fraction of uninfected HIV individuals with the least amount of new HIV infections and cost associated with the control measures.

The contribution of Rodrigues [7] entitled *Isogeometric Analysis for Fluid Shear Stress in Cancer Cells* constitutes the seventh and last paper of this book. In this article, the author considers the modelling of a cancer cell using non-uniform rational b-splines (NURBS) and uses isogeometric analysis to model the fluid-generated forces that tumor cells are exposed to, in the vascular and tumor microenvironments, during the metastatic process. The aim of the article is focused on the geometrical sensitivities to the shear stress exhibition of the cell membrane.

At this point, as editors of this book, we would like to express our deep gratitude for the opportunity to publish with MDPI. This acknowledgment is deservedly extensive to the MCA Editorial Office and more particularly to Mr. Everett Zhu, who has permanently supported us in this process. It was a great pleasure to work in such conditions. We look forward to collaborating with MCA in the future.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Bridging Symbolic Computation and Economics: A Dynamic and Interactive Tool to Analyze the Price Elasticity of Supply**

#### **Jorge M. Andraz 1,2, Renato Candeias <sup>2</sup> and Ana C. Conceição 3,\***


Received: 1 August 2019; Accepted: 9 October 2019; Published: 10 October 2019

**Abstract:** It is not possible to achieve the objectives and skills of a program in economics, at the secondary and undergraduate levels, without resorting to graphic illustrations. In this way, the use of educational software has been increasingly recognized as a useful tool to promote students' motivation to deal with, and understand, new economic concepts. Current digital technology allows students to work with a large number and variety of graphics in an interactive way, complementing the theoretical results and the so often used paper and pencil calculations. The computer algebra system *Mathematica* is a very powerful software that allows the implementation of many interactive visual applications. Thanks to the symbolic and numerical capabilities of *Mathematica*, these applications allow the user to interact with the graphical and analytical information in real time. However, *Mathematica* is a commercially distributed application which makes it difficult for teachers and students to access. The main goal of this paper is to present a new dynamic and interactive tool, created with *Mathematica* and available in the Computable Document Format. This format allows anyone with a computer to use, at no cost, the PES(Linear)-Tool, even without an active Wolfram *Mathematica* license. The PES(Linear)-Tool can be used as an active learning tool to promote better student activity and engagement in the learning process, among students enrolled in socio-economic programs. This tool is very intuitive to use which makes it suitable for less experienced users.

**Keywords:** symbolic computation; dynamic and interactive tool; socio-economic sciences; F-Tool concept; PES(Linear)-Tool; Wolfram *Mathematica*; computable document format

#### **1. Introduction**

The use of educational software has been increasingly recognized as a useful tool to promote students' motivation to deal with, and understand, new concepts in different study fields (see, for instance, [1–12]). In fact, educational software tools have a great potential of applicability, particularly at the university level, where the knowledge of various areas by different careers is required [8]. Current digital technology allows students to work with a large number and variety of graphics, in an interactive way, complementing the theoretical results and the so often used paper and pencil calculations. Obviously, calculations with this kind of support do not replace paper and pencil calculations, and they should be properly combined with other methods of calculation, including mental calculation. Some studies conclude that students using computer algebra systems are at least as good in "pencil and paper" skills as their traditional counterparts [13]. This aspect is not

of minor relevance. Although the "pencil and paper" work can be done by computers, students should learn how calculations are made and also should learn how the computer algebra systems work [14] (we thank an anonymous referee for this observation). Also, the use of technology in the classroom can lead to advances in conceptualization, contributing thereby to students' engagements and motivation [15]. According to [16], one of the reasons for students to use computer algebra systems is their belief that these tools help their understanding of new concepts.

The computer algebra system *Mathematica*, conceived by Stephen Wolfram, and developed by Wolfram Research, is a very powerful software that allows the implementation of many interactive visual applications. Thanks to the symbolic and numerical capabilities of Mathematica, these applications are eminently dynamic tools, where the user can interact with the graphical and analytical information in real time. More importantly, the graphics are taken out of the textbook and they are placed under the user's control, where the user can manipulate, investigate, and explore their characteristics. Students who have used *Mathematica* for at least one year identified this kind of visualization as one of the significant benefits they found with the use of *Mathematica* [16].

Graphics are always helpful in the learning process, but [16] states that it makes a difference whether the students' interaction with graphic visualization is active or passive. As reported by [17], academics in higher education institutions should not only worry about the contents, but also give attention to the learning environment as they face students with different motivations and different levels of involvement. Such differences will likely affect the teaching and learning process. Moreover, teachers can expect that, in any classroom, some students prefer to be receivers (observers or listeners), while others prefer to be active participants. In fact, there are students with a more active attitude, who, even in a more traditional class, theorize, apply and relate, and there are those who exhibit more passive behavior. Clearly, these students require different orientation and teaching methods so that they are able to fully engage in the classroom activities as agents of a truly active learning process. This type of learning denotes a style of teaching that provides opportunities for students to talk, to listen to, and to reflect on what they have learned, as they participate in a variety of learning activities [18,19]. We should note that teachers who employ active learning strategies in their classrooms are unlikely to please all students all the time [20], but neither is a teacher who relies regularly on traditional lectures. The active learning also aims to improve the students' performance and develop the skills they need, for example, to obtain a better grade in a specific curricular unit [19]. In many cases, active learning can be employed without increased costs and with only a modest change in current teaching practices with a reduced risk and a high return [20]. Unfortunately, there are gaps between teaching and learning, between teaching and testing, and between educational research and practice in higher education institutions [21]. A serious gap also exists between how faculty members typically teach (i.e., relying largely on the "lecture method") and how they know they should teach (i.e., employing active learning strategies to develop intellectual skills, and to shape personal attitudes and values). Moreover, teachers see few incentives to change mainly because the use of educational software in classrooms is time-consuming. In fact, any faculty member who has ever attempted to lead a true one hour class discussion, in which students talk and respond to one another, knows how difficult it is to have control over the discussion.

Notwithstanding the above, the importance of using educational software in mathematics, as an efficient tool to help students grasp with hard-to-understand concepts and to more quickly gain a deeper understanding of the materials being taught firsthand, is acknowledged (see, for instance, [1,2,10,16,22]) and thereby such software can help to promote an active learning environment inside the classroom.

Although it is recognized that some economic concepts can be more easily understood when the students work with a large number and variety of graphics in an interactive way, with the support of the appropriate technology, the use of computer algebra systems is rare and under-studied in economics education (we thank an anonymous referee for this observation). In fact, the use of educational software in economics has been limited to some specific economic concepts (see, for instance, [23–25]). According to [23,24], there are automatic algebraic simplifiers, but simplicity is often in the eye of

the beholder and such tools are sparingly used by economic theorists. Furthermore, computers have already been used to generate numerical examples, providing only approximate, rather than exact, results. This gap opens a window of opportunity for the development of new educational tools directed to socio-economic science students. In a previous work [26], it was shown how some dynamic and interactive mathematical tools, created with *Mathematica*, can be used to promote better student activity and engagement in the learning process. Another work [27] discusses some teaching possibilities offered by the F-Tool concept that can provide an active learning environment in socio-economic science subjects.

The current paper intends to present a new interactive and dynamic mathematical tool for the study of the price elasticity of supply concept, the new PES(Linear)-Tool (see Supplementary Materials), which allows students to change a function's parameter values and get the analytical and graphical results in real time. Furthermore, the interactive and dynamic features of this tool make it suitable to promote an active learning environment and it is available, at no cost, in the Computable Document Format. This format allows the use of the PES(Linear)-Tool, even without an active Wolfram *Mathematica* license (additional information about how to work with the CDF format can be found at http://www.wolfram.com/cdf-player/). The potentialities of the PES(Linear)-Tool will be exhaustively explored to introduce and deal with multiple features of the price elasticity of supply, a central concept in economics. In our opinion, its use in classrooms can promote better student activity and engagement in the learning process, among students enrolled in socio-economic programs.

This paper is structured as follows. After this brief introduction section, Section 2 introduces some basic economic concepts which frame the application of the new tool. Section 3 details the F-tool concept and its application. Section 4 presents the design of the PES(Linear)-Tool. Section 5 is dedicated to some final remarks.

#### **2. Basic Economic Concepts**

This section introduces some basic concepts related to the price elasticity of supply.

#### *2.1. The Market Supply Curve and the Market Supply Function of a Good*

The producers in a given industry will supply a certain quantity of a produced good at a given price. At this price, the sum of all units gives the total market supply of that good. This corresponds to a point on a curve for the commodity. Continuously changing the price and summing individual supply across all suppliers, we can trace out the market supply curve for the good. That is, a market supply curve of a good shows the total units of that good that are supplied at different prices. More specifically, the short-run market supply curve is the horizontal summation of the individual producers' supply curves, that is:

$$Q(P) = \sum\_{i=1}^{n} q\_i(P),\tag{1}$$

where *n* represents the total number of producers in the industry and *qi*(*P*) represents the producer *i*'s supply function.

#### The Linear Case

Considering a linear specification, the market supply function can be written in the general form

$$Q(P) = \mathfrak{a}P + \beta\_{\prime} \tag{2}$$

with *<sup>α</sup>*, *<sup>β</sup>* <sup>∈</sup> <sup>R</sup>, *<sup>α</sup>* - 0 and *P* max -−*β <sup>α</sup>* , 0 . These restrictions are according to the economic theory.

In this paper, we consider the market supply inverse function (when *α* > 0), which can be expressed as

$$P(Q) = aQ + b\_{\prime} \tag{3}$$

with *a* = <sup>1</sup> *<sup>α</sup>* and *<sup>b</sup>* <sup>=</sup> <sup>−</sup>*<sup>β</sup> α* .

According to the above restrictions, *Q* max(*β*, 0), that is,

$$Q \geqslant \max\left(-\frac{b}{a'}, 0\right). \tag{4}$$

#### *2.2. Measurement and Interpretation of Price Elasticity of Supply*

The price elasticity of supply (PES) is a measure used in economics to show the responsiveness of the quantity supplied of a good or service to a change in its price. The elasticity, in a numerical form, is defined as the percentage change in the quantity supplied divided by the percentage change in price, that is,

$$PES(Q) = \lim\_{\Delta P \to 0} \frac{Percentage\ change\ in\ quantity\ supplied}{Percentage\ change\ in\ price}.\tag{5}$$

Given that we consider the linear case with *α* > 0 and *a* = 1/*α*, algebraically, the price elasticity of supply is given by the following expression:

$$PES(Q\_0) = \lim\_{\Delta P \to 0} \frac{\frac{Q\_0 - Q\_0}{Q\_0}}{\frac{P\_0 - P\_0}{P\_0}} = \lim\_{\Delta P \to 0} \frac{\frac{\triangle Q}{Q\_0}}{\frac{\triangle P}{P\_0}} = \frac{1}{\lim\_{\Delta Q \to 0} \frac{\triangle P}{\triangle Q}} \frac{P\_0}{Q\_0} \,\tag{6}$$

where *Q*<sup>0</sup> is the (positive) quantity supplied and *P* represents the price.

So, the expression for the price elasticity of supply can be expressed through the derivative of the function defined by (3) as

$$PES(Q\_0) = \frac{1}{\frac{dP}{dQ}(Q\_0)} \frac{P\_0}{Q\_0} = \frac{1}{P'(Q\_0)} \frac{P\_0}{Q\_0}.\tag{7}$$

Obviously, the price elasticity of supply takes only non-negative values. Relatively large values of the PES imply that market supply is responsive to price changes, whereas low values indicate that the supply is not very reactive to price changes.

The elasticity takes the value of zero if the quantity does not react to price changes. In this case, the supply is said to be perfectly inelastic. The elasticity takes a value between 0 and 1 if a price change causes a lower change in the quantity supplied. In this case, the supply is said to be inelastic or rigid with respect to the price. The elasticity takes the value of 1 if a price change causes identical change in the quantity supplied. In this case, the supply is said to have a unitary elasticity. Finally, the elasticity takes a value above 1 if a price change causes a higher change in the quantity supplied. In this case, the supply is said to be elastic with respect to price. The limit case occurs when the elasticity is infinite. In this case, the supply is said to be perfectly elastic.

The Linear Case

Considering a linear specification of the market supply function (2) we get the following expression:

$$PES(Q\_0) = \frac{1}{a} \frac{P\_0}{Q\_0} \,\tag{8}$$

that is,

$$PES(Q\_0) = 1 + \frac{b}{a}Q\_0^{-1} \,\text{.}\tag{9}$$

and the following situations must be considered in the design of a dynamic and interactive tool.

Perfectly elastic supply: The limit case occurs when *a* = 0 and *b* > 0. This corresponds to an infinite PES (see Figure 7).

Remark: In this case the market supply function (2) is not defined since the function (3) is not an invertible function.

Elastic supply: This case occurs when *a* > 0 and *b* > 0. This corresponds to a PES above 1 (see Figure 8).

Unit elastic supply: This case occurs when *a* > 0 and *b* = 0. This corresponds to a PES equal to 1 (see Figures 10–12).

Inelastic supply: This case occurs when *a* > 0 and *b* < 0. This corresponds to a PES below 1 (see Figure 9).

Perfectly inelastic supply: This limit case occurs when *<sup>b</sup> a Q*−<sup>1</sup> <sup>0</sup> = −1. This corresponds to a PES equal zero (see Figure 14).

Remark: In this case *α* = 0 in the market supply function (2). So, (2) is not an invertible function.

#### **3. Dynamic and Interactive Tools**

Faculty members who regularly use strategies to promote active learning typically find several ways to ensure that students learn the assigned content: promoting the dialog and reflection, promoting the acquisition of new knowledge and the transmission of the acquired knowledge, and doing short-assessments every week.

Currently, several software applications can be (free of charge or for a cost) downloaded from the World Wide Web. In particular, there are many dynamic and interactive tools dealing with some specific economic concepts implemented with the computer algebra system *Mathematica*, which is already available in the Wolfram Demonstrations Project website. In this project (http://demonstrations. wolfram.com) the creators of *Mathematica* promote and divulge globally the innovations designed by its users. Some of these applications provide only analytical information (the *Inflation-Adjusted Yield* tool, available at http://demonstrations.wolfram.com/InflationAdjustedYield/, illustrates how one's investment life planning turns on the net of nominal investment yield and inflation, according to its author). Several other tools provide only graphical information (the *Short-Run Cost Curves* tool, available at http://demonstrations.wolfram.com/ShortRunCostCurves/, provides graphical information about the cubic cost function and its average and marginal cost curves; the *Monopoly Profit and Loss* tool, available at http://demonstrations.wolfram.com/MonopolyProfitAndLoss/, provides graphical information about the marginal cost and the average cost curves). In particular, for the elasticity of demand concept there are tools that provide non-rigorous analytical information such as *The Price Elasticity of Demand* tool (available at http://demonstrations.wolfram.com/ ThePriceElasticityOfDemand/) which shows two ways to calculate the price elasticity of demand), and tools that provide only graphical information (the *Constant Price Elasticity of Demand* tool, available at http://demonstrations.wolfram.com/ThePriceElasticityOfDemand/ illustrates the price elasticity of demand for a specific inverse demand function). However, none of these applications provide all rigorous and exhaustive required information for a global and deep understanding of economic concepts introduced at undergraduate levels, in higher education institutions. Furthermore, these existing materials can hardly be adapted to explain specific concepts in socio-economic sciences or they would require additional resources from both the teacher and the students. This is a gap that the new educational tool described in this paper intends to fulfill since it is adapted to specific training programs to meet educational goals. It allows the design of tasks for independent work and the analysis of individual special cases that are important to recent graduate economists.

#### *3.1. The F-Tool Concept*

The F-Tool concept, which was first presented in the 1st National Conference on Symbolic Computation in Education and Research (Portugal 2012), where it was distinguished with the *Timberlake Award for Best Article by a Young Researcher*, was created as an interactive *Mathematica* notebook, specifically to explore the concept of real functions and their graphics, by analyzing the effects caused

by changing the values of the parameters of general analytical expressions [28]. Each F-Tool allows the study of a typical class of functions. For each class, a set of parameters is considered such that the class is fully determined by the corresponding analytical expression. This means that each F-Tool provides graphical and rigorous analytical information for all the functions within the corresponding class. In fact, unlike the other tools available in the *Wolfram Demonstrations Project* website, all the tools created under the F-Tool concept provide all the graphical and analytical information desired by the user. Additionally, the user can get exact or approximate analytical results. Finally, the new PES(Linear)-Tool has a very intuitive interface that allows even the most inexperienced user, with no previous knowledge in educational software, to start using all its features in an efficient and autonomous way.

The existing F-Tool are available, free of charge, in the Computable Document Format and the corresponding CDF files can be downloaded for free at https://sapientia.ualg.pt. This format allows anyone with a computer to fully use it, even without an active Wolfram *Mathematica* license.

The F-Tool's framework is composed by three main panels (see Figure 1):

**Figure 1.** A general example of the price elasticity of supply (PES)(Linear)-Tool: How to get the market supply function in terms of the variable *Q*.

In the left panel, the user can set the parameters' values, and choose which functions related with the main function are to be displayed.

In the middle panel, all the functions are plotted, according to the options defined in the left panel.

In the right panel, all the analytical information is displayed in accordance with the options chosen by the user in the left panel.

In summary, all the controls and options for all functionalities are located in the left panel. As the user interacts dynamically with the tool, all the graphical and analytical results are displayed in real time in the middle and right panels, respectively. When choosing the option , the user will then see the corresponding graphics moving continuously and the analytical information changing accordingly. It is through this kind of dynamic interaction that "computer algebra systems present new opportunities for teaching and learning" [29].

The use of the F-Tool concept in the classroom allows a dynamic approach to various concepts related to the study of functions and promotes new ways of reasoning/thinking, evaluating, teaching, and learning. The F-Tool concept was conceived as an active learning tool, that is, its adequate use provides a context of teaching and learning where students and teachers are both invited to fully participate [30]. Through dynamic changes of the parameters values, it is possible to obtain rigorous analytical information, presented in exact or approximate arithmetic, as well as static and non-static visual information [22]. Although it is a dynamic and interactive educational software, the F-Tool can also be used in the construction of multiple choice and open response evaluation questions [1].

#### *3.2. The F-Tool Concept Adapted to the Socio-economic Sciences*

Taking into account our experience of using dynamic and interactive mathematical tools [1] as active learning tools in natural science courses, we decided to adapt this type of approach to some economic concepts. The idea is to focus the teaching process on the students, stimulating their participation and motivating those with a level of math knowledge, often insufficient, to obtain new knowledge in a solid way. In this way, it becomes possible to teach new concepts in a solid and consistent way.

The most common way for faculty members to engage students in active learning is by stimulating the discussion [20]. A variety of materials and techniques can be used to trigger the discussion and each teacher can provide several experiences that will stimulate the discussion among students. Demonstrations during a lecture can be used to stimulate the students' curiosity and to improve their understanding of conceptual material and processes [31], particularly when the demonstration invites students to participate in research activities through the use of questions such as "What would happen if we change dynamically the parameter *b*? Would the price elasticity of supply change? And what would happen if the parameter *a* changes dynamically?" (see Figures 8–10). So, the faculty member can encourage the discussion, dialogue, and reflection in the classroom, proposing stimulating exercises that lead to a supervised constructive debate among the students.

In Section 4 we present the new dynamic and interactive economic tool, called the PES(Linear)-Tool, created under the F-Tool concept. The usefulness of this tool is illustrated by introducing the price elasticity of supply concept in a microeconomics class, as well as all the analytical and graphical information involved with the analysis of this concept.

#### **4. Designing the New PES(Linear)-Tool**

The use of the symbolic computation capabilities of *Mathematica*, and its own programming language (along with the pretty-print functionality that allows one to write mathematical expressions on the computer using the traditional notation, as on paper), enables us to implement on a computer, and in a rather straightforward manner, all the ideas that go into the F-Tool concept.

The PES(Linear)-Tool was created as an interactive *Mathematica* notebook and it is available online, in the Computable Document Format, as a supplement to this article. It allows the exploration of concepts related to a market supply function (3), where *<sup>a</sup>*, *<sup>b</sup>* <sup>∈</sup> <sup>R</sup>, *<sup>a</sup>* <sup>&</sup>gt; 0 and *<sup>Q</sup>* <sup>&</sup>gt; max - −*b <sup>a</sup>* , 0 . It should be noted that the particular case of *a* = 0 was also included to exemplify the perfectly elastic supply (when *Q* > 0) (see Figure 7). In terms of implementation and in spite of their mathematic simplicity, constant functions should be dealt with separately because they have no inverse function (see Figure 7). This means that the constant case has to be coded separately, in order to generate the correct analytical information for those functions. The PES(Linear)-Tool provides all graphical and analytical information of the inverse function of *P*(*Q*) (that is, the market supply function). As students often confuse the concepts of elasticity and derivative, the tool provides the option "Derivative" on the left panel (see Figures 10–12). The PES(Linear)-Tool displays graphical information on the value of the PES(*Q*0) whenever this option is selected. This allows the user to visualize the change from an economic model with an elastic supply to a model with an inelastic supply (going through a unitary elastic supply). As in the F-Tool, the user can interact with this information in real time.

As an illustration of this tool, let us to consider the plot of the inverse function (3) as depicted in Figure 1, and the market supply function (in terms of the variable *Q*). In this case, the exact analytical expressions of the function and its inverse are displayed, once the exact arithmetic option has been selected. The dashed line displayed on the plot is described by the equation *y* = *x* and corresponds to the symmetry axis of the inverse transformation.

The PES(Linear)-Tool is essentially created by a single Manipulate command (see Figure 2), whose output is not just a static result but a running program that we can interact with. In fact, the code consists of some initial definitions followed by the single command Manipulate. This command is responsible for creating the interactive object that contains the three panels. In particular, the command Manipulate generates all the functional controls, such as the sliders for the parameters' values and checkboxes for the plots' options. Through dynamic changes of the parameters' values, it is possible to obtain approximate or exact analytical information, as well as static and non-static visual information [28].


**Figure 2.** General code structure of the PES(Linear)-Tool.

To create the PES(Linear)-Tool we used part of the code of the educational software F-Tool. Obviously, to provide all the graphical and analytical information for the price elasticity of supply, several adaptations were performed and new fields related to this socio-economic concept were added. Figure 3 displays the code block that generates the value of the price elasticity of supply at a given quantity *<sup>Q</sup>*0. It should be noticed that the cases of *<sup>a</sup>* 0 and/or *<sup>Q</sup>*<sup>0</sup> max - −*b <sup>a</sup>* , 0 (see Figures 5, 6, 13) should be considered separately.

**Figure 3.** Code snippet of the PES(Linear)-Tool. This is part of the code that generates the analytical information about the price elasticity of supply.

#### *4.1. Parameters a and b*

In order to create a consistent tool that considers all the mathematical possibilities for which the economic model makes sense, several situations concerning the values of the parameters *a* and *b* should be implemented. Given the function (3), only non-negative values for the parameter *a* are considered in the code (the range of values that run through the slider, see Figure 4).

**Figure 4.** Code snippet of the PES(Linear)-Tool. This is part of the code that generates the range of values for the parameter *a* that run through the slider.

The user can also introduce directly the parameters' values. However, for certain values of the parameters, the correspondent market supply function is not defined and therefore, the PES(Linear)-Tool will exhibit the following message: "Does not make sense to analyze PES(Q0)!", whenever the PES button is selected (see Figure 5). Consequently, all options will be unavailable until acceptable parameter values are considered. This situation occurs when the user chooses a non-positive value for the parameter *a*, and/or the user chooses a non-positive value for the variable *Q* (see Figures 6 and 13).

**Figure 5.** An example of the information obtained when a negative value for the parameter *a* is chosen.

Depending on the values considered for the parameters *a* and/or *b*, *Q* can assume values in different numeric sets. So, the values of *Q*<sup>0</sup> that can be considered depend on the values of *a* and *b*. The PES(Linear)-Tool can be used to improve the students' understanding of this conceptualization because it enables students to analyze the relationship between the null value of the function *P*(*Q*) and the range of acceptable values for *Q*<sup>0</sup> (see Figure 6).

**Figure 6.** Example of a non economic model due the fact that *<sup>Q</sup>*<sup>0</sup> max - −*b <sup>a</sup>* , 0 .

#### *4.2. Perfectly Elastic Supply*

This subsection illustrates how the PES(Linear)-Tool can be used to improve the students' understanding of the price elasticity of supply concept (see Figure 7). In the classroom the faculty can explain that this is a limit case that occurs when the market supply function is not defined. The teacher may ask the students if a change in the parameter *b* causes any change in the type of price elasticity. Then it can be asked about the effects of a possible change in parameter *a*.

**Figure 7.** Example of a perfectly elastic supply.

#### *4.3. Elastic and Inelastic Supplies*

This subsection describes how the PES(Linear)-Tool can be used to improve the students' understanding of the price elasticity of supply concept.

By using the PES(Linear)-Tool the faculty can ask the students to interpret the value PES(*Q*0) depending on the values of *a*, *b*, and *Q*0. The faculty can start with an example of an elastic supply and ask the students to identify the parameter to be changed in order to get an inelastic supply and how the value of *Q*<sup>0</sup> affects the elasticity's value (see Figures 8–10).

**Figure 9.** Example of an inelastic supply.

#### *4.4. Unit Elastic Supply*

It is generally acknowledged that there is often a confusion between the concepts of elasticity and derivative among students. In order to illustrate the contribution of the PES(Linear)-Tool to distinguish such concepts, this subsection presents some examples of unit elastic supply functions associated to different derivatives' values. Figures 10–12 present examples of unitary elasticity supply functions associated to derivatives' values above, below and equal to 1, respectively.

**Figure 10.** Example of an unit elastic supply with a derivative value above 1.

**Figure 11.** Example of an unit elastic supply with a derivative value below 1.

**Figure 12.** Example of an unit elastic supply with a derivative value equal to 1.

The graphical and analytical information, reported by the tool, confirm that despite the existence of a relationship between the two concepts (elasticity and derivative), their values are not directly connected.

Finally, Figure <sup>13</sup> exhibits a non economic model in which *<sup>Q</sup>*<sup>0</sup> max -−*b <sup>a</sup>* , 0 (if any positive value of *Q*<sup>0</sup> is considered, the economic model would have a unitary elasticity).

**Figure 13.** Example of a non economic model due the fact that *<sup>Q</sup>*<sup>0</sup> max -−*b <sup>a</sup>* , 0 (if any positive value of *Q*<sup>0</sup> is considered the economic model would have a unitary elasticity).

#### *4.5. Perfectly Inelastic Supply*

Although the PES(Linear)-Tool cannot fully illustrate the perfectly inelastic supply case, it can be used to make this case easier for students to understand. Once the information that this limit case occurs when *<sup>b</sup> a Q*−<sup>1</sup> <sup>0</sup> = −1 has been transmitted to the students, the immediate conclusion is that the corresponding PES is zero. In this case, the teacher should state that *α*, in (2), is also null (or question why) and therefore (2) is not an invertible function. This case is depicted in Figure 14.

**Figure 14.** Example of an almost perfectly inelastic economic model.

#### **5. Final Remarks**

This paper presents a new dynamic and interactive tool created with the computer algebra system *Mathematica*, the PES(Linear)-Tool, designed to be applied in economics education, a domain where the use of computer algebra has been particularly limited. Although there are several free of charge applications available at the Wolfram Demonstrations Project website, none of those provide all the graphical and analytical information necessary for a good understanding of the concepts introduced in socio-economic undergraduate courses in universities. Moreover, these applications provide either graphical or analytical information, but not both, and/or only for some particular cases, and they can hardly be adapted to explain specific concepts in social economic sciences, or that adaptation would require additional resources from both the teacher and the students.

The above mentioned issues constitute several gaps which the new tool intends to fulfill. The PES(Linear)-Tool is a computer algebra tool directed to the study of one of the most widely used concepts in socio-economics courses—the price elasticity of supply. Starting with the specification of the market supply function, the design, functionalities and capabilities of the PES(Linear)-Tool are exhaustively explored in this paper to analyze the price elasticity of supply, accounting for all the mathematical possibilities for which the economic model makes sense. This tool also differs from other existing tools in that it can be downloaded at no cost and allows the complete analysis of multiple situations involving the study of the price elasticity supply in a dynamic and interactive way. Specifically, the tool offers the students the possibility of changing the parameters' values in the economic model and getting both the analytical and graphical effects in real time. The new PES(Linear)-Tool has a very intuitive interface that allows even the most inexperienced user, with

no previous knowledge in educational software, to start using all the features in an efficient and autonomous way.

Given the recognition in the literature that some economic concepts can be more easily understood when students work with a large number and variety of graphics in an interactive way, with the support of the appropriate technology, we believe that the use of the PES(Linear)-Tool in the classroom can promote new ways of reasoning/thinking, evaluating, teaching, and learning in a context where students and teachers are invited to contribute. In this way, this tool promotes the active learning in classrooms and simultaneously students' autonomous work, by allowing the design of challenging problems based on dynamic and interactive exercises using the CDF format, which students can work on and then send in their results to the faculty by email.

The design of the PES(Linear)-Tool can be generalized to other economic models that can be studied through other classes of functions, and also opens the possibility for the development of other interactive tools associated with other economic concepts.

Going forward, a statistically rigorous study in loco to assess the students' perception when using the PES(Linear)-Tool, and therefore to estimate the tool's pedagogical value is of extremely importance. We believe that this study can be an important help for the future development of these kind of educational tools.

**Supplementary Materials:** The PES(Linear)-Tool is available online at http://www.mdpi.com/2297-8747/24/4/ 87/s1.

**Author Contributions:** The dynamic and interactive tool presented in this paper was designed by A.C.C. (under the F-Tool concept created by A.C.C., José C. Pereira, Cátia M. Silva, and Cristina R. Simão). All authors contributed to the improvement of the tool. The implementation with the computer algebra system *Mathematica* was made by R.C. and A.C.C. The conceptualization and methodology was performed by J.M.A. and A.C.C. The paper was written by J.M.A. and A.C.C. All authors reviewed the manuscript.

**Funding:** This research was funded by Fundação para a CiênciaeaTecnologia within the project UID/ECO/04007/2019.

**Acknowledgments:** The authors thank the contributions and suggestions of two anonymous referees.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **The Invariant Two-Parameter Function of Algebras** *ψ***¯**

#### **José María Escobar 1, Juan Núñez-Valdés 1,\* and Pedro Pérez-Fernández 2,3**


Received: 12 September 2019; Accepted: 11 October 2019; Published: 14 October 2019

**Abstract:** At present, the research on invariant functions for algebras is very extended since Hrivnák and Novotný defined in 2007 the invariant functions *ψ* and *ϕ* as a tool to study the Inönü–Wigner contractions (IW-contractions), previously introduced by those authors in 1953. In this paper, we introduce a new invariant two-parameter function of algebras, which we call *ψ*¯, as a tool which makes easier the computations and allows researchers to deal with contractions of algebras. Our study of this new function is mainly focused in Malcev algebras of the type Lie, although it can also be used with any other types of algebras. The main goal of the paper is to prove, by means of this function, that the five-dimensional classical-mechanical model built upon certain types of five-dimensional Lie algebras cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra. As an example of other applications of the new function obtained, its computation in the case of the Lie algebra induced by the Lorentz group *SO*(3, 1) is shown and some open physical problems related to contractions are also formulated.

**Keywords:** invariant functions; contractions of algebras; Lie algebras; Malcev algebras; Heisenberg algebras

#### **1. Introduction**

Regarding the concept of *limit process* between physical theories in terms of contractions of their associated symmetry groups, formulated by Erdal Inönü and Eugene Wigner [1,2], these authors introduced the so-called *Inönü–Wigner contractions (IW-contractions)* in 1953. Later, other extensions of these IW-contractions have also been addressed, for instance the *generalized Inönü–Wigner contractions*, introduced by Melsheiner [3], the *parametric degenerations* [4–6], widely used in the Algebraic Invariants Theory, and the *singular contractions* [2]. To study these contractions, Hrivnák and Novotný introduced the invariant functions *ψ* and *ϕ* as a tool in 2007 [7]. These invariant functions depend on one parameter.

Continuing with this topic, the main goal of this paper is to introduce a new invariant function, in this case depending on two parameters, which we call *the two-parameter invariant function ψ*¯, to get some advances on this research. Indeed, the objective is to prove, by means of this function, that the five-dimensional classical-mechanical model built upon certain types of five-dimensional Lie algebras cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra.

Indeed, the study of this function is mainly focused in the frame of the Malcev algebras of the type Lie. Thus, this paper can be considered as the natural continuation of a previous one dealing with Lie algebras [8]. We try to generalize the properties obtained on that to the case of Malcev algebras.

The structure of the paper is as follows. In Section 2, we recall some preliminaries on the mathematical objects dealt with in this paper, Lie algebras and Malcev algebras. Section 3 is devoted to introducing and proving the main properties of the two-parameter invariant function *ψ*¯. For computations, we used the SAGE symbolic computation package and in this section we prove that this new function is different from others previously defined, which are used as a tool to study contractions of algebras. We also prove the main result of the paper: no proper contraction between a fifth Heisenberg algebra and a filiform Lie algebra of dimension 5 exists. It implies that the five-dimensional classical-mechanical model built upon a five-dimensional filiform Lie algebra cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra. In this way, the new function allows us to step forward in the research on contractions. In Section 4, we show some of our discussion and conclusions regarding the research done. Finally, in Section 5, we give some comments on the materials and methods used in such a research.

#### **2. Preliminaries**

We show in this section some preliminaries on Lie algebras, Malcev algebras and on Heisenberg algebras, which are the main mathematical objects used in the paper.

#### *2.1. Preliminaries on Lie Algebras*

In this subsection, we show some preliminaries on Lie algebras. For a further review on this topic, the reader can consult [9].

An *<sup>n</sup>*-dimensional *Lie algebra* g over a field *<sup>K</sup>* is an *<sup>n</sup>*-dimensional vector space over *<sup>K</sup>* endowed with a second inner law, named *bracket product*, which is bilinear and anti-commutative and satisfies the *Jacobi identity*

$$f(u, \upsilon, w) = [u, [\upsilon, w]] + [\upsilon, [w, u]] + [w, [u, \upsilon]] = 0,\text{ for all } u, \upsilon, w \in \mathfrak{g}.\tag{1}$$

The law of the *<sup>n</sup>*-dimensional Lie algebra g is determined by the products

$$[e\_{i\prime}e\_j] = \sum\_{k=1}^n c\_{ij}^k e\_{k\prime} \quad \text{for} \quad 1 \le i < j \le n\_{\prime}$$

where *c<sup>k</sup> <sup>i</sup>*,*<sup>j</sup>* <sup>∈</sup> *<sup>K</sup>* are called *structure constants* of <sup>g</sup>. If all these constants are zero, then the Lie algebra is called *abelian*.

Two Lie algebras g and h are *isomorphic* if there exists a vector space isomorphism *<sup>f</sup>* between them such that *<sup>f</sup>*([*u*, *<sup>v</sup>*]) = [ *<sup>f</sup>*(*u*), *<sup>f</sup>*(*v*)], for all *<sup>u</sup>*, *<sup>v</sup>* <sup>∈</sup> g.

A mapping *<sup>d</sup>* : g −→ g is a *derivation* of g if *<sup>d</sup>*([*u*, *<sup>v</sup>*]) = [*d*(*u*), *<sup>v</sup>*]+[*u*, *<sup>d</sup>*(*v*)], for all *<sup>u</sup>*, *<sup>v</sup>* <sup>∈</sup> g. The set of derivations of g is denoted by *Der*g.

The *lower central series* of a Lie algebra g is defined as g<sup>1</sup> <sup>=</sup> <sup>g</sup>, <sup>g</sup><sup>2</sup> = [g1, <sup>g</sup>], ..., <sup>g</sup>*<sup>k</sup>* = [g*k*−1, <sup>g</sup>], ...

If there exists *<sup>m</sup>* <sup>∈</sup> <sup>N</sup> such that g*<sup>m</sup>* <sup>≡</sup> 0, then g is called *nilpotent*. The *nilpotency class* of g is the smallest natural *<sup>c</sup>* such that g*c*+<sup>1</sup> <sup>≡</sup> 0.

An *<sup>n</sup>*-dimensional nilpotent Lie algebra g is said to be *filiform* if it is verified that dim g*<sup>k</sup>* <sup>=</sup> *<sup>n</sup>* <sup>−</sup> *k*, for all *k* ∈ {2, ... , *n*}. Filiform Lie algebras were introduced by Vergne in her Ph.D. Thesis, in 1966, later published in [10] in 1970.

The only *n*-dimensional filiform Lie algebra for *n* < 3 is the abelian. For *n* ≥ 3, it is always possible to find an *adapted basis* {*e*1, ... ,*en*} of g such that [*e*1,*e*2] = 0, [*e*1,*ej*] = *ej*−1, for all *<sup>j</sup>* ∈ {3, ... , *<sup>n</sup>*} and [*e*2,*ej*]=[*e*3,*ej*] = 0, for all *j* ∈ {3, . . . , *n*}.

From the condition of filiformity and the Jacobi identity in Equation (1), the bracket product of g is determined by

$$[\mathfrak{e}\_{i\prime}, \mathfrak{e}\_{j}] = \sum\_{k=2}^{\min\{i-1, n-2\}} \mathfrak{e}\_{ij}^{k} \mathfrak{e}\_{k\prime} \qquad \text{for} \qquad 4 \le i < j \le n\prime$$

where *c<sup>k</sup> <sup>i</sup>*,*<sup>j</sup>* <sup>∈</sup> *<sup>K</sup>* are called *structure constants* of <sup>g</sup>. If all these constants are zero, then the filiform Lie algebra <sup>g</sup> is called *model*. The model algebra is not isomorphic to any other algebra of the same dimension and every *<sup>n</sup>*-dimensional filiform Lie algebra <sup>g</sup> having an adapted basis {*e*1, ... ,*en*} verifies that g<sup>2</sup> <sup>=</sup> *<sup>e</sup>*2,...,*en*−<sup>1</sup>, g<sup>3</sup> <sup>=</sup> *<sup>e</sup>*2,...,*en*−<sup>2</sup>,..., g*n*−<sup>1</sup> <sup>=</sup> *<sup>e</sup>*<sup>2</sup>, g*<sup>n</sup>* <sup>=</sup> 0.

#### *2.2. Preliminaries on Malcev Algebras*

Now, we recall some preliminary concepts on Malcev algebras, taking into account that a general overview can be consulted in [11]. From here on, we only consider finite-dimensional Malcev algebras over the complex number field C.

A *Malcev algebra* M is a vector space with a second bilinear inner composition law ([·, ·]) called the *bracket product* or *commutator*, which satisfies: (a) [*u*, *v*] = −[*v*, *u*], ∀*u*, *v* ∈ M; and (b) [[*u*, *v*], [*u*, *w*]] = [[[*u*, *v*], *w*], *u*] + [[[*v*, *w*], *u*], *u*] + [[[*w*, *u*], *u*], *v*], ∀*u*, *v*, *w* ∈ M. Condition (b) is named *Malcev identity* and we use the notation *M*(*u*, *v*, *w*) = [[*u*, *v*], [*u*, *w*]] − [[[*u*, *v*], *w*], *u*] − [[[*v*, *w*], *u*], *u*] − [[[*w*, *u*], *u*], *v*].

Given a basis {*ei*}*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> of a *<sup>n</sup>*-dimensional Malcev algebra M, the *structure constants <sup>c</sup><sup>h</sup> <sup>i</sup>*,*<sup>j</sup>* are defined as [*ei*,*ej*] = ∑*<sup>n</sup> <sup>h</sup>*=<sup>1</sup> *<sup>c</sup><sup>h</sup> i*,*j eh*, for 1 ≤ *i*, *j* ≤ *n*.

It is immediate to see that Malcev algebras and Lie algebras are not disjoint sets. Indeed, every Lie algebra is a Malcev algebra, but the converse is not true. Therefore, we can distinguish between Malcev algebras of the type Lie and Malcev algebras of the type non-Lie. Obviously, those Malcev algebras which are of the type Lie verify both identities: Jacobi and Malcev.

If the Jacobi identity does not hold, then the Malcev algebra is said to have a *Jacobi anomaly*. In quantum mechanics, the existence of Jacobi anomalies in the underlying non-associative algebraic structure related to the coordinates and momenta of a quantum non-Hamiltonian dissipative system was already claimed by Dirac [12] in the process of taking Poisson brackets. In string theory, for instance, one such anomaly is involved by the non-associative algebraic structure that is defined by coordinates (*x*) and velocities or momenta (*v*) of an electron moving in the field of a constant magnetic charge distribution, at the position of the location of the magnetic monopole [13]. In particular, *<sup>J</sup>*(*v*1, *<sup>v</sup>*2, *<sup>v</sup>*3) = −∇ ◦ *<sup>B</sup>*(*x*), where ∇ ◦ *B*(*x*) denotes the divergence of the magnetic field *B*(*x*). The underlying algebraic structure constitutes a non-Lie Malcev algebra [14], with the commutation relations [*xa*, *xb*] = 0, [*xa*, *vb*] = *i δab* and [*va*, *vb*] = *i εabc Bc*(*x*), where *a*, *b*, *c* ∈ {1, 2, 3}, *δab* denotes the Kronecker delta and *εabc* denotes the Levi–Civita symbol. If the magnetic field is proportional to the coordinates, the latter can be normalized and *Bc*(*x*) can then be supposed to coincide with *xc*. The resulting algebra is then called magnetic [15]. A generalization to electric charges has recently been considered [15] by defining the products [*xa*, *xb*] = <sup>−</sup>*<sup>i</sup> <sup>ε</sup>abc Ec*(*x*,*v*), where the electric field *E* as well as the magnetic field *B* must depend not only on coordinates but also on velocities. It is worth remarking that both magnetic and electric algebras constitute magma algebras (see [16] for this last concept).

If g is a Malcev algebra of the type Lie and *<sup>D</sup>* <sup>∈</sup> *Der*g a derivation of g, then, according to the anti-commutative property of g and the Jacobi identity in Equation (1) of Lie algebras, we get that

$$d\left[d[\mathbf{x},\mathbf{y}]\_\prime,[\mathbf{x},\mathbf{z}]\right] + \left[[\mathbf{x},\mathbf{y}]\_\prime,d[\mathbf{x},\mathbf{z}]\right] = d\left[\left[[\mathbf{x},\mathbf{z}]\_\prime,\mathbf{y}\right]\_\prime,\mathbf{x}\right] + d\left[\left[[\mathbf{z},\mathbf{x}]\_\prime,\mathbf{x}\right]\_\prime,\mathbf{y}\right] \quad \forall \mathbf{x},\mathbf{y},\mathbf{z} \in \mathfrak{g}$$

Starting from here and due to reasons of length, only Malcev algebras of type Lie, that is to say, actually Lie algebras, are used in this paper. Malcev algebras of type non-Lie will be dealt with in future work.

#### *2.3. Preliminaries on Heisenberg Algebras*

Let *n* be a non-negative integer or infinity. The *n*th Heisenberg algebra (so-called after Werner Karl Heisenberg) is the Lie algebra with basis B = {*p*1, ... , *pn*, *q*1, ... , *qn*, *z*} with the following relations, known as *canonical commutation relations*

1. [*pi*, *qj*] = *cij z*, 1 ≤ *i*, *j* ≤ *n*.

2. [*pi*, *z*]=[*qi*, *z*]=[*pi*, *pj*]=[*qi*, *qj*] = 0, 1 ≤ *i*, *j* ≤ *n*.

Note that the dimension of an *n*th Heisenberg algebra is not *n*, but 2*n* + 1. In fact, the *n* in the above definition is called the *rank* of the Heisenberg algebra, although it is not, however, a rank in any of the usual meanings that this word has in the theory of Lie algebras. Thus, this Lie algebra is also known as the *Heisenberg algebra of rank n*.

In any case, from here on and to avoid confusions we designate under the notation fifth Heisenberg algebras to those Heisenberg algebras generated by five generators.

#### **3. Results**

In this section, which is divided by subheadings, we provide a concise and precise description of our experimental results. They are the following.

#### *3.1. Introducing a New Invariant Function*

Let g = (*V*, [ , ]) be a Lie algebra. *End* <sup>g</sup> denotes the vector space of all linear operators of <sup>g</sup> over *<sup>V</sup>*.

**Definition 1.** *Let* g *be a Lie algebra. The set*

$$\operatorname{Der}\_{(a,\emptyset,\gamma,\tau)}\mathfrak{g} = \{d \in \operatorname{End}\mathfrak{g} : a[d[\mathbf{x},y],[\mathbf{x},z]] + \beta[[\mathbf{x},y],d[\mathbf{x},z]] = \gamma d[[[\mathbf{x},z],y]\mathbf{x}] + \tau d[[[[z,\mathbf{x}],\mathbf{x}],y]]\}$$

<sup>∀</sup>(*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup>4*, is called* the set of the (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*)-derivations *of the algebra* g. *It is denoted by Der*(*α*,*β*,*γ*,*τ*)g.

It is obvious that dim *Der*(1,1,1,1)g = dim *Der*g . Then, as dim *Der*g is an invariant of g, it follows that dim *Der*(1,1,1,1)g is an invariant of g. This leads the following result.

**Proposition 1.** *If* <sup>g</sup> *is a Lie algebra, then dim*(1,1,1,1)g *is an algebraic invariant of* g*.*

**Theorem 1.** *Let* g *and* g¯ *be two Malcev algebras of the type Lie and let <sup>f</sup>* : g <sup>→</sup> g¯ *be an isomorphism. Then, the mapping <sup>ρ</sup>* : *End* g <sup>→</sup> *End* g¯, *defined by <sup>D</sup>* −→ *fDf* <sup>−</sup>1*, is an isomorphism between the vector spaces Der*(*α*,*β*,*γ*,*τ*)g *and Der*(*α*,*β*,*γ*,*τ*)g¯, <sup>∀</sup>(*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup>4.*,*

**Proof.** Let g = (*V*, ·) and g¯ = (*V*¯ , <sup>∗</sup>) be two Malcev algebras of the type Lie and let us consider *<sup>D</sup>* <sup>∈</sup> *Der*(*α*,*β*,*γ*,*τ*)g, for any (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup><sup>4</sup> and for all *<sup>x</sup>*, *<sup>y</sup>*, *<sup>z</sup>* <sup>∈</sup> g¯. Then,

$$\begin{split} &\alpha D\Big(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{y})\Big)\cdot\Big(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\Big) + \beta \Big(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{y})\Big)\cdot D\Big(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\Big) = \\ &\gamma D\Big(\Big(\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\right)\cdot f^{-1}(\mathbf{y})\Big)\cdot f^{-1}(\mathbf{x})\Big) + \tau D\Big(\Big(\left(f^{-1}(\mathbf{z})\cdot f^{-1}(\mathbf{x})\right)\cdot f^{-1}(\mathbf{x})\Big)\cdot f^{-1}(\mathbf{y})\Big). \end{split}$$

It is deduced that

$$\gamma D\left(\left(\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\right)\cdot f^{-1}(\mathbf{y})\right)\cdot f^{-1}(\mathbf{x})\right) = \gamma D\left(\left(f^{-1}(\mathbf{x}\ast\mathbf{z})\cdot f^{-1}(\mathbf{y})\right)\cdot f^{-1}(\mathbf{x})\right) = 0$$

$$\gamma Df^{-1}((\mathbf{x}\ast\mathbf{z})\ast g)\cdot f^{-1}(\mathbf{x})) = \gamma Df^{-1}(((\mathbf{x}\ast\mathbf{z})\ast g)\ast \mathbf{x}),$$

and, similarly,

$$
\tau D\left(\left(\left(f^{-1}(\mathbf{z})\cdot f^{-1}(\mathbf{x})\right)\cdot f^{-1}(\mathbf{x})\right)\cdot f^{-1}(\mathbf{y})\right) = \tau Df^{-1}(\left(\left(\mathbf{z}\ast\mathbf{x}\right)\ast\mathbf{y}\right)\cdot\mathbf{z})
$$

$$
aD\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{y})\right)\cdot\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\right) = aDf^{-1}(\mathbf{x}\ast\mathbf{y})\cdot f^{-1}(\mathbf{x}\ast\mathbf{z})
$$

$$
\beta\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{y})\right)\cdot D\left(f^{-1}(\mathbf{x})\cdot f^{-1}(\mathbf{z})\right) = \beta f^{-1}(\mathbf{x}\ast\mathbf{y})\cdot Df^{-1}(\mathbf{x}\ast\mathbf{z}).
$$

Thus,

*<sup>α</sup>D f* <sup>−</sup>1(*<sup>x</sup>* <sup>∗</sup> *<sup>y</sup>*)· *<sup>f</sup>* <sup>−</sup>1(*<sup>x</sup>* <sup>∗</sup> *<sup>z</sup>*) + *<sup>β</sup> <sup>f</sup>* <sup>−</sup>1(*<sup>x</sup>* <sup>∗</sup> *<sup>y</sup>*)· *D f* <sup>−</sup>1(*<sup>x</sup>* <sup>∗</sup> *<sup>z</sup>*) = *<sup>γ</sup>D f* <sup>−</sup>1(((*<sup>x</sup>* <sup>∗</sup> *<sup>z</sup>*) <sup>∗</sup> *<sup>y</sup>*) <sup>∗</sup> *<sup>x</sup>*) + *<sup>τ</sup>D f* <sup>−</sup>1(((*<sup>z</sup>* <sup>∗</sup> *<sup>x</sup>*) <sup>∗</sup> *<sup>x</sup>*) <sup>∗</sup> *<sup>y</sup>*).

Now, the result of applying *f* to the previous expression is

$$\mu \left( f \mathcal{D} f^{-1} \right) (\mathbf{x} \ast y) \ast (\mathbf{x} \ast z) + \beta (\mathbf{x} \ast y) \ast \left( f \mathcal{D} f^{-1} \right) (\mathbf{x} \ast z) = \gamma \left( f \mathcal{D} f^{-1} \right) (((\mathbf{x} \ast z) \ast y) \ast \mathbf{x}) + \tau \left( f \mathcal{D} f^{-1} \right) (((\mathbf{z} \ast \mathbf{x}) \ast \mathbf{x}) \ast y) \ast \mathbf{x}$$

Thus, *fDf* <sup>−</sup><sup>1</sup> <sup>∈</sup> *Der*(*α*,*β*,*γ*,*τ*)g¯, which concludes the proof.

An immediate consequence of this result is the following.

**Corollary 1.** *Let* <sup>g</sup> *be a Lie algebra. The dimension of the vector space Der*(*α*,*β*,*γ*,*τ*)<sup>g</sup> *is an invariant of the algebra, for all* (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup>4.

**Lemma 1.** *(Technical Lemma) Let d be a derivation of a Lie algebra* g. *The following expressions are verified*


**Proof.** All expressions are immediate consequences of the properties of the derivations (see Section 2).

**Lemma 2.** *Let* g = (*V*, [, ]) *be a Lie algebra. Then,*

$$Der\_{(\mathfrak{a},\mathfrak{b},\gamma,\tau)}\mathfrak{g} = Der\_{(\mathfrak{a}+\mathfrak{f},\mathfrak{a}+\mathfrak{f},2\gamma,2\tau)}\mathfrak{g} \cap Der\_{(\mathfrak{a}-\mathfrak{f},\mathfrak{b}-\mathfrak{a},0,0)}\mathfrak{g}$$

**Proof.** Suppose *<sup>D</sup>* <sup>∈</sup> *Der*(*α*,*β*,*γ*,*τ*)g. Then, for all (*x*, *<sup>y</sup>*, *<sup>z</sup>*) <sup>∈</sup> g, we have

$$d[d[\mathbf{x}, y], [\mathbf{x}, z]] + \beta[[\mathbf{x}, y], d[\mathbf{x}, z]] = \gamma d[[[\mathbf{x}, z], y], \mathbf{x}] + \tau d[[[z, \mathbf{x}], \mathbf{x}], y].$$

Charging now *y* and *z* between themselves, we have

$$d a[d[\mathbf{x}, z], [\mathbf{x}, y]] + \beta[[\mathbf{x}, z], d[\mathbf{x}, y]] = \gamma d[[[\mathbf{x}, y], z], \mathbf{x}] + \tau d[[[y, \mathbf{x}], \mathbf{x}], z]$$

and by adding the two first expressions of Lemma 1 and taking the anti-skew property of the Lie bracket into consideration, we have

$$\begin{aligned} &(\boldsymbol{a}-\boldsymbol{\beta})[d[\boldsymbol{x},\boldsymbol{y}],[\boldsymbol{x},\boldsymbol{z}]]+(\boldsymbol{\beta}-\boldsymbol{\alpha})[[\boldsymbol{x},\boldsymbol{y}],d[\boldsymbol{x},\boldsymbol{z}]] = \\ &\boldsymbol{\gamma}(d[[[\boldsymbol{x},\boldsymbol{z}],\boldsymbol{y}]\boldsymbol{x}]+d[[[\boldsymbol{x},\boldsymbol{y}],\boldsymbol{z}]\boldsymbol{x}])+\boldsymbol{\tau}(d[[[\boldsymbol{z},\boldsymbol{x}],\boldsymbol{x}],\boldsymbol{y}]+d[[[[\boldsymbol{y},\boldsymbol{x}],\boldsymbol{x}],\boldsymbol{z}])).\end{aligned}$$

Similarly, starting from the two last expressions of Lemma 1, we obtain

$$d\left[\left[\left[z,x\right],x\right],y\right] + d\left[\left[\left[y,x\right],x\right],z\right] = 0$$

and by repeating the same procedure we obtain

$$d\left[\left[\left[x,z\right],y\right],x\right] + d\left[\left[\left[x,y\right],z\right],x\right] = 0.$$

Now, starting from both expressions, we have

$$(\alpha - \beta)[d[\mathbf{x}, y], [\mathbf{x}, z]] + (\beta - \alpha)[[\mathbf{x}, y], d[\mathbf{x}, z]] = 0.$$

Therefore, *<sup>D</sup>* <sup>∈</sup> *Der*(*α*−*β*,*β*−*α*,0,0)g.

Now, by subtracting the two first expressions of the proof and taking into account the anti-skew property, we have (*α* + *β*)[*d*[*x*, *y*], [*x*, *z*]] + (*β* + *α*)[[*x*, *y*], *d*[*x*, *z*]] = *γ*(*d*[[[*x*, *z*], *y*], *x*] − *d*[[[*x*, *y*], *z*], *x*]) + *τ*(*d*[[[*z*, *x*], *x*], *y*] − *d*[[[*y*, *x*], *x*], *z*]).

We use now in the previous equality the two expressions *d*[[[*y*, *x*], *x*], *z*] = −*d*[[[*z*, *x*], *x*], *y*] and *d*[[[*x*, *y*], *z*], *x*] = −*d*[[[*x*, *z*], *y*], *x*], respectively, obtained from previous expressions.

We have that (*α* + *β*)[*d*[*x*, *y*], [*x*, *z*]] + (*α* + *β*)[[*x*, *y*], *d*[*x*, *z*]] = 2*γd*[[[*x*, *z*], *y*], *x*] + 2*τd*[[[*z*, *x*], *x*], *y*]. It involves that *<sup>D</sup>* <sup>∈</sup> *Der*(*α*+*β*,*α*+*β*,2*γ*,2*τ*)g. Therefore, it is verified that *Der*(*α*,*β*,*γ*,*τ*)<sup>g</sup> <sup>⊂</sup> *Der*(*α*+*β*,*α*+*β*,2*γ*,2*τ*)<sup>g</sup> <sup>∩</sup> *Der*(*α*−*β*,*β*−*α*,0,0)g.

If *<sup>D</sup>* <sup>∈</sup> *Der*(*α*+*β*,*α*+*β*,2*γ*,2*τ*)g <sup>∩</sup> *Der*(*α*−*β*,*β*−*α*,0,0)g, then *<sup>D</sup>* verifies both equations (*<sup>α</sup>* <sup>+</sup> *<sup>β</sup>*)[*d*[*x*, *<sup>y</sup>*], [*x*, *<sup>z</sup>*]] + (*α* + *β*)[[*x*, *y*], *d*[*x*, *z*]] = 2*γd*[[[*x*, *z*], *y*], *x*] + 2*τd*[[[*z*, *x*], *x*], *y*] and (*α* − *β*)[*d*[*x*, *y*], [*x*, *z*]] + (*β* − *α*)[[*x*, *y*], *d*[*x*, *z*]] = 0.

Then, by adding these last equations and simplifying, we observe that *D* verifies

$$d[d[\mathbf{x}, y], [\mathbf{x}, z]] + \beta[[\mathbf{x}, y], d[\mathbf{x}, z]] = \gamma d[[[\mathbf{x}, z], y], \mathbf{x}] + \tau d[[[z, \mathbf{x}], \mathbf{x}], y].$$

Thus, *<sup>D</sup>* <sup>∈</sup> *Der*(*α*,*β*,*γ*,*τ*)<sup>g</sup> <sup>=</sup> *Der*(*α*+*β*,*α*+*β*,2*γ*,2*τ*)<sup>g</sup> <sup>∩</sup> *Der*(*α*−*β*,*β*−*α*,0,0)g. Therefore, *Der*(*α*,*β*,*γ*,*τ*)<sup>g</sup> <sup>=</sup> *Der*(*α*+*β*,*α*+*β*,2*γ*,2*τ*)<sup>g</sup> <sup>∩</sup> *Der*(*α*−*β*,*β*−*α*,0,0)g, which completes the proof.

**Theorem 2.** *Let* g *be a Lie algebra. Then, for all* (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup>4, *it exists* (*λ*1, *<sup>λ</sup>*2) <sup>∈</sup> <sup>C</sup><sup>2</sup> *such that Der*(*α*,*β*,*γ*,*τ*)g <sup>⊂</sup> <sup>C</sup><sup>2</sup> *is one of the following four sets: Der*(0,0,*λ*1,*λ*2)g; *Der*(1,−1,*λ*1,*λ*2)g; *Der*(1,0,*λ*1,*λ*2)g; *or Der*(1,1,*λ*1,*λ*2)g.

**Proof.** Consider (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*) <sup>∈</sup> <sup>C</sup>4. We distinguish the following cases

Case 1: *α* + *β* = 0. We distinguish now the following two subcases:


$$\operatorname{Der}\_{(a,\mathfrak{E},\gamma,\tau)}\mathfrak{g} = \operatorname{Der}\_{(0,0,2\gamma,2\tau)}\mathfrak{g} \cap \operatorname{Der}\_{(-2\beta,2\beta,0,0)}\mathfrak{g} = \operatorname{Der}\_{(0,0,\gamma,\tau)}\mathfrak{g} \cap \operatorname{Der}\_{(-1,1,0,0)}\mathfrak{g}.$$

Apart from that, it is also verified that

$$\operatorname{Der}\_{(-1,1,\gamma,\tau)}\mathfrak{g} = \operatorname{Der}\_{(0,0,2\gamma,2\tau)}\mathfrak{g} \cap \operatorname{Der}\_{(-2,2,0,0)}\mathfrak{g} = \operatorname{Der}\_{(0,0,\gamma,\tau)}\mathfrak{g} \cap \operatorname{Der}\_{(-1,1,0,0)}\mathfrak{g}.$$

Therefore, *Der*(*α*,*β*,*γ*,*τ*)<sup>g</sup> <sup>=</sup> *Der*(−1,1,*γ*,*τ*)g. It involves that *<sup>λ</sup>*<sup>1</sup> <sup>=</sup> *<sup>γ</sup>* and *<sup>λ</sup>*<sup>2</sup> <sup>=</sup> *<sup>τ</sup>*.

Case 2: *α* + *β* = 0. Two subcases are also considered:


$$\mathbb{D}$$

These two two-parameter sets *Der*(1,0,*λ*1,*λ*2)g and *Der*(1,1,*λ*1,*λ*2)g previously defined allow us to define the following invariant two-parameter functions of Lie algebras.

**Definition 2.** *The functions <sup>ψ</sup>*¯g, *<sup>ψ</sup>*¯0 g : <sup>C</sup><sup>2</sup> → <sup>N</sup> *defined, respectively, as* (*ψ*¯g)(*α*, *<sup>β</sup>*) = dim *Der*(1,1,*α*,*β*)<sup>g</sup> *and* (*ψ*¯0 g)(*α*, *<sup>β</sup>*) = dim *Der*(1,0,*α*,*β*)<sup>g</sup> *are called <sup>ψ</sup>*¯<sup>g</sup> *and <sup>ψ</sup>*¯0 g invariant functions *corresponding to the* (*α*, *<sup>β</sup>*, *<sup>γ</sup>*, *<sup>τ</sup>*)*-derivations of* g*.*

**Corollary 2.** *If two Malcev algebras of the type Lie* <sup>g</sup> *and* <sup>f</sup> *are isormorphic, then <sup>ψ</sup>*¯<sup>g</sup> <sup>=</sup> *<sup>ψ</sup>*¯f *and <sup>ψ</sup>*¯0 g <sup>=</sup> *<sup>ψ</sup>*¯0 f .

Note that the function *ψ*¯ is a two-parameter function, whereas the function *ψ* by Novotný and Hrivnák [7] is one-parameter. It implies that both functions are structurally different. However, it can be thought that *ψ* could be obtained as a particular case of *ψ*¯ by simply taking one of the parameters as a constant. The following counter-example shows that it is not possible.

Indeed, we now compare the function *ψ*¯ with the invariant function *ψ* and prove that both functions are totally different. To do this, we compute both functions for a same Lie algebra, in the particular case of being *α* = 1. Concretely, we use the Lie algebra induced by the Lorentz group *SO*(3, 1), which we denote by g6.

## **Computing** *<sup>ψ</sup>*g<sup>6</sup> , **for** *<sup>α</sup>* <sup>=</sup> <sup>1</sup>

Let us recall that Minkowski defined the spacetime as a four-dimensional manifold with the metric *ds*<sup>2</sup> = −*c*2*dt*<sup>2</sup> + *dx*<sup>2</sup> + *dy*<sup>2</sup> + *dz*2. We introduce the metric tensor

$$
\eta = \begin{bmatrix} -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}.
$$

If we rename (*ct*, *<sup>x</sup>*, *<sup>y</sup>*, *<sup>z</sup>*) <sup>→</sup> (*x*0, *<sup>x</sup>*1, *<sup>x</sup>*2, *<sup>x</sup>*3), then the expression *ds*<sup>2</sup> can be written as *ds*<sup>2</sup> <sup>=</sup> *ημγdxμdx<sup>γ</sup>* (summed over *μ* and *γ*). Recall that this distance is invariant under the following type of transformations *<sup>x</sup><sup>μ</sup>* <sup>→</sup> *<sup>λ</sup><sup>μ</sup> <sup>γ</sup>x<sup>γ</sup>* such that the coefficients *<sup>λ</sup><sup>μ</sup> <sup>γ</sup>* are the elements of a matrix Λ (which is called *Lorentz transformations*) that satisfies Λ*<sup>t</sup> η*Λ = *η*. Since the metric in the three-dimensional Euclidean space corresponds to the identity matrix, if *R* is the matrix of a rotation, then *R<sup>t</sup> 1R* = 1 and comparing this expression with Λ*<sup>t</sup> η*Λ = *η* it is possible to say that the Lorentz transformations are rotations in the Minkowski space. These transformations form a group called the *Lorentz group SO*(3, 1).

Now, we focus our study on the *infinitesimal* Lorentz transformations. A Lorentz transformation matrix can be written as Λ*<sup>μ</sup> <sup>γ</sup>* = *δ μ <sup>γ</sup>* <sup>+</sup> *<sup>λ</sup><sup>μ</sup> <sup>γ</sup>*, where the parameters *<sup>λ</sup><sup>μ</sup> <sup>γ</sup>* are infinitesimal and verify that *<sup>λ</sup><sup>μ</sup> <sup>γ</sup>* <sup>=</sup> <sup>−</sup>*λ<sup>γ</sup> μ* so that the Lorentz transformation is valid. The action of this transformation on the coordinates *x<sup>μ</sup>* in the Minkowski space can be written as *δx<sup>μ</sup>* = Λ*<sup>μ</sup> <sup>γ</sup>xγ*.

If we define *<sup>A</sup>ρσ* such that <sup>Λ</sup>*<sup>μ</sup> <sup>γ</sup>* <sup>=</sup> <sup>1</sup> 2 *λρσ*(*Aρσ*) *μ <sup>γ</sup>*, we can write the above action as *<sup>δ</sup>x<sup>μ</sup>* <sup>=</sup> <sup>1</sup> 2 *λρσ*(*Aρσ*) *μ <sup>γ</sup>xγ*. Then, it is easily proved that (*Aρσ*) *μ <sup>γ</sup>* = *δ μ <sup>ρ</sup> ησγ* − *δ μ <sup>σ</sup> ηργ*.

Explicitly,

$$A\_{10} = \begin{pmatrix} 0 & -1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \quad A\_{20} = \begin{pmatrix} 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \quad A\_{30} = \begin{pmatrix} 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 \end{pmatrix}$$

$$A\_{12} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \qquad A\_{23} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -1 & 0 \end{pmatrix} \qquad A\_{31} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{pmatrix}$$

Now, by defining the Lie product as the usual commutator [*Aij*, *Ahk*] = *Aij* · *Ahk* − *Ahk* · *Aij*, *<sup>A</sup>*10, *<sup>A</sup>*20, *<sup>A</sup>*30, *<sup>A</sup>*12, *<sup>A</sup>*<sup>23</sup> and *<sup>A</sup>*<sup>31</sup> generate a Lie algebra, which we denote by <sup>g</sup>6.

Let us consider *<sup>d</sup>* <sup>∈</sup> *Der*(1,1,1,1)g6 and let *<sup>A</sup>* = (*aij*), 1 <sup>≤</sup> *<sup>i</sup>*, *<sup>j</sup>* <sup>≤</sup> 6 be the 6 <sup>×</sup> 6 square matrix associated with the endomorphism *d*.

To obtain the elements of this matrix, for the pair of generators (*ei*,*ej*), with *i* < *j*, the derivation *d* satisfies *d* [*ei*,*ej*] = [*d*(*ei*),*ej*]+[*ei*, *d*(*ej*)] and *d*(*ei*) = ∑<sup>6</sup> *<sup>h</sup>*=<sup>1</sup> *aih eh*. In this way, the following conditions are obtained. This can be seen in the following Table 1.


**Table 1.** Condition obtained.

From these conditions on *aij* and <sup>∀</sup> *<sup>a</sup>*41, *<sup>a</sup>*42, *<sup>a</sup>*44, *<sup>a</sup>*46, *<sup>a</sup>*55, *<sup>a</sup>*61, *<sup>a</sup>*65, *<sup>a</sup>*<sup>66</sup> <sup>∈</sup> <sup>C</sup>, we have the following conditions shown in Table 2.

**Table 2.** Conditions obtained.


This implies that *<sup>ψ</sup>*g<sup>6</sup> (1) = *dim Der*(1,1,1,1)g<sup>6</sup> = 8.

## **Computing** *<sup>ψ</sup>*¯g<sup>6</sup> , **for** *<sup>α</sup>* <sup>=</sup> <sup>1</sup>

Let us consider *<sup>d</sup>* <sup>∈</sup> *Der*(1,1,1,1)g6. Then, [*d*[*u*, *<sup>v</sup>*], [*u*, *<sup>w</sup>*]] + [[*u*, *<sup>v</sup>*], *<sup>d</sup>*[*u*, *<sup>w</sup>*]] = *<sup>d</sup>*[[[*u*, *<sup>w</sup>*], *<sup>v</sup>*], *<sup>u</sup>*] + *<sup>d</sup>*[[[*w*, *<sup>u</sup>*], *<sup>u</sup>*], *<sup>v</sup>*], <sup>∀</sup>*u*, *<sup>v</sup>*, *<sup>w</sup>* <sup>∈</sup> g6.

To obtain the elements *aij* of the corresponding 6 × 6 square matrix associated with *d*, we see that for each triplets of generators (*ei*,*ej*,*ek*) of the algebra, the previous expression is written as

$$d\left[d[e\_i, e\_j], [e\_i, e\_k]\right] + \left[[e\_i, e\_j], d[e\_i, e\_k]\right] = d\left[\left([e\_i, e\_k], e\_j\right), e\_i\right] + d\left[\left([e\_k, e\_i], e\_i\right), e\_j\right].$$

Starting from it, we obtain the following conditions shown in Table 3.


**Table 3.** Conditions obtained.

It follows from these conditions for *aij* that *aij* = 0, ∀ *i*, *j* ∈ {1,2,3,4,5,6 }. This implies that *<sup>ψ</sup>*¯g<sup>6</sup> (1, 1) = *dim Der*(1,1,1,1)g6 = 0, which proves that *<sup>ψ</sup>* = *<sup>ψ</sup>*¯ in general.

*3.2. The Quantum-Mechanical Model Based on a* 5*th Heisenberg Algebra*

In this section, and by using the invariant function previously introduced *ψ*¯, we prove the following result.

#### **Theorem 3.** Main Theorem

*The five-dimensional classical-mechanical model built upon certain types of five-dimensional Lie algebras cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra.*

**Proof.** Let <sup>H</sup><sup>5</sup> be the fifth Heisenberg algebra generated by {*e*1, ... ,*e*5} and defined by the brackets [*e*1,*e*3] = *e*<sup>5</sup> and [*e*2,*e*4] = *e*5.

Let us consider *<sup>d</sup>* <sup>∈</sup> *Der*(1,1,1,1)H5. Then, [*d*[*u*, *<sup>v</sup>*], [*u*, *<sup>w</sup>*]] + [[*u*, *<sup>v</sup>*], *<sup>d</sup>*[*u*, *<sup>w</sup>*]] = *<sup>d</sup>*[[[*u*, *<sup>w</sup>*], *<sup>v</sup>*], *<sup>u</sup>*] + *<sup>d</sup>*[[[*w*, *<sup>u</sup>*], *<sup>u</sup>*], *<sup>v</sup>*], <sup>∀</sup>*u*, *<sup>v</sup>*, *<sup>w</sup>* <sup>∈</sup> <sup>H</sup>5.

To obtain the elements *aij* of the corresponding 5 × 5 square matrix associated with *d*, we see that for each triplet of generators (*ei*,*ej*,*ek*) of the algebra, the previous expression is written as

[*d*[*ei*,*ej*], [*ei*,*ek*]] + [[*ei*,*ej*], *d*[*ei*,*ek*]] = *d*[[[*ei*,*ek*],*ej*],*ei*] + *d*[[[*ek*,*ei*],*ei*],*ej*].

Note that, in this case, there is no restriction on the elements of the matrix associated with *d* and, thus, *ψ*¯H<sup>5</sup> (1, 1) = *dim Der*(1,1,1,1)H<sup>5</sup> = 25.

For another part, let <sup>f</sup><sup>5</sup> be the five-dimensional filiform Lie algebra, defined by [*e*1,*e*3] = *<sup>e</sup>*2, [*e*1,*e*4] = *<sup>e</sup>*<sup>3</sup> and [*e*1,*e*5] = *e*4.

Let us consider *<sup>d</sup>* <sup>∈</sup> *Der*(1,1,1,1)f5. Then, it is verified that [*d*[*u*, *<sup>v</sup>*], [*u*, *<sup>w</sup>*]] + [[*u*, *<sup>v</sup>*], *<sup>d</sup>*[*u*, *<sup>w</sup>*]] = *<sup>d</sup>*[[[*u*, *<sup>w</sup>*], *<sup>v</sup>*], *<sup>u</sup>*] + *<sup>d</sup>*[[[*w*, *<sup>u</sup>*], *<sup>u</sup>*], *<sup>v</sup>*], <sup>∀</sup>*u*, *<sup>v</sup>*, *<sup>w</sup>* <sup>∈</sup> f5.

Similar to the previous case, to obtain the elements *aij* of the corresponding 5 × 5 square matrix associated with *d*, we see that, for each triplet of generators (*ei*,*ej*,*ek*) of the algebra, the previous expression is written as

$$[d[e\_i, e\_j], [e\_i, e\_k]] + [[e\_i, e\_j], d[e\_i, e\_k]] = d[[[[e\_i, e\_k], e\_j], e\_i] + d[[[[e\_k, e\_i], e\_i], e\_j]].$$

In this case, the restrictions of the matrix associated with *d* are *a*<sup>21</sup> = 0, obtained from the bracket (*e*1,*e*3,*e*5) and *a*<sup>31</sup> = *a*<sup>41</sup> from (*e*1,*e*4,*e*5), therefore *ψ*¯(1, 1) = 23.

Next, we use the highly non-trivial result, which was originally proved by Borel [17]: *If* <sup>g</sup><sup>0</sup> *is a proper contraction of a complex Lie algebra* g*, then it holds:* dim *Der*g < dim *Der*g<sup>0</sup> .

Indeed, according to Proposition 1 we obtain that

$$\psi\_{\mathbb{H}\_{\mathbb{S}}}(1,1) = \dim\left(Der\_{(1,1,1)}\mathbb{H}\_{\mathbb{S}}\right) = \dim\left(Der\left(\mathbb{H}\_{\mathbb{S}}\right)\right) = 25$$

and

$$\psi\_{\mathfrak{f}\mathfrak{s}}(1,1) = \dim\left(Der\_{(1,1,1,1)}\mathfrak{f}\_{\mathfrak{f}}\right) = \dim\left(Der\left(\mathfrak{f}\_{\mathfrak{f}}\right)\right) = 23.$$

It implies that no proper contraction transforming the Heisemberg algebra H<sup>5</sup> into the filiform Lie algebra <sup>f</sup><sup>5</sup> exists. Thus, since both algebras are not isomorphic, the five-dimensional classical-mechanical model built upon a five-dimensional filiform Lie algebra cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra.

#### **4. Discussion and Conclusions**

In this paper, we introduce an invariant two-parameter function of algebras, *ψ*¯, and we have used it as a tool to study contractions of certain particular types of algebras.

Indeed, by means of this function, we have proved that there is no proper contraction between a fifth Heisenberg algebra and a filiform Lie algebra of dimension 5. It implies, as a main result, that the five-dimensional classical-mechanical model built upon a five-dimensional filiform Lie algebra cannot be obtained as a limit process of a quantum-mechanical model based on a fifth Heisenberg algebra.

We have also computed this function in the case of other types of algebras, for instance, Malcev algebras of the type Lie and the Lie algebra induced by the Lorentz group *SO*(3, 1).

Apart from continuing this study with with higher-dimensional algebras, we indicate next some open problems to be dealt with in future work, most of them with the objective of trying to find some possible interesting physical applications for the filiform Lie algebras. They are the following


To perform this task, it is necessary to write the boson operators involved in the Hamiltonian in term of new ones that fulfill the commutation relations for a given filiform Lie algebra. However, at that point, we find the difficulty that we should employ a tensorial product of two filiform Lie algebras in order to describe the system properly. That means that an isomorphism between the semi-simple Lie algebra of the original hamiltonian and the filiform Lie algebra proposed to describe the physical system should exist. Fortunately, it seems that we have obtained a theorem that can confirm that kind of isomorphism.

Now, the advantage that we gain employing a filiform Lie algebra instead of a semi-simple Lie algebra is that we could map a non-linear problem such as the problem described by a system with up to two-body interactions onto a linear problem with just one-body interactions. On the other hand, once we have described the system in terms of the filiform Lie algebra, it is necessary to define the branching rules, that is to find the irreducible representations of an algebra g contained in a given representation of g. Since the representations are interpreted as quantum mechanical states, it is necessary to provide a complete set of quantum numbers (labels) to characterize uniquely the basis of the system. This is a non-trivial task that it may even lead to a further research.

4. Another possible physical applications of the present topic is to study phase spaces by using filiform Lie algebras as a tool.

In this respect, Arzano and Nettel [18] in 2016 introduced a general framework for describing deformed phase spaces with group valued momenta. Using techniques from the theory of Poisson–Lie groups and Lie bialgebras, they developed tools for constructing Poisson structures on the deformed phase space starting from the minimal input of the algebraic structure of the generators of the momentum Lie group. These tools developed are used to derive Poisson structures on examples of group momentum space much studied in the literature such as the *n*-dimensional generalization of the *κ*-deformed momentum space and the *SL*(2, *R*) momentum space in three space-time dimensions. They also discussed classical momentum observables associated to multiparticle systems and argued that these combined according the usual four-vector addition despite the non-Abelian group structure of momentum space (see [18] for further information).

In that paper, the authors work with a phase space Γ = *T* × *G*, given by the Cartesian product of a *n*-dimensional Lie group configuration space *T* and a *n*-dimensional Lie group momentum space *<sup>G</sup>*. Since *<sup>T</sup>* and *<sup>G</sup>* are Lie groups, we can consider their associated Lie algebras t and g so that we can define a Lie–Poisson algebra, which can endow a mathematical structure to the phase space Γ. Indeed, Arzano and Nettel considered a phase space Γ in which the component related to momentum is an *n*-dimensional Lie sub-group of the (*n* + 2)-dimensional Lorentz group *SO*(*n* + 1, 1), denoted as *AN*(*n*).

Taking into consideration this paper, we have tried to construct a phase space similar to the one by those authors, although we have taken the (*n* + 2)-dimensional Lorentz group *SO*(*n* + 1, 2) as the Lie group related to momentum.

We began our research on this subject considering the Lie group *SO*(2, 2) and using the same procedure as Arzano and Nettel did. However, we realized that that attempt was going to be very complicated because of the great dimensions of the matrices involved (in the computations, a 49 × 49 *r*−matrix appeared).

Therefore, the fact of finding a Poisson structure that allows us to endow the phase space Γ = *T* × *SO*(*n* + 1, 2) with a mathematical structure is another problem, which we consider open.

5. Finally, semi-invariant functions of algebras could also be considered to study contractions of Lie Algebras (see [19], for instance).

We will dedicate our efforts to these objectives in future work.

#### **5. Materials and Methods**

Since this is a work on pure and applied mathematics, no type of materials different from the usual ones in a theoretical investigation was needed. Indeed, on the one hand, only the existing bibliography on the subject and, on the other hand, a suitable symbolic computation package were used. In the same way, with regard to the methodology used for the writing of the manuscript, it was also the usual one in research work of this nature, namely, based on already established hypotheses and known results.

We used the SAGE symbolic computation package for computations. SageMath, which is a free open-source mathematics software system licensed under the GPL, builds on top of many existing open-source packages, such as matplotlib, Sympy, Maxima, GAP, R and many more (see [20], for instance).

**Author Contributions:** All of the authors wrote the paper and J.N.-V. and P.P.-F. checked the proofs.

**Funding:** This research was funded by the Spanish Ministerio de Ciencia e Innovación and Junta de Andalucía via grants No. MTM2013-40455-P and No. FQM-326 (J.N.-V.) and No. FQM-160 (P.P.-F.).

**Acknowledgments:** The authors gratefully acknowledge the financial support above mentioned.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Inönü, E.; Wigner, E. On the contraction of groups and their representations. *Proc. Nat. Acad. Sci. USA* **1953**, *39*, 510–524. [CrossRef] [PubMed]


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Solving Nonholonomic Systems with the Tau Method**

**Alexandra Gavina 1,\*, José M. A. Matos 1,2 and Paulo B. Vasconcelos 2,3**


Received: 27 September 2019; Accepted: 17 October 2019; Published: 19 October 2019

**Abstract:** A numerical procedure based on the spectral Tau method to solve nonholonomic systems is provided. Nonholonomic systems are characterized as systems with constraints imposed on the motion. The dynamics is described by a system of differential equations involving control functions and several problems that arise from nonholonomic systems can be formulated as optimal control problems. Applying the Pontryagins maximum principle, the necessary optimality conditions along with the transversality condition, a boundary value problem is obtained. Finally, a numerical approach to tackle the boundary value problem is required. Here we propose the Lanczos spectral Tau method to obtain an approximate solution of these problems exploiting the Tau toolbox software library, which allows for ease of use as well as accurate results.

**Keywords:** Tau method; nonholonomic systems

#### **1. Introduction**

Nonholonomic systems are a class of nonlinear systems that cannot be stabilized by a continuous time-invariant feedback, i.e., at a certain time or state there are constraints imposed on the motion (nonholonomic constraints). These systems are controllable but they cannot move instantaneously in certain directions. They belong to a class of nonlinear differential systems with nonintegrable constraints imposed on the motion [1].

Nonholonomic control systems, which result from formulations of nonholonomic systems that include control inputs, are nonlinear control problems requiring nonlinear treatment. There is ample literature on the formulation of the equations of motion and on the dynamics of nonholonomic systems, being [2] an excellent survey for examples. Nonholonomic control systems have been studied in the context of robot manipulation, mobile robots, wheeled vehicles, and space robotics, just to mention a few. In the case of wheeled vehicles, the kinematics and dynamics can be modeled based on the assumption that the wheels are ideally rolling. Typical constraints of wheeled vehicles are rolling contact, like rolling between the wheels and the ground without slipping, or sliding contact such the sliding of skates.

The solution of nonholonomic optimal control problems can be obtained following a standard procedure, which consists of applying Pontryagin's maximum/minimum principle to obtain set of equations along with initial and terminal conditions resulting into a two-point boundary value problem (BVP) [3].

In this work we propose the use of the spectral Tau method to obtain approximate solutions of nonholonomic optimal control problems through the associated BVP.

The spectral Tau method produces a polynomial approximation of the solution of the differential problem. It is based on solving a system of linear algebraic equations, obtained by imposing that all conditions are verified exactly and the residual is orthogonal to the first elements of an orthogonal polynomial basis [4].

The paper is organized as follows. Section 2 describes the system model of the nonholonomic wheeled vehicle and Section 3 explains the optimal control formulation. A brief description of the Tau method is presented in Section 4. An illustrative example with numerical results is provided in Section 5 and some conclusions are drawn in Section 6.

#### **2. Nonholonomic Wheeled Vehicle Model**

Vehicle models are usually described by a set of ordinary differential equations that define the dynamics of the vehicle and the relationship between the state variables and control input. The kinematic model of a wheeled vehicle can be defined by the following differential equations

$$\begin{array}{l} \dot{\mathfrak{x}} = f\_1(\mathfrak{x}, \mathfrak{y}, \mathfrak{v}, \theta, \phi) \\ \dot{\mathfrak{y}} = f\_2(\mathfrak{x}, \mathfrak{y}, \mathfrak{v}, \theta, \phi) \\ \dot{\theta} = f\_3(\mathfrak{x}, \mathfrak{y}, \mathfrak{v}, \theta, \phi) \end{array}$$

where (*x*, *<sup>y</sup>*) <sup>∈</sup> <sup>R</sup><sup>2</sup> is the robot's position in space, *<sup>θ</sup>* is the angle with respect to the *<sup>x</sup>*-axis, *<sup>φ</sup>* is the steering wheel's angle with respect to the robot's longitudinal axis, and *v* is the velocity. (see Figure 1).

**Figure 1.** Car-like robot model.

A nonholonomic car-like robot is a car model which rolls without slipping between the wheels and the ground. This constraint is expressed by the equation [5]

$$
\dot{x}\sin(\theta) - \dot{y}\cos(\theta) = 0.\tag{1}
$$

The simplest model corresponds to a robot with a single wheel: the unicycle model. In this model the wheel rolls on a plane while keeping its body vertical. It is an unrestricted model since it can rotate freely while standing in its position (*x*, *y*). Furthermore, the dynamics are characterized by

$$\begin{cases} \dot{\mathbf{x}} = v \cos(\theta) \\ \dot{\mathbf{y}} = v \sin(\theta) \\ \dot{\theta} = \phi. \end{cases} \tag{2}$$

The kinematic model of a car-like robot has the same state variables as the unicycle model and its dynamic is represented by

$$\begin{cases} \dot{\mathbf{x}} = v \cos(\theta) \\ \dot{\mathbf{y}} = v \sin(\theta) \\ \dot{\theta} = vu\_\prime \end{cases} \tag{3}$$

where *<sup>u</sup>* <sup>∈</sup> [−<sup>1</sup> *r* , 1 *<sup>r</sup>* ] stands for the curvature and *r* for the turning radius of the robot that corresponds the maximum curvature [6].

#### **3. Optimal Control Problems**

An optimal control problem (OCP) can be formulated as

$$\text{Minimize } J(\mathbf{x}(t), \mathbf{u}(t)) = \int\_{t\_0}^{t\_f} F(\mathbf{x}(t), \mathbf{u}(t)) \, dt + G(\mathbf{x}(t\_f), t\_f)$$

subject to

$$\begin{aligned} \dot{\mathbf{x}}(t) &= f(\mathbf{x}(t), \mathbf{u}(t)), \quad t \in [t\_0, t\_f] \\ \mathbf{x}(t\_0) &= \mathbf{x}\_0 \\ \mathbf{x}(t\_f) &= \mathbf{x}\_f \\ \mathbf{x}(\cdot) &\in X \\ \mathbf{u}(\cdot) &\in U \\ t\_f &\in [t\_{0'} + \infty[\_{\prime}] \end{aligned} \tag{4}$$

where *J* is the cost function, x is the state vector representing the dynamics, u is the control vector, x<sup>0</sup> is the initial configuration and x*<sup>f</sup>* is the final configuration.

The solution of the OCP can be obtained following a standard procedure, which consists in applying Pontryagin's maximum principle and obtaining the necessary optimality conditions along with the transversality condition resulting into a two-point boundary value problem (BVP).

#### *Pontryagin Maximum Principle*

Considering the Hamiltonian function

$$H(\mathbf{x}(t), \mathbf{u}(t), \lambda) = F(\mathbf{x}(t), \mathbf{u}(t)) + \lambda f(\mathbf{x}(t), \mathbf{u}(t)),\tag{5}$$

where *F* and *f* are the functions described above and *λ* = [*λ*1(*t*), *λ*2(*t*), ... , *λn*(*t*)] is a vector of co-state variables, and considering that (x, u∗) is a controlled trajectory defined over the interval [*t*0, *tf* ] then (x, u∗) is optimal, for all admissible controls u, if the Pontryagin's maximum principle holds, i.e.,

$$H(\mathbf{x}, \mathbf{u}^\*, \lambda) \ge H(\mathbf{x}, \mathbf{u}, \lambda).$$

The Pontryagin maximum principle guarantees that if (x, u∗) is an optimal pair, a solution of the problem (4), then the first order necessary conditions

$$H\_{\mathbf{x}} = -\lambda \tag{6}$$

together with the stationary conditions

$$H\_{\mathfrak{u}} = 0 \tag{7}$$

satisfies the Hamiltonian maximization with transversality conditions given by [3]

$$\begin{aligned} \lambda(t\_f) &= 0 \qquad \text{if} \qquad t\_f = \infty \text{ or } G(.) = 0\\ \text{or} \\ \lambda(t\_f) &= \left. \frac{\partial G}{\partial \mathbf{x}} \right|\_{t = t\_f} . \end{aligned} \tag{8}$$

This reduces the constrained problem (4) to an unconstrained differential equations system (6)–(8). Usually nonlinear, this system of differential equations can be approximately solved by a numerical method. The next section is devoted to introducing the Tau method, which will be used to numerically tackle the problem.

#### **4. Tau Method**

The spectral Tau method produces a polynomial approximation, *yn*(*x*), of the solution, *y*(*x*), of a given differential problem *Dy*(*x*) = *f*(*x*), satisfying a set of conditions defined on an interval ]*a*, *b*[. Introduced by Lanczos in 1938 [7] to compute approximate solutions of linear differential problems with polynomial coefficients and right-hand side, the Tau method solves a tuned system of linear algebraic equations obtained by imposing that the conditions are verified exactly and that the residual is minimized in a quadrature sense, i.e, is orthogonal to the first elements of an orthogonal polynomial basis. It can be applied, indifferently, to initial, boundary or mixed value problems and it can be implemented with any orthogonal basis. We begin by introducing the method for the original case and then shed some light on how to extend it for the solution of nonlinear problems with non-polynomial coefficients.

#### *4.1. Preliminaries and Notation*

Let <sup>P</sup> be the space of all algebraic polynomials and let *<sup>D</sup>* : <sup>P</sup> <sup>→</sup> <sup>P</sup> be a linear differential operator of order *ν* ≥ 1 with polynomial coefficients represented by

$$D \equiv \sum\_{r=1}^{\nu} p\_r \frac{d^r}{d\mathbf{x}^r} \tag{9}$$

and *y* be the exact solution of the differential problem

$$\begin{cases} \, \, Dy(\mathbf{x}) = f(\mathbf{x}), & a < \mathbf{x} < b \\\, \, g\_j(y) = \sigma\_{j\prime} & j = 1, \dots, \nu\_{\prime} \end{cases} \tag{10}$$

where *<sup>f</sup>* <sup>∈</sup> <sup>P</sup> is a polynomial or a convenient polynomial approximation and *gj* are *<sup>ν</sup>* linear functionals, acting on *Cν*[*a*, *b*], representing the (supplementary) conditions.

The main idea of the Tau method is to approximate *y* by the polynomial *yn*, solution of the perturbed problem

$$\begin{cases} \, \, D y\_n(\mathbf{x}) = f(\mathbf{x}) + \tau\_n(\mathbf{x}), & a < \mathbf{x} < b \\\, \, g\_j(y\_n) = \sigma\_{j\prime} & j = 1, \dots, \nu \end{cases} \tag{11}$$

where *τ<sup>n</sup>* is a polynomial perturbation close to zero in ]*a*, *b*[. Choosing an orthogonal polynomial basis **P** = [*P*0, *P*<sup>1</sup> ...], then the coefficients of *yn* are determined imposing that *τ<sup>n</sup>* is orthogonal to *Pi*, *i* = 0, 1, . . . , *n* − *ν*.

The original idea of the Tau method [8] is based on the minimax property of Chebyshev polynomials and on the fact that the solution *yn* of (11) depends continuously on the residual *τn*. Later generalized to more general bases [4], the method looks for a residual *τ<sup>n</sup>* that minimizes the

weighted norm ||.||*<sup>w</sup>* associated to the sequence **P**. Indeed, **P** is orthogonal with respect to the weight function *w*(*x*) on [*a*, *b*]

$$\langle P\_{\dot{\imath}\prime}P\_{\dot{\jmath}}\rangle = \int\_{a}^{b} w(\mathfrak{x})P\_{\dot{\imath}}(\mathfrak{x})P\_{\dot{\jmath}}(\mathfrak{x})d\mathfrak{x} = \,w\_{\dot{\imath}}\delta\_{\dot{\imath}\prime}$$

where *wi* = *Pi*, *Pi* = ||*Pi*||<sup>2</sup> *<sup>w</sup>* and *δij* is the Kronecker delta, and (11) is achieved by imposing

$$\langle \pi\_{\mathfrak{n}}, P\_{\mathfrak{j}} \rangle = 0, \mathfrak{j} = 0, 1, \dots, n - \nu\_{\mathfrak{n}}$$

that is, *τ<sup>n</sup>* = O(*Pn*−*ν*).

Using suitable matrices, see for example [4,9,10], the differential problem (11) is translated into an algebraic problem. Such matrices must be computed with criteria in order to ensure stable computations [10].

#### *4.2. Nonlinear Problems*

Nonlinear differential problems are solved iteratively by first linearizing the problem and then applying the Tau method to the linear inner problem.

Let

$$F(y) = 0, \ x \in ]a, b[\tag{12}$$

be a differential equation, where *F* is a differential operator that could be nonlinear in *y* and on its derivatives.

From *ym*, an approximation of the exact solution *y*, we take the first order Taylor polynomial of *F* centered in *ym* to approximate *F*

$$D\_{\mathbf{m}}y = F(y\_{\mathbf{m}}) + \sum\_{k=0}^{\nu} (y^{(k)} - y\_{\mathbf{m}}^{(k)}) \left. \frac{\partial F}{\partial y^{(k)}} \right|\_{y\_{\mathbf{m}}}.\tag{13}$$

*F* can be replaced by *Dm* in (12) to solve

$$\left. \sum\_{k=0}^{\nu} y^{(k)} \left. \frac{\partial F}{\partial y^{(k)}} \right|\_{y\_m} = -F(y\_m) + \sum\_{k=0}^{\nu} y\_m^{(k)} \left. \frac{\partial F}{\partial y^{(k)}} \right|\_{y\_m} \right. \tag{14}$$

Applying the Tau method to the linear differential Equation (14), an iterative process is implemented to get increasingly better approximations for the differential problem

$$D\_m y\_{m+1} = 0, \; m = 0, 1, \dots \tag{15}$$

For additional details on the use of the Tau method for nonlinear problems the reader is invited to read [11].

#### **5. Numerical Experiments**

The proposed example is based on the work of [12] and it will be tackled by the Tau method, described in Section 4.

The problem at hand is an optimal control problem of the form (4) with *F* = <sup>1</sup> <sup>2</sup> (*u*<sup>2</sup> <sup>1</sup> + *<sup>u</sup>*<sup>2</sup> <sup>2</sup>) and *G* = <sup>1</sup> <sup>2</sup> **<sup>x</sup>***T*(*tf*)*Q***x**(*tf*) where **<sup>x</sup>** is the vector of state variables and *<sup>Q</sup>* is the weighting matrix.

#### *5.1. System Model*

Consider an automated vehicle that moves on a horizontal plane, the contact of each wheel with the floor is assumed to satisfy the rolling without slipping condition and the control inputs are the torque generated by two motors mounted on the wheels. For a fixed final time, it is desired to find the control inputs that minimizes the energy of the final state. This system can be modeled by

$$J = \min \int\_{t\_0}^{t\_f} \frac{1}{2} (u\_1^2 + u\_2^2) dt + \frac{1}{2} \mathbf{x}^T(t\_f) \mathbf{Q} \mathbf{x}(t\_f) \tag{16}$$

subject to

$$
\lambda m \ddot{x} = \frac{\cos(\theta)}{R} (\mu\_1 + \mu\_2) - \lambda \sin(\theta) \tag{17}
$$

$$m\ddot{y} = \frac{\sin(\theta)}{R}(\mu\_1 + \mu\_2) + \lambda\cos(\theta) \tag{18}$$

$$I\ddot{\theta} = \frac{L}{R}(\mu\_1 - \mu\_2) \tag{19}$$

with nonholonomic constraint

$$
\dot{x}\sin(\theta) - \dot{y}\cos(\theta) = 0,
$$

where the position coordinates (*x*, *y*), and the heading angle *θ*, define the system configuration. The mass *m*, the inertia *I*, the wheels' radius *R*, the half-length of the axis *L*, are parameters of the system and *μ*<sup>1</sup> and *μ*<sup>2</sup> are torques generated by the motors.

Defining the control inputs **u** = [*u*<sup>1</sup> *u*2] *<sup>T</sup>* as

$$
\mu\_1 = \frac{1}{mR}(\mu\_1 + \mu\_2), \quad \mu\_2 = \frac{L}{IR}(\mu\_1 - \mu\_2)
$$

and the state variables **x** = [*x*<sup>1</sup> *x*<sup>2</sup> *x*<sup>3</sup> *x*<sup>4</sup> *x*5] *<sup>T</sup>* as

$$\begin{aligned} x\_1 &= x \cos(\theta) + y \sin(\theta) \\ x\_2 &= \theta \\ x\_3 &= -x \sin(\theta) + y \cos(\theta) \\ x\_4 &= \dot{x} \cos(\theta) + \dot{y} \sin(\theta) - \dot{\theta} (x \sin(\theta) - y \cos(\theta)) \\ x\_5 &= \theta. \end{aligned}$$

Equations (17)–(19) can be reduced to the following system of differential equations:

$$\begin{cases} \dot{\mathbf{x}}\_1 = \mathbf{x}\_4 \\ \dot{\mathbf{x}}\_2 = \mathbf{x}\_5 \\ \dot{\mathbf{x}}\_3 = -\mathbf{x}\_1 \mathbf{x}\_5 \\ \dot{\mathbf{x}}\_4 = -\mathbf{x}\_1 \mathbf{x}\_5^2 + \boldsymbol{\mu}\_1 + \boldsymbol{\mu}\_2 \mathbf{x}\_3 \\ \dot{\mathbf{x}}\_5 = \boldsymbol{\mu}\_2. \end{cases}$$

Applying the Hamiltonian *H*, defined in (5),

$$H(\mathbf{x}, \mathbf{u}, \lambda) = \frac{1}{2} (\mu\_1^2 + \mu\_2^2) + \lambda \dot{\mathbf{x}}^2$$

with *λ* = [*λ*<sup>1</sup> *λ*<sup>2</sup> *λ*<sup>3</sup> *λ*<sup>4</sup> *λ*5], and calculating the necessary conditions (6), the stationary conditions (7) and the transversality conditions (8), the following second order system of differential equations is obtained:

$$\begin{cases} \mathbf{x}\_1 \dot{\mathbf{x}}\_2 + \dot{\mathbf{x}}\_3 &= 0 \\ \ddot{\mathbf{x}}\_1 - \dot{\mathbf{x}}\_2 \dot{\mathbf{x}}\_3 - \ddot{\mathbf{x}}\_2 \mathbf{x}\_3 + \lambda\_4 &= 0 \\ \ddot{\mathbf{x}}\_2 + \mathbf{x}\_3 \lambda\_4 + \lambda\_5 &= 0 \\ -\dot{\mathbf{x}}\_2 \lambda\_3 - \dot{\mathbf{x}}\_2^2 \lambda\_4 + \dot{\lambda}\_4 &= 0 \\ \ddot{\mathbf{x}}\_2 \lambda\_4 + \dot{\lambda}\_3 &= 0 \\ -\mathbf{x}\_1 \lambda\_3 + 2 \dot{\mathbf{x}}\_3 \lambda\_4 + \lambda\_2 + \dot{\lambda}\_5 &= 0 \end{cases}$$

for *<sup>x</sup>*˙1 <sup>=</sup> *<sup>x</sup>*4, *<sup>x</sup>*˙2 <sup>=</sup> *<sup>x</sup>*5, *<sup>λ</sup>*˙ <sup>4</sup> <sup>=</sup> <sup>−</sup>*λ*<sup>1</sup> and *<sup>λ</sup>*˙ <sup>2</sup> <sup>=</sup> 0, with initial and transversality conditions given by

$$\begin{cases} \begin{aligned} \mathbf{x}\_{i}(t\_{0}) &= \mathbf{x}\_{i0}, \qquad i \in \{1,2,3\} \\ \dot{\mathbf{x}}\_{1}(t\_{0}) &= \mathbf{x}\_{40} \\ \dot{\mathbf{x}}\_{2}(t\_{0}) &= \mathbf{x}\_{50} \\ \lambda\_{i}(t\_{f}) &= \lambda\_{if}, \qquad i \in \{3,4,5\} \\ \lambda\_{4}(t\_{f}) &= -\lambda\_{1f}. \end{aligned} \tag{20}$$

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Since this is a nonlinear differential system, in order to implement the Tau method, differential equations need to be linearized. Expressions of the form *uv* and *uv*<sup>2</sup> will be replaced, respectively, by

$$\begin{aligned} \mu \upsilon &\approx \upsilon\_m \iota + \iota\_m \upsilon - \iota\_m \upsilon\_m \\ \mu \upsilon^2 &\approx \upsilon\_m^2 \iota + 2 \iota\_m \upsilon\_m \upsilon - 2 \iota\_m \upsilon\_m^2 \iota \end{aligned}$$

Thus, the differential system becomes

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *x*˙2,*mx*<sup>1</sup> + *x*1,*mx*˙2 + *x*˙3 = *x*1,*mx*˙2,*<sup>m</sup> x*¨1 − *x*˙3,*mx*˙2 − *x*3,*mx*¨2 − *x*¨2,*mx*<sup>3</sup> − *x*˙2,*mx*˙3 + *λ*<sup>4</sup> = −*x*˙2,*mx*˙3,*<sup>m</sup>* − *x*¨2,*mx*3,*<sup>m</sup> x*¨2 + *λ*4,*mx*<sup>3</sup> + *x*3,*mλ*<sup>4</sup> + *λ*<sup>5</sup> = *x*3,*mλ*4,*<sup>m</sup>* (*λ*3,*<sup>m</sup>* + 2*x*˙2,*mλ*4,*m*)*x*˙2 + *x*˙2,*mλ*<sup>3</sup> + *x*˙ 2 2,*mλ*<sup>4</sup> <sup>+</sup> *<sup>λ</sup>*¨ <sup>4</sup> <sup>=</sup> *<sup>x</sup>*˙2,*mλ*3,*<sup>m</sup>* <sup>+</sup> <sup>2</sup>*x*˙ 2 2,*mλ*4,*<sup>m</sup> λ*4,*mx*¨2 + *λ*˙ <sup>3</sup> + *x*¨2,*mλ*<sup>4</sup> = *x*¨2,*mλ*4,*<sup>m</sup> <sup>λ</sup>*3,*mx*<sup>1</sup> <sup>−</sup> <sup>2</sup>*λ*4,*mx*˙3 <sup>+</sup> *<sup>x</sup>*1,*mλ*<sup>3</sup> <sup>−</sup> <sup>2</sup>*x*˙3,*mλ*<sup>4</sup> <sup>−</sup> *<sup>λ</sup>*˙ <sup>5</sup> <sup>=</sup> *<sup>x</sup>*1,*mλ*3,*<sup>m</sup>* <sup>−</sup> <sup>2</sup>*x*˙3,*mλ*4,*<sup>m</sup>* <sup>+</sup> *<sup>λ</sup>*2,0 (21)

where *xi*,*m*, *i* = 1, 2, 3 and *λi*,*m*, *i* = 3, 4, 5 are approximations to *x*1, *x*2, *x*<sup>3</sup> and *λ*3, *λ*4, *λ*5.

The matrix representation of the differential problem (21), together with the conditions (20), is of the form *Ta* = *b*, where

$$T = \begin{bmatrix} C \\ D \end{bmatrix} \quad \text{and} \quad b = \begin{bmatrix} s \\ f \end{bmatrix}$$

$$\mathbf{C} = \begin{bmatrix} \mathbf{P}(t\_0) & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbf{P}(t\_0) & 0 & 0 & 0 & 0 \\ 0 & 0 & \mathbf{P}(t\_0) & 0 & 0 & 0 \\ \mathbf{P}'(t\_0) & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbf{P}'(t\_0) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \mathbf{P}(t\_f) & 0 & 0 \\ 0 & 0 & 0 & 0 & \mathbf{P}(t\_f) & 0 \\ 0 & 0 & 0 & 0 & 0 & \mathbf{P}(t\_f) \\ 0 & 0 & 0 & 0 & \mathbf{P}'(t\_f) & 0 \end{bmatrix}$$

with

where **P**(*ti*)=[*P*0(*ti*), *P*1(*ti*), ...] represents the action of boundary conditions (20) over the polynomial base elements and **P'**(*ti*) its derivatives. Matrix *D* represents the differential operator

$$D = \begin{bmatrix} \dot{\mathbf{x}}\_{2,\mathsf{m}}(M) & \mathbf{x}\_{1,\mathsf{m}}(M)N & N & 0 & 0 & 0\\ N^2 & A\_{2,2} & A\_{2,3} & 0 & I & 0\\ 0 & N^2 & \lambda\_{4,\mathsf{m}}(M) & 0 & \mathbf{x}\_{3,\mathsf{m}}(M) & I\\ 0 & A\_{4,2} & 0 & -\dot{\mathbf{x}}\_{2,\mathsf{m}}(M) & A\_{4,5} & 0\\ 0 & \lambda\_{4,\mathsf{m}}(M)N^2 & 0 & N & \ddot{\mathbf{x}}\_{2,\mathsf{m}}(M) & 0\\ \lambda\_{3,\mathsf{m}}(M) & 0 & -2\lambda\_{4,\mathsf{m}}(M)N & \mathbf{x}\_{1,\mathsf{m}}(M) & -2\dot{\mathbf{x}}\_{3,\mathsf{m}}(M) & N \end{bmatrix}.$$

Finally,

$$s = \left[ \mathbf{x}\_{10\prime} \mathbf{x}\_{20\prime} \mathbf{x}\_{30\prime} \mathbf{x}\_{40\prime} \mathbf{x}\_{50\prime} \boldsymbol{\lambda}\_{3f\prime} \boldsymbol{\lambda}\_{4f\prime} \boldsymbol{\lambda}\_{5f\prime} - \boldsymbol{\lambda}\_{1f} \right]^2$$

$$f = \begin{bmatrix} \begin{array}{c} \chi\_{1,\text{m}} \dot{\chi}\_{2,\text{m}} \\ -\dot{\chi}\_{2,\text{m}} \dot{\chi}\_{3,\text{m}} - \ddot{\chi}\_{2,\text{m}} \chi\_{3,\text{m}} \end{array} \\ \begin{array}{c} \chi\_{3,\text{m}} \lambda\_{4,\text{m}} \\ \dot{\chi}\_{2,\text{m}} \lambda\_{3,\text{m}} + 2 \dot{\chi}\_{2,\text{m}}^{2} \lambda\_{4,\text{m}} \\ \chi\_{1,\text{m}} \lambda\_{3,\text{m}} - 2 \dot{\chi}\_{3,\text{m}} \lambda\_{4,\text{m}} + \lambda\_{20} \end{array} \end{array}' $$

where

$$\begin{cases} A\_{2,2} = -\dot{\mathfrak{x}}\_{3,\mathfrak{m}}(M)N - \mathfrak{x}\_{3,\mathfrak{m}}(M)N^2 \\ A\_{2,3} = -\dot{\mathfrak{x}}\_{2,\mathfrak{m}}(M) - \dot{\mathfrak{x}}\_{2,\mathfrak{m}}(M)N \\ A\_{4,2} = \left(\lambda\_{3,\mathfrak{m}}(M) + 2\dot{\mathfrak{x}}\_{2,\mathfrak{m}}\lambda\_{4,\mathfrak{m}}(M)\right)N \\ A\_{4,5} = -\dot{\mathfrak{x}}\_{2,\mathfrak{m}}^2(M) + N^2 \end{cases}$$

and *M* and *N* are matrices described in [4] representing, respectively, the multiplication and the differentiation operator.

#### *5.2. Numerical Results*

In this section we report the numerical results for the example described in Section 5.1 with initial positions (*x*, *y*)=(10, 3) and heading angle 0◦. The time interval is [*t*0, *f <sup>f</sup>* ]=[0, 5] and the system parameters used in the simulation are *<sup>m</sup>* = 10 kg, *<sup>I</sup>* = 1.2 kg·m2, *<sup>R</sup>* = 0.05 m and *<sup>L</sup>* = 0.1 m. The weighting matrix is *Q* = 10*I*, where *I* stands for the identity matrix.

The simulation results were obtained using the Tau Toolbox with Chebyshev polynomials.

The state trajectories *x*1, *x*<sup>2</sup> and *x*<sup>3</sup> and the optimal trajectory for the position (*x*, *y*) are illustrated in Figures 2 and 3. The trajectories were obtained with 5th order Chebyshev polynomials.

**Figure 2.** State trajectories *x*1, *x*<sup>2</sup> and *x*3.

**Figure 3.** State trajectories *x*1, *x*<sup>2</sup> and *x*3.

These approximate solutions only required *m* = 8 iterations to satisfy the stopping criterion ||*x<sup>m</sup> <sup>i</sup>* <sup>−</sup> *<sup>x</sup>m*−<sup>1</sup> *<sup>i</sup>* || ≤ <sup>10</sup><sup>−</sup>14. For larger degree polynomials (10 and 15) machine precision can be achieved.

The residual produced by the Tau method in the iteration *m* is *rm* = *Dym*, where *ym* is the approximating solution of the system of homogeneous differential equations *Dy* = 0. Figure 4 plots the residual *r*1,*m*, *r*2,*<sup>m</sup>* and *r*1,*<sup>m</sup>* produced by the Tau method for the state variables *x*1, *x*<sup>2</sup> and *x*3, respectively.

**Figure 4.** Residual, *r*<sup>9</sup> = *Dy*9, produced by the Tau method for the state variables *x*1, *x*<sup>2</sup> and *x*3.

Table 1 presents the values for the functional *J* defined in (16) using polynomials of degree *<sup>n</sup>* <sup>=</sup> 5, 10 and 15. Since *<sup>u</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>*λ*<sup>4</sup> and *<sup>u</sup>*<sup>2</sup> <sup>=</sup> *<sup>x</sup>*˙5, the integral <sup>5</sup> 0 1 2 (*u*<sup>2</sup> <sup>1</sup> + *<sup>u</sup>*<sup>2</sup> <sup>2</sup>)*dt* and **x**(5)*TQ***x**(5) can be calculated using the the approximate solutions for the state variables *xi*, *i* = 1, ... , 5 and the co-state variable *λ*4.

As the polynomial multiplication and differentiation, the integration can be set into algebraic operation as well, using a suitable matrix [4].

**Table 1.** *J* values for several polynomial degree approximations.


#### **6. Conclusions**

The Lanczos spectral Tau method was used to compute approximate polynomial solutions for nonholonomic systems. A detailed illustration on the approximation procedure is offered. The Tau toolbox provides the appropriate environment to solve systems of ordinary differential problems while allowing for accurate solutions, whenever the sought solution is regular. Numerical results for this dynamical optimization problem confirm both aspects: ease of use and accuracy of approximation.

**Author Contributions:** A.G. and J.M.A.M. conceived and designed the experiments; J.M.A.M. and P.B.V. performed the experiments; A.G. and J.M.A.M. analyzed the data; A.G. and P.B.V. wrote the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Almost Exact Computation of Eigenvalues in Approximate Differential Problems**

#### **José M. A. Matos 1,2,\* and Maria João Rodrigues 2,3**


Received: 29 September 2019; Accepted: 13 November 2019; Published: 14 November 2019

**Abstract:** Differential eigenvalue problems arise in many fields of Mathematics and Physics, often arriving, as auxiliary problems, when solving partial differential equations. In this work, we present a method for eigenvalues computation following the Tau method philosophy and using Tau Toolbox tools. This Matlab toolbox was recently presented and here we explore its potential use and suitability for this problem. The first step is to translate the eigenvalue differential problem into an algebraic approximated eigenvalues problem. In a second step, making use of symbolic computations, we arrive at the exact polynomial expression of the determinant of the algebraic problem matrix, allowing us to get high accuracy approximations of differential eigenvalues.

**Keywords:** eigenvalue differential problems; spectral methods; Sturm–Liouville problems

**MSC:** 34B09; 34L15; 65L15

#### **1. Introduction**

Finding eigenfunctions of differential problems can be a hard task, at least for some classical problems. Among others, we can find literature in Sturm–Liouville problems, in Mathieu problems or in Orr–Sommerfeld problems describing the difficulties involved in the resolution of those problems [1–8]. The first difficulty consists of finding accurate numerical approximations for the respective eigenvalues.

In this work, we present a procedure based on the Ortiz and Samara's operational approach to the Tau method described in [9], where the differential problem is translated into an algebraic problem. This is achieved using the called operational matrices that represent the action of differential operators in a function. We have deduced explicit formulae for the elements of these matrices [10,11] obtained by performing operations on the bases of orthogonal polynomials and, for some families, we have exact formulae, which enables the construction of very accurate operational matrices. The Tau method has already been used for these kinds of problems [5,9,12,13]; however, our work on matrix calculation formulas adds efficiency and precision to the method.

Our main purpose is to use the Tau Toolbox, a Matlab numerical library that is being developed by our research group [14–16]. This library allows a stable implementation of the Tau method for the construction of accurate approximate solutions for integro-differential problems. In particular, the construction of the operational matrices is done automatically. These facts led us to think that the Tau Toolbox seems to be useful for these kinds of problems.

Finally, operating with symbolic variables, we define the determinant of those matrices as polynomials and use its roots as eigenvalues' approximations.

We present some examples showing that, using this technique in the Tau Toolbox, we are able to obtain results comparable with those reported in the literature and sometimes even better.

#### **2. The Tau Method**

Let <sup>D</sup> : <sup>E</sup> → <sup>F</sup> be an order *<sup>ν</sup>* differential operator, where <sup>E</sup> and <sup>F</sup> are some function spaces, and let *gi* : <sup>E</sup> → <sup>R</sup>, *<sup>i</sup>* <sup>=</sup> 1, . . . , *<sup>ν</sup>* be *<sup>ν</sup>* functionals representing boundary conditions, so that

$$\begin{cases} \mathcal{D}y = f, \quad f \in \mathbb{F}, \\ g\_i(y) = \phi\_i, \quad i = 1, \dots, \nu \end{cases} \tag{1}$$

is a well posed differential problem.

#### *2.1. The Tau Method Principle*

A particular implementation of the Tau method depends on the choice of an orthogonal basis for <sup>F</sup>. A sequence of orthogonal polynomials {*Pn*(*x*)}<sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> with respect to the weight function *w*(*x*) on a given interval of orthogonality [*a*, *b*] satisfies

$$\langle P\_{\dot{\imath}}, P\_{\dot{\jmath}} \rangle \, = \int\_{a}^{b} w(\mathfrak{x}) P\_{\dot{\imath}}(\mathfrak{x}) P\_{\dot{\jmath}}(\mathfrak{x}) d\mathfrak{x} = \, w\_{\dot{\imath}} \delta\_{\dot{\imath}\acute{\jmath}\prime}$$

where *wi* = *Pi*, *Pi* and *δij* is the Kronecker delta [17].

Let P be the space of algebraic polynomials of any degree and let us suppose that P is dense in F; then, the solution *y* of (1) has a series representation *y* ∼ ∑*j*≥<sup>1</sup> *ajPj*−1. A polynomial approximation of degree *<sup>n</sup>* <sup>−</sup> <sup>1</sup> <sup>∈</sup> <sup>N</sup> is achieved by

$$y\_n = \sum\_{j=1}^n a\_{n,j} P\_{j-1} = \mathcal{P}\_n \mathbf{a}\_{n\prime} \tag{2}$$

where <sup>P</sup>*<sup>n</sup>* = [*P*0, *<sup>P</sup>*1,..., *Pn*−1] and <sup>a</sup>*<sup>n</sup>* = [*an*,1,..., <sup>a</sup>*n*,*n*] *T*.

In the Tau method sense, *yn* is a polynomial satisfying the boundary conditions in (1) and solving the differential equation with a residual *τ<sup>n</sup>* = D*yn* − *f* of maximal order. Thus, the differential problem is reduced to an algebraic one of finding the *n* coefficients *an*,*j*, *j* = 1, . . . , *n* in (2) such that

$$\begin{cases} g\_i(y\_n) = \phi\_{i\prime} \text{ i } = 1, \dots, \nu\_{\prime} \\ \langle P\_{i-1\prime} \mathcal{D} y\_n - f \rangle = 0, \; i = 1, \dots, n - \nu\_{\prime} \end{cases} \tag{3}$$

and so the residual *τ<sup>n</sup>* = O(*Pn*−*ν*).

#### *2.2. Operational Formulation*

For a given *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>, *<sup>n</sup>* <sup>&</sup>gt; *<sup>ν</sup>*, we define the matrix

$$\mathbf{T}\_{n} = \begin{bmatrix} \mathbf{B}\_{\nu \times n} \\ \mathbf{D}\_{(n-\nu)\times n} \end{bmatrix} = (t\_{i,j})\_{i,j=1}^{n} : t\_{i,j} = \begin{cases} g\_{i}(P\_{j-1})\_{\times} \, \, i = 1, \ldots, \nu \\ \left\langle \frac{P\_{i-\nu-1} \, \mathcal{D}P\_{j-1}}{w\_{i-\nu-1}} \right\rangle \, \, i = \nu+1, \ldots, n \end{cases} \tag{4}$$

and the vector

$$\mathbf{b}\_{\mathsf{il}} = \left[ \begin{array}{c} \mathbf{g}\_{\mathsf{v}} \\ \mathbf{f}\_{\mathsf{il}-\mathsf{v}} \end{array} \right] = (b\_{\mathsf{i}})\_{\mathsf{i}=1}^{\mathsf{n}} : \ b\_{\mathsf{i}} = \begin{cases} \boldsymbol{\phi}\_{\mathsf{i}} \, \mathsf{i} = 1, \dots, \mathsf{v} \\ \frac{\langle P\_{\mathsf{i}-\mathsf{v}-1}, f \rangle}{w\_{\mathsf{i}-\mathsf{v}-1}} \, \mathsf{i} = \mathsf{v} + 1, \dots, \mathsf{n}. \end{cases}$$

If, in problem (1), the differential operator D is linear and the *gj* are *ν* linear functionals, then problem (3) can be put in matrix form as

$$
\mathsf{T}\_n \mathsf{a}\_n = \mathsf{b}\_n.
$$

The matrix T*n*, called the Tau matrix, can be evaluated from operational matrices, that is, matrices translating into coefficients vectors the action of a differential operator D in a function *y*.

**Proposition 1.** *Let* <sup>P</sup> = [*P*0, *<sup>P</sup>*1, ...] *be an orthogonal polynomial basis, <sup>y</sup>* <sup>=</sup> <sup>P</sup>a *and* M, N *infinite matrices such that*

$$\mathbf{x}\mathcal{P} = \mathcal{P}\mathbf{M} \quad \text{and} \quad \frac{d}{d\mathbf{x}}\mathcal{P} = \mathcal{P}\mathbf{N}.$$

*Then, for each k* <sup>∈</sup> <sup>N</sup>0*,*

$$\mathbf{x}^k y = \mathcal{P}\mathbb{M}^k \mathbf{a} \quad \text{and} \quad \frac{d^k}{d\mathbf{x}^k} y = \mathcal{P}\mathbb{M}^k \mathbf{a}. \tag{5}$$

**Proof.** For *<sup>k</sup>* <sup>=</sup> 1, the result is true by hypothesis. Now, supposing that (5) is true for a *<sup>k</sup>* <sup>∈</sup> <sup>N</sup>, then

$$\mathbf{x}^{k+1}y = \mathbf{x}(\mathbf{x}^k y) = (\mathbf{x}\mathcal{P})\mathbb{M}^k \mathbf{a} = \mathcal{P}\mathbb{M}^{k+1} \mathbf{a}$$

and

$$\frac{d^{k+1}}{dx^{k+1}}y = \frac{d}{dx}(\frac{d^k}{dx^k}y) = (\frac{d}{dx}\mathcal{P})\mathbb{N}^k\\\mathbf{a} = \mathcal{P}\mathbb{N}^{k+1}\mathbf{a}.$$

ending the proof by induction.

The following result generalizes the algebraic representation from the previous proposition to differential operators.

**Corollary 1.** *Let* <sup>D</sup> : <sup>P</sup>*<sup>n</sup>* → <sup>P</sup>*n*+*<sup>h</sup> be a linear differential operator with polynomial coefficients*

$$\mathcal{D} = \sum\_{i=0}^{\nu} p\_i \frac{d^i}{d\mathbf{x}^{i'}} , \; p\_i \in \mathbb{P}\_{n\_i} \tag{6}$$

*and let h* = max*i*=0,...,*ν*{*ni* − *i*}*.*

*If yn* <sup>=</sup> <sup>P</sup>*n*a*n, then* <sup>D</sup>*yn* <sup>=</sup> <sup>P</sup>*n*+*h*D(*n*+*h*)×*n*a*<sup>n</sup> with*

$$\mathsf{D}\_{(n+h)\times n} = \sum\_{i=0}^{\nu} p\_i(\mathsf{M}) \mathsf{N}^{i}\_{(n-i)\times n\prime}$$

*where pi*(M) = *ni* ∑ *k*=0 *pik*M*<sup>k</sup>* (*n*+*h*)×(*n*−*i*) *when pi* <sup>=</sup> *ni* ∑ *k*=0 *pikx<sup>k</sup> , and* A*<sup>k</sup> <sup>m</sup>*×*<sup>n</sup> denotes the main <sup>m</sup>* × *<sup>n</sup> block of the matrix* (A*p*)*<sup>k</sup> with p* <sup>=</sup> max{*m*, *<sup>n</sup>*}*.*

In [9], the authors discussed the application of this operational formulation of the Tau method to the numerical approximation of eigenvalues defined by differential equations. They proved that, for a differential eigenvalue problem, where in (1)

$$D = \sum\_{r=0}^{t} \lambda^r D\_r$$

and *<sup>λ</sup>* is a parameter, the zeros of det(T*n*(*λ*)) approach the eigenvalues of (1).

#### *2.3. Tau Matrices' Properties*

Given that we are dealing with a general orthogonal polynomial basis, instead of particular cases like Chebyshev or Legendre, we can only make assumptions about general properties of Tau matrices *Tn*. Anyway, we can't expect to have symmetric matrices and, in general, they can be considered sparse but with a low level of sparsity.

Since <sup>P</sup> in Proposition <sup>1</sup> is an orthogonal basis, then M is the tridiagonal matrix with the coefficients of its three term recurrence relation. Therefore, for problems with polynomial coefficients, matrices *pi*(M) of Corollary <sup>1</sup> are banded matrices, with all non-zero elements between the ±*ni* diagonals.

Matrices <sup>N</sup> are always strictly upper triangular and so *pi*(*M*)*N<sup>i</sup>* are *ni* <sup>−</sup> *<sup>i</sup>* upper Hessenberg matrices. The resulting D(*n*−*ν*)×*<sup>n</sup>* block of *Tn* defined in (4) is a general *h* upper Hessenberg matrix.

Moreover, one advantage of the Tau method is its ability to deal with boundary conditions, allowing the treatment of any linear combination of values of *y* and of its derivatives for *gi* in (1). Thus, the *ν* × *n* block *Bν*×*<sup>n</sup>* in *Tn* is usually dense, with its entries *gi*(*Pj*), made by linear combinations of orthogonal polynomial values *Pj*(*xk*) and of its derivatives *<sup>P</sup>*(*l*) *<sup>j</sup>* (*xk*), in prescribed abscissas *xk*.

Assembling those blocks *<sup>B</sup>ν*×*<sup>n</sup>* and D(*n*−*ν*)×*<sup>n</sup>* in *Tn*, we get an *<sup>ν</sup>* + *<sup>h</sup>* upper Hessenberg matrix.

For some problems, whose dependence on the eigenvalues *λ* is verified only in the differential equation, we can use Schur complements to reduce matrix sizes. Considering matrix *Tn* in (4) partitioned as

$$\mathbf{T}\_n = \begin{bmatrix} \mathbf{B}\_1 & \mathbf{B}\_2 \\ \mathbf{D}\_1 & \mathbf{D}\_2 \end{bmatrix}$$

where <sup>B</sup><sup>1</sup> is *<sup>ν</sup>* <sup>×</sup> *<sup>ν</sup>* and the other blocks are partitioned accordingly. If <sup>B</sup><sup>1</sup> is non-singular, then

$$\det(\mathsf{T}\_n) = \det(\mathsf{B}\_1) \det(\mathsf{D}\_2 - \mathsf{D}\_1 \mathsf{B}\_1^{-1} \mathsf{B}\_2)$$

and the problem is reduced to solve

$$\det(\mathbb{C}\_{\mathbb{N}}) = 0, \; \mathbb{C}\_{\mathbb{N}} = \mathbb{D}\_2 - \mathbb{D}\_1 \mathbb{B}\_1^{-1} \mathbb{B}\_{2\prime} \tag{7}$$

reducing to *n* − *ν* the problem dimension. In the worst case, when *B*<sup>1</sup> is singular, we have to work with the *n* × *n* matrix *Tn*.

In the following sections, we illustrate the application of the Tau method to approximate eigenvalues in some classical problems.

#### **3. Problems with Polynomial Coefficients**

Sturm–Liouville problems arise from vibration problems in continuum mechanics. The general form of a fourth order Sturm–Liouville equation is

$$(p(\mathbf{x})y''(\mathbf{x}))'' - (q(\mathbf{x})y'(\mathbf{x}))' + r(\mathbf{x})y(\mathbf{x}) = \lambda \mu(\mathbf{x})y(\mathbf{x}), \ a < \mathbf{x} < b \tag{8}$$

with appropriate initial and boundary conditions, where *<sup>a</sup>* <sup>&</sup>lt; *<sup>b</sup>* <sup>∈</sup> <sup>R</sup>, *<sup>p</sup>*, *<sup>q</sup>*,*r*, and *<sup>μ</sup>* are given piecewise continuous functions, with *p*(*x*) > 0 and *μ*(*x*) > 0. These conditions mean that (8) has an infinite sequence of real eigenvalues, bounded from above, and each one has multiplicity of at most 2 [1].

If *p* and *q* are differentiable functions, it is an elementary task to give (8) the form

$$p(\mathbf{x})y^{(4)}(\mathbf{x}) + 2p'(\mathbf{x})y'''(\mathbf{x}) + (p''(\mathbf{x}) - q(\mathbf{x}))y''(\mathbf{x}) - q'(\mathbf{x})y'(\mathbf{x}) + (r(\mathbf{x}) - \lambda \mu(\mathbf{x}))y(\mathbf{x}) = 0.$$

From this equation, we derive the operational matrix for the general form of the fourth order Sturm–Liouville differential operator associated with (8)

$$\mathbf{D} = p(\mathsf{M})\mathsf{N}^{4} + 2p'(\mathsf{M})\mathsf{N}^{3} + (p''(\mathsf{M}) - q(\mathsf{M}))\mathsf{N}^{2} - q'(\mathsf{M})\mathsf{N} + r(\mathsf{M}) - \lambda\mu(\mathsf{M}).\tag{9}$$

Assuming that coefficients *p*, *q*,*r*, and *μ* are polynomials, or convenient polynomial approximations of the coefficient functions, then the height of this differential operator is well defined as

$$h = \max\{\deg(p) - 4, \deg(q) - 2, \deg(r), \deg(\mu)\}\_{\ell}$$

where deg(.) is the polynomial degree. One consequence of Corollary 1 is that to evaluate the block <sup>D</sup>(*n*−*ν*)×*<sup>n</sup>* in (4) we have to apply (9) with <sup>M</sup> and <sup>N</sup> truncated to its first *<sup>n</sup>* <sup>+</sup> *<sup>h</sup>* lines and columns.

The Tau matrix of a fourth order Sturm–Liouville problem is the *n* × *n* matrix *Tn* = [B4×*n*; <sup>D</sup>(*n*−4)×*n*], where <sup>B</sup>4×*<sup>n</sup>* is the 4 <sup>×</sup> *<sup>n</sup>* matrix representing boundary conditions and <sup>D</sup>(*n*−4)×*<sup>n</sup>* is the first (*<sup>n</sup>* <sup>−</sup> <sup>4</sup>) <sup>×</sup> *<sup>n</sup>* main block of D.

**Example 1.** *Consider the Sturm–Liouville boundary value problem*

$$\begin{cases} y^{(4)}(\mathbf{x}) = \lambda y(\mathbf{x}), \quad 0 < \mathbf{x} < 1, \\ y(0) = y(1) = y'(0) = y''(1) = 0, \end{cases} \tag{10}$$

*whose exact eigenvalues satisfy [1,2]*

$$
\tanh(\sqrt[4]{\lambda}) - \tan(\sqrt[4]{\lambda}) = 0.\tag{11}
$$

In that case D <sup>=</sup> N<sup>4</sup> <sup>−</sup> *<sup>λ</sup>I*, where *<sup>I</sup>* is the identity matrix, and the boundary conditions can be represented by <sup>B</sup>*<sup>n</sup>* = [*v*0; *<sup>v</sup>*1; *<sup>v</sup>*0N; *<sup>v</sup>*1N2], where *<sup>v</sup>*<sup>0</sup> = [1, <sup>−</sup>1, ··· ,(−1)*n*−1] and *<sup>v</sup>*<sup>1</sup> = [1, 1, ··· , 1] are length *n* line vectors with the polynomial base values in the boundary domain.

For each *<sup>n</sup>* <sup>&</sup>gt; 4, <sup>C</sup>*<sup>n</sup>* in (7) is an *<sup>n</sup>* <sup>−</sup> 4 square matrix and its determinant an *<sup>n</sup>* <sup>−</sup> 4 degree polynomial. We use the Matlab function *roots* to find its zeros and we inspect their accuracy by testing if they satisfy relation (11).

In Figure 1, we present | tanh( <sup>4</sup> *λn*,*k*) − tan( <sup>4</sup> *λn*,*k*)| for *k* = 1, ... , 10 the first 10 eigenvalues approximations obtained with *n* = 21, 28, 35 and *n* = 42, with Chebyshev and Legendre bases shifted to [0, 1].

**Figure 1.** | tanh( <sup>4</sup> *λn*,*k*) − tan( <sup>4</sup> *λn*,*k*)|, *k* = 1, 2, . . . , 10, *λn*,*<sup>k</sup>* being the roots of det(T*n*) in Example 1.

**Example 2.** *A very similar problem, presented as the clamped rod problem in [12,18], is*

$$\begin{cases} y^{(4)}(\mathbf{x}) = \lambda y(\mathbf{x}), \quad -1 < \mathbf{x} < 1, \\ y(\pm 1) = y'(\pm 1) = 0. \end{cases} \tag{12}$$

In that case, and whenever we have a symmetric problem and a symmetric base, the matrix *Cn* in (7) has zeros intercalating all non-zero elements. We can reduce the problem dimension, defining two matrices *CO* and *CE* with respectively the odd and the even entries of *Cn*, then det(*Cn*) = det(*CO*) det(*CE*). In that case, since *Cn* is an 4-upper Hessenberg matrix, *CO* and *CE* are 2-upper Hessenberg matrices. The sparsity pattern of those two matrices, in Legendre basis and with *n* = 52, are showed in Figure 2.

**Figure 2.** Sparsity pattern of *CO* and *CE* with *n* = 48 in Legendre base, for Example 2.

The first 14 eigenvalues, evaluated with an 16 × 16 matrix, are presented in [12]. In Table 1, we compare those values with our results in the Legendre basis and with *n* = 16 and *n* = 52. We present values of *λ*52,*k*/*k*4, which allows us to verify that our estimates satisfy the property that the *k*th eigenvalue is proportional to *k*<sup>4</sup> [18].

**Table 1.** Eigenvalues of Example 2 presented in [12] and *λn*,*<sup>k</sup>* with *n* = 20 and *n* = 52 in Legendre basis.


**Example 3.** *Consider the following Sturm–Liouville problem with non-null q and r coefficients*

$$\begin{cases} y^{(4)}(\mathbf{x}) - (\alpha x^2 y'(\mathbf{x}))' + (\beta x^4 - a)y(\mathbf{x}) = \lambda y(\mathbf{x}), & 0 < \mathbf{x} < \mathbf{5}, \\ y(0) = y''(0) = y(\mathbf{5}) = y''(\mathbf{5}) = 0, \end{cases} \tag{13}$$

*with constants <sup>α</sup>*, *<sup>β</sup>* <sup>∈</sup> <sup>R</sup>*.*

The operational matrix (9) for this case is

$$\mathcal{D} = \mathsf{N}^4 - \alpha \mathsf{M}^2 \mathsf{N}^2 - 2\alpha \mathsf{M} \mathsf{N} + \beta \mathsf{M}^4 - (\alpha + \lambda)I$$

and <sup>B</sup>*<sup>n</sup>* = [*v*0; *<sup>v</sup>*0N2; *<sup>v</sup>*1; *<sup>v</sup>*1N2], with the same *<sup>v</sup>*<sup>0</sup> and *<sup>v</sup>*<sup>1</sup> vectors of the previous example.

If *<sup>λ</sup>n*,*<sup>k</sup>* is the *<sup>k</sup>*th root of det(T*n*), and considering ˜ *<sup>δ</sup>n*,*<sup>k</sup>* <sup>=</sup> *<sup>λ</sup>n*,*k*−*λn*−1,*<sup>k</sup> <sup>λ</sup>n*,*<sup>k</sup>* as an estimative of the relative error in *<sup>λ</sup>n*−1,*k*, then *<sup>δ</sup><sup>n</sup>* <sup>=</sup> max*k*=1,...,*<sup>m</sup>* <sup>|</sup> ˜ *δn*,*k*| is an estimative of the maximum relative error in the first *m* eigenvalues of the problem. In Figure 3 left, we present *δn*, with *m* = 6 and with *m* = 8 for Example 3 with *<sup>α</sup>* <sup>=</sup> 0.02 and *<sup>β</sup>* <sup>=</sup> 0.0001 for *<sup>n</sup>* <sup>=</sup> 16, ... , 45. In Figure <sup>3</sup> right, the absolute relative error <sup>|</sup> ˜ *δn*,1| in the lowest eigenvalue is presented for the same *n* values.

**Figure 3.** *δ<sup>n</sup>* = max*k*=1,...,*<sup>m</sup>* ˜ *δn*,*k*, *m* = 6 and *m* = 8, (left) and *δn*,1 (right), in Example 3.

In Table 2, we compare our results with those of [2] for the first six eigenvalues, and of [8] for the first 4, obtained with values *α* = 0.02 and *β* = 0.0001.

**Table 2.** Eigenvalues of Example 3 presented in [8] and [2] and *λn*,*<sup>k</sup>* with *n* = 35.


**Example 4.** *Now, we consider the Orr–Sommerfeld problem*

$$\begin{cases} y^{(4)}(\mathbf{x}) - 2a^2 y'' + a^4 y = iaR[(\mathcal{U} - \lambda)(y'' - a^2 y) - \mathcal{U}'' y], \quad -1 < \mathbf{x} < 1, \\ y(\pm 1) = y'(\pm 1) = 0, \end{cases} \tag{14}$$

*with fixed constants α*, *R and function U.*

The particular case *<sup>U</sup>* = <sup>1</sup> − *<sup>x</sup>*<sup>2</sup> is the Poiseuille flow and, with *<sup>α</sup>* = 1 and *<sup>R</sup>* = 10000 was treated in [3–5,12]. The operational matrix in that case is

$$\mathbf{D} = \mathbb{N}^4 - [(2a^2 + (1 - \lambda)iaR)I - iaR\mathbb{M}^2]\mathbb{N}^2 + (a^4 + (1 - \lambda)ia^3R - 2iaR)I - ia^3R\mathbb{M}^2.$$

Like in Example 2, this is an upper Hessenberg matrix with zeros intercalating its non-zero elements and we can reduce the problem dimension, splitting in two matrices the Schur complement <sup>C</sup>*<sup>n</sup>* of the resulting Tau matrix T*n*. Choosing Chebyshev basis, this is the operational version of the Tau procedure of [5], where the author was confined to eigenvalues associated with symmetric eigenfunctions, which is equivalent to finding the eigenvalues of C*E*.

In [5], the author obtained for *λ*<sup>1</sup> = 0.23752649 + 0.00373967*i* as an 8 decimal places exact value for the most unstable mode of this problem. Working with double-precision arithmetic, we obtain *λ*<sup>1</sup> = 0.2375264889204038 + 0.003739670740170985*i*. This value results with *n* = 58 that is an 29 × 29 matrix C*E*, the same dimension used in [5].

In addition, with *Rc* = 5772.22, the smallest value of *R* for which an unstable eigenmode exists [5], and *α<sup>c</sup>* = 1.02056, we get the results presented in Table 3, together with those of [5].

**Table 3.** Values of *λ*<sup>1</sup> of Example 4 with critical values *Rc* = 5772.22 and *α* = 1.02056.


#### **4. Non-Polynomial Coefficients**

In the previous section, we solved problems in the conditions of Corollary 1, i.e., with differential operators acting in polynomial spaces. In a more general situation, if some of the coefficients *pi* in (6) are non-polynomial functions, then the corresponding matrices *pi*(M) are functions of M instead of polynomial expressions.

If a non-polynomial function *pi* in (6) can be defined implicitly by a differential problem, with polynomial coefficients, then we can first of all use the Tau method to find a polynomial approximation *<sup>p</sup>*˜*<sup>i</sup>* to *pi* and use *<sup>p</sup>*˜*i*(M) to approximate the matrix *pi*(M).

**Example 5.** *Mathieu's equation appears related to wave equations in elliptic cylinders [19]. For an arbitrary parameter q, the problem is to find the values of λ for which non-trivial solutions of*

$$y''(x) + (\lambda - 2q \cos(2x))y(x) = 0\tag{15}$$

*exist with prescribed boundary conditions.*

It can be shown that there exists a countably infinite set of eigenvalues *ar* associated with even periodic eigenfunctions and a countably infinite set of eigenvalues *br* associated with odd periodic eigenfunctions [19]. We are interested in reproducing some of those values given in there.

The operational matrix for problem (15) is

$$
\mathbb{D} = \mathbb{N}^2 + \lambda I - 2q \cos(2\mathbb{M})\dots
$$

Our first step to approximate Mathieu's eigenvalues is to approximate matrix cos(2M). This can be done by, firstly, considering the function *z*(*x*) = cos(2*x*) as the solution of a differential problem, using Tau method to get a polynomial approximation *zn*(*x*) ≈ *z*(*x*). In a second step, the operational matrix D is approximated by

$$
\bar{\mathsf{D}} = \mathsf{N}^2 + \lambda I - 2qz\_n(\mathsf{M}),
$$

and, finally, the last step consists in building the Tau matrix <sup>T</sup>*<sup>n</sup>* and evaluating the zeros of its determinant.

We take integer values *q* = 0, 1, ... , 16 and boundary conditions *y* (−1) = *y* (*π*/2) = 0 to get *ar*(*q*) for even *r*, *y* (−1) = *y*(*π*/2) = 0 for odd *r*, and *y*(−1) = *y*(*π*/2) = 0 to get *br*(*q*) for odd *r* and *y*(−1) = *y* (*π*/2) = 0 for even *r*.

In Figure 4, we show Mathieu eigenvalues *ar*(*q*), *r* = 0, ... , 5 and *br*(*q*), *r* = 1, ... , 6. Values were obtained with a 18th degree polynomial approximation *<sup>z</sup>*<sup>18</sup> <sup>≈</sup> cos(2*x*) and a 36 <sup>×</sup> 36 Tau matrix <sup>T</sup><sup>36</sup> in Chebyshev polynomials.

We can observe, as pointed out in [19], that, for a fixed *q* > 0, we have *a*<sup>0</sup> < *b*<sup>1</sup> < *a*<sup>1</sup> < *b*2 < ··· and that *ar*(*q*), *br*(*q*) approach *r*<sup>2</sup> as *q* approaches zero.

**Figure 4.** Mathieu eigenvalues *ar*(*q*), *r* = 0, ... , 5 and *br*(*q*), *r* = 1, ... , 6, for *q* = 0, ... , 16 in Example 5.

**Example 6.** *Mathieu's equation also appears coupled with a modified Mathieu's equation in systems of differential equations as multi parameter eigenvalues problems. The particular case*

$$\begin{cases} y\_1''(\mathbf{x}\_1) + (\lambda - 2q \cos(2\mathbf{x}\_1)) y\_1(\mathbf{x}\_1) = 0, \ 0 < \mathbf{x}\_1 < \frac{\pi}{2}, \\ y\_1'(0) = y\_1'(\frac{\pi}{2}) = 0, \\ y\_2''(\mathbf{x}\_2) - (\lambda - 2q \cosh(2\mathbf{x}\_2)) y\_2(\mathbf{x}\_2) = 0, \ 0 < \mathbf{x}\_2 < 2, \\ y\_2'(0) = y\_2(2) = 0 \end{cases} \tag{16}$$

*is studied in [6,7] associated with the eigenfrequencies of an elliptic membrane with semi axes α* = cosh(2) *and β* = sinh(2)*.*

To approximate eigenvalues for this problem, we first have to approximate cos(2*x*) and cosh(2*x*) by polynomials. Considering, as in the previous example, *z*<sup>16</sup> ≈ cos(2*x*) the 16th degree Tau solution of

$$\begin{cases} z'''(\mathbf{x}) + 4z(\mathbf{x}) = 0, \; 0 < \mathbf{x} < \frac{\pi}{2},\\ z(0) = 1, \; z(\frac{\pi}{2}) = -1, \end{cases} \tag{17}$$

and *w*<sup>16</sup> ≈ cosh(2*x*) as the same degree Tau solution of

$$\begin{cases} w''(\mathbf{x}) - 4w(\mathbf{x}) = 0, \; 0 < \mathbf{x} < 2, \\ w(0) = 1, \; w'(0) = 0, \end{cases} \tag{18}$$

then

$$\mathbb{D}\_1 = \mathbb{N}^2 + \lambda I - 2qz\_{16}(\mathbb{M})^2$$

and

$$
\tilde{\mathbb{D}}\_2 = \mathbb{N}^2 - \lambda I + 2qw\_{16}(\mathbb{M}),
$$

are matrices approximating the operational matrices associated with differential equations (16).

For each fixed *q*, we define matrices Tau *T*<sup>1</sup> and *T*2, representing Mathieu and modified Mathieu equations, respectively. Defining *an*(*q*) the nth eigenvalue of *T*1, in ascending order, and *Am*(*q*) the mth eigenvalue of *T*2, in descending order, in [6], it was proved that *an*(*q*) and *Am*(*q*) are analytical functions of *q*. Moreover, for each pair (*m*, *n*), it was proved the existence and uniqueness of an intersection point of curves *an*(*q*) and *Am*(*q*). Those intersections identify the eigenmodes of the elliptic membrane.

In Figure 5, we recover, and extend, figures presented in [6] and in [7]. Intersection points of *an*(*q*), the almost vertical curves, and *Am*(*q*), the oblique curves, are the eigenpairs (*a*, *q*) of (16).

**Figure 5.** Mathieu eigenvalues *an*(*q*) and *Am*(*q*), for *q* = 0, ... , 12 in Example 5. Only values *a*(*q*), *A*(*q*) ∈ [−15, 85] are presented.

#### **5. Nonlinear Eigenvalues Problem**

In some differential problems, the eigenvalues can arise in a nonlinear relation with eigenfunctions. Let us consider the following second order problem, related to Weber's equation:

**Example 7.**

$$\begin{cases} y''(\mathbf{x}) + (\lambda + \lambda^2 \mathbf{x}^2) y(\mathbf{x}) = 0, \ -1 < \mathbf{x} < 1, \\ y(-1) = y(1) = 0. \end{cases} \tag{19}$$

The operational matrix corresponding to the differential equation is

$$\mathbb{D} = \mathbb{N}^2 + \lambda I + \lambda^2 \mathbb{M}^2$$

and so det(T*n*) is a polynomial with degree 2(*<sup>n</sup>* <sup>−</sup> <sup>2</sup>) in *<sup>λ</sup>*.

In Table 4, we present the 10 eigenvalues closest to zero, obtained with *n* = 30 and with *n* = 31. We can verify that max*k*=−5,...,5 <sup>|</sup>*λk*(30) <sup>−</sup> *<sup>λ</sup>k*(31)<sup>|</sup> <sup>&</sup>lt; <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>9.

**Table 4.** Eigenvalues of Example 7 with *n* = 30 and *n* = 31. Decimal places presented are those that coincide, until to the first distinct two.


#### **6. Conclusions**

Since the pioneering works of Orzag [5] and Ortiz and Samara [9], the Tau method has been scarcely used to solve differential eigenvalues' problems. With our work, we conclude that the Tau method is a competitive one if we want to evaluate with high accuracy the first eigenvalues, in a large kind of differential problem.

**Author Contributions:** Formal analysis, J.M.A.M. and M.J.R.; Investigation, J.M.A.M. and M.J.R.; Writing—original draft, J.M.A.M. and M.J.R.; Writing—review and editing, J.M.A.M. and M.J.R.

**Funding:** This research was partially supported by CMUP (UID/ MAT/ 00144/ 2019), which is funded by FCT (Portugal) with national (MEC) and European structural funds through the programs FEDER, under the partnership agreement PT2020.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Factors for Marketing Innovation in Portuguese Firms CIS 2014**

#### **Patrícia Monteiro \*, Aldina Correia \* and Vítor Braga \***

School of Management and Technology, Polytechnic of Porto, Rua do Curral, Casa do Curral, Margaride, 4610-156 Felgueiras, Portugal

**\*** Correspondence: 8150194@estg.ipp.pt (P.M.); aic@estg.ipp.pt (A.C.); vbraga@estg.ipp.pt (V.B.)

Received: 28 September 2019; Accepted: 4 November 2019; Published: 22 November 2019

**Abstract:** Globalization, radical and frequent changes as well as the increasing importance of applying knowledge through the efficient implementation of innovation is critical under the current circumstances. Innovation has been the source of businesses competitive advantage, but it is not restricted to technological innovations, and thus marketing innovation also plays a central role. This is a significant topic in the marketing field and not yet deeply analysed in academic research. The main objective of this study is to understand what factors influence marketing innovation and to establish a business profile of firms that innovate or do not in marketing. We used multivariate statistical techniques, such as, multiple linear regression (with the Marketing Innovation Index as dependent variable) and discriminant analysis where the dependent variable is a dummy variable indicating if the firm innovates or not in marketing. The results suggest that there are several factors explaining marketing innovation, although in this study, we find that the factors contributing the most for marketing innovation are: the Organizational Innovation Index, customer and/or user suggestions, and intellectual property rights and licensing (IPRL). Most of the literature has studied these factors separately. This research studied such factors together, and it is clear that both organizational innovation and IPRL play an important role that drives firms to innovate in marketing, which differs from some literature; customer suggestions help in the process of marketing innovation, as some authors argue that customers do not always know what they want until they have it. In parallel, this study proved to be useful in understanding that the different values for the Marketing Innovation Index display no influence on the results, since they were equivalent when a dummy variable (innovated/not innovated in marketing) was used as a dependent variable. In practice, we realize that the factors are useful to clarify what Portuguese firms innovate or not in marketing, with no different results when we the four marketing innovation levels (design, distribution, advertising and price) are considered.

**Keywords:** marketing innovation; CIS 2014; multiple linear regression; discriminant analysis

#### **1. Introduction**

The era of globalization brought radical and frequent changes, as well as a higher recognition of the importance of knowledge through the successful implementation of innovation.

In fact, the changes are constant, and appear in different ways and at an increasing speed. These changes become a challenge for firms which need to, first, identify trends, through well-defined marketing strategies, and subsequently innovate. Innovation, according to the Organization for Economic Cooperation and Development (OECD) and Eurostat [1], requires the implementation of a new or significantly improved product (good or service), or a process, or a new marketing method, or a new organizational method in business practices, within the organization or external relations.

The role of marketing in an organization is very important since it allows increased sales by establishing a long-term relationship with customers. In fact, in addition to financial issues, marketing allows a better understanding of the customer profile leading to co-creation of value.

In order to become more competitive, firms must design new marketing approaches. Marketing innovation is considered by the literature as a non-technological innovation that lacks the same importance as technological innovations (example: product innovation) [2]. According to Mendonça et al. [3] non-technological innovation is an important factor in competitiveness and productivity growth in the economy, specifically in the service industry.

The OECD and Eurostat [1] define marketing innovation as the implementation of a new marketing concept or strategy that differs significantly from existing ones and that has not been previously used.

This study aims to gain a clearer understanding of the role of marketing innovation in Portuguese firms. First, one needs to understand which factors influence and/or impact and secondly to establish a profile of firms regarding marketing innovation.

Marketing innovation is a recent approach with a significant number of publications from 2009 [4]. Therefore, exploring what factors mostly influence Innovation in Marketing is pertinent since the literature contains limited approaches in this regard. According to Correia et al. [5], to achieve the benefits of innovation in terms of economic growth and business competitiveness, it is important to understand its determinants.

Our paper starts with a literature review, that supports the study, followed by the identification of the goals, assumptions and variables used. Subsequently, multivariate analysis of the sample taken from the CIS database (Community Innovation Survey) 2014 was performed and, finally, a connection between the literature and the results of the two statistical techniques: multiple linear regression and discriminant analysis using the SPSS (Statistical Package for the Social Sciences) are assessed.

The results suggest that there are several factors explaining marketing innovation, although, in this study, we find that the factors with higher contribution to marketing innovation are: The Organizational Innovation Index, customer and/or user suggestions and intellectual property rights and licensing (IPRL). In fact, IPRL increase the capacity of marketing innovation in the sense that firms feel more confident in sharing knowledge since they are protected [6]. In turn, the positive contribution of organizational innovation can be explained by the fact that firms increasingly apply improvements in organizational management through innovative marketing measures [7]. Finally, the contribution of customer suggestions and/or users may be related to the fact that they are the consumers of the innovations implemented through products and/or services, so they perceive of what they want to buy [8].

In parallel, this study proved to be useful for understanding that the different values for the Marketing Innovation Index display no influence on our results, since they were equivalent when using a dummy variable (innovated/not innovated in marketing) as dependent variable. In practice, we realize that such factors are useful to classify Portuguese firms that conduct marketing innovation or not, with no different results when one takes the 4 marketing innovation levels (design, distribution, advertising and price) into account.

#### **2. Literature Review and Research Hypothesis**

The main objective of this study is to identify the main factors that influence marketing innovation. Therefore, a survey of scientific production was conducted. Firstly, a literature review was carried out aiming to deepen the knowledge about the subject, promoting ideas for research, identifying gaps in the literature and later reviewing it, considering the methodological purpose of this study.

#### *2.1. Marketing, Innovation and Marketing Innovation Concepts*

Marketing is one of the most important business areas, in addition to promoting the brand of the firm, accelerating sales and business, it involves customers in the dynamics of the firm allowing a better understanding of the value proposition in a creative way. Modern consumers value the experience the brand can provide through marketing dynamics, in contrast to the price of the product and/or service. As a result, the objective of firms is to establish a lasting relationship, giving importance to the client's opinion and involving them in the business [9].

There are numerous definitions of marketing but one of the most relevant is from the American Marketing Association [10] that defines marketing as "the activity, set of institutions, and processes for creating, communicating, delivering and exchanging offerings that have value for customers, clients, partners, and society in general". In turn, Kotler and Armstrong [11] argue that marketing is a social and management process by which individuals and organizations obtain what they need by creating and exchanging value with each other. In a restricted business context, marketing involves building profitable and valuable trading relationships with customers. Thus, the authors conceptualize marketing as the process by which firms create value and build strong relationships with customers, aiming to return this value to them [11].

Both definitions, regardless of the temporal emergence, point **customers** as focus of the firm and, consequently, marketing practices.

Dantas and Moreira [12] point out that is through innovating that one can design irreverent advertising that captivates **customers**, it allows low-price traps by competitors, namely innovation should be part of the DNA of competitive organizations. They also argue that not to innovate does not mean dying but it means being vulnerable to the most direct competitors, showing the importance of innovation to organizations.

**So, what does innovation mean?** In a Yesple way, according to the same authors, *"Innovating is creating new things, doing things di*ff*erently."* The concept of innovation has been approached by several authors and it depends on its application. Table 1 points out some of existing perspectives:


**Table 1.** Innovation definitions | Source: Own Elaboration.

In fact, these definitions are based around 3 main areas: the product (new or improved), processes and organizations (organizational innovation, management or marketing).

The OECD and Eurostat [1] present a structure (Figure 1) that shows innovation as a system and entails the different types of innovation within a firm, the connection of the firm with other organizations and the market demand.

**Figure 1.** The structure of innovation | Source: [1].

The term innovation has been subject of different adjustments due to its importance in the competitive advantage of firms, thus encompassing fields beyond technological improvements, such as marketing management [18].

In fact, marketing and innovation coexist (Figure 2) and Martin [18] argues that successful modern firms are those that successfully combine innovation and marketing. For example, it is essential to firstly identify trends so that innovation can take place at a subsequent stage, considering what the market and customers need. Indeed, in recent years, new ways of collecting information about consumers through innovative marketing programs have allowed firms to reach their target audience more efficiently by using price strategies that were previously not viable [19].

**Figure 2.** Marketing innovation | Source: Own Elaboration.

According to Hume et al. [20], Marketing Innovation develops the marketing philosophy throughout the entire innovation process that goes from the emergence of the idea (based on what the customer needs and meets their needs) to the control of the results associated to the launch of the innovation.

On the other hand, the OECD and Eurostat [1] conceptualize marketing innovation as corresponding to the implementation of a new concept or marketing strategy that differs significantly from the existing ones and that has not been previously used by firms. It requires significant changes in appearance/aesthetic or packaging, placement/distribution, promotion or on product pricing policies. It excludes seasonal changes, and regular or other routine changes in marketing methods. This definition is used throughout the present study to support our dependent variable: "Marketing Innovation Index".

#### *2.2. Marketing Innovation and Product Innovation (Good or Service)*

According to the OECD and Eurostat [1], product innovation corresponds to the introduction of new goods or services or significantly improved ones in the market, about their abilities or inborn abilities, ease of use, components or subsystems.

Currently, the business community strategically uses different types of innovation; one example is marketing and product innovation. The synergy between both seems to be intuitive, but there are few studies in this area. According to Gupta et al. [21], in their research on the relationship between product innovation and marketing, firms operating product innovation tend to rely on marketing as they face uncertainty about how the product will be understood by consumers. On the other hand, Junge et al. [22] concluded that firms that innovate in the product in parallel with marketing achieve a higher productivity growth. In the same line of thought, Ganzer et al. [23] tried to verify the relationship between skilled labour, turnover and number of employees with the amount of investment in product innovation, innovation process, marketing innovation and organizational innovation and concluded that firms that invest in new products or the improvement of existing products tend to innovate in marketing and, consequently, in the management of the firm. Consequently, our hypothesis is:

#### **Hypothesis 1.** *Product innovation contributes positively to marketing innovation.*

Instead, Rebane's [8] study shows different results, since complementarity between product innovation and marketing innovation could not be verified. However, for the services sector the results were different because service providers, when implementing innovation in services and marketing, display greater productivity. Considering these results, the following hypothesis is presented:

**Hypothesis 2.** *Innovation in services contributes positively to marketing innovation.*

#### *2.3. Marketing Innovation and Organizational Innovation*

The OECD and Eurostat [1] show that organizational innovation corresponds to the introduction of a new organizational method in business practices (including knowledge management), in the organization or in the firm's external relations. Higgins [24] mentions that organizational innovation is essential for firms willing to pursue strategic challenges, as they result in improvements in the management of the organization.

The relationship between Organizational Innovation and Marketing Innovation is poorly explored in the literature, but Fleacă et al. [7] studied the extent to which a marketing research process is essential in Organizational Innovation. Their article aimed to understand the importance of using well-defined processes and innovative marketing research, linking the organization's stakeholders to improve work and the overall results of the business.

Marketing research is a sub-process of marketing included in the core processes of a firm, since an effective model of market research allows an organization to more directly and economically commercialize its innovative products, according to current market trends.

The modeling marketing research workflow has drawn valuable results from the APQC (the business process classification framework that allows firms to compare their business processes with other firms [24]). Process classification frameworks developed by the worldwide leader organization in business practice, benchmarking and knowledge management [7].

In this way, a process analyst may be able to structure the necessary steps, such as research objectives, collection, methods and data analysis techniques and information to communicate their findings and implications to those responsible for organizational **decision-making** [7].

Conversely, Ganzer et al. [23] studied the relationship between: product innovation, process, marketing and organization of the knitting industry and concluded that there is a moderate positive correlation between the amount invested in product innovation with the value invested in marketing and organizational innovation. Our hypothesis is:

**Hypothesis 3.** *Innovative changes in organizational forms contribute to the innovation of marketing techniques.*

*2.4. Marketing Innovation and Suggestions of Clients and*/*or Users in Their Innovation Activities and in the Production of Their Innovative Goods or Services*

Clients play a key role in creating and promoting the essential conditions for an innovation project as they allow firms to better understand their needs and desires [25]. Truly, customers are often the consumers of innovations implemented through products and/or services, so they provide important insights about what the market is looking for [8].

Figure 3, proposed by Kilinc et al. [25], reinforces the literature, demonstrating the role of customers in the different stages of the innovation value chain and the impact of the primary roles customers play in the major innovation variables.

**Figure 3.** Role of clients | Prepared by: [7].

In contrast, Cabigiosu and Campagnolo [26] report that customers are a source of relevant knowledge, but cannot be used as the main or exclusive source because (i) on one hand, to develop solutions that address the specific needs of customers, there may be a limited match probability of such solutions to other market opportunities and, (ii) on the other hand, according to Tauber [27], customers often do not realize that they need certain innovative products until they are available in the market.

In fact, cooperation with customers may have a positive effect on firms; however, there are still many costs associated with cooperation with customers and negative aspects to introduce radical or revolutionary changes [8].

The literature points to the importance of customer suggestions in innovations and this article aims to understand, in addition to other factors, how customer suggestions contribute to a non-technological innovation [1], such as marketing innovation. Consequently, the following hypothesis is considered:

#### **Hypothesis 4.** *Customer suggestions contribute to marketing innovation.*

#### *2.5. Marketing Innovation and Intellectual Property Rights and Licensing*

The Oslo manual considers IPRL as requests by firms for patents, European utility models, industrial design rights and trademark registrations [1].

The connection between, for example, registration of brands and product innovation is relatively straightforward and clear, since the marketing of new products is, sometimes, associated with the creation of a new brand to communicate such innovation [3]. As far as marketing innovation is concerned, the connection between them is more complex. According to Mendonça et al. [3] amongst the four types of Innovation in the Oslo manual, only innovation in the promotion of products is not registered, all others can be registered, for example:


Indeed, given the competitiveness of the market, the construction of strong brands may demand marketing innovation, in order to differ from the competition.

In their study, Olaisen and Revang [6] concluded that IPRL increases the innovation capacity, since when IPRL are in place firms feel more confident in sharing knowledge. Also, in this study it was observed that IPRL has no impact on the innovative design of the products. Therefore, the hypothesis for our study is:

#### **Hypothesis 5.** *Firms with intellectual property rights and licensing contribute to marketing innovation.*

#### *2.6. Marketing Innovation and Socioeconomic Characteristics of the Firm*

The success of innovation can be influenced by the type of organization as well as by the characteristics of its employees [21]. The success of marketing practices depends on the creation of an effective multifunctional team that works as a unit creating value for customers [28]. Consequently, the literature points out that firms involved in product innovation and marketing have qualified employees and with the adequate skills [22,29]. This indication of the literature leads to the hypothesis:

#### **Hypothesis 6.** *The academic degree of employees is relevant for marketing innovation.*

Employees are providers of competitive advantage for organizations, and together with turnover they define the size of businesses, i.e., whether the firm is micro, small, medium and/or large. The role of size of the firm is addressed in many studies on Innovation, since it is important to learn about their influence on marketing innovation. Sok et al. [30] state that it is essential, especially, for small and medium enterprises (SMEs) to guarantee the supply of new products, new forms/channels of distribution, to ensure customer satisfaction. The same authors further state that the Yesultaneous implementation of product innovation and marketing combined with qualified employees allow SMEs to be more competitive and achieve better results.

Another aspect leading SMEs to innovate in marketing are circumstantial austerity measures, which do not allow a more permanent support to firms. Therefore, it is imperative that SMEs maximize their internal resources and engage in marketing innovation to better understand the market [31]. In this way the following hypothesis was formulated:

#### **Hypothesis 7.** *The business size has an impact on marketing innovation.*

Larger firms are more likely to innovate in marketing techniques than SMEs due to the investment pressure they experience [32]. Notwithstanding the importance of the size of firms, it's also crucial to study the markets in which they operate. The market action defines the strategic path of the firm, so their decisions consider the type of market in which they choose to operate. This factor may contribute to marketing innovation since firms are currently operating in a globalized environment, which forces them, in competitive terms, to modernize and follow the market trends [33]. Thus, we can consider the hypothesis:

**Hypothesis 8.** *Geographic markets are relevant to firms that innovate in marketing.*

Moreira [33], in his doctoral thesis on the determinants of marketing innovation, conclude that international markets display greater propensity to innovate in Marketing, however, a variable "emerge in national markets" also has a positive and significant effect on innovation marketing. Thus, we can propose the hypothesis:

**Hypothesis 9.** *Internationalization may explain marketing innovation.*

To achieve a broader explanation for the phenomenon of Marketing Innovation, we will try to understand the synergy between firms that belong to the same innovation group in marketing practices. The literature reports that the effects of synergy between firms of the same group and innovation should be treated with caution due to several factors [34].

However, through the study of Entezarkheir and Moshiri [35] it can be understood that mergers can improve incentives for innovation, promoting economies of scale, increasing the capacity to deal with uncertainty, among other things. It was also concluded that mergers are positively and significantly correlated with firm innovation. Therefore, we try to confirm that:

**Hypothesis 10.** *Cooperation between firms of the same group is conducive to an innovative marketing environment.*

Figure 4 and Table 2 summarize the research hypothesis, pointed by literature review and considered in this work.

**Figure 4.** Proposed explanatory model | Source: Own Elaboration.


**Table 2.** Hypothesis synthesis and theoretical support | Source: Own Elaboration.

#### **3. Methodology**

The Community Innovation Survey (CIS) 2014 database was used for the study of Marketing Innovation. The CIS is a notation of the National Statistical System regulated by the European Union aiming to measure and characterize innovation activities in European firms. CIS 2014 covers four types of innovation: product innovation, organizational innovation, process innovation and marketing innovation, being this last innovation the focus of this study. This questionnaire is based on Eurostat guidelines and on the principles of the Oslo manual. In fact, this study, in the literature review, tried to approach the definitions contained in the manual whenever possible.

#### *3.1. Population, Sample and Data Collection*

The data from CIS 2014 database was the basis of our analysis. Our population was all firms located in Portugal over a period of three years, in which the sample initially consists of 8736 firms and after correction by 7083 valid firms. CIS 2014 collected data on the four types of innovation over the period 2012–2014. The database initially contained 187 variables.

#### *3.2. Exploratory Analysis of Data and Study Variables*

Table 3 (Frequency tables and charts in attach) presents a synthesis of the sample used in our study, which was aimed to represent and characterize the data contained in the database. Effectively, it is essential to understand our data before proceeding to multivariate statistical techniques. It can be concluded, from the analysis of Table 3, that most of the Portuguese firms in the sample did not

innovate in product, organization and marketing. Within the firms that innovate in marketing, the most frequent innovation is the innovation in the appearance/aesthetic or in the packaging of the products.


**Table 3.** Exploratory data analysis.

**Table 3.** *Cont.*


For this study 5 variables were created by the authors in order to investigate the validity of the research hypothesis and to provide the interpretation of the results.

Therefore, the following innovation measures were defined:


policies (MKTPRI) considering their sum, i.e., Inov\_Mark = MKTDGP + MKTPDP + MKTPDL + MKTPRI, with values between 0 (no item selected) and 4 (all items selected) [5].

Subsequently, the following variables were also created:

Ö Customer and/or User Suggestions, calculated considering the sum:

Sug\_User = market studies (CLUFEED) + consumer groups (CLUMKT) + discussion groups and interviews (CLUSUR) + surveys of user needs (CLUFOR) + development forums (CLUADA) + development of new goods or services by customers and/or users and that the firm has produced and introduced to the market (CLUDEV)

Intellectual Property Rights and Licensing calculated considering the sum: Prop\_Intellectual = acquired a patent (PROPAT) + required a European utility model (PROEUM) + registered a design right industry (PRODSG) + registered a trademark (PROTM), with values between 0 (no item selected) and 4 (all items selected).

Ö Geographic Markets: M\_GEO = geographic market to sell its goods and/or services the local/regional market (MARLOC) + national market (MARNAT) + market to the European market (MAREUR) + market to other countries not associated with the European Union (MAROTH), with values between 0 (no item selected) and 4 (all items selected).

#### *3.3. Explanatory Variables*

Considering the data analysis and the literature review, a database was built with the variables that could allow a better understanding of Marketing Innovation. Thus, the independent variables pointed out for this multivariate study are summarized in Table 4:


**Table 4.** Explanatory Variables, Expected Signal and Theoretical Support | Source: Own Elaboration.

#### **4. Factors that Influence Marketing Innovation**

Multiple linear regression was used for predicting the value of a variable based on the value of two or more variables [40]. The dependent variable was "Marketing Innovation Index". The variables used to predict the value of the dependent variable are the independent variables: GP—"Belonging to a Group of Firms", Inov\_Org—"Organizational Innovation Index", M\_GEO—"Geographic Markets", Prop\_Intellectual—"Intellectual Property Rights and Licensing", Sug\_Users—"Customer and/or User Suggestions", INPDGD—"Goods Innovation", INPDSV—"Service

Innovation", EMPUD—"% Of Employees with Higher Education", SLO14—"Internationalization" and SIZE14\_COD—"Business Size".

Firstly, we used the forward method in which variables are introduced one by one. The first variable to be introduced is the one with the highest correlation coefficient with the dependent variable Marketing Innovation Index. Subsequently, the variables with the highest coefficient of partial correlation are introduced sequentially [41]. Once the forward analysis was performed it was concluded that the EMPUD, SIZE14\_COD, M\_GEO and GP variables at a significance level of 5% are not significant for the model (Appendix A Table A12). Consequently, hypothesis H6, H7, H8 and H10 are rejected, i.e., the academic level of employees, the business size, the geographic markets and the probability of belonging to a group of firms do not contribute to explain the Index of Marketing Innovation.

After this, linear regression by the stepwise method was conducted in order to eliminate these variables from the model. By the Stepwise method of the 10 independent variables initially considered, only 6 variables were used for the estimation of the model, and the EMPUD, SIZE14\_COD, M\_GEO and GP variables were eliminated as expected (Appendix A Table A13).

Analyzing the summary of the multiple linear regression model (Table 5) we conclude that *Ra*<sup>2</sup> = 0.204 so, approximately 20.4% of the Marketing Innovation Index is explained by the independent variables.


According to the analysis of the ANOVA test (Table 6), *p*-value ≈ 0.000 so, H0 is reject, then we are faced with a highly significant model in which at least one independent variable has a considerable effect on the variation of the dependent variable of marketing innovation.

**Table 6.** Analysis of variance (ANOVA) test | linear regression.


<sup>g</sup> Predictors: (Constant), Inov\_Org, Sug\_Users, Prop\_Intellectual, INPDSV, INPDGD, SLO14.

The variables Organizational Innovation Index (with a standardized coefficient of 0.258), customer suggestions and/or users (with a standardized coefficient of 0.173) and intellectual property and licensing (with a standardized coefficient of 0.147) **are those that contribute** more to explain the Index of Marketing Innovation (Table 7).



<sup>a</sup> Dependent Variable: Inov\_Mark.

$$\text{Then, the adjusted model is: Inov\\_Mark} = 0.314 + 0.258 \text{ Inov\\_Org} + 0.173 \text{ Sug\\_Users} + 0.147$$

$$\text{Prop\\_Intilectual} + 0.088 \text{ INPDSV} + 0.091 \text{ INPDGD} - 0.070 \text{ SLO14}$$

These results are, to some extent, contradictory to the literature review insofar a positive sign was expected for all independent variables (Table 4).

Contrary to expectations (Table 8—NS represent non-significant in the regression model) hypothesis H1 and H10 are rejected, then there is no statistical evidence to consider Product and Organizational Innovation as a factor to Marketing Innovation, as well as H6, H7, H8 and H10.


**Table 8.** Explanatory Variables, Obtained Signal and Theoretical Support | Source: Own Elaboration.

As expected, organizational innovation, customer and/or user suggestions and intellectual property rights and licensing are proved to be important for increasing marketing innovation, as pointed by literature, as well as internationalization, but the latter with opposite sign to the expected. Thus, taking into account our data, the factors promoting marketing innovation are organizational innovation, customer and/or user suggestions and intellectual property rights and licensing, and internationalization are an obstacle to innovate in marketing.

#### *4.1. Testing the Assumptions of Multiple Linear Regression Analysis*

In order to validate the assumptions of the Multiple Linear Regression model, a residual analysis was developed. We analyzed if the residuals follow a normal distribution and had a constant variance (using KS test and dispersion diagrams) and to understand if the residuals are independent, we used the Durbin–Watson test).

Table 9 shows the summary of the multiple linear regression model and the overall adjustment statistics. The Durbin–Watson returned a value of d = 2.002, (approximate to 2) and thus the residuals are not correlated [42]. Consequently, one could proceed with multiple linear regression.


**Table 9.** Durbin–Watson test | linear regression.

<sup>f</sup> Predictors: (Constant), Inov\_Org, Sug\_Users, Prop\_Intellectual, INPDSV, INPDGD, SLO14. <sup>g</sup> Dependent Variable: Inov\_Mark.

The standard predicted and residual values show approximate maximum and minimum values but are not proportional (Table 10).


**Table 10.** Residuals statistics | linear regression.

<sup>a</sup> Dependent Variable: Inov\_Mark.

Through the normal P-P plot of the regression standardized residual in Figure 5, one can conclude that some points are distant from the diagonal. This may indicate that the residuals do not follow a normal distribution.

**Figure 5.** Normal P-P plot of regression standardized residual | linear regression.

In turn, the Scatterplot (Figure 6) presents horizontal lines due to the rounding errors of the values predicted by the regression model for the values of a discrete variable [28].

**Figure 6.** Scatterplot | linear regression.

There is an absence of correlation between independent variables (absence of multicollinearity). Another assumption for linear regression is that none or few collinearities are present. Collinearity occurs when two independent variables are highly correlated [43].

Table 11 shows that no independent variable presents multicollinearity problems since the T is not adjacent to 0 and the Variance Inflation Factor (VIF) displays values below 5.


**Table 11.** Collinearity statistics | linear regression.

In the diagnosis of collinearity (Table 12), it follows that the values of the condition index are not close to 30 and the values themselves are distant from 0.

**Table 12.** Collinearity diagnosis | linear regression.


<sup>a</sup> Dependent Variable: Inov\_Mark.

Most of the proportions of variance, except for a few, show values that are distant from 50%, and may not indicate a multicollinearity problem.

Thus, generically the model meets the multiple linear regression model assumptions, and our model is significant, and there is statistical evidence in the data to consider the conclusions valid.

#### *4.2. Features that Distinguish Firms that Innovate in Marketing*

According to Maroco [44], discriminant analysis is "a dependent multivariate technique used to investigate, evaluate differences between groups and classify entities within groups, based on known discretionary variables." In fact, it is used to discriminate between groups, using a categorical dependent variable and independent interval scale variables [45].

As the discriminant analysis aims to discover the characteristics that distinguish the members of one group from members of a different one, the characteristics of a new individual allows predicting the group it belongs to [45]. We aimed to study which are the characteristics of firms that do not innovate in marketing and those that innovate in marketing. In particular we are interested in comparing the results with the previous analysis, where we consider marketing innovation as an index ranging between 0 (no item selected) and 3 (all items selected). For the analysis we considered Marketing Innovation as a dummy variable being 0 for non-innovative in marketing firms and 1 for innovative in marketing firms.

The non-metric dependent variable marketing innovation consists of 2 mutually exclusive categories. The independent metric variables were selected taking the literature into account. Continuously, Figure 7 presents the metric and non-metric variables under study.

**Figure 7.** Variables metrics and non-metrics | discriminant analysis | Source: Own Elaboration.

The discriminant analysis requires the verification of the following assumptions:


Considering the assumptions, the following tests were performed in order to understand if the discriminant analysis could be performed.

#### 4.2.1. Multivariate Normality

In relation to the first assumption, a K-S test was previously developed and H0 rejected, indicating that the variables do not follow a normal distribution. In order to overcome this problem, we used the central limit theorem which indicates that the larger the size of a sample, the distribution of the mean will be closer to a normal distribution. In this case, the sample contains more than 30 cases so the distribution of the mean can be satisfactorily approximated by a normal distribution [44,46]. The remaining assumptions will then be verified in the output of the discriminant analysis.

4.2.2. Analysis of Variance (ANOVA): Analysis of Differences between Groups

Hypothesis to be tested:


Looking at the test of equality of the groups means, it can be concluded that the Wilks' λ is generally approximate to 1 suggesting that the groups means are equal (Table 13).


**Table 13.** Tests of equality of group means | discriminant analysis.

Concerning the F-test, a small value indicates that when independent variables are considered individually, they do not differ between groups. In turn, the variable Inov\_Org presents a high F suggesting being a variable that is able to differentiate the groups.

For the significance levels, most variables display a *p*-value < 0.05, thus rejecting the null hypothesis, i.e., the means in the two groups, of innovative and non-innovative firms, for all variables are equal.

In contrast with the others, the GP and SIZE14\_COD present *p*-values above 5% (Table 13), indicating that these variables probably do not contribute to the model since the null hypothesis cannot be rejected.

#### 4.2.3. Multivariate Homoskedasticity—Box's M test

Hypothesis to be tested:

H0: Equivalent matrices of variance–covariance for the two groups

H1: Different matrices of variance–covariance for the two groups

Analyzing the Box's M test (Table 14), it is verified that the *p*-value ≈ 0.000 < 0.05 then rejecting H0, i.e., the variance-covariance matrices are the same for the two groups. Therefore, instead of presenting homoscedasticity, required by the analysis, the data shows heteroscedasticity, becoming a limitation of the analysis.


**Table 14.** Box's M test | discriminant analysis.

#### 4.2.4. Absence of Multicollinearity

One of the assumptions of the discriminant analysis is that there is no multicollinearity. Table 15 shows no multicollinearity, i.e., there is no high correlation between the variables, since the values are smaller than 50% presenting, in this case, levels of correlation between variables generally weak [47].

**Table 15.** Pooled within-groups matrices | discriminant analysis.


#### 4.2.5. Stepwise Method

Since, previously, it was verified that the variables GP and SIZE14\_COD do not show significant discriminant power, we used the Stepwise method. This method selects the variables with discriminative capacity, so that the analysis is only done with such variables. In fact, the stepwise method starts without variables and in the following steps variables are added or removed, depending on their discriminative ability [44].

In this analysis the method used for the inclusion/removal of variables was Wilk's λ. Consequently, the variables by this method are included (or removed) according to their inclusion, it greatly decreases (or not) the lambda value [44].

By the stepwise method, only 7 (out of 10) independent variables were considered for the model estimation, and the variables GP, EMPUD and M\_GEO were eliminated (Table 16).


**Table 16.** Variables not considered in the analysis | discriminant analysis.

Table 17 shows that as variables were introduced, the Wilks' λ decreased. Considering that a variable with little tolerance contributes little to the model, Internationalization (SLO14), shows the smallest tolerance (0.866). Prop\_Intellectual and Inov\_Org are two variables that present high tolerance values which indicates that they are the ones that most contribute to the model. However, all variables present high tolerance values, thus showing their relevance to the model and the absence of multicollinearity as they approach 1.


**Table 17.** Variables in the analysis | discriminant analysis.

Hypothesis to be tested:

H0: The group averages are equal

H1: The group averages are different

To understand if the functions are discriminant the Wilks' λ test (Table 18) was performed and it was concluded that one must to reject H0, since the test shows a *p*-value below 5%, i.e., the means of the groups in the function are not equal. Therefore, the functions are discriminant.

**Table 18.** Wilks' λ | discriminant analysis.


#### **To estimate the coe**ffi**cients of the discriminant function, assuring the significance of the functions**

Table 19 shows that there is 1 discriminant function and the eigenvalue attributed to function 1 is 0.134 and represents 100% of the explained variance.

**Table 19.** Eigenvalues | discriminant analysis.


<sup>a</sup> First 1 canonical discriminant functions were used in the analysis.

Regarding canonical correlation, function 1 presents a canonical correlation (0.343)2 corresponding to 0.117649 so, approximately 11.8% of the variance of the groups is explained by the discriminant function 1.

#### **Find the contribution of the variables to the function**

Table 20 allows us to understand which variables contribute to the discriminant function. This indicates that for function 1 the variables that most contribute to distinguish innovative from non-innovative firms are Inov\_Org, Prop\_Intellectual and Sug\_Users. On a different perspective, Organizational Innovation Index, intellectual property and licensing and suggestions of clients and/or users display a positive contribution to be classified in the group of firms that innovate in marketing.


**Table 20.** Standardized canonical discriminant functional coefficients | discriminant analysis.

In turn, the structured matrix (Table 21) allows examining the contribution (ordered by the absolute value) of each variable to the discriminant function, without the effect of collinearity. In this way, organizational innovation is the factor that most positively contributes to function 1, followed by the intellectual property rights and licensing and customer and/or user suggestions. The results, without the effect of collinearity, remained almost equal to Table 20 since the collinearity test resulted negative for the discriminant analysis.

**Table 21.** Structured matrix | discriminant analysis.


<sup>a</sup> This variable not used in the analysis.

#### **Classify cases**

Table 22 allows observing coefficients by Fisher function which, in turn, allow the classification of cases into groups. Thus, it follows that the classification models:

D0 (Don't Innovate in Marketing) = ƺ3.858 + 0.839 \* INPDGD + 0.439 \* INPDSV ƺ 0.203 \* SLO14 + 3.488 SIZE14\_COD + 0.237 \* Inov\_Org ƺ 0.096 \* Prop\_Intellectual + 0.155 \* Sug\_U

D1 (Innovate in at least 1 item of Marketing Innovation) = ƺ4.435 + 1.109 \* INPDGD + 0.681 \* INPDSV ƺ 0.627 \* SLO14 + 3.309 SIZE14\_COD + 0.629 \* Inov\_Org + 0.496 \* Prop\_Intellectual + 0.222 \* Sug\_U


**Table 22.** Classification function coefficients.

#### **Interpretation of the results of discrimination and validation**

Considering Table 23, 63.7% of the cases were correctly classified. In cross-validation, the percentage is almost the same (63.5%) of the original classification.



<sup>a</sup> 63.7% of original grouped cases correctly classified; <sup>b</sup> Cross validation is done only for those cases in the analysis. In cross validation, each case is classified by the functions derived from all cases other than that case; <sup>c</sup> 63.5% of cross-validated grouped cases correctly classified.

In Table 24 a comparison between the linear regression model (MLR) and discriminant analysis (DA) results is presented.

As in MLR, organizational innovation, customer and/or user suggestions and intellectual property rights and licensing are proved to be important for differentiating positively innovative firms from non-innovative ones, as pointed by literature. Internationalization, Yesilar to the MLR analysis, proved to be an obstacle to the promotion of marketing innovation, as much as the business size. In addition, with this DA, product innovation and organizational innovation, proved to be differentiators for distinguish marketing innovative from non-innovative firms, although not significant in differentiating the level of marketing innovation in MLR.


**Table 24.** Linear regression model vs. discriminant analysis.

#### **5. Conclusions**

This study explored marketing innovation in Portuguese firms between 2012 and 2014. Two multivariate statistical techniques were performed to confirm the hypothesis resulting from the literature review, namely: multiple linear regression and discriminant analysis. Both had different objectives. In the first one, it was aimed to understand which factors contributed more to explain the Marketing Innovation Index or marketing innovation level of firms and in the second one it was aimed to define a profile of the firms that do not innovate and innovate in marketing.

Regarding multiple linear regression, it was concluded that the model is significant, and it explains 20.4% of the Marketing Innovation Index. Organizational Innovation Index, customer suggestions and/or users and IPRL were the variables with the greatest contribution to explain Marketing Innovation Index. In fact, about the contribution of IPRL, they can increase the capacity of Marketing Innovation in the sense that firms feel more confident in sharing knowledge because they are protected [6]. In turn, the positive contribution of organizational innovation can be explained by the fact that firms increasingly apply improvements in organizational management through innovative marketing measures [7]. Finally, the contribution of customer suggestions and/or users may be related to the fact that they are the consumers of the innovations implemented through products and/or services, so they have a good perception of what they want to buy [8].

Discriminant analysis reinforced the results obtained through multiple linear regression and proved useful to understand that the different indices of Marketing Innovation display no influence on the results, since they were equivalent when used a dummy variable (innovated/not innovated in marketing). In order to summarise the results of the discriminant analysis, the variables show little discriminative power, however, most of the 7,083 cases (both in the original classification and in the cross validation) were correctly classified. Product Innovation and Organizational Innovation, proved to be important to distinguish innovative from non-innovative in marketing firms, but not relevant to explain the increase of the level of marketing innovation.

Geographic markets, a higher academic level of the employees and belonging to a group of firms do not contribute to explain the Marketing Innovation, thus rejecting the hypothesis initially placed: H6, H8 and H10.

Internationalization, proved to be an obstacle to promotion of marketing innovation, as much as the business size, thus H7 and H9 are verified but with a sign different from expected in the literature.

This study offers some difficulties and limitations, namely that most of the existing literature on Marketing Innovation is considerably recent and that, since 2014 (date of the CIS database) to present, behavioral changes may occur in firms regarding the importance of marketing and innovation itself.

The results show that there is still room for exploring the factors explaining marketing innovation, and this study took some steps in this direction. In fact, a future study may consider other variables, such as cooperation, marketing activities and/or public financial support [33], since the variables used in this study although relevant, are insufficient to fully explain Marketing Innovation. In parallel, it would be relevant to obtain more recent data through primary data, for example, firm surveys, in order to enrich and complement this study or to expand the research.

**Author Contributions:** Cooperative work throughout the manuscript. All authors read and approved the final manuscript.

**Acknowledgments:** We are grateful for ESTG—P. PORTO and CIISESI for the support in the preparation of this manuscript and in the participation in SYMCOMP 2019—4th International Conference on Numerical and Symbolic Computation. Developments and Applications.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Appendix A.1. Frequency Tables*


**Table A1.** Business size | frequency table.


**Table A3.** Geographic markets | frequency table.



**Table A4.** Higher education of employees | frequency table.

**Table A5.** Intellectual property rights and licensing | frequency table.


**Table A6.** Marketing Innovation Index | frequency table.





**Table A8.** Organizational Innovation Index | frequency table.


**Table A10.** Goods innovation | frequency table.




#### *Appendix A.2. Graphics*

**Figure A1.** Internationalization | bar chart.

#### *Appendix A.3. Multiple Linear Regression* | *Tables*


<sup>a</sup> Dependent Variable: Inov\_Mark.

**Table A13.** Variables entered/removed | Stepwise.


<sup>a</sup> Dependent Variable: Inov\_Mark.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Optimal Control of HIV Transmission in Octave/MATLAB**

**Carlos Campos 1,†, Cristiana J. Silva 2,\*,† and Delfim F. M. Torres 2,†**


Received: 1 October 2019; Accepted: 17 December 2019; Published: 19 December 2019

**Abstract:** We provide easy and readable GNU Octave/MATLAB code for the simulation of mathematical models described by ordinary differential equations and for the solution of optimal control problems through Pontryagin's maximum principle. For that, we consider a normalized HIV/AIDS transmission dynamics model based on the one proposed in our recent contribution (Silva, C.J.; Torres, D.F.M. A SICA compartmental model in epidemiology with application to HIV/AIDS in Cape Verde. *Ecol. Complex.* **2017**, *30*, 70–75), given by a system of four ordinary differential equations. An HIV initial value problem is solved numerically using the ode45 GNU Octave function and three standard methods implemented by us in Octave/MATLAB: Euler method and second-order and fourth-order Runge–Kutta methods. Afterwards, a control function is introduced into the normalized HIV model and an optimal control problem is formulated, where the goal is to find the optimal HIV prevention strategy that maximizes the fraction of uninfected HIV individuals with the least HIV new infections and cost associated with the control measures. The optimal control problem is characterized analytically using the Pontryagin Maximum Principle, and the extremals are computed numerically by implementing a forward-backward fourth-order Runge–Kutta method. Complete algorithms, for both uncontrolled initial value and optimal control problems, developed under the free GNU Octave software and compatible with MATLAB are provided along the article.

**Keywords:** numerical algorithms; optimal control; HIV/AIDS model; GNU Octave; open source code for optimal control through Pontryagin Maximum Principle

**MSC:** 34K28; 49N90; 92D30

#### **1. Introduction**

In recent years, mathematical modeling of processes in biology and medicine, in particular in epidemiology, has led to significant scientific advances both in mathematics and biosciences [1,2]. Applications of mathematics in biology are completely opening new pathways of interactions, and this is certainly true in the area of optimal control: a branch of applied mathematics that deals with finding control laws for dynamical systems over a period of time such that an objective functional is optimized [3,4]. It has numerous applications in both biology and medicine [5–8].

To find the best possible control for taking a dynamical system from one state to another, one uses, in optimal control theory, the celebrated Pontryagin's maximum principle (PMP), formulated in 1956 by the Russian mathematician Lev Pontryagin and his collaborators [3]. Roughly speaking, PMP states that it is necessary for any optimal control along with the optimal state trajectory to satisfy the so-called Hamiltonian system, which is a two-point boundary value problem, plus a maximality condition on the Hamiltonian. Although a classical result, PMP is usually not taught to biologists

and mathematicians working on mathematical biology. Here, we show how such scientists can easily implement the necessary optimality conditions given by the PMP, numerically, and can benefit from the power of optimal control theory. For that, we consider a mathematical model for HIV.

HIV modeling and optimal control is a subject under strong current research: see, e.g., Reference [9] and the references therein. Here, we consider the SICA epidemic model for HIV transmission proposed in References [10,11], formulate an optimal control problem with the goal to find the optimal HIV prevention strategy that maximizes the fraction of uninfected HIV individuals with least HIV new infections and cost associated with the control measures, and give complete algorithms in GNU Octave to solve the considered problems. We trust that our work, by providing the algorithms in an open programming language, contributes to reducing the so-called "replication crisis" (an ongoing methodological crisis in which it has been found that many scientific studies are difficult or impossible to replicate or reproduce [12]) in the area of optimal biomedical research. We trust our current work will be very useful to a practitioner from the disease control area and will become a reference in the field of epidemiology for those interested to include an optimal control component in their work.

#### **2. A Normalized SICA HIV/AIDS Model**

We consider the SICA epidemic model for HIV transmission proposed in References [10,11], which is given by the following system of ordinary differential equations:

$$\begin{cases} S'(t) = bN(t) - \lambda(t)S(t) - \mu S(t) \\ I'(t) = \lambda(t)S(t) - (\rho + \phi + \mu)I(t) + aA(t) + \omega \mathcal{C}(t) \\ \mathcal{C}'(t) = \phi I(t) - (\omega + \mu)\mathcal{C}(t) \\ A'(t) = \rho \, I(t) - (a + \mu + d)A(t) \end{cases} \tag{1}$$

The model in Equation (1) subdivides human population into four mutually exclusive compartments: susceptible individuals (*S*); HIV-infected individuals with no clinical symptoms of AIDS (the virus is living or developing in the individuals but without producing symptoms or only mild ones) but able to transmit HIV to other individuals (*I*); HIV-infected individuals under ART treatment (the so-called chronic stage) with a viral load remaining low (*C*); and HIV-infected individuals with AIDS clinical symptoms (*A*). The total population at time *t*, denoted by *N*(*t*), is given by *N*(*t*) = *S*(*t*) + *I*(*t*) + *C*(*t*) + *A*(*t*). Effective contact with people infected with HIV is at a rate *λ*(*t*), given by

$$\lambda(t) = \frac{\beta}{N(t)} \left( I(t) + \eta\_C C(t) + \eta\_A A(t) \right) \tau$$

where *β* is the effective contact rate for HIV transmission. The modification parameter *η<sup>A</sup>* ≥ 1 accounts for the relative infectiousness of individuals with AIDS symptoms in comparison to those infected with HIV with no AIDS symptoms (individuals with AIDS symptoms are more infectious than HIV-infected individuals—pre-AIDS). On the other hand, *η<sup>C</sup>* ≤ 1 translates the partial restoration of immune function of individuals with HIV infection that use ART correctly [10]. All individuals suffer from natural death at a constant rate *μ*. Both HIV-infected individuals with and without AIDS symptoms have access to ART treatment: HIV-infected individuals with no AIDS symptoms *I* progress to the class of individuals with HIV infection under ART treatment *C* at a rate *φ*, and HIV-infected individuals with AIDS symptoms are treated for HIV at rate *α*. An HIV-infected individual with AIDS symptoms *A* that starts treatment moves to the class of HIV-infected individuals *I* and after, if the treatment is maintained, will be transferred to the chronic class *C*. Individuals in the class *C* leave to the class *I* at a rate *ω* due to a default treatment. HIV-infected individuals with no AIDS symptoms *I* that do not take ART treatment progress to the AIDS class *A* at rate *ρ*. Finally, HIV-infected individuals with AIDS symptoms *A* suffer from an AIDS-induced death at a rate *d*.

In the situation where the total population size *N*(*t*) is not constant, it is often convenient to consider the proportions of each compartment of individuals in the population, namely

$$s = \text{S/N}, \quad \text{i} = I/N, \quad \text{c} = \text{C/N}, \quad r = \text{R/N} \dots$$

The state variables *s*, *i*, *c*, and *a* satisfy the following system of differential equations:

$$\begin{cases} \mathbf{s}'(t) = b(1 - \mathbf{s}(t)) - \beta(\dot{\mathbf{i}}(t) + \eta\_{\mathbf{C}}\mathbf{c}(t) + \eta\_{A}a(t))\mathbf{s}(t) + d\,a(t)\,\mathbf{s}(t) \\\\ \mathbf{i}'(t) = \beta\left(\dot{\mathbf{i}}(t) + \eta\_{\mathbf{C}}\mathbf{c}(t) + \eta\_{A}a(t)\right)\mathbf{s}(t) - (\rho + \phi + b)\dot{\mathbf{i}}(t) + aa(t) + \omega c(t) + d\,a(t)\,\mathbf{i}(t) \\\\ \mathbf{c}'(t) = \phi\mathbf{i}(t) - (\omega + b)c(t) + d\,a(t)\,\mathbf{c}(t) \\\\ \mathbf{a}'(t) = \rho\,\mathbf{i}(t) - (a + b + d)a(t) + d\,a^2(t) \end{cases} \tag{2}$$

with *s*(*t*) + *i*(*t*) + *c*(*t*) + *a*(*t*) = 1 for all *t* ∈ [0, *T*].

#### **3. Numerical Solution of the SICA HIV/AIDS Model**

In this section, we consider Equation (2) subject to the initial conditions given by

$$s(0) = 0.6, \quad i(0) = 0.2, \quad c(0) = 0.1, \quad a(0) = 0.1,\tag{3}$$

by the fixed parameter values from Table 1, and by the final time value of *T* = 20 (years).



All our algorithms, developed to solve numerically the initial value problems in Equations (2) and (3), are developed under the free GNU Octave software (version 5.1.0), a high-level programming language primarily intended for numerical computations and is the major free and mostly compatible alternative to MATLAB [13]. We implement three standard basic numerical techniques: Euler, second-order Runge–Kutta, and fourth-order Runge–Kutta. We compare the obtained solutions with the one obtained using the ode45 GNU Octave function.

#### *3.1. Default* ode45 *Routine of GNU Octave*

Using the provided ode45 function of GNU Octave, which solves a set of non-stiff ordinary differential equations with the well known explicit Dormand–Prince method of order 4, one can solve the initial value problems in Equations (2) and (3) as follows:

```
function dy = odeHIVsystem(t,y)
 % Parameters of the model
 mi = 1.0 / 69.54; b = 2.1 * mi; beta = 1.6;
 etaC = 0.015; etaA = 1.3; fi = 1.0; ro = 0.1;
```

```
alfa = 0.33; omega = 0.09; d = 1.0;
% Differential equations of the model
dy = zeros(4,1);
aux1 = beta * (y(2) + etaC * y(3) + etaA * y(4)) * y(1);
aux2 = d * y(4);
dy(1) = b * (1 - y(1)) - aux1 + aux2 * y(1);
dy(2) = aux1 - (ro + fi + b - aux2) * y(2) + alfa * y(4) + omega * y(3);
dy(3) = fi * y(2) - (omega + b - aux2) * y(3);
dy(4) = ro * y(2) - (alfa + b+d- aux2) * y(4);
```
On the GNU Octave interface, one should then type the following instructions:

>> T = 20; N = 100; >> [vT,vY] = ode45(@odeHIVsystem,[0:T/N:T],[0.6 0.2 0.1 0.1]);

Next, we show how such approach compares with standard numerical techniques.

*3.2. Euler's Method*

Given a well-posed initial-value problem

$$\frac{dy}{dt} = f\left(t, y\right) \quad \text{with} \quad y\left(a\right) = a \quad \text{and} \quad a \le t \le b,$$

Euler's method constructs a sequence of approximation points (*t*, *w*) ≈ (*t*, *y* (*t*)) to the exact solution of the ordinary differential equation by *ti*+<sup>1</sup> = *ti* + *h* and *wi*+<sup>1</sup> = *wi* + *h f* (*ti*, *wi*), *i* = 0, 1, ... , *N* − 1, where *t*<sup>0</sup> = *a*, *w*<sup>0</sup> = *α*, and *h* = (*b* − *a*) /*N*. Let us apply Euler's method to approximate each one of the four state variables of the system of ordinary differential equations (Equation (2)). Our odeEuler GNU Octave implementation is as follows:

```
function dy = odeEuler(T)
 % Parameters of the model
 mi = 1.0 / 69.54; b = 2.1 * mi; beta = 1.6;
 etaC = 0.015; etaA = 1.3; fi = 1.0; ro = 0.1;
 alfa = 0.33; omega = 0.09; d = 1.0;
 % Parameters of the Euler method
 test = -1; deltaError = 0.001; M = 100;
 t = linspace(0,T,M+1); h = T / M;
 S = zeros(1,M+1); I = zeros(1,M+1);
 C = zeros(1,M+1); A = zeros(1,M+1);
 % Initial conditions of the model
 S(1) = 0.6; I(1) = 0.2; C(1) = 0.1; A(1) = 0.1;
 % Iterations of the method
 while(test < 0)
   oldS = S; oldI = I; oldC = C; oldA = A;
   for i = 1:M
      % Differential equations of the model
      aux1 = beta * (I(i) + etaC * C(i) + etaA * A(i)) * S(i);
      aux2 = d * A(i);
      auxS = b * (1 - S(i)) - aux1 + aux2 * S(i);
      auxI = aux1 - (ro + fi + b - aux2) * I(i) + alfa * A(i) + omega * C(i);
      auxC = fi * I(i) - (omega + b - aux2) * C(i);
```

```
auxA = ro * I(i) - (alfa + b+d- aux2) * A(i);
    % Euler new approximation
    S(i+1) = S(i)+h* auxS;
    I(i+1) = I(i)+h* auxI;
    C(i+1) = C(i)+h* auxC;
    A(i+1) = A(i)+h* auxA;
 end
 % Absolute error for convergence
 temp1 = deltaError * sum(abs(S)) - sum(abs(oldS - S));
 temp2 = deltaError * sum(abs(I)) - sum(abs(oldI - I));
 temp3 = deltaError * sum(abs(C)) - sum(abs(oldC - C));
 temp4 = deltaError * sum(abs(A)) - sum(abs(oldA - A));
 test = min(temp1,min(temp2,min(temp3,temp4)));
end
dy(1,:) = t; dy(2,:) = S; dy(3,:) = I;
dy(4,:) = C; dy(5,:) = A;
```
Figure 1 shows the solution of the system of ordinary differential equations (Equation (2)) with the initial conditions (Equation (3)), computed by the ode45 GNU Octave function (dashed line) *versus* the implemented Euler's method (solid line). As depicted, Euler's method, although being the simplest method, gives a very good approximation to the behaviour of each of the four system variables. Both implementations use the same discretization knots in the interval [0, *T*] with a step size given by *h* = *T*/100.

**Figure 1.** HIV/AIDS system (Equation (2)) behaviour: GNU Octave versus Euler's method.

Euler's method has a global error (total accumulated error) of *O* (*h*), and therefore, the error bound depends linearly on the step size *h*, which implies that the error is expected to grow in no worse than a linear manner. Consequently, diminishing the step size should give correspondingly greater accuracy to the approximations. Table 2 lists the norm of the difference vector, where each component of this vector is the absolute difference between the results obtained by the ode45 GNU Octave function and our implementation of Euler's method, calculated by the vector norms 1, 2, and ∞.


**Table 2.** Norms 1, 2, and ∞ of the difference vector between ode45 GNU Octave and Euler's results.

#### *3.3. Runge–Kutta of Order Two*

Given a well-posed initial-value problem, the Runge–Kutta method of order two constructs a sequence of approximation points (*t*, *w*) ≈ (*t*, *y* (*t*)) to the exact solution of the ordinary differential equation by *ti*+<sup>1</sup> = *ti* + *h*, *K*<sup>1</sup> = *f* (*ti*, *wi*), *K*<sup>2</sup> = *f* (*ti*+1, *wi* + *hK*1), and *wi*+<sup>1</sup> = *wi* + *h K*<sup>1</sup> + *K*<sup>2</sup> <sup>2</sup> , for each *i* = 0, 1, ... , *N* − 1, where *t*<sup>0</sup> = *a*, *w*<sup>0</sup> = *α*, and *h* = (*b* − *a*) /*N*. Our GNU Octave implementation of the Runge–Kutta method of order two applies the above formulation to approximate each of the four variables of the system in Equation (2). We implement the odeRungeKutta\_order2 function through the following GNU Octave instructions:

```
function dy = odeRungeKutta_order2(T)
 % Parameters of the model
 mi = 1.0 / 69.54; b = 2.1 * mi; beta = 1.6;
 etaC = 0.015; etaA = 1.3; fi = 1.0; ro = 0.1;
 alfa = 0.33; omega = 0.09; d = 1.0;
 % Parameters of the Runge-Kutta (2nd order) method
 test = -1; deltaError = 0.001; M = 100;
 t = linspace(0,T,M+1); h = T / M; h2 = h / 2;
 S = zeros(1,M+1); I = zeros(1,M+1);
 C = zeros(1,M+1); A = zeros(1,M+1);
 % Initial conditions of the model
 S(1) = 0.6; I(1) = 0.2; C(1) = 0.1; A(1) = 0.1;
 % Iterations of the method
 while(test < 0)
   oldS = S; oldI = I; oldC = C; oldA = A;
   for i = 1:M
      % Differential equations of the model
      % First Runge-Kutta parameter
      aux1 = beta * (I(i) + etaC * C(i) + etaA * A(i)) * S(i);
      aux2 = d * A(i);
      auxS1 = b * (1 - S(i)) - aux1 + aux2 * S(i);
      auxI1 = aux1 - (ro + fi + b - aux2) * I(i) + alfa * A(i) + omega * C(i);
      auxC1 = fi * I(i) - (omega + b - aux2) * C(i);
      auxA1 = ro * I(i) - (alfa + b+d- aux2) * A(i);
```

```
% Second Runge-Kutta parameter
     auxS = S(i) + h * auxS1; auxI = I(i)+h* auxI1;
     auxC = C(i) + h * auxC1; auxA = A(i)+h* auxA1;
     aux1 = beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
     aux2 = d * auxA;
     auxS2 = b * (1 - auxS) - aux1 + aux2 * auxS;
     auxI2 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
     auxC2 = fi * auxI - (omega + b - aux2) * auxC;
     auxA2 = ro * auxI - (alfa + b+d- aux2) * auxA;
     % Runge-Kutta new approximation
     S(i+1) = S(i) + h2 * (auxS1 + auxS2);
     I(i+1) = I(i) + h2 * (auxI1 + auxI2);
     C(i+1) = C(i) + h2 * (auxC1 + auxC2);
     A(i+1) = A(i) + h2 * (auxA1 + auxA2);
  end
  % Absolute error for convergence
  temp1 = deltaError * sum(abs(S)) - sum(abs(oldS - S));
  temp2 = deltaError * sum(abs(I)) - sum(abs(oldI - I));
  temp3 = deltaError * sum(abs(C)) - sum(abs(oldC - C));
  temp4 = deltaError * sum(abs(A)) - sum(abs(oldA - A));
  test = min(temp1,min(temp2,min(temp3,temp4)));
end
dy(1,:) = t; dy(2,:) = S; dy(3,:) = I; dy(4,:) = C; dy(5,:) = A;
```
Figure 2 shows the solution of the system of Equation (2) with the initial value conditions in Equation (3) computed by the ode45 GNU Octave function (dashed line) *versus* our implementation of the Runge–Kutta method of order two (solid line). As we can see, Runge–Kutta's method produces a better approximation than Euler's method, since both curves in each plot of Figure 2 are indistinguishable.

**Figure 2.** HIV/AIDS system (Equation (2)): GNU Octave *versus* Runge–Kutta's method of order two.

Runge–Kutta's method of order two (RK2) has a global truncation error of order *O h*2 , and as it is known, this truncation error at a specified step measures the amount by which the exact solution to the differential equation fails to satisfy the difference equation being used for the approximation at that step. This might seems like an unlikely way to compare the error of various methods, since we really want to know how well the approximations generated by the methods satisfy the differential equation, not the other way around. However, we do not know the exact solution, so we cannot generally determine this, and the truncation error will serve quite well to determine not only the error of a method but also the actual approximation error. Table 3 lists the norm of the difference vector between the results from ode45 routine and Runge–Kutta's method of order two results.


**Table 3.** Norms 1, 2, and ∞ of the difference vector between ode45 GNU Octave and RK2 results.

#### *3.4. Runge–Kutta of Order Four*

The Runge–Kutta method of order four (RK4) constructs a sequence of approximation points (*t*, *w*) ≈ (*t*, *y* (*t*)) to the exact solution of an ordinary differential equation by *ti*+<sup>1</sup> = *ti* + *h*, *K*<sup>1</sup> = *f* (*ti*, *wi*), *K*<sup>2</sup> = *f* - *ti* + *<sup>h</sup>* <sup>2</sup> , *wi* <sup>+</sup> *<sup>h</sup>* <sup>2</sup>*K*<sup>1</sup> , *K*<sup>3</sup> = *f* - *ti* + *<sup>h</sup>* <sup>2</sup> , *wi* <sup>+</sup> *<sup>h</sup>* <sup>2</sup>*K*<sup>2</sup> , *K*<sup>4</sup> = *f* (*ti*+1, *wi* + *hK*3), and *wi*+<sup>1</sup> = *wi* + *h* 6 (*K*<sup>1</sup> + 2*K*<sup>2</sup> + 2*K*<sup>3</sup> + *K*4), for each *i* = 0, 1, ... , *N* − 1, where *t*<sup>0</sup> = *a*, *w*<sup>0</sup> = *α*, and *h* = (*b* − *a*) /*N*. Our GNU Octave implementation of the Runge–Kutta method of order four applies the above formulation to approximate the solution of the system in Equation (2) with the initial conditions of Equation (3) through the following instructions:

```
function dy = odeRungeKutta_order4(T)
 % Parameters of the model
 mi = 1.0 / 69.54; b = 2.1 * mi; beta = 1.6;
 etaC = 0.015; etaA = 1.3; fi = 1.0; ro = 0.1;
 alfa = 0.33; omega = 0.09; d = 1.0;
 % Parameters of the Runge-Kutta (4th order) method
 test = -1; deltaError = 0.001; M = 100;
 t = linspace(0,T,M+1);
 h = T / M; h2 = h / 2; h6 = h / 6;
 S = zeros(1,M+1); I = zeros(1,M+1);
 C = zeros(1,M+1); A = zeros(1,M+1);
 % Initial conditions of the model
 S(1) = 0.6; I(1) = 0.2; C(1) = 0.1; A(1) = 0.1;
 % Iterations of the method
 while(test < 0)
   oldS = S; oldI = I; oldC = C; oldA = A;
   for i = 1:M
      % Differential equations of the model
      % First Runge-Kutta parameter
      aux1 = beta * (I(i) + etaC * C(i) + etaA * A(i)) * S(i);
      aux2 = d * A(i);
      auxS1 = b * (1 - S(i)) - aux1 + aux2 * S(i);
      auxI1 = aux1 - (ro + fi + b - aux2) * I(i) + alfa * A(i) + omega * C(i);
```
end

```
auxC1 = fi * I(i) - (omega + b - aux2) * C(i);
     auxA1 = ro * I(i) - (alfa + b+d- aux2) * A(i);
     % Second Runge-Kutta parameter
     auxS = S(i) + h2 * auxS1; auxI = I(i) + h2 * auxI1;
     auxC = C(i) + h2 * auxC1; auxA = A(i) + h2 * auxA1;
     aux1 = beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
     aux2 = d * auxA;
     auxS2 = b * (1 - auxS) - aux1 + aux2 * auxS;
     auxI2 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
     auxC2 = fi * auxI - (omega + b - aux2) * auxC;
     auxA2 = ro * auxI - (alfa + b+d- aux2) * auxA;
     % Fird Runge-Kutta parameter
     auxS = S(i) + h2 * auxS2; auxI = I(i) + h2 * auxI2;
     auxC = C(i) + h2 * auxC2; auxA = A(i) + h2 * auxA2;
     aux1 = beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
     aux2 = d * auxA;
     auxS3 = b * (1 - auxS) - aux1 + aux2 * auxS;
     auxI3 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
     auxC3 = fi * auxI - (omega + b - aux2) * auxC;
     auxA3 = ro * auxI - (alfa + b+d- aux2) * auxA;
     % Fourth Runge-Kutta parameter
     auxS = S(i) + h * auxS3; auxI = I(i)+h* auxI3;
     auxC = C(i) + h * auxC3; auxA = A(i)+h* auxA3;
     aux1 = beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
     aux2 = d * auxA;
     auxS4 = b * (1 - auxS) - aux1 + aux2 * auxS;
     auxI4 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
     auxC4 = fi * auxI - (omega + b - aux2) * auxC;
     auxA4 = ro * auxI - (alfa + b+d- aux2) * auxA;
     % Runge-Kutta new approximation
     S(i+1) = S(i) + h6 * (auxS1 + 2 * (auxS2 + auxS3) + auxS4);
     I(i+1) = I(i) + h6 * (auxI1 + 2 * (auxI2 + auxI3) + auxI4);
     C(i+1) = C(i) + h6 * (auxC1 + 2 * (auxC2 + auxC3) + auxC4);
     A(i+1) = A(i) + h6 * (auxA1 + 2 * (auxA2 + auxA3) + auxA4);
  end
  % Absolute error for convergence
  temp1 = deltaError * sum(abs(S)) - sum(abs(oldS - S));
  temp2 = deltaError * sum(abs(I)) - sum(abs(oldI - I));
  temp3 = deltaError * sum(abs(C)) - sum(abs(oldC - C));
  temp4 = deltaError * sum(abs(A)) - sum(abs(oldA - A));
  test = min(temp1,min(temp2,min(temp3,temp4)));
dy(1,:) = t; dy(2,:) = S; dy(3,:) = I;
dy(4,:) = C; dy(5,:) = A;
```
Figure 3 shows the solution of the initial value problem in Equations (2) and (3) computed by the ode45 GNU Octave function (dashed line) *versus* our implementation of the Runge–Kutta method

of order four (solid line). The results of the Runge–Kutta method of order four are extremely good. Moreover, this method requires four evaluations per step and its global truncation error is *O h*4 .

**Figure 3.** HIV/AIDS system (Equation (2)): GNU Octave versus Runge–Kutta's method of order four.

Table 4 lists the norm of the difference vector between results obtained by the Octave routine ode45 and the 4th order Runge–Kutta method.


**Table 4.** Norms 1, 2, and ∞ of the difference vector between ode45 GNU Octave and RK4 results.

#### **4. Optimal Control of HIV Transmission**

In this section, we propose an optimal control problem that will be solved numerically in Octave/MATLAB in Section 4.2. We introduce a control function *u*(·) in the model of Equation (2), which represents the effort on HIV prevention measures, such as condom use (used consistently and correctly during every sex act) or oral pre-exposure prophylasis (PrEP). The control system is given by

$$\begin{cases} \dot{s}'(t) = b(1 - s(t)) - (1 - u(t))\beta(i(t) + \eta\_C c(t) + \eta\_A a(t))s(t) + d \, a(t) \, s(t) \\\\ \dot{i}'(t) = (1 - u(t))\beta\left(i(t) + \eta\_C c(t) + \eta\_A a(t)\right) s(t) - (\rho + \phi + b)i(t) + a a(t) + \omega c(t) + d \, a(t) \, i(t) \\\\ c'(t) = \phi i(t) - (\omega + b)c(t) + d \, a(t) \, c(t) \\\\ a'(t) = \rho \, i(t) - (a + b + d)a(t) + d \, a^2(t), \end{cases} \tag{4}$$

where the control *u*(·) is bounded between 0 and *u*max, with *u*max < 1. When the control vanishes, no extra preventive measure for HIV transmission is being used by susceptible individuals. We assume that *umax* is never equal to 1, since it makes the model more realistic from a medical point of view.

The goal is to find the optimal value *u*∗ of the control *u* along time, such that the associated state trajectories *s*∗, *i* ∗, *c*∗, and *a*∗ are solutions of the system in Equation (4) in the time interval [0, *T*] with the following initial given conditions:

$$s(0) \ge 0, \quad i(0) \ge 0, \quad c(0) \ge 0, \quad a(0) \ge 0,\tag{5}$$

and *u*∗(·) maximizes the objective functional given by

$$J(\boldsymbol{u}(\cdot)) = \int\_0^T \left(\boldsymbol{s}(t) - \boldsymbol{i}(t) - \boldsymbol{u}^2(t)\right) \, dt,\tag{6}$$

which considers the fraction of susceptible individuals (*s*) and HIV-infected individuals without AIDS symptoms (*i*) and the cost associated with the support of HIV transmission measures (*u*).

The control system in Equation (4) of ordinary differential equations in R<sup>4</sup> is considered with the set of admissible control functions given by

$$\Omega \cap \Omega = \{ u(\cdot) \in L^{\infty}(0, T) \, | \, 0 \le u(t) \le u\_{\text{max}} \, \forall t \in [0, T] \}. \tag{7}$$

We consider the optimal control problem of determining (*s*∗(·), *i* <sup>∗</sup>(·), *c*∗(·), *a*∗(·)) associated to an admissible control *u*∗(·) ∈ Ω on the time interval [0, *T*], satisfying Equation (4) and the initial conditions of Equation (5) and maximizing the cost functional of Equation (6):

$$J(\mu^\*(\cdot)) = \max\_{\Omega} J(\mu(\cdot))\,. \tag{8}$$

Note that we are considering a *L*2-cost function: the integrand of the cost functional *J* is concave with respect to the control *u*. Moreover, the control system of Equation (4) is Lipschitz with respect to the state variables (*s*, *i*, *c*, *a*). These properties ensure the existence of an optimal control *u*∗(·) of the optimal control problem in Equations (4)–(8) (see, e.g., Reference [14]).

To solve optimal control problems, two approaches are possible: direct and indirect. Direct methods consist in the discretization of the optimal control problem, reducing it to a nonlinear programming problem [15,16]. For such an approach, one only needs to use the Octave/MATLAB fmincon routine. Indirect methods are more sound because they are based on Pontryagin's Maximum Principle but less widespread since they are not immediately available in Octave/MATLAB. Here, we show how one can use Octave/MATLAB to solve optimal control problems through Pontryagin's Maximum Principle, reducing the optimal control problem to the solution of a boundary value problem. *Math. Comput. Appl.* **2020**, *25*, 1

#### *4.1. Pontryagin's Maximum Principle*

According to celebrated Pontryagin's Maximum Principle (see, e.g., Reference [3]), if *u*∗(·) is optimal for Equations (4)–(8) with fixed final time *T*, then there exists a nontrivial absolutely continuous mapping <sup>Λ</sup> : [0, *<sup>T</sup>*] <sup>→</sup> <sup>R</sup>4, <sup>Λ</sup>(*t*) = (*λ*1(*t*), *<sup>λ</sup>*2(*t*), *<sup>λ</sup>*3(*t*), *<sup>λ</sup>*4(*t*)), called the *adjoint vector*, such that

$$s' = \frac{\partial H}{\partial \lambda\_1}, \; l' = \frac{\partial H}{\partial \lambda\_2}, \; c' = \frac{\partial H}{\partial \lambda\_3}, \; a' = \frac{\partial H}{\partial \lambda\_4}, \; \lambda\_1' = -\frac{\partial H}{\partial \mathbf{s}}, \; \lambda\_2' = -\frac{\partial H}{\partial i}, \; \lambda\_3' = -\frac{\partial H}{\partial \mathbf{c}}, \; \lambda\_4' = -\frac{\partial H}{\partial a}, \; \lambda\_4' = -\frac{\partial H}{\partial a}, \; \lambda\_4' = -\frac{\partial H}{\partial b}, \; \lambda\_5' = -\frac{\partial H}{\partial c}, \; \lambda\_6' = -\frac{\partial H}{\partial b}, \; \lambda\_7' = -\frac{\partial H}{\partial c}, \; \lambda\_8' = -\frac{\partial H}{\partial c}, \; \lambda\_9' = -\frac{\partial H}{\partial c}$$

where

$$\begin{aligned} H &= H\left(s(t), i(t), c(t), a(t), \Lambda(t), u(t)\right) = s(t) - i(t) - u^2(t) \\ &+ \lambda\_1(t) \left(b(1 - s(t)) - (1 - u(t))\beta(i(t) + \eta\_{\mathbb{C}}c(t) + \eta\_A a(t))s(t) + d \, a(t) \, s(t)\right) \\ &+ \lambda\_2(t) \left((1 - u(t))\beta\left(i(t) + \eta\_{\mathbb{C}}c(t) + \eta\_A a(t)\right)s(t) - \left(\rho + \phi + b\right)i(t) + a a(t) + \omega c(t) + d \, a(t)\, i(t)\right) \\ &+ \lambda\_3(t) \left(\phi i(t) - (\omega + b)c(t) + d \, a(t)\, c(t)\right) \\ &+ \lambda\_4(t) \left(\rho \, i(t) - (a + b + d)a(t) + d \, a^2(t)\right) \end{aligned}$$

is called the *Hamiltonian* and the maximality condition

$$H(\mathbf{s}^\*(t), \mathbf{i}^\*(t), \mathbf{c}^\*(t), \mathbf{a}^\*(t), \Lambda(t), \mathbf{u}^\*(t)) = \max\_{0 \le u \le u\_{\text{max}}} H(\mathbf{s}^\*(t), \mathbf{i}^\*(t), \mathbf{c}^\*(t), \mathbf{a}^\*(t), \Lambda(t), \mathbf{u})$$

holds almost everywhere on [0, *T*]. Moreover, the transversality conditions

$$
\lambda\_i(T) = 0, \qquad i = 1, \ldots, 4, \quad 
$$

hold. Applying the Pontryagin maximum principle to the optimal control problem in Equations (4)–(8), the following theorem follows.

**Theorem 1.** *The optimal control problem of Equations* (4)*–*(8) *with fixed final time T admits a unique optimal solution* (*s*∗(·), *i* <sup>∗</sup>(·), *c*∗(·), *a*∗(·)) *associated to the optimal control u*∗(·) *on* [0, *T*] *described by*

$$u^\*(t) = \min\left\{ \max\left\{ 0, \frac{\beta\left(i^\*(t) + \eta\_C c^\*(t) + \eta\_A a^\*(t)\right) s^\*(t) \left(\lambda\_1(t) - \lambda\_2(t)\right)}{2} \right\}, u\_{\max} \right\},\tag{9}$$

*where the adjoint functions satisfy*

$$\begin{cases} \lambda\_1'(t) = -1 + \lambda\_1(t) \left( b + \left( 1 - u^\*(t) \right) \beta \left( i^\*(t) + \eta\_C c^\*(t) + \eta\_A a^\*(t) \right) - d a^\*(t) \right), \\ \qquad - \lambda\_2(t) \left( 1 - u^\*(t) \right) \beta \left( i^\*(t) + \eta\_C c^\*(t) + \eta\_A a^\*(t) \right) \\ \qquad \lambda\_2'(t) = 1 + \lambda\_1^\*(t) \left( 1 - u^\*(t) \right) \beta \kappa^\*(t) - \lambda\_2(t) \left( \left( 1 - u^\*(t) \right) \beta \kappa^\*(t) - \left( \rho + \phi + s^\*(t) \right) + d a^\*(t) \right) \\ \qquad - \lambda\_3(t) \phi - \lambda\_4(t) \rho, \\ \lambda\_3'(t) = \lambda\_1(t) \left( 1 - u^\*(t) \right) \beta \eta\_C s^\*(t) - \lambda\_2(t) \left( \left( 1 - u^\*(t) \right) \beta \eta\_C s^\*(t) + \omega \right) + \lambda\_3(t) \left( \omega + b - d a^\*(t) \right), \\ \lambda\_4'(t) = \lambda\_1(t) \left( \left( 1 - u^\*(t) \right) \beta \eta\_A s^\*(t) + d s^\*(t) \right) - \lambda\_2(t) \left( \left( 1 - u^\*(t) \right) \beta \eta\_A s^\*(t) + a + d i^\*(t) \right) \\ \qquad - \lambda\_3(t) d c^\*(t) + \lambda\_4(t) \left( a + b + d - 2d a^\*(t) \right), \end{cases} \tag{10}$$

*subject to the transversality conditions λi*(*T*) = 0*, i* = 1, . . . , 4*.*

**Remark 1.** *The uniqueness of the optimal control u*∗ *is due to the boundedness of the state and adjoint functions and the Lipschitz property of the systems in Equations* (4) *and* (10) *(see References [17,18] and references cited therein).*

We implement Theorem 1 numerically in Octave/MATLAB in Section 4.2, and the optimal solution (*s*∗(·), *i* <sup>∗</sup>(·), *c*∗(·), *a*∗(·)) associated to the optimal control *u*∗(·) is computed for given parameter values and initial conditions.

#### *4.2. Numerical Solution of the HIV Optimal Control Problem*

The extremal given by Theorem 1 is now computed numerically by implementing a forwardbackward fourth-order Runge–Kutta method (see, e.g., Reference [19]). This iterative method consists in solving the system in Equation (4) with a guess for the controls over the time interval [0, *T*] using a forward fourth-order Runge–Kutta scheme and the transversality conditions *λi*(*T*) = 0, *i* = 1, ... , 4. Then, the adjoint system in Equation (10) is solved by a backward fourth-order Runge–Kutta scheme using the current iteration solution of Equation (4). The controls are updated by using a convex combination of the previous controls and the values from Equation (9). The iteration is stopped if the values of unknowns at the previous iteration are very close to the ones at the present iteration. Our odeRungeKutta\_order4\_WithControl function is implemented by the following GNU Octave instructions:

```
function dy = odeRungeKutta_order4_WithControl(T)
 % Parameters of the model
 mi = 1.0 / 69.54; b = 2.1 * mi; beta = 1.6;
 etaC = 0.015; etaA = 1.3; fi = 1.0; ro = 0.1;
 alfa = 0.33; omega = 0.09; d = 1.0;
 % Parameters of the Runge-Kutta (4th order) method
 test = -1; deltaError = 0.001; M = 1000;
 t = linspace(0,T,M+1);
 h = T / M; h2 = h / 2; h6 = h / 6;
 S = zeros(1,M+1); I = zeros(1,M+1);
 C = zeros(1,M+1); A = zeros(1,M+1);
 % Initial conditions of the model
 S(1) = 0.6; I(1) = 0.2; C(1) = 0.1; A(1) = 0.1;
 %Vectors for system restrictions and control
 Lambda1 = zeros(1,M+1); Lambda2 = zeros(1,M+1);
 Lambda3 = zeros(1,M+1); Lambda4 = zeros(1,M+1);
 U = zeros(1,M+1);
 % Iterations of the method
 while(test < 0)
   oldS = S; oldI = I; oldC = C; oldA = A;
   oldLambda1 = Lambda1; oldLambda2 = Lambda2;
   oldLambda3 = Lambda3; oldLambda4 = Lambda4;
   oldU = U;
   %Forward Runge-Kutta iterations
   for i = 1:M
     % Differential equations of the model
     % First Runge-Kutta parameter
     aux1 = (1 - U(i)) * beta * (I(i) + etaC * C(i) + etaA * A(i)) * S(i);
     aux2 = d * A(i);
     auxS1 = b * (1 - S(i)) - aux1 + aux2 * S(i);
     auxI1 = aux1 - (ro + fi + b - aux2) * I(i) + alfa * A(i) + omega * C(i);
     auxC1 = fi * I(i) - (omega + b - aux2) * C(i);
     auxA1 = ro * I(i) - (alfa + b+d- aux2) * A(i);
```

```
% Second Runge-Kutta parameter
  auxU = 0.5 * (U(i) + U(i+1));
  auxS = S(i) + h2 * auxS1; auxI = I(i) + h2 * auxI1;
  auxC = C(i) + h2 * auxC1; auxA = A(i) + h2 * auxA1;
  aux1 = (1 - auxU) * beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
  aux2 = d * auxA;
  auxS2 = b * (1 - auxS) - aux1 + aux2 * auxS;
  auxI2 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
  auxC2 = fi * auxI - (omega + b - aux2) * auxC;
  auxA2 = ro * auxI - (alfa + b+d- aux2) * auxA;
  % Third Runge-Kutta parameter
  auxS = S(i) + h2 * auxS2; auxI = I(i) + h2 * auxI2;
  auxC = C(i) + h2 * auxC2; auxA = A(i) + h2 * auxA2;
  aux1 = (1 - auxU) * beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
  aux2 = d * auxA;
  auxS3 = b * (1 - auxS) - aux1 + aux2 * auxS;
  auxI3 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
  auxC3 = fi * auxI - (omega + b - aux2) * auxC;
  auxA3 = ro * auxI - (alfa + b+d- aux2) * auxA;
  % Fourth Runge-Kutta parameter
  auxS = S(i) + h * auxS3; auxI = I(i)+h* auxI3;
  auxC = C(i) + h * auxC3; auxA = A(i)+h* auxA3;
  aux1 = (1 - U(i+1)) * beta * (auxI + etaC * auxC + etaA * auxA) * auxS;
  aux2 = d * auxA;
  auxS4 = b * (1 - auxS) - aux1 + aux2 * auxS;
  auxI4 = aux1 - (ro + fi + b - aux2) * auxI + alfa * auxA + omega * auxC;
  auxC4 = fi * auxI - (omega + b - aux2) * auxC;
  auxA4 = ro * auxI - (alfa + b+d- aux2) * auxA;
  % Runge-Kutta new approximation
  S(i+1) = S(i) + h6 * (auxS1 + 2 * (auxS2 + auxS3) + auxS4);
  I(i+1) = I(i) + h6 * (auxI1 + 2 * (auxI2 + auxI3) + auxI4);
  C(i+1) = C(i) + h6 * (auxC1 + 2 * (auxC2 + auxC3) + auxC4);
  A(i+1) = A(i) + h6 * (auxA1 + 2 * (auxA2 + auxA3) + auxA4);
end
%Backward Runge-Kutta iterations
for i = 1:M
  j = M + 2 - i;
  % Differential equations of the model
  % First Runge-Kutta parameter
  auxU = 1 - U(j);
  aux1 = auxU * beta * (I(j) + etaC * C(j) + etaA * A(j));
  aux2 = d * A(j);
  auxLambda11 = -1 + Lambda1(j) * (b + aux1 - aux2) - Lambda2(j) * aux1;
  aux1 = auxU * beta * S(j);
  auxLambda21 = 1 + Lambda1(j) * aux1 - Lambda2(j) * (aux1 - (ro + fi + b) + ...
                + aux2) - Lambda3(j) * fi - Lambda4(j) * ro;
```

```
aux1 = auxU * beta * etaC * S(j);
auxLambda31 = Lambda1(j) * aux1 - Lambda2(j) * (aux1 + ...
              + omega) + Lambda3(j) * (omega + b - aux2);
aux1 = auxU * beta * etaA * S(j);
auxLambda41 = Lambda1(j) * (aux1 + d * S(j)) ...
              - Lambda2(j) * (aux1 + alfa + ...
              +d*I(j)) - Lambda3(j) * d * C(j) + ...
              + Lambda4(j) * (alfa+b+d-2* aux2);
% Second Runge-Kutta parameter
auxU = 1 - 0.5 * (U(j) + U(j-1));
auxS = 0.5 * (S(j) + S(j-1));
auxI = 0.5 * (I(j) + I(j-1));
auxC = 0.5 * (C(j) + C(j-1));
auxA = 0.5 * (A(j) + A(j-1));
aux1 = auxU * beta * (auxI + etaC * auxC + etaA * auxA);
aux2 = d * auxA;
auxLambda1 = Lambda1(j) - h2 * auxLambda11;
auxLambda2 = Lambda2(j) - h2 * auxLambda21;
auxLambda3 = Lambda3(j) - h2 * auxLambda31;
auxLambda4 = Lambda4(j) - h2 * auxLambda41;
auxLambda12 = -1 + auxLambda1 * (b + aux1 - aux2) - auxLambda2 * aux1;
aux1 = auxU * beta * auxS;
auxLambda22 = 1 + auxLambda1 * aux1 - auxLambda2 * (aux1 - (ro + fi + b) + ...
              + aux2) - auxLambda3 * fi - auxLambda4 * ro;
aux1 = auxU * beta * etaC * auxS;
auxLambda32 = auxLambda1 * aux1 - auxLambda2 * (aux1 + ...
              + omega) + auxLambda3 * (omega + b - aux2);
aux1 = auxU * beta * etaA * auxS;
auxLambda42 = auxLambda1 * (aux1 + d * auxS) ...
              - auxLambda2 * (aux1 + alfa + ...
              +d*auxI) - auxLambda3 * d * auxC + ...
              + auxLambda4 * (alfa+b+d-2* aux2);
% Third Runge-Kutta parameter
aux1 = auxU * beta * (auxI + etaC * auxC + etaA * auxA);
auxLambda1 = Lambda1(j) - h2 * auxLambda12;
auxLambda2 = Lambda2(j) - h2 * auxLambda22;
auxLambda3 = Lambda3(j) - h2 * auxLambda32;
auxLambda4 = Lambda4(j) - h2 * auxLambda42;
auxLambda13 = -1 + auxLambda1 * (b + aux1 - aux2) - auxLambda2 * aux1;
aux1 = auxU * beta * auxS;
auxLambda23 = 1 + auxLambda1 * aux1 ...
              - auxLambda2 * (aux1 - (ro + fi + b) + ...
              + aux2) - auxLambda3 * fi - auxLambda4 * ro;
aux1 = auxU * beta * etaC * auxS;
auxLambda33 = auxLambda1 * aux1 - auxLambda2 * (aux1 + ...
              + omega) + auxLambda3 * (omega + b - aux2);
aux1 = auxU * beta * etaA * auxS;
auxLambda43 = auxLambda1 * (aux1 + d * auxS) ...
        - auxLambda2 * (aux1 + alfa + d * auxI) - auxLambda3*d* auxC + ...
        + auxLambda4 * (alfa+b+d-2* aux2);
```

```
% Fourth Runge-Kutta parameter
  auxU = 1 - U(j-1); auxS = S(j-1);
  auxI = I(j-1); auxC = C(j-1); auxA = A(j-1);
  aux1 = auxU * beta * (auxI + etaC * auxC + etaA * auxA);
  aux2 = d * auxA;
  auxLambda1 = Lambda1(j) - h * auxLambda13;
  auxLambda2 = Lambda2(j) - h * auxLambda23;
  auxLambda3 = Lambda3(j) - h * auxLambda33;
  auxLambda4 = Lambda4(j) - h * auxLambda43;
  auxLambda14 = -1 + auxLambda1 * (b + aux1 - aux2) ...
                - auxLambda2 * aux1;
  aux1 = auxU * beta * auxS;
  auxLambda24 = 1 + auxLambda1 * aux1 ...
                - auxLambda2 * (aux1 - (ro + fi + b) + ...
                + aux2) - auxLambda3 * fi - auxLambda4 * ro;
  aux1 = auxU * beta * etaC * auxS;
  auxLambda34 = auxLambda1 * aux1 - auxLambda2 * (aux1 + ...
                + omega) + auxLambda3 * (omega + b - aux2);
  aux1 = auxU * beta * etaA * auxS;
  auxLambda44 = auxLambda1 * (aux1 + d * auxS) ...
                - auxLambda2 * (aux1 + alfa + ...
                +d*auxI) - auxLambda3 * d * auxC + ...
                + auxLambda4 * (alfa+b+d-2* aux2);
  % Runge-Kutta new approximation
  Lambda1(j-1) = Lambda1(j) - h6 * (auxLambda11 + ...
                 + 2 * (auxLambda12 + auxLambda13) + auxLambda14);
  Lambda2(j-1) = Lambda2(j) - h6 * (auxLambda21 + ...
                 + 2 * (auxLambda22 + auxLambda23) + auxLambda24);
  Lambda3(j-1) = Lambda3(j) - h6 * (auxLambda31 + ...
                 + 2 * (auxLambda32 + auxLambda33) + auxLambda34);
  Lambda4(j-1) = Lambda4(j) - h6 * (auxLambda41 + ...
                 + 2 * (auxLambda42 + auxLambda43) + auxLambda44);
end
% New vector control
for i = 1:M+1
  vAux(i) = 0.5 * beta * (I(i) + etaC * C(i) + ...
            + etaA * A(i)) * S(i) * (Lambda1(i) - Lambda2(i));
  auxU = min([max([0.0 vAux(i)]) 0.5]);
  U(i) = 0.5 * (auxU + oldU(i));
end
% Absolute error for convergence
temp1 = deltaError * sum(abs(S)) - sum(abs(oldS - S));
temp2 = deltaError * sum(abs(I)) - sum(abs(oldI - I));
temp3 = deltaError * sum(abs(C)) - sum(abs(oldC - C));
temp4 = deltaError * sum(abs(A)) - sum(abs(oldA - A));
temp5 = deltaError * sum(abs(U)) - sum(abs(oldU - U));
temp6 = deltaError * sum(abs(Lambda1)) - sum(abs(oldLambda1 - Lambda1));
temp7 = deltaError * sum(abs(Lambda2)) - sum(abs(oldLambda2 - Lambda2));
temp8 = deltaError * sum(abs(Lambda3)) - sum(abs(oldLambda3 - Lambda3));
```

```
temp9 = deltaError * sum(abs(Lambda4)) - sum(abs(oldLambda4 - Lambda4));
  test = min(temp1,min(temp2,min(temp3,min(temp4, ...
         min(temp5,min(temp6,min(temp7,min(temp8,temp9))))))));
end
dy(1,:) = t; dy(2,:) = S; dy(3,:) = I;
dy(4,:) = C; dy(5,:) = A; dy(6,:) = U;
disp("Value of LAMBDA at FINAL TIME");
disp([Lambda1(M+1) Lambda2(M+1) Lambda3(M+1) Lambda4(M+1)]);
```
For the numerical simulations, we consider *u*max = 0.5, representing a lack of resources or misuse of the preventive HIV measures *u*(·), that is, the set of admissible controls is given by

$$\Omega = \{ u(\cdot) \in L^{\infty}(0, T) \, | \, 0 \le u(t) \le 0.5, \,\forall t \in [0, T] \}\tag{11}$$

with *T* = 20 (years). Figures 4 and 5 show the numerical solution to the optimal control problem of Equations (4)–(8) with the initial conditions of Equation (3) and the admissible control set in Equation (11) computed by our odeRungeKutta\_order4\_WithControl function. Figure 6 depicts the extremal control behaviour of *u*∗.

**Figure 4.** Optimal state variables for the control problem in Equations (4)–(8) subject to the initial conditions in Equation (3) and the admissible control set in Equation (11) *versus* trajectories without control measures.

**Figure 5.** Comparison: solutions to the initial value problem in Equations (2)–(3) *versus* solutions to the optimal control problem in Equations (4)–(8) subject to initial conditions in Equation (3) and the admissible control set in Equation (11).

**Figure 6.** Optimal control *u*∗ for the HIV optimal control problem in Equations (4)–(8) subject to the initial conditions in Equation (3) and the admissible control set in Equation (11).

#### **5. Conclusions**

The paper provides a study on numerical methods to deal with modelling and optimal control of epidemic problems. Simple but effective Octave/MATLAB code is fully provided for a recent model proposed in Reference [10]. The given numerical procedures are robust with respect to the parameters: we have used the same values as the ones in Reference [10], but the code is valid for other values of the parameters and easily modified to other models. The results show the effectiveness of optimal control theory in medicine and the usefulness of a scientific computing system such as GNU Octave: using the control measure as predicted by Pontryagin's maximum principle and numerically computed by our Octave code, one sees that the number of HIV/AIDS-infected and -chronic individuals diminish and, as a consequence, the number of susceptible (not ill) individuals increase. We trust our paper will be very useful to a practitioner from the disease control area. Indeed, this work has been motivated by many emails we received and continue to receive, asking us to provide the software code associated to our research papers on applications of optimal control theory in epidemiology, e.g., References [20,21].

**Author Contributions:** Each author equally contributed to this paper, read and approved the final manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by the Portuguese Foundation for Science and Technology (FCT) within projects UID/MAT/04106/2019 (CIDMA) and PTDC/EEI-AUT/2933/2014 (TOCCATTA) and was cofunded by FEDER funds through COMPETE2020—Programa Operacional Competitividade e Internacionalização (POCI) and by national funds (FCT). Silva is also supported by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in numbers 4–6 of article 23 of the Decree-Law 57/2016 of August 29, changed by Law 57/2017 of July 19.

**Acknowledgments:** The authors are sincerely grateful to four anonymous reviewers for several useful comments, suggestions, and criticisms, which helped them to improve the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Isogeometric Analysis for Fluid Shear Stress in Cancer Cells**

#### **José A. Rodrigues**

Centro de Investigação em Modelação e Optimização de Sistemas Multifuncionais (CIMOSM), Instituto Superior de Engenharia de Lisboa (ISEL), Instituto Politécnico de Lisboa, Rua Conselheiro Emídio Navarro, 1959-007 Lisboa, Portugal; jose.rodrigues@isel.pt

Received: 20 December 2019; Accepted: 2 April 2020; Published: 3 April 2020

**Abstract:** The microenvironment of the tumor is a key factor regulating tumor cell invasion and metastasis. The effects of physical factors in tumorigenesis is unclear. Shear stress, induced by liquid flow, plays a key role in proliferation, apoptosis, invasion, and metastasis of tumor cells. The mathematical models have the potential to elucidate the metastatic behavior of the cells' membrane exposed to these microenvironment forces. Due to the shape configuration of the cancer cells, Non-uniform Rational B-splines (NURBS) lines are very adequate to define its geometric model. The Isogeometric Analysis allows a simplified transition of exact CAD models into the analysis avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques. In this work, we use an isogeometric analysis to model the fluid-generated forces that tumor cells are exposed to in the vascular and tumor microenvironments, in the metastatic process. Using information provided by experimental tests in vitro, we present a suite of numerical experiments which indicate, for standard configurations, the metastatic behavior of cells exposed to such forces. The focus of this paper is strictly on geometrical sensitivities to the shear stress' exhibition for the cell membrane, this being its innovation.

**Keywords:** Darcy; Brinkman; incompressible; isogeometric analysis; shear stress; interstitial flow; cancer; NURBS

#### **1. Introduction**

The formation of a secondary tumor at a site distant from the the primary tumor is known as metastasis by a cancerous tumor. To initiate the metastatic spread of cancer through the bloodstream, tumor cells must transit through microenvironments of dramatically varying physical forces. Cancer cells are able to migrate through the stroma, intravasate through the endothelium into the blood or lymphatic vessels, to flow within the vessels, to extravasate from the vessel through the endothelium and colonize in tissue at a secondary site [1].

In soft tissues, cancer cells are exposed to mechanical forces due to fluid shear stress, hydrostatic pressure, tension and compression forces. Fluid shear stress is one of the most important forces that cells are exposed to, and its effects on blood cells, endothelial cells, smooth muscle cells, and others have been extensively studied. However, much less is known about fluid shear stress effects on tumor cells. Cancer cells experience two main kinds of fluid shear stress: stresses generated by blood flow in the vascular microenvironment, and those generated by interstitial flows in the tumor microenvironment [2]. Stresses generated by interstitial and blood flows could contribute to the metastatic process by enhancing tumor cell invasion and circulating tumor cell adhesion to blood vessels, respectively. However, it is difficult to predict tumor cell behavior to such forces and it is difficult to experimentally measure such flows in the tumor microenvironment [3].

In this work, we apply Mathematical and Mechanical processes to analyze the metastatic behavior of the cells' membrane exposed to microenvironement forces. Providing an Isogeometric Analysis (IGA) [4] computational method to model and predict how cancer cells respond to such forces, we allow for new insight and new decision tools for medical problems.

The biological fluid study at microenvironments, in vivo or in vitro, is a recent topic following the microfluidic technology and mechanical methodology methods developments [5–8], with the ethical and budget restrictions, increasing the evidence that fluid shear stress is an essential factor affecting fluid mechanics [1].

The current work serves as an introduction to this line of research on mathematical modeling to help us understand the omics data produced by experimental techniques and to bridge the gap between the developments of technologies and systemic modelling of the biological process in cancer research.

Although it is accepted that the effectiveness of IGA methods is well-established for problems with complex geometries, its effect in biologic models has not been extensively studied. Indeed, at present, general research in this subject is still in this infancy for the case studies presented [9,10].

First, we introduce and describe the methodology. We then provide a mathematical and numerical analysis overview for the interstitial flow governing model and in order to estimate the fluid shear stress on cells in tissues.

Standard shape for cancer cells is very irregular, for this reason Non-uniform Rational B-splines (NURBS) lines are very adequate to define its geometric model. The Isogeometric Analysis allows for simplified transition of exact CAD models into the analysis, avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques. For the exhibit models, NURBS provides a high-quality geometric model, quite analogous to a physical model, and also defines the basis of the discrete space in which the partial differential equation solution is approximated with great accuracy per degree of freedom, as referred to in Section 3.

We conclude with simulations to predict the effect of fluid shear stress on cancer cells for several scenarios of microenvironment tumor implantation.

In this paper, we focus on the effects of shear stress, induced by liquid flow, on apoptosis, invasion, and metastasis of tumor cells, following the theoretical [11] and numerical [12] work for a Stokes flow.

#### **2. Materials and Methods**

#### *2.1. Model Geometry*

Tumor cells have a large variety in shape and size. Therefore, several different geometries will be considered in this study, always considering these as standard cells. We will consider bidimensional projections (as shown in Figure 24) of the model shown in Figure 1 to inspire the geometrical model for the presented cases.

**Figure 1.** Tridimensional tumor's geometrical model [7]. Bidimensional projections of this model will be considered for the study cases.

For these geometrical models, Non-uniform Rational B-splines (NURBS) lines are very adequate to define it due to simplified transition of exact CAD models into the analysis, avoiding the geometrical discontinuities of the traditional Galerkin traditional techniques.

#### *2.2. Governing Equations*

Average equations describing viscous flow through a porous region are of great theoretical and practical interest. At the microscale level, the Stokes equations apply and provide a complete description of the entire flow field. However, Darcy's law, formally derived by performing appropriate volume averages of Stokes equations, is applicable. The qualitative difference between these two flow descriptions motivate Brinkman to suggest a general equation that interpolates between the Stokes equation and Darcy's law [13]. We introduce and describe the mathematical model in the next section.

#### **3. Mathematical Formulation**

#### *3.1. The Darcy–Brinkman Equation*

Darcy's law is a phenomenologically derived constitutive equation that describes the flow of a fluid through a porous medium [9,10]. This law, as an expression of conservation of momentum, was determined experimentally by Darcy. It has since been derived from the Navier–Stokes equations via homogenization. It is analogous to Fourier's law in the field of heat conduction.

$$
\vec{u} = -\frac{k}{\nu} \nabla p \tag{1}
$$

where *k* is the permeability of the medium, ∇*p* is the pressure gradient vector, *ν* is the viscosity of the fluid and *u* is the fluid's average velocity through a plane region represented by Ω. Region Ω is supposed regular enough to ensure the later theoretical regularity for *u* and *p*.

To account for interstitial flows between boundaries, Brinkman has developed a second-order term, taking into account no-slip boundary conditions cells' membrane (Figure 2).

**Figure 2.** Outline of cells membrane interstitial zone and no-slip boundary conditions representation.

For governing interstitial flow between boundaries we will consider the Darcy–Brinkman equation (2)

$$-\nu\Delta\vec{u} + \frac{\nu}{k}\vec{u} + \nabla p = 0\tag{2}$$

The mass balance equation for a steady state incompressible fluid is that the divergence of the fluid is zero

$$\nabla \cdot \vec{u} = 0 \tag{3}$$

To apply the boundary conditions we decompose the boundary of the region Ω in two non-overlaping regions, as shown in Figure 3, *∂*Ω = Γ*<sup>T</sup>* ∪ Γ*I*, the boundary of the tumor region and the boundary of the interstitial region considered, respectively. We consider homogeneous Dirichlet boundary condition for the velocity over the cells membrane and non-homogeneous on the interstitial region boundaries, i.e., *u* = (*u*<sup>0</sup> *<sup>x</sup>*, *u*<sup>0</sup> *<sup>y</sup>*) and parameters *ν* and *k* given by Table 1, for the cases under study.

**Figure 3.** Outline of boundary decomposition *∂*Ω = Γ*<sup>T</sup>* ∪ Γ*I*.

To facilitate understanding, we write the component-wise formulation in two dimensions. Let *u* = (*ux*, *uy*). Then Equations (2) and (3) consists of three equations

$$\begin{array}{rcl} -\nu\Delta u\_{\overline{x}} + \frac{\nu}{k}u\_{\overline{x}} + \frac{\partial p}{\partial x} & = & 0\\ -\nu\Delta u\_{\overline{y}} + \frac{\nu}{k}u\_{\overline{y}} + \frac{\partial p}{\partial y} & = & 0 \end{array} \tag{4}$$

$$\begin{array}{rcl} \frac{\partial u\_{\overline{x}}}{\partial y} + \frac{\partial u\_{\overline{y}}}{\partial y} & = & 0 \end{array}$$

We are interested to find *u* = (*ux*, *uy*) and *p* solution of the problem defined by Equation (4) and the Dirichlet boundary conditions showed at Table 1.

With the solution, *u* = (*ux*, *uy*) we evaluate the fluid shear stress. For bidimensional Newtonian fluid the Stress is written [14]

$$T = -\begin{bmatrix} p & 0\\ 0 & p \end{bmatrix} + 2\nu \begin{bmatrix} \frac{\partial u\_x}{\partial x} & \frac{1}{2} \left( \frac{\partial u\_x}{\partial y} + \frac{\partial u\_y}{\partial x} \right) \\\ \frac{1}{2} \left( \frac{\partial u\_y}{\partial x} + \frac{\partial u\_x}{\partial y} \right) & \frac{\partial u\_y}{\partial y} \end{bmatrix} \tag{5}$$

and the Shear Stress is described by

$$T\_{xy} = \nu \frac{1}{2} \left( \frac{\partial u\_x}{\partial y} + \frac{\partial u\_y}{\partial x} \right) . \tag{6}$$

#### *3.2. Isogeometric Discretization*

In order to represent complex shapes, the use of polynomials or rational segments may often be inadequate or imprecise. On the other hand, B-spline and NURBS functions enjoy some major advantages that make them extremely convenient for complex geometrical representations.

The main idea behind the isogeometric approach [4] is to discretize the unknowns of the problem with the same set of basis functions that CAD employs for the construction of geometries. Let *p* be the prescribed degree and *n* control points, we define by

$$\Xi = \{t\_1, \dots, t\_n, t\_{n+p+1}\} \tag{7}$$

the knots vector, with *ti* ∈ [0, 1], *i* ∈ {1, ··· , *n* + *p* + 1}. Cox-De Boor's formula [15] defines *n* one-dimensional B-spline basis functions recursively as

$$B\_{i,0}\left(t\right) = \begin{cases} 1 & t\_i \le t < t\_{i+1} \\ 0 & \text{otherwise} \end{cases} \tag{8}$$

$$B\_{i,p}\left(t\right) = \frac{t - t\_i}{t\_{i+p} - t\_i} B\_{i,p-1}\left(t\right) + \frac{t\_{i+p+1} - t}{t\_{i+p+1} - t\_{i+1}} B\_{i+1,p-1}\left(t\right)$$

for *i* ∈ {1, ··· , *n*}

As referred to in [12] the support of a B-spline of degree *p* is always *p* + 1 knot spans and, as a consequence, each *p*-th degree function has *p* − 1 continuous derivatives across the element boundaries, or across the knots, if they are not repeated. Repetition of knots can be exploited to prescribe the regularity.

NURBS of degree *p* are defined as rational B-splines, associating to each B-spline function a weight *wi*

$$N\_{i,p}\left(t\right) = \frac{w\_i B\_{i,p}\left(t\right)}{\sum\_{j} w\_j B\_{j,p}\left(t\right)}\tag{9}$$

with *wi* called the weight parameter. Geometries in the projective space may be described by using the concept of homogeneous coordinates, which are frequently denoted as weights *wi*. A weighted polynomial B-spline geometry of IR*d*+<sup>1</sup> is obtained by first multiplying its control point data with the homogeneous coordinates. For values *wi* > 1, the object moves toward the control polygon, whereas for weights smaller than one, the influence of the control point on the geometry decreases. Control points with *wi* = 0 do not affect the geometric object at all. If all *wi* are equal to one, the NURBS basis simplifies to the polynomial B-spline basis [15] and allows for the partition of unity property.

As in [12] we define bidimensional B-splines and NURBS using a tensor product approach. Considering **Ξ** = Ξ<sup>1</sup> × Ξ<sup>2</sup> the knot vectors, **p** the degrees and **n** the number of basis functions the bivariate B-spline is given by

$$N\_{\mathbf{i},\mathbf{p}}\left(\mathbf{t}\right) = N\_{\mathbf{i}\_1,p\_1}\left(t\_1\right)N\_{\mathbf{i}\_2,p\_2}\left(t\_2\right) \tag{10}$$

where **t** = (*t*1, *t*2), **p** = (*p*1, *p*2), **n** = (*n*1, *n*2) and **i** = (*i*1, *i*2) is a multi-index in the set *i*<sup>1</sup> ∈ {1, ··· , *n*1} and *i*<sup>2</sup> ∈ {1, ··· , *n*2}.

It is straightforward to notice that there is a parametric Cartesian mesh T*<sup>h</sup>* associated with **Ξ**. The knot vectors partitioning the parametric domain [0, 1] <sup>2</sup> into parallelograms. For each element *Q* ∈ T*<sup>h</sup>* we associate a parametric mesh size *hQ* = diam(*Q*), and *h* = max *hQ*, *Q* ∈ T*<sup>h</sup>* ! .

In the following, we refer to the basis functions indicating the global index, and we will denote by *Sh* (**Ξ**) the bivariate B-spline space spanned by the basis functions **Ni**, 1 ≤ *i* ≤ *n*. For convenience we also use the notation *<sup>S</sup>p*1,*p*<sup>2</sup> *<sup>α</sup>*1,*α*<sup>2</sup> (**Ξ**) to designate the associated space of splines of order *<sup>p</sup>*<sup>1</sup> in the *<sup>x</sup>* direction, *p*<sup>2</sup> in the *y* direction, and smoothness *α*<sup>1</sup> and *α*<sup>2</sup> respectively. Notice that *α<sup>i</sup>* is determined by the knot vector Ξ*i*, spanned by the basis functions of degree *pd* and regularity *α<sup>d</sup>* for each direction *d* ∈ {1, 2}.

A NURBS curve is defined by a set of control points **P** which act as weights for the linear combination of the basis functions, giving the mapping to the physical space. In particular, given *<sup>n</sup>* one-dimensional basis functions *Ni*,*<sup>p</sup>* and n control points **<sup>P</sup>***<sup>i</sup>* <sup>∈</sup> IR2, *<sup>i</sup>* ∈ {1, ··· , *<sup>n</sup>*} a curve parameterization is given by:

$$\mathbf{C}\left(t\right) = \sum\_{i=1}^{n} \mathbf{P}\_{i} N\_{i}\left(t\right). \tag{11}$$

The control points define the so-called control mesh but this does not, in general, conform to the actual geometry (cf Figure 5b). In many real-world applications, the computational domain may be too complicated to be represented by a single NURBS mapping from the reference domain to the physical space. This could be due to topological reasons (or to the presence of different materials). In these cases it is common practice to resort to the so-called multipatch approach [4,12] (cf Figure 4). Here, the physical domain is split into simpler subdomains Ω*<sup>k</sup>* such as Ω = ∪*k*Ω*<sup>k</sup>* and Ω*<sup>i</sup>* ∩ Ω*<sup>j</sup>* = ∅, or a point, or an edge, for *i* = *j*.

**Figure 4.** Case 1: the multipatch domain.

#### *3.3. Isogeometric Conforming Spaces*

Analogously to classical Finite Element Methods, IGA is based on a Galerkin approach: the equations are written in their variational formulation, and the solution is sought in a finite-dimensional space with the correct approximation properties. In IGA the basis function space is inherited from the one used to parametrize the geometry.

Let us now consider a domain <sup>Ω</sup> <sup>⊂</sup> *<sup>V</sup>* that can be exactly parametrized with a mapping *F*

$$\mathbb{P}' \colon [0,1]^2 \to \Omega \tag{12}$$

For a multipatch approach we will consider a map for each Ω*k*. Then, the discrete space in the physical domain is defined applying the isoparametric concept as

$$Q\_h \quad = \quad \left\{ q\_h := \eta \circ \mathbb{I}^{t-1}, \quad \eta \in S^{p\_1, p\_2}\_{a\_1, a\_2} \left( \Xi \right) \right\} \tag{13}$$

$$\mathbf{V}\_h \quad = \quad \left\{ \vec{v}\_h := \vec{\varphi} \circ \vec{F}^{-1}, \quad \vec{\varphi} \in \left( S^{p\_1 + 1, p\_2 + 1}\_{\mathfrak{a}\_1, \mathfrak{a}\_2} \left( \mathbf{II} \right) \right)^2 \right\}, \tag{14}$$

where each **Π** has the same knots **Ξ** but the multiplicity has been increased by one, that means the velocity components' space have the same continuity as the pressure space [16]. We recall the important property of B-spline: at a knot of multiplicity *m*, basis function *Ni*,*<sup>p</sup>* is *Cp*−*<sup>m</sup>* = *C<sup>α</sup>* continuous.

**Proposition 1.** *With the notation and assumptions above, the space*

$$M\_{\rm h} = \left\{ \eta\_{\rm h} = \nabla \cdot \vec{v}\_{\rm h} \,\,\,\vec{v}\_{\rm h} \in \mathbf{V}\_{\rm h} \right\}.\tag{15}$$

*is subspace of Qh.*

**Proof of Proposition 1.** For *w* ∈ - *<sup>S</sup>p*1+1,*p*2+<sup>1</sup> *<sup>α</sup>*1,*α*<sup>2</sup> (**Π**) 2 we obtain

$$\nabla \cdot \vec{w} \in S^{p\_1, p\_2}\_{a\_1 - 1, a\_2 - 1} \left( \Pi \right) = S^{p\_1, p\_2}\_{a\_1, a\_2} \left( \Xi \right) \tag{7}$$

The subspace **V**<sup>0</sup> *<sup>h</sup>* ⊂ **V***<sup>h</sup>* h is the space of discrete functions that vanish on the boundary of Ω.

#### *3.4. The Variational Formulation Discretization*

We consider a classic mixed variational discretization of problem (4) in primitive variables [16,17], in which an approximation (*uh*, *ph*) to the exact solution (*u*, *p*) of (4) is obtained by solving the problem

$$\begin{aligned} a\left(\vec{u}\_{h\prime}\vec{v}\_{h}\right) + b\left(\vec{v}\_{h\prime}p\_{h}\right) &=& 0 \quad \forall \vec{v}\_{h} \in \mathsf{V}^{0}\_{h} \\ b\left(\vec{u}\_{h\prime}q\_{h}\right) &=& 0 \quad \forall q\_{h} \in \mathsf{M}\_{h\prime} \end{aligned} \tag{16}$$

where (**V***h*, *Mh*) are couples of finite-dimensional spaces parameterized with *h*, as introduced above, and with the bilinear forms *a* and *b* defined as

$$\begin{aligned} a \left( \vec{w}\_{\text{h}}, \vec{v}\_{\text{h}} \right) &= \quad \nu \int\_{\Omega} \nabla \vec{w}\_{\text{h}} : \nabla \vec{v}\_{\text{h}} \, d\Omega + \frac{\nu}{k} \int\_{\Omega} \vec{w}\_{\text{h}} \cdot \vec{v}\_{\text{h}} \, d\Omega, \\\\ b \left( \vec{v}\_{\text{h}}, q\_{\text{h}} \right) &= \quad - \int\_{\Omega} \nabla \cdot \vec{v}\_{\text{h}} \, q\_{\text{h}} \, d\Omega. \end{aligned}$$

The conditions for the well posedness of a saddle point system is known as inf-sup conditions or Ladyzhenskaya–Babuška–Breezi (LBB) condition:

$$\inf\_{q \in M\_h, q\_h \neq 0} \sup\_{\vec{v}\_h \in \mathbf{V}\_h, \vec{v}\_h \neq \vec{0}} \frac{b \left(\vec{v}\_h, q\_h\right)}{||q\_h||\_{L^2} ||\vec{v}\_h||\_{H^1}} \ge \mathcal{C}\_\prime \tag{17}$$

with *C* independent of the discritization parameter *h*, ·*L*<sup>2</sup> ·*H*<sup>1</sup> are the classic norms defined over the spaces *L*<sup>2</sup> (Ω) and *H*<sup>1</sup> (Ω), respectively, as in [11]. The subspace choice *Mh* ensures the discrete velocity *uh* solution of (16) is divergence-free.

In the context of IGA, we will use spline-based spaces **V**0 *<sup>h</sup>*, *Mh* , which satisfy the inf-sup condition (17). The elements considered can be seen as spline generalization of well-known finite elements spaces, namely Taylor–Hood elements [11]. Thanks to the high interelement regularity, which is the main feature of splines, the proposed discretizations are conforming, i.e., produce globally continuous, discrete velocities. Moreover, the spline generalization of Taylor–Hood elements also enjoys property [18] and thus provides divergence-free discrete solutions (15).

#### **4. Numerical Results**

In this section, we present some numerical experiments to predict the effect of fluid shear stress on cancer cells localized in the interstitial region to assess its metastatic behavior. For all cases, we use experimental data obtained in the works of [1,5–8].

We use GeoPDEs, an open source and free package for isogeometric analysis in Matlab [12,19], to achieved the numerical implementation of the discretized problem (16).

As a peculiarity of IGA is to allow for high degree and high regularity discretization spaces, most of the tests below are performed for degree *p* = 5 and regularity *α* = 4.

We use multipatch geometries (Figure 4) and these discretizations are *C*<sup>0</sup> between the patches. When applied to incompressible flows, these discretizations produce pointwise divergence-free velocity fields and hence exactly satisfy mass conservation. We enforce the Dirichlet boundary conditions weakly by Nitsche's method, allowing this method to default to a compatible discretization of Darcy flow in the limit of vanishing viscosity [12].

Using the data from Table 1 we perform different experiments with different choices of region configuration and tumor size.


**Table 1.** Boundary conditions and parameters values.

#### *4.1. Case 1*

Let us begin by considering a simple circular case which will be used as a pattern for the others ones. We consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular boundary, with radius *r* = 2(*μ*).

For this case we present the physical initial mesh, with control points, at Figure 5, an *h*-refinement (by multiple knot insertion) sequence, Figures 6 and 7, and the correspondent results for velocity and shear stress.

**Figure 5.** Case 1: physical domain, initial mesh and control points. (**a**) The physical mesh; (**b**) the control point for the NURBS surface.

**Figure 6.** Case 1: physical domain, *h*1-refinement mesh and correspondent control points. (**a**) The physical mesh; (**b**) the control point for the NURBS surface.

**Figure 7.** Case 1: physical domain, *h*2-refinement mesh and correspondent control points. (**a**) The physical mesh; (**b**) the control point for the NURBS surface.

With this discretization the velocity magnitude and stream lines obtained is represented in Figure 8. The flow-generate shear stress is between <sup>−</sup><sup>7</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> and 7 <sup>×</sup> <sup>10</sup>−<sup>2</sup> dyn/cm2 (we recall 1 dyn/cm2 = 0.1 Pa) over the whole domain, as shown in Figure 9, and maximum is achieved on the cell surface, this means on the membrane as expected.

**Figure 8.** Case 1: the velocity representation with stream lines.

**Figure 9.** Case 1: the fluid shear stress representation, which is between <sup>−</sup><sup>7</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> and 7 <sup>×</sup> 10−<sup>2</sup> dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

Next, for this case, we will repeat the simulation with a *h*-refinement on the domain. First we consider at patches interfaces the new knots array Ξ*h*<sup>1</sup> = [0.8, 0.95] .

With this discretization the flow-generate shear stress is between <sup>−</sup>1.5 and 1.5 dyn/cm2 over the whole domain, as shown in Figure 10, and velocity magnitude is between 0 and 2 μs<sup>−</sup>1, Figure 11.

**Figure 10.** Case 1: the fluid shear stress representation, which is between <sup>−</sup>1.5 and 1.5 dyn/cm2, for *h*1-refinement. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 11.** Case 1: the velocity representation with stream lines for domain *h*1-refinement.

Continuing with the *h*-refinement, we consider at paches interfaces the new knots array Ξ*h*<sup>2</sup> = [0.4, 0.6, 0.75, 0.85, 0.95].

With this discretization the flow-generate shear stress is between <sup>−</sup>1.5 and 1.5 dyn/cm2 over the whole domain, as shown in Figure 12, and velocity magnitude is between 0 and 2 μs<sup>−</sup>1, Figure 13, as in the previous refinement case.

**Figure 12.** Case 1: the fluid shear stress representation, which is between <sup>−</sup>1.5 and 1.5 dyn/cm2, for *h*2-refinement. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 13.** Case 1: the velocity representation with stream lines for domain *h*2-refinement.

These numerical results show that mesh refinement does not affect results for T*h*<sup>2</sup> . At this point, the model and its results are independent of the mesh.

Finally, we point to the velocity magnitude computed at any point *P* ∈ Ω using discrete <sup>2</sup> norm

$$\|\|\vec{u}(P)\|\| = \sqrt{u\_x(P)^2 + u\_y(P)^2}.\tag{18}$$

*4.2. Case 2*

Case 2 is similarly to case 1, we consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular irregular boundary, as shown at Figure 14. Whereas the irregular circular line is obtained from perturbations on the line of the pattern case 1, we quantify the boundary irregularities by comparison between the cells' perimeter, introducing the parameter I (the irregularity parameter indexed with the case study number) defined as

$$\mathcal{Z}\_2 = \frac{|d\_2 - d\_1|}{d\_1}.\tag{19}$$

with *d*<sup>1</sup> = 4*π* = 12.5664 and *d*<sup>2</sup> = 13.3937. For Case 2 we obtain I<sup>2</sup> = 0.0658.

**Figure 14.** Case 2: physical domain and mesh with *h*2-refinement.

With this cell configuration, I<sup>2</sup> = 0.0658, the flow-generate shear stress is between −1.7 and 1.7 dyn/cm2 over the whole domain, as shown in Figure 15, and velocity magnitude is between 0 and 2.2 μs<sup>−</sup>1, Figure 16.

**Figure 15.** Case 2: the fluid shear stress representation, which is between <sup>−</sup>1.7 and 1.7 dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 16.** Case 2: the velocity representation with stream lines.

#### *4.3. Case 3*

As in case 2, we consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular irregular boundary. Here we increase the irregularity parameter to <sup>I</sup><sup>3</sup> <sup>=</sup> <sup>|</sup>13.8401−12.5664<sup>|</sup> 12.5664 = 0.1014, as shown in Figure 17.

**Figure 17.** Case 3: physical domain and mesh with *h*2-refinement.

With this cell configuration, I<sup>3</sup> = 0.1014, the flow-generate shear stress is between −2.2 and 2.2 dyn/cm2 over the whole domain, as shown in Figure 18, and velocity magnitude is between 0 and 2.5 μs<sup>−</sup>1, Figure 19.

**Figure 18.** Case 3: the fluid shear stress representation, which is between <sup>−</sup>2.2 and 2.2 dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 19.** Case 3: the velocity representation with stream lines.

#### *4.4. Case 4*

Once again we consider Ω = [−10, 10] × [−10, 10] and a centred tumor region with a circular irregular boundary and the irregularity parameter <sup>I</sup><sup>4</sup> <sup>=</sup> <sup>|</sup>20.2692−12.5664<sup>|</sup> 12.5664 = 0.613, as shown at Figure 20.

**Figure 20.** Case 4: physical domain and mesh with *h*2-refinement.

With this cell configuration, I<sup>4</sup> = 0.613, the flow-generate shear stress is between −3 and 3 dyn/cm2 over the whole domain, as shown in Figure 21, and velocity magnitude is between 0 and 2.9 μs<sup>−</sup>1, Figure 22.

**Figure 21.** Case 4: the fluid shear stress representation, which is between <sup>−</sup>3 and 3 dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 22.** Case 4: the velocity representation with stream lines.

With the results from 1 to 4 we can define an evolution graphic in Figure 23 for shear stress maximum, the blue line, and velocity magnitude maximum, the red line. Regarding cell membrane resistance, at microenvironment context this evolution is quite important. We observed a linear dependence between the irregularity parameter and shear stress. The Figure 23 shows the dependence of shear stress on velocity and the irregularity parameter, given that the two lines have different growth rates.

**Figure 23.** Evolution for shear stress maximum and velocity magnitude maximum by the irregularity parameter.

To conclude we present two more cases quite near to reality. Inspired by Figure 1 we define the outline shown in Figure 24. The following cases are determined from the outline 3D by cut planes for case 5, the green one, and case 6, the blue one.

**Figure 24.** Tumor's 3D outline and cut planes for case 5, the green one, and case 6, the blue one.

#### *4.5. Case 5*

With this cell configuration, domain in Figure 25 (corresponding to the green section in Figure 24) we get <sup>I</sup><sup>5</sup> <sup>=</sup> 0.613, the flow-generate shear stress is between <sup>−</sup>2.8 and 2.5 dyn/cm2 over the whole domain, as shown in Figure 26, and velocity magnitude is between 0 and 2.8 μs<sup>−</sup>1, Figure 27. For this value of I<sup>5</sup> we obtain the expected from comparison with the above cases. We point to the geometric asymmetric effect on the maximum value of shear stress.

**Figure 25.** Case 5: physical domain and mesh.

**Figure 26.** Case 5: the fluid shear stress representation, which is between <sup>−</sup>2.8 and 2.5 dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone.

**Figure 27.** Case 5: the velocity representation with stream lines.

#### *4.6. Case 6*

Finally we consider for the cell configuration the blue section of the model in Figure 24, as shown in Figure 28. The flow-generate shear stress is between <sup>−</sup>4.1 and 3.3 dyn/cm2 over the whole domain, as shown in Figure 29, and velocity magnitude is between 0 and 1.9 μs<sup>−</sup>1, Figure 30.

With this result, we can see that the zone with the highest shear stress value is over the "body" of the cell and not in its branches. In fact, it is the branches we should focus on because they promote the displacement or duplication of the cell.

**Figure 29.** Case 6: the fluid shear stress representation, which is between <sup>−</sup>4.1 and 3.3 dyn/cm2. (**a**) The fluid shear stress overview; (**b**) zoom over the higher shear stress zone; (**c**) zoom over the higher shear stress zone; (**d**) zoom over the higher shear stress zone.

**Figure 30.** Case 6: the velocity representation with stream lines.

#### **5. Conclusions**

From this preliminary study, it is possible to establish the adequacy of the isogeometric analysis and NURBS to model irregular tumor cells and evaluate the shear forces at tumor microenvironments.

With some numerical examples, using different cell configurations, we have identified conditions that allow the increase or decrease of the fluid shear stress, which could contribute to the metastatic process by enhancing tumor cell invasion and circulating tumor cells. The cancer cell invasive potential is significantly reduced, as much as 92%, upon exposure to 0.55 dyn/cm2 fluid shear stress [21]. This work contributed to establish a relationship between a quantification of cell membrane irregularity and the maximum value of shear stress evaluated on the cell membrane, by the effect of interstitial fluid in the microenvironment. In follow-up works, we must understand how the cell reacts to these forces.

Mathematical models of fluid shear stress effects coupled with in vitro and in vivo experimental validation, may better predict cell behavior in such dynamic microenvironments, and potentially provide novel approaches for the prevention of metastasis.

A particular feature of cancer modeling revolves around the idea that a tumors' size and shape changes over time and its resistance to the fluid shear stress is also affected by the presence of other substances in the interstitial region meaning, for this reason, is the motivation for further work.

**Acknowledgments:** The author acknowledges the support of IPL IDI and CA through the Project IGACFC 2018. **Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Mathematical and Computational Applications* Editorial Office E-mail: mca@mdpi.com www.mdpi.com/journal/mca

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18