How to Make Systems of Nonlinear Autonomous ODEs with Attractor-Behavior, by First Making the General Solutions: Part Two

Magne Stensland

Moldjord, Norway.

DOI: 10.4236/jamp.2023.111009                

Abstract

This paper is presenting a new method for making first-order systems of nonlinear autonomous ODEs that exhibit limit cycles with a specific geometric shape in two and three dimensions, or systems of ODEs where surfaces in three dimensions have attractor behavior. The method is to make the general solutions first by using the exponential function, sine, and cosine. We are building up the general solutions bit for bit according to constant terms that contain the formula of the desired limit cycle, and differentiating them. In Part One, we used only formulas for closed curves where all parts of the formula were of the same degree. In order to use many other formulas for closed curves, the method in this paper is to introduce an additional variable, and we will get an additional ODE. We will choose the part of the formula with the highest degree and multiply the other parts with an extra variable, so that all parts of the formula have the same degree, creating a constant term containing this new formula. We will place it under the fraction line in the solutions, building up the rest of the solutions according to this constant term and differentiating. Keeping this extra variable constant, we will achieve almost the desired result. Using the methods described in this paper, it is possible to make some systems of nonlinear ODEs that are exhibiting limit cycles with a distinct geometric shape in two or three dimensions and some surfaces having attractor behavior, where not all parts of the formulas are the same degree. The pictures show the result.

Keywords

System of Nonlinear ODEs, Limit Cycle, General Solution, Attractor

1. Introduction

In the introduction to [1], we can read: Nonlinear dynamical systems exhibiting limit cycles are found in a large variety of fields including biology, chemistry, mechanics, and electronics.

A limit cycle is a closed trajectory in phase space having the property that at least one other trajectory spirals into it, either as time approaches to infinity or as time approaches to negative infinity. In other words, the limit cycle is an isolated trajectory (isolated in the sense that neighboring trajectories are not closed); it spirals either toward or away from the limit cycle. If all neighboring trajectories approach the limit cycle, we say the limit cycle is stable or attractive, that is, in the case where all the neighboring trajectories approach the limit cycle as time approaches to infinity [2] [3].

The limit cycle is an oscillation peculiar to nonlinear systems. The oscillatory behavior, unexplainable in terms of linear theory, is characterized by a constant amplitude and frequency determined by the nonlinear properties of the system. Limit cycles are distinguishable from linear oscillations in that their amplitude of oscillation is independent of initial conditions. For instance, if a system has a stable limit cycle, the system will tend to fall into the limit cycle, with the output approaching the amplitude of that limit cycle regardless of the initial conditions and forcing function [4].

A limit cycle is said asymptotic stable if all trajectories in the vicinity of the limit cycle converge to it as $t\to \infty$. Otherwise, the limit cycle is said semi-stable or unstable [5].

Roughly speaking, an attracting set for a dynamical system is a closed subset A of its phase space such that for “many” choices of initial point the system will evolve towards A [6]. An attractor of a dynamical system is a subset of the state space to which orbits originating from typical initial conditions tend as time increases. It is very common for dynamical systems to have more than one attractor [7]. An attractor is defined as the smallest unit which cannot be itself decomposed into two or more attractors with distinct basins of attraction. This restriction is necessary since a dynamic system may have multiple attractors, each with its own basin of attraction.

Until the 1960s, attractors were thought of as being simple geometric subsets of the phase space, like points, lines, and surfaces.

If somebody wants to make a system of ODEs that is exhibiting a Limit Cycle (LC) with a certain geometric shape, or a distinct surface that has attractor behavior, how can we do it?

1) For example, the graph of ${\left(x^2+y^2\right)}^4-P^2{\left(x^3-3x{y}^2\right)}^2=f^8$ as a limit cycle.

2) Or making a system of ODEs exhibiting a LC with the three-dimensional shape determined by the projection of the closed curve in a) onto a sphere.

3) Or making a system of three ODEs where the closed surface ${\left(x^2+y^2+z^2\right)}^4-P{\left(z^3-3x{y}^2\right)}^2=f^8$ has attractor-behavior, which is attracting the solution curves from both inside and outside this closed surface, and then following the surface as the variable t goes to infinity. P and f are parameters.

This paper is a continuation of Part One [8], and we will continue to describe some techniques I have developed for making nonlinear first-order systems with some special behavior. For the most part, we are going to work with systems of ODEs that are exhibiting limit cycles in two and three dimensions.

We repeat from Part One how to write the general solutions in order to make them easier to work with, and also the theory behind LC.

The method that we will use in this paper is new to the literature, at least to my knowledge. My background for the investigations is the textbook Differential Equations [9]. The number I have given each system of equations is the number they have in my collection.

2. How to Write the Solutions

First, we will look at the way I have chosen to write the solutions, in order to “simplify” the derivations.

The differential equation

$\frac{\text{d}x}{\text{d}t}=x-x^2$ (1)

has the general solution

$x\left(t\right)=\frac{1}{1-\frac{1}{x_0}\left(x_0-1\right){\text{e}}^{-t}}$ (2)

This is how we find the solution in the textbooks.

Multiply the numerator and the denominator with $x_0{\text{e}}^t$ :

$x\left(t\right)=\frac{x_0{\text{e}}^t}{x_0{\text{e}}^t-x_0+1}$ (3)

The drawback with this last expression is that we have got the same exponential function both in the numerator and the denominator. I prefer to use (3) because:

1) The exponent in the exponential function has the same sign than the coefficient in the differential equation. It is easier to see the connection between the solution and the equation.

2) It is easier to differentiate the solution when it is written on the form (3).

I have consequently chosen to write the solutions on the form (3). But when we want to analyze the solutions, we must use (2).

Derivation of fraction:

$\frac{\text{d}}{\text{d}x}\frac{f\left(x\right)}{g\left(x\right)}=\frac{f^{\prime }\left(x\right)g\left(x\right)-g^{\prime }\left(x\right)f(x)}{g{\left(x\right)}^2}=\frac{f^{\prime }\left(x\right)}{g\left(x\right)}-\frac{f\left(x\right)g^{\prime }\left(x\right)}{g{\left(x\right)}^2}$ (4)

We will use the last expression when derivation of solutions.

We can demonstrate by differentiating (3):

$\frac{\text{d}x}{\text{d}t}=\frac{x_0{\text{e}}^t}{x_0{\text{e}}^t-x_0+1}-\frac{x_0^2{\text{e}}^{2t}}{{\left(x_0{\text{e}}^t-x_0+1\right)}^2}=x-x^2$ (5)

Here, we see how simple the derivation becomes, and how easy it is to see the connection between the derivative and the solution $x\left(t\right)$, when we write the solution on the form (3).

The general expression of solution (3) is:

$x\left(t\right)=\frac{f{x}_0{\text{e}}^{at}}{x_0{\text{e}}^{at}-x_0+f}$ (6)

$\frac{\text{d}x}{\text{d}t}=ax-\frac{a}{f}{x}^2$ (7)

In order to practice the method that we are going to use later, it is crucial necessary that we write the solutions, so that $x\left(0\right)=x_0$, $y\left(0\right)=y_0$ and $z\left(0\right)=z_0$. If we put $t=0$ in the solution (6), we will get $x\left(0\right)=x_0$.

Another solution we will use is when the denominator is square root:

$x\left(t\right)=\frac{f{x}_0{\text{e}}^{at}}{\sqrt{x_0^2{\text{e}}^{2at}-x_0^2+f^2}}$ (8)

$x\left(0\right)=x_0$

$\begin{array}{c}\frac{\text{d}x}{\text{d}t}=\frac{af{x}_0{\text{e}}^{at}}{\sqrt{x_0^2{\text{e}}^{2at}-x_0^2+f^2}}-\frac{f{x}_0{\text{e}}^{at}\frac{1}{2}2a{x}_0^2{\text{e}}^{2at}}{\sqrt{x_0^2{\text{e}}^{2at}-x_0^2+f^2}\left(x_0^2{\text{e}}^{2at}-x_0^2+f^2\right)}\\ =ax-\frac{a}{f^2}{x}^3\end{array}$ (9)

And once again, we can see how easy it is to see the connection between the derivative $\frac{\text{d}x}{\text{d}t}$ and the solution $x\left(t\right)$ when we write the solution in this way. Notice that all parts in the square root are of same degree, namely the second degree. If the denominator is fourth root, all parts inside the root must be of fourth degree, and so on. This is important when we shall build up the solutions bit for bit, unless it will be wrong.

3. How to Make Limit Cycles (LCs)

3.1. The Theory behind LC

A spiral contains two types of motions: One rotating motion according to the function sine and cosine, and one outwards or inwards motion according to the exponential function. Remember that we are using only these three functions in the solutions described in this paper. In order to make a stable limit cycle, the equilibrium point must be spiral source. If we choose initial values $\left(x_0,y_0\right)$ inside the LC, the solution curve will make an outwards rotating spiral ending up at the LC. And if we choose initial values outside the LC, the solution curve will make an inwards rotating spiral ending up at the LC. When the solution curves arrive the LC, either from outside or inside the LC, the outwards or inwards motion will stop. We have only rotating motion left. What has happened?

Regarding the general solutions this means that the exponential functions have disappeared, and the solutions contain only the functions sine and cosine. If we choose initial values $\left(x_0,y_0\right)$ wherever on a stable LC, we will observe only rotation along the same closed curve over and over when the variable t goes to infinity.

This tells us that we must make the general solutions, so that if we choose initial values wherever on the LC, all exponential functions are gone. And if we choose initial values wherever else, the general solutions will contain exponential functions. The general solutions must contain a constant term made, so that if we choose initial values on the LC, this constant term become zero.

3.2. LC in Two Dimensions

The general solutions to a system of linear differential equations making spiral are:

$x\left(t\right)={\text{e}}^{at}\left(x_0\cos \left(bt\right)+y_0\sin \left(bt\right)\right)$ (10)

$y\left(t\right)={\text{e}}^{at}\left(y_0\cos \left(bt\right)-x_0\sin \left(bt\right)\right)$ (11)

In Part One, we made the general solutions to a system of ODEs exhibiting an elliptic LC, where we placed the formula of an ellipse in a constant term under the fraction line:

$-\left(\frac{x_0^2}{g^2}+\frac{y_0^2}{h^2}\right)+1$ (12)

Choosing initial values wherever on the ellipse $\frac{x^2}{g^2}+\frac{y^2}{h^2}=1$, this constant term will become zero, and the exponential functions will disappear, as we can see in the solutions below:

$x\left(t\right)=\frac{{\text{e}}^{at}\left(x_0\cos \left(bt\right)+y_0\sin \left(bt\right)\right)}{\sqrt{\frac{1}{g^2}{\text{e}}^{2at}{\left(x_0\cos \left(bt\right)+y_0\sin \left(bt\right)\right)}^2+\frac{1}{h^2}{\text{e}}^{2at}{\left(y_0\cos \left(bt\right)-x_0\sin \left(bt\right)\right)}^2-\left(\frac{x_0^2}{g^2}+\frac{y_0^2}{h^2}\right)+1}}$ (13)

$y\left(t\right)=\frac{{\text{e}}^{at}\left(y_0\cos \left(bt\right)-x_0\sin \left(bt\right)\right)}{\sqrt{\frac{1}{g^2}{\text{e}}^{2at}{\left(x_0\cos \left(bt\right)+y_0\sin \left(bt\right)\right)}^2+\frac{1}{h^2}{\text{e}}^{2at}{\left(y_0\cos \left(bt\right)-x_0\sin \left(bt\right)\right)}^2-\left(\frac{x_0^2}{g^2}+\frac{y_0^2}{h^2}\right)+1}}$ (14)

Notice that all parts inside the square root are of same degree, namely second degree. This is import, unless these methods will not work, and we will not obtain an ODE that is possible to express with the solutions $x\left(t\right)$ and $y\left(t\right)$.

If we want to make a system of ODEs exhibiting a LC expressed by the formula $x^4+y^4-4x{y}^2=1$, how can we do it? Two parts in the formula have degree four and one part has degree three. In order to give all parts in the formula the same degree, we can introduce an extra variable and name it w. $x^4+y^4-4x{y}^{2}w=1$. Now, all parts in the formula are of same degree, namely fourth degree.

Making this formula a little more general, the constant term under the fraction line must be:

$-\left(x_0^4+y_0^4-2P{x}_0{y}_0^2{w}_0\right)+f^4$ (15)

P and f are parameters.

Since we now have three variables, we must make three general solutions. After the derivations, we will have three ODEs. The fraction line will be too long for this page, so we put $x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$ and $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$ :

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[4]{{\text{e}}^{4at}{\left(A\right)}^4+{\text{e}}^{4at}{\left(B\right)}^4-2P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left(x_0^4+y_0^4-2P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (16)

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[4]{{\text{e}}^{4at}{\left(A\right)}^4+{\text{e}}^{4at}{\left(B\right)}^4-2P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left(x_0^4+y_0^4-2P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (17)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[4]{{\text{e}}^{4at}{\left(A\right)}^4+{\text{e}}^{4at}{\left(B\right)}^4-2P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left(x_0^4+y_0^4-2P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (18)

Notice that when $x_0^4+y_0^4-2P{w}_0{x}_0{y}_0^2=f^4$, the exponential function in all three solutions will disappear, and only the functions sine and cosine are left.

All three solutions contain the same exponential function ${\text{e}}^{at}$, and the denominator is the same in the three solutions. When differentiating these solutions, we obtain a system of three ODEs, system (1358):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^4}\left[a\left(x^4+y^4-2Pwx{y}^2\right)+bxy\left(x^2-y^2\right)-\frac{b}{2}P{y}^{3}w+bPw{x}^{2}y\right]$ (19)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^4}\left[a\left(x^4+y^4-2Pwx{y}^2\right)+bxy\left(x^2-y^2\right)-\frac{b}{2}P{y}^{3}w+bPw{x}^{2}y\right]$ (20)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^4}\left[a\left(x^4+y^4-2Pwx{y}^2\right)+bxy\left(x^2-y^2\right)-\frac{b}{2}P{y}^{3}w+bPw{x}^{2}y\right]$ (21)

Notice that when $x^4+y^4-2Pwx{y}^2=f^4$, all parts containing the parameter a are gone. The parameter a belongs to the only exponential function in the three ODEs.

The secret now is to keep the variable w constant, for example $w=1$. Then we have obtained the system that we wanted to make:

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^4}\left[a\left(x^4+y^4-2Px{y}^2\right)+bxy\left(x^2-y^2\right)-\frac{b}{2}P{y}^3+bP{x}^{2}y\right]$ (22)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^4}\left[a\left(x^4+y^4-2Px{y}^2\right)+bxy\left(x^2-y^2\right)-\frac{b}{2}P{y}^3+bP{x}^{2}y\right]$ (23)

I don’t have solutions to this system, but that don’t matter. We have obtained what we wanted to do, making a system of two ODEs with the desired behavior.

See Figure 1 where $a=4,b=1,f=1,P=2$ and compare with Figure 2.

My experience is that choosing greater values of the parameter a will give a better similarity with the desired geometric shape of the LC.

Figure 1. A two-dimensional LC.

Figure 2. The graph of $x^4+y^4-4x{y}^2=1$.

Inspired of these results, let us try to give the closed curve ${\left(x^2+y^2\right)}^2-2x^3+6x{y}^2=1$ attractor behavior. Multiply two parts in the formula with an extra variable, and make it a little more general:

${\left(x^2+y^2\right)}^2-Pw{x}^3+3Pxw{y}^2=f^4$ (24)

P and f are parameters.

Now, all parts in the formula have degree four, and then we have the constant term:

$-\left({\left(x_0^2+y_0^2\right)}^2-P{w}_0{x}_0^3+3P{x}_0{w}_0{y}_0^2\right)+f^4$ (25)

Build up the general solutions $x\left(t\right),y\left(t\right)$ and $w\left(t\right)$ according to this constant term, and define

$x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$, $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$ in order to compress the solutions.

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[4]{{\text{e}}^{4at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^2-P{w}_0{\text{e}}^{4at}{\left(A\right)}^3+3P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left({\left(x_0^2+y_0^2\right)}^2-P{w}_0{x}_0^3+3P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (26)

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[4]{{\text{e}}^{4at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^2-P{w}_0{\text{e}}^{4at}{\left(A\right)}^3+3P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left({\left(x_0^2+y_0^2\right)}^2-P{w}_0{x}_0^3+3P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (27)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[4]{{\text{e}}^{4at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^2-P{w}_0{\text{e}}^{4at}{\left(A\right)}^3+3P{\text{e}}^{4at}{w}_0\left(A\right){\left(B\right)}^2-\left({\left(x_0^2+y_0^2\right)}^2-P{w}_0{x}_0^3+3P{w}_0{x}_0{y}_0^2\right)+f^4}}$ (28)

Differentiating these solutions give system (1356):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^4}\left[a{\left(x^2+y^2\right)}^2-aP{x}^{3}w+3aPx{y}^{2}w-\frac{9}{4}Pb{x}^{2}yw+\frac{3}{4}Pb{y}^{3}w\right]$ (29)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^4}\left[a{\left(x^2+y^2\right)}^2-aP{x}^{3}w+3aPx{y}^{2}w-\frac{9}{4}Pb{x}^{2}yw+\frac{3}{4}Pb{y}^{3}w\right]$ (30)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^4}\left[a{\left(x^2+y^2\right)}^2-aP{x}^{3}w+3aPx{y}^{2}w-\frac{9}{4}Pb{x}^{2}w+\frac{3}{4}Pb{y}^{3}w\right]$ (31)

We keep the extra variable w constant, choose $w=1$, and obtain a system of two ODEs with the desired behavior:

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^4}\left[a{\left(x^2+y^2\right)}^2-aP{x}^3+3aPx{y}^2-\frac{9}{4}Pb{x}^{2}y+\frac{3}{4}Pb{y}^3\right]$ (32)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^4}\left[a{\left(x^2+y^2\right)}^2-aP{x}^3+3aPx{y}^2-\frac{9}{4}Pb{x}^{2}y+\frac{3}{4}Pb{y}^3\right]$ (33)

See Figure 3 and compare with Figure 4.

In Figure 3, $a=6,b=1,w=1,f=1,P=2$.

Level curve in height $c=\left\{\left(x,y\right)\in R^2|f\left(x,y\right)=c\right\}$ [10].

In this case $w=f\left(x,y\right)=\frac{f^4-{\left(x^2+y^2\right)}^2}{2Px{y}^2-P{x}^3}=c$, a constant. (34)

The LC in Figure 3 is a level curve of $f\left(x,y\right)$ in height $c=1$.

When ${\left(x^2+y^2\right)}^2-P{x}^3+2Px{y}^2=f^4$ in Equations (32) and (33), all parts of the system that contain the parameter a are gone. Along this closed curve all exponential functions in the solutions have disappeared, and we have a limit cycle as we can see in the phase portrait.

Let us now try to make a LC that reminds a bit of the quadrifolium, the closed curve:

${\left(x^2+y^2\right)}^3-4P^2{x}^2{y}^2=f^6$ (35)

Figure 3. A LC that reminds a bit of the trifolium.

Figure 4. The graph of ${\left(x^2+y^2\right)}^2-2x^3+6x{y}^2=1$.

Use an extra variable w to give all parts in the formula the same degree:

${\left(x^2+y^2\right)}^3-4P^2{x}^2{y}^2{w}^2=f^6$ (36)

The constant term under the fraction line:

$-\left({\left(x_0^2+y_0^2\right)}^3-4P^2{x}_0^2{y}_0^2{w}_0^2\right)+f^6$ (37)

And then build up the solutions $x\left(t\right),y\left(t\right)$ and $w\left(t\right)$ according to this constant term. We put $x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$, $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$ in order to compress the solutions.

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^3-4P^2{\text{e}}^{6at}{\left(A\right)}^2{\left(B\right)}^2{w}_0^2-\left({\left(x_0^2+y_0^2\right)}^3-4P^2{x}_0^2{y}_0^2{w}_0^2\right)+f^6}}$ (38)

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^3-4P^2{e}^{6at}{\left(A\right)}^2{\left(B\right)}^2{w}_0^2-\left({\left(x_0^2+y_0^2\right)}^3-4P^2{x}_0^2{y}_0^2{w}_0^2\right)+f^6}}$ (39)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^3-4P^2{\text{e}}^{6at}{\left(A\right)}^2{\left(B\right)}^2{w}_0^2-\left({\left(x_0^2+y_0^2\right)}^3-4P^2{x}_0^2{y}_0^2{w}_0^2\right)+f^6}}$ (40)

Notice the sixth root of the denominator and that all parts inside the root are of sixth degree, and that when ${\left(x_0^2+y_0^2\right)}^3-4P^2{x}_0^2{y}_0^2{w}_0^2=f^6$, all exponential functions are gone.

Differentiating these three solutions give system (1354):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^6}\left[a{\left(x^2+y^2\right)}^3-4a{P}^2{x}^2{y}^2{w}^2-\frac{4}{3}b{P}^{2}x{y}^3{w}^2+\frac{4}{3}b{P}^2{x}^{3}y{w}^2\right]$ (41)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^6}\left[a{\left(x^2+y^2\right)}^3-4a{P}^2{x}^2{y}^2{w}^2-\frac{4}{3}b{P}^{2}x{y}^3{w}^2+\frac{4}{3}b{P}^2{x}^{3}y{w}^2\right]$ (42)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^6}\left[a{\left(x^2+y^2\right)}^3-4a{P}^2{x}^2{y}^2{w}^2-\frac{4}{3}b{P}^{2}x{y}^3{w}^2+\frac{4}{3}b{P}^2{x}^{3}y{w}^2\right]$ (43)

When ${\left(x^2+y^2\right)}^3-4P^2{x}^2{y}^2{w}^2=f^6$ all parts containing the parameter a are gone. The parameter a belongs to the exponential functions in the solutions.

Keep the variable w constant, and we have obtained a system of 2 ODEs with the desired behavior. See Figure 5 and compare with Figure 6. In Figure 5, $a=8,b=1,f=2,P=4,w^2=2$

Level curve at height $c=\left\{\left(x,y\right)\in R^2|f\left(x,y\right)=c\right\}$.

In this case

$f\left(x,y\right)=w=\pm \frac{\sqrt{{\left(x^2+y^2\right)}^3-f^6}}{2Pxy}=c$ (44)

For $w^2=2$, as in Figure 5, we have two equal limit cycles as level curves at the heights $c=\pm \sqrt{2}$.

Now, we are able to solve problem a) in the Introduction: Make a system of ODEs exhibiting a limit cycle with the geometric shape ${\left(x^2+y^2\right)}^4-P^2{\left(x^3-3x{y}^2\right)}^2=f^8$. The constant term under the fraction line must be:

$-\left({\left(x_0^2+y_0^2\right)}^4-P^2{w}_0^2{\left(x_0^3-3x_0{y}_0^2\right)}^2\right)=f^8$ (45)

We build up the solutions according to this constant term, and put $x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$, $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$ and ${\left(x_0^2+y_0^2\right)}^4-P^2{w}_0^2{\left(x_0^3-3x_0{y}_0^2\right)}^2=C$

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^4-P^2{\text{e}}^{8at}{w}_0^2{\left({\left(A\right)}^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (46)

Figure 5. A LC that reminds a bit of the quadrifolium.

Figure 6. The graph of ${\left(x^2+y^2\right)}^3-128x^2{y}^2=64$.

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^4-P^2{\text{e}}^{8at}{w}_0^2{\left({\left(A\right)}^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (47)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2\right)}^4-P^2{\text{e}}^{8at}{w}_0^2{\left({\left(A\right)}^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (48)

Notice that we have taken the eighth root of the denominator and all parts inside the root have degree eight.

Differentiating these solutions give system (1357):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^8}\left[a{\left(x^2+y^2\right)}^4-a{P}^2{w}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2{w}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right]$ (49)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^8}\left[a{\left(x^2+y^2\right)}^4-a{P}^2{w}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2{w}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right]$ (50)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^8}\left[a{\left(x^2+y^2\right)}^4-a{P}^2{w}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2{w}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right]$ (51)

We keep the variable w constant, for example $w^2=1$, and we have a system with the desired behavior:

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^8}\left[a{\left(x^2+y^2\right)}^4-a{P}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right]$ (52)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^8}\left[a{\left(x^2+y^2\right)}^4-a{P}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right]$ (53)

I don’t have solutions to the system (52) and (53), but that is no problem since this system gives us a LC with the desired geometric shape. See Figure 7 and compare with Figure 8. In Figure 7, $a=6,b=1,P=2,f=1,w=1$.

My experience is that choosing higher values of the parameter a gives a better similarity between the shape of the LC and the graph of the closed curve. In Figure 7 the parameter $a=6$, and is giving a good similarity, but the curves are not absolutely similar.

Level curve in height $c=\left\{(x,y)\in R^2|f\left(x,y\right)=c\right\}$

In this case $f\left(x,y\right)=w=\pm \frac{\sqrt{{\left(x^2+y^2\right)}^4-f^8}}{P\left(x^3-3x{y}^2\right)}=c$

For $w^2=1$, as in Figure 7, we have two equal limit cycles as level curves at height $c=\pm 1$.

Figure 7. A LC that reminds a bit of the six-folium.

Figure 8. The graph of ${\left(x^2+y^2\right)}^4-4{\left(x^3-3x{y}^2\right)}^2=1$.

3.3. LC in Three Dimensions

Let us try to give the LC in Figure 7 a three-dimensional shape by projecting it on surfaces that we know the formula to, so that the LC will lie on these surfaces. The solution curves will then follow the surface, and repeat the same closed curve over and over when the variable t goes to infinity.

We must make the solution $z\left(t\right)$, so that the exponential function belonging to $z\left(t\right)$ has disappeared when the solution curves have arrived this surface. We are going to make a constant term in a such way that when we choose initial values on this surface, this constant term will become zero, and for all other initial values is this constant term not zero.

For example, the paraboloid or hyperbolic paraboloid $z=c{x}^2+d{y}^2$. The solution $z\left(t\right)$ contain two constant terms: One for the closed curve (45), as we earlier have placed under the fraction line in the solutions $x\left(t\right),y\left(t\right)$ and $w\left(t\right)$. And one constant term containing the formula for the paraboloid or hyperbolic paraboloid $+\left(z_0-c{x}_0^2-d{y}_0^2\right)$, and multiply this last constant term with an exponential function. When both constant terms are zero, are absolutely all exponential functions gone, and the general solutions contain only the functions sine and cosine. We have then obtained a three-dimensional LC.

$z\left(t\right)=cx{\left(t\right)}^2+dy{\left(t\right)}^2+\left(z_0-c{x}_0^2-d{y}_0^2\right){\text{e}}^{kt}$ (54)

$\frac{\text{d}z}{\text{d}t}=c2x\frac{\text{d}x}{\text{d}t}+d2y\frac{\text{d}y}{\text{d}t}+k\left(z_0-c{x}_0^2-d{y}_0^2\right){\text{e}}^{kt}$ (55)

Subtract kz from both sides and put in the expressions for $\frac{\text{d}y}{\text{d}t}$ (49) and for $\frac{\text{d}y}{\text{d}t}$ (50).

$\begin{array}{c}\frac{\text{d}z}{\text{d}t}=kz+\left(2a-k\right)\left(c{x}^2+d{y}^2\right)+2bxy\left(c-d\right)-\frac{1}{f^8}\left(c{x}^2+d{y}^2\right)[2a{\left(x^2+y^2\right)}^4\\ -2a{P}^2{w}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{2}b{P}^2{w}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)]\end{array}$ (56)

Choose the extra variable $w^2=1$. Together with (52) and (53) we have the system (1425)

See Figure 9, where $a=6,b=1,c=1,d=1,P=2,k=-1,w=1$.

Now, we are able to solve problem b) in the Introduction: Make a system of ODEs exhibiting a LC with the three-dimensional shape determined by the projection of the closed curve in Figure 8 onto a sphere. We have already made the solutions $x\left(t\right)$ (46), $y\left(t\right)$ (47) and $w\left(t\right)$ (48). We will obtain two LC, one on the upper hemisphere and one on the lower hemisphere.

The solution $z\left(t\right)$ must contain the formula for a sphere as a constant term.

$x^2+y^2+z^2=R^2$

$z=\pm \sqrt{R^2-x^2-y^2}$

$z\left(t\right)=\pm \sqrt{R^2-x{\left(t\right)}^2-y{\left(t\right)}^2}+\left(z_0\mp \sqrt{R^2-x_0^2-y_0^2}\right){\text{e}}^{kt}$ (57)

$\frac{\text{d}z}{\text{d}t}=\pm \frac{-2x\frac{\text{d}x}{\text{d}t}-2y\frac{\text{d}y}{\text{d}t}}{2\sqrt{R^2-x^2-y^2}}+k\left(z_0\mp \sqrt{R^2-x_0^2-y_0^2}\right){\text{e}}^{kt}$ (58)

Subtract kz from both sides.

$\begin{array}{l}\frac{\text{d}z}{\text{d}t}=kz\mp k\sqrt{R^2-x^2-y^2}\pm \frac{1}{\sqrt{R^2-x^2-y^2}}[-a{x}^2-a{y}^2+\frac{x^2+y^2}{f^8}\\ \times \left(a{\left(x^2+y^2\right)}^4-a{P}^2{w}^2{\left(x^3-3x{y}^2\right)}^2-\frac{3}{4}b{P}^2{w}^2\left(x^3-3x{y}^2\right)\left(3x^{2}y-y^3\right)\right)]\end{array}$ (59)

Figure 9. The LC in Figure 7 on the paraboloid $z=x^2+y^2$.

Keep the extra variable constant, $w=1$. Together with Equations (52) and (53), we have the system (1429). See Figure 10 with a three-dimensional LC on the upper hemisphere, and Figure 11 that shows a three-dimensional LC on the lower hemisphere.

In Figure 10 and Figure 11, $a=6,b=1,P=2,f=1,R=3,w=1,k=-1$.

The initial values are the same in both figures.

Many more systems of ODEs exhibiting a three-dimensional LC are possible to make by projecting whatever closed curve that you know the formula to, onto whatever surface that you know the formula to.

Figure 10. A LC on the upper hemisphere.

Figure 11. A LC on the lower hemisphere.

4. Surfaces as Attractors

In Part One [2], we made systems of ODEs where a distinct closed surface had attractor behavior, that is attracting the solution curves from both inside and outside the closed surface, and following the surface over and over in the same orbit as the parameter t goes to infinity [5]. The orbit is determined by the initial values. It is a necessary condition that the equilibrium point is spiral source. All exponential functions have disappeared along the whole surface. This is what characterizes a stable three-dimensional attractor.

In this paper, we will use formulas to closed surfaces where not all parts of the formulas are of same degree. As in the previous sections, we are going to use an extra variable w. We will get an extra ODE, and achieve a system of four ODEs with four variables. Keeping this extra variable constant, the problem is solved. The denominator is the same in all four solutions, and they contain the same exponential function.

Consider the formula to a closed surface:

${\left(x^2+y^2+z^2\right)}^3+P^2{x}^2{y}^3=f^6$ (60)

P and f are parameters. The constant term under the fraction line must be:

$-\left({\left(x_0^2+y_0^2+z_0^2\right)}^3+P^2{x}_0^2{y}_0^3{w}_0\right)+f^6$ (61)

Build up the solutions $x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right)$ according to this constant term. In order to compress the solutions, put $x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$, $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^3+P^2{\text{e}}^{6at}{w}_0{\left(A\right)}^2{\left(B\right)}^3-\left({\left(x_0^2+y_0^2+z_0^2\right)}^3+P^2{x}_0^2{y}_0^3{w}_0\right)+f^6}}$ (62)

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^3+P^2{\text{e}}^{6at}{w}_0{\left(A\right)}^2{\left(B\right)}^3-\left({\left(x_0^2+y_0^2+z_0^2\right)}^3+P^2{x}_0^2{y}_0^3{w}_0\right)+f^6}}$ (63)

$z\left(t\right)=\frac{f{z}_0{\text{e}}^{at}}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^3+P^2{\text{e}}^{6at}{w}_0{\left(A\right)}^2{\left(B\right)}^3-\left({\left(x_0^2+y_0^2+z_0^2\right)}^3+P^2{x}_0^2{y}_0^3{w}_0\right)+f^6}}$ (64)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[6]{{\text{e}}^{6at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^3+P^2{\text{e}}^{6at}{w}_0{\left(A\right)}^2{\left(B\right)}^3-\left({\left(x_0^2+y_0^2+z_0^2\right)}^3+P^2{x}_0^2{y}_0^3{w}_0\right)+f^6}}$ (65)

Differentiating these four solutions give system (1686):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^6}\left[a{\left(x^2+y^2+z^2\right)}^3+a{P}^2{x}^2{y}^{3}w+\frac{b}{3}{P}^{2}x{y}^{4}w-\frac{b}{2}{P}^2{x}^3{y}^{2}w\right]$ (66)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^6}\left[a{\left(x^2+y^2+z^2\right)}^3+a{P}^2{x}^2{y}^{3}w+\frac{b}{3}{P}^{2}x{y}^{4}w-\frac{b}{2}{P}^2{x}^3{y}^{2}w\right]$ (67)

$\frac{\text{d}z}{\text{d}t}=az-\frac{z}{f^6}\left[a{\left(x^2+y^2+z^2\right)}^3+a{P}^2{x}^2{y}^{3}w+\frac{b}{3}{P}^{2}x{y}^{4}w-\frac{b}{2}{P}^2{x}^3{y}^{2}w\right]$ (68)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^6}\left[a{\left(x^2+y^2+z^2\right)}^3+a{P}^2{x}^2{y}^{3}w+\frac{b}{3}{P}^{2}x{y}^{4}w-\frac{b}{2}{P}^2{x}^3{y}^{2}w\right]$ (69)

Keeping the variable w constant, we have made a system of three ODEs with the desired behavior.

Notice that when ${\left(x^2+y^2+z^2\right)}^3+P^2{x}^2{y}^{3}w=f^6$ all parts containing the parameter a have disappeared, that means all exponential functions are gone. When $a>0$, the solution curves will follow the closed surface as the variable t goes to infinity.

See Figure 12, where $a=8,b=1,f=3,P=2,w=6$.

We see 9 limit cycles on this closed surface, one for each initial value.

Level set at height $c=\left\{(x,y,z)\in R^3|f\left(x,y,z\right)=c\right\}$.

In this case $w=f\left(x,y,z\right)=\frac{f^6-{\left(x^2+y^2+z^2\right)}^3}{P^2{x}^2{y}^3}=c$, a constant. (70)

Now, we are able to solve problem c) in the Introduction: Make a system of ODEs where the closed surface ${\left(x^2+y^2+z^2\right)}^4-P{\left(z^3-3x{y}^2\right)}^2=f^8$ has attractor-behavior.

The constant term under the fraction line to all the solutions $x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right)$ must be:

$-\left({\left(x_0^2+y_0^2+z_0^2\right)}^4-P{w}_0^2{\left(z_0^3-3x_0{y}_0^2\right)}^2\right)+f^8$ (71)

In order to compress the solutions, put $x_0\cos \left(bt\right)+y_0\sin \left(bt\right)=A$, $y_0\cos \left(bt\right)-x_0\sin \left(bt\right)=B$, ${\left(x_0^2+y_0^2+z_0^2\right)}^4-P{w}_0^2{\left(z_0^3-3x_0{y}_0^2\right)}^2=C$ The fraction line will be too long for this page.

$x\left(t\right)=\frac{f{\text{e}}^{at}\left(A\right)}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^4-P{\text{e}}^{8at}{w}_0^2{\left(z_0^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (72)

Figure 12. 9 LC on a closed surface.

$y\left(t\right)=\frac{f{\text{e}}^{at}\left(B\right)}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^4-P{\text{e}}^{8at}{w}_0^2{\left(z_0^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (73)

$z\left(t\right)=\frac{f{z}_0{\text{e}}^{at}}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^4-P{\text{e}}^{8at}{w}_0^2{\left(z_0^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (74)

$w\left(t\right)=\frac{f{w}_0{\text{e}}^{at}}{\sqrt[8]{{\text{e}}^{8at}{\left({\left(A\right)}^2+{\left(B\right)}^2+z_0^2\right)}^4-P{\text{e}}^{8at}{w}_0^2{\left(z_0^3-3\left(A\right){\left(B\right)}^2\right)}^2-\left(C\right)+f^8}}$ (75)

Notice that when $\left(C\right)=f^8$, the exponential functions in all four solutions are gone.

Differentiating these four solutions give system (1762):

$\frac{\text{d}x}{\text{d}t}=ax+by-\frac{x}{f^8}\left[a{\left(x^2+y^2+z^2\right)}^4-aP{w}^2{\left(z^3-3x{y}^2\right)}^2-\frac{3}{4}Pb{w}^2\left(z^3-3x{y}^2\right)\left(2x^{2}y-y^3\right)\right]$ (76)

$\frac{\text{d}y}{\text{d}t}=ay-bx-\frac{y}{f^8}\left[a{\left(x^2+y^2+z^2\right)}^4-aP{w}^2{\left(z^3-3x{y}^2\right)}^2-\frac{3}{4}Pb{w}^2\left(z^3-3x{y}^2\right)\left(2x^{2}y-y^3\right)\right]$ (77)

$\frac{\text{d}z}{\text{d}t}=az-\frac{z}{f^8}\left[a{\left(x^2+y^2+z^2\right)}^4-aP{w}^2{\left(z^3-3x{y}^2\right)}^2-\frac{3}{4}Pb{w}^2\left(z^3-3x{y}^2\right)\left(2x^{2}y-y^3\right)\right]$ (78)

$\frac{\text{d}w}{\text{d}t}=aw-\frac{w}{f^8}\left[a{\left(x^2+y^2+z^2\right)}^4-aP{w}^2{\left(z^3-3x{y}^2\right)}^2-\frac{3}{4}Pb{w}^2\left(z^3-3x{y}^2\right)\left(2x^{2}y-y^3\right)\right]$ (79)

Keeping the extra variable w constant, the problem is solved. We have general solutions to a system of four ODEs (1762), but we don’t have solutions to the system of three ODEs in level height $w=c$, a constant.

See Figure 13, where $a=8,b=1,P=6,f=3,w^2=4$.

Figure 13. 7 LC on a closed surface.

Notice that when ${\left(x^2+y^2+z^2\right)}^4-P{w}^2{\left(z^3-3x{y}^2\right)}^2=f^8$, all parts in the ODEs containing the parameter a have disappeared. That means all exponential functions in the solutions are gone along this surface. When $a>0$, this surface is an attractor.

Level set at height $c=\left\{\left(x,y,z\right)\in R^3|f\left(x,y,z\right)=c\right\}$.

In this case $w=f\left(x,y,z\right)=\pm \sqrt{\frac{{\left(x^2+y^2+z^2\right)}^4-f^8}{P{\left(z^3-3x{y}^2\right)}^2}}=\pm c$, a constant.

We have two equal solutions in level height $w=c=\pm 2$.

5. Conclusions

It is possible to make systems of nonlinear ODEs that are exhibiting limit cycles by making a system of general solutions first. A short summary of the techniques that we have used in this paper:

1) Choose a formula for a closed curve.

2) If not all parts of the formula have the same degree, choose the part with the highest degree and multiply the other parts with an extra variable, so that all parts of the formula have the same degree.

3) Create a constant term containing this formula, so that if you choose initial values on the closed curve, this constant term will become zero.

4) Place this constant term under the fraction line and build up the solutions according to the constant term. The denominator must be the same in all solutions.

5) Differentiate all solutions and keep the extra variable constant, and you have obtained the desired system of ODEs.

6) If you want to give this two-dimensional limit cycle a three-dimensional shape, make a constant term containing the formula for a surface, no matter how complicated this formula is, and multiply it with an exponential function. Make the constant term, so that it will become zero when you choose initial values wherever on the surface.

7) Differentiate the solution and keep the extra variable constant with the same value as in the other solutions, and you have obtained a system of three ODEs with the desired behavior.

8) If you want to make a system of ODEs where a closed surface has attractor-behavior and not all parts of the formula of the closed surface are of the same degree, choose the part with the highest degree and multiply the other parts with an extra variable, so that all parts of the formula have the same degree.

9) Create a constant term containing this formula, so that if you choose initial values wherever on the surface, this constant term will become zero.

10) Place this constant term under the fraction line in all solutions and build up the solutions according to the constant term. The denominator must be the same in all solutions and they must contain the same exponential function.

11) Differentiate the four solutions. Keep the extra variable constant and you have obtained a system of three ODEs with the desired behavior.

You can place as many limit cycles on a closed surface, and some open surfaces, as you want to. Choosing greater values of the parameter a and using a better computer, I am sure you will achieve a better similarity with the desired geometric shape.

Many more limit cycles in two and three dimensions are easy to make using these methods. Only the imagination is the limit. The possibilities for a lot of mathematical fun seem unlimited.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Maliki, O.S. and Sesan, O. (2018) On Existence of Periodic Solutions of Certain Second Order Nonlinear ODEs via Phase Portrait Analysis. Applied Mathematics, 9, 1225-1237.
[2]	Espinosa-Paredes, G. (2021) Chapter 7—Nonlinear BWR Dynamics with a Fractional Reduced Order Model. Fractional-Order Models for Nuclear Reactor Analysis, 247-295. https://doi.org/10.1016/B978-0-12-823665-9.00007-9
[3]	Prieto-Guerrero, A. and Espinosa-Paredes, G. (2019) 2—Description of Boiling Water Reactors. Linear and Non-Linear Stability Analysis in Boiling Water Reactors, 25-55. https://doi.org/10.1016/B978-0-08-102445-4.00002-3
[4]	Van Valkenberg, M.E. and Middleton, W.M. (Ed.) (2002) Reference Data for Engineers: Radio, Electronics, Computers and Communications. Elsevier, Amsterdam.
[5]	Salehi Fathabadi, H. (2012) Behavior of Limit Cycles in Nonlinear Systems. International Journal of Information and Electronics Engineering, 24, 523-526. https://doi.org/10.7763/IJIEE.2012.V2.152
[6]	Milnor, J.W. (2006) Attractor. Scholarpedia, 1, Article No.1815.
[7]	Ott, E. (2006) Basin of Attraction. Scholarpedia, 1, Article No. 1701. https://doi.org/10.4249/scholarpedia.1701
[8]	Stensland, M. (2022) How to Make Systems of Nonlinear Autonomous ODEs with Attractor-Behavior, by First Making the General Solutions: Part One. Journal of Applied Mathematics and Physics, 10, 3814-3835. https://doi.org/10.4236/jamp.2022.1012253
[9]	Blanchard, P., Devaney, R. L. and Hall, G.R. (2002) Differential Equations. 2nd Edition, BROOKS/COLE, Pacific Grove, 97, 197, 438.
[10]	Colley, S.J. (2002) Vector Calculus. Prentice Hall, Upper Saddle River, NJ, 90-97.