# Preface

This section studies some first order nonlinear ordinary differential equations describing the time evolution (or “motion”) of those hamiltonian systems provided with a first integral linking implicitly both variables to a motion constant. An application has been performed on the Lotka--Volterra predator-prey system, turning to a strongly nonlinear differential equation in the phase variables.

Introduction to Linear Algebra with Mathematica

The EquationTrekker package is a great package for plotting and exploring phase space

# Limit Cycles

Consider a system of first order autonomous ordinary differential equations in normal form
$\frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x} ) ,$
where
${\bf x} (t) = \begin{bmatrix} x_1 (t) \\ x_2 (t) \\ \vdots \\ x_n (t) \end{bmatrix} , \qquad {\bf f} (x_1 , x_2 , \ldots , x_n , t) = \begin{bmatrix} f_1 (x_1 , x_2 , \ldots , x_n ) \\ f_2 (x_1 , x_2 , \ldots , x_2 ) \\ \vdots \\ f_n (x_1 , x_2 , \ldots , x_n ) \end{bmatrix}$
are n-column vectors. In engineering and physics, it is a custom to denote a derivative with respect to time variable t by dot: $$\dot{\bf x} = {\text d}{\bf x} / {\text d} t.$$

If an initial position of the vector $${\bf x} (t)$$ is known, we get an initial value problem:

$\dot{\bf x} \equiv \frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x} ) , \qquad {\bf x} (t_0 ) = {\bf x}_0 ,$
where $${\bf x}_0$$ is a given column vector.

# 2.3.5. Limited Circles

Liénard's Theorem: A Liénard system $$\ddot{x} +f(x)\,\dot{x} + g(x) =0$$ has a unique and stable limit cycle surrounding the origin if it satisfies the following additional properties:
• g(x) > 0 for all x > 0;
• $$\lim_{x\to \infty} F(x) \equiv \lim_{x\to \infty} \int_0^x f(\xi )\,{\text d}\xi = \infty ;$$
• F(x) has exactly one positive root at some value p, where F(x) < 0 for 0 < x < p and F(x) > 0 and monotonic for x > p.

Example 1: Consider the undamped Duffing equation
$\ddot{x} +x +0.1\,x^3 =0,$
which we transfer into the first-order system of equations
$v= \dot{\bf x} \qquad \dot{v} = -x-0.1\,x^3 .$
dirfield = VectorPlot[{v, -x - 0.1*x^3}, {x, -4, 4}, {v, -5, 5}, Prolog -> {Thickness[0.001]},
Axes -> True, AspectRatio -> 1, VectorPoints -> 15, AxesLabel -> {x, v}] There is a strong suggestion of periodic motion, in the way that the arrows circle around. Now we will construct some solutions and then display them with the direction field. We define the equation first.
duff[x0_, v0_] := {x'[t] == v[t], v'[t] == -x[t] - 0.1*(x[t])^3, x == x0, v == v0}
As before we define a function which will produce, but not display, the phase plane plot for given initial position $$(x_0, y_0 ) .$$
graph[x0_, v0_] :=
Module[{xans, vans},
{xans, vans} = {x[t], v[t]} /. Flatten[NDSolve[duff[x0, v0], {x[t], v[t]}, {t, 0, 8}]];
ParametricPlot[{xans, vans}, {t, 0, 8}, PlotRange -> {{-5.1, 5.1}, {-5, 5}},
AxesLabel -> {"x", "v"}, AspectRatio -> 1, DisplayFunction -> Identity,
PlotStyle -> RGBColor[1, 0, 0]]];
We try this for one initial condition
Show[graph[3.5, 0], DisplayFunction -> $DisplayFunction] We get a closed curve, which is the signature of a periodic solution. Now we construct a family of solutions, using the list of initial conditions given below. initlist = {{0.5, 0}, {1, 0}, {1.5, 0}, {2, 0}, {2.5, 0}, {3, 0}, {3.5, 0}, {4, 0}}; Now we produce all of the graphs without showing them individually Module[{i, x0temp, v0temp, pair, newgraph}, graphlist = {}; Do[pair = initlist[[i]]; x0temp = First[pair]; v0temp = Last[pair]; newgraph = graph[x0temp, v0temp]; graphlist = Append[graphlist, newgraph], {i, 1, 8}]] Now we plot them. solgraph = Show[graphlist, DisplayFunction ->$DisplayFunction] Finally we combine the solutions and the direction field.
Show[dirfield, solgraph] Our graph shows that this interesting system has the property that every initial condition gives rise to a periodic solution. Of course if we remember that the equations describe an undamped spring mass system, such behavior is not so surprising. The equilibrium point at the origin is enclosed by the periodic solutions, but the system does not go to the equilibrium point. If we introduced some damping, all of the closed curves would change to spirals which would asymptotically approach the equilibrium point.     ■

Example 2: Consider the system of nonlinear differential equations
$\begin{split} \dot{x} &= -y - x^2 , \\ \dot{y} &= 2x - y^3 . \end{split}$
According to Mathematica, the iven system has two critical points: the origin and
$\left( - 2^{1/5} , - 2^{2/5} \right) \approx \left( -1.1487 , - 1.31951 \right) .$
Solve[2*x + x^6 == 0, x]
{{x -> 0}, {x -> (-2)^( 1/5)}, {x -> -2^(1/5)}, {x -> -(-1)^(2/5) 2^(1/5)}, {x -> (-1)^( 3/5) 2^(1/5)}, {x -> -(-1)^(4/5) 2^(1/5)}}
NSolve[2*x + x^6 == 0, x]
{{x -> -1.1487}, {x -> -0.354967 - 1.09248 I}, {x -> -0.354967 + 1.09248 I}, {x -> 0.}, {x -> 0.929316 - 0.675188 I}, {x -> 0.929316 + 0.675188 I}}
Upon plotting the direction field, we see that the origin is the center, but (-1.1487, -1.3195) is an unstable saddle point. StreamPlot[{-y - x^2, 2*x - y^3}, {x, -2.5, 2.5}, {y, -2.5, 2.5}, StreamScale -> Large]
Then we solve the system numerically and plot x versus y to identify the limit cycle. A = StreamPlot[{-y - x^2, 2*x - y^3}, {x, -2.5, 2.5}, {y, -2.5, 2.5}, StreamScale -> Large, PlotRange -> {{3, -3}, {3, -3}}]; sol = NDSolve[{x'[t] == -y[t] - (x[t])^2, y'[t] == 2*x[t] - (y[t])^3, x == 1, y == 1}, {x, y}, {t, 0, 150}]; B = ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, -.1, 100}, PlotStyle -> {Blue, Thick}, PlotRange -> {{3, -3}, {3, -3}}] ; B2 = ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 95, 100}, PlotStyle -> {Red, Thickness[.012]}, PlotRange -> {{3, -3}, {3, -3}}]; Show[A, B, B2, PlotRange -> {{3, -3}, {3, -3}}] Show[A, B, B2, PlotRange -> {{.5, -.5}, {.5, -.5}}]
■

Example 3: Consider the Rabinovich-- Fabrikant system
\begin{equation} \label{EqChaos.1} \begin{split} \dot{x} &= y \left( z-1 + x^2 \right) + \gamma \,x , \\ \dot{y} &= x \left( 3z + 1 - x^2 \right) + \gamma \, y , \\ \dot{z} &= -2z \left( \alpha + xy \right) , \end{split} \end{equation}
where α and γ are constants that control the evolution of the system. We plot solutions to RF system when α = 1.14 and γ = 0.1 with the initial conditions $x(0) = -1, \quad y(0) = 0, \quad z(0) = 0.5.$ The corresponding limit circle (in red color) is clearly visible. a = 1.14; gamma = 0.1; sol = NDSolve[{x'[t] == y[t]*(z[t] - 1 + (x[t])^2) + gamma*x[t], y'[t] == x[t]*(3*z[t] + 1 - (x[t])^2) + gamma*y[t], z'[t] == -2*z[t]*(a + x[t]*y[t]), x == -1, y == 0, z == 0.5}, {x, y, z}, {t, 0, 200}]; pp = ParametricPlot3D[ Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 100}, PlotRange -> All, PlotTheme -> "Monochrome"]; lim = ParametricPlot3D[ Evaluate[{x[t], y[t], z[t]} /. sol], {t, 35, 100}, PlotRange -> All, PlotTheme -> "Monochrome", PlotStyle -> {Red, Thickness[.005]}]; Show[lim, pp] Rabinovich–Fabrikant system Mathematica code

■

# 2.3.6. Applications

F[x_, y_] := -y^3;
G[x_, y_] := x^3;
sol = NDSolve[{x'[t] == F[x[t], y[t]], y'[t] == G[x[t], y[t]],
x == 1, y == 0}, {x, y}, {t, 0, 3*Pi},
WorkingPrecision -> 20]
ParametricPlot[Evaluate[{x[t], y[t]}] /. sol, {t, 0, 3*Pi}]

Out=
{{x -> InterpolatingFunction[{{0, 9.4247779607693797154}}, <>],
y -> InterpolatingFunction[{{0, 9.4247779607693797154}}, <>]}} X[t_] := Evaluate[x[t] /. sol]
Y[t_] := Evaluate[y[t] /. sol]
fns[t_] := {X[t], Y[t]};
len := Length[fns[t]];
Plot[Evaluate[fns[t]], {t, 0, 3*Pi}]

## V. Pendulum Equation

fx = Function[{x, y}, y]; fy = Function[{x, y}, -Sin[x] - 0.25*y]
portrait =
StreamPlot[{fx[x, y], fy[x, y]}, {x, -6, 6}, {y, -3, 3},
AspectRatio -> Automatic]
solution =
Function[point,
Function @@ {t, ({x[t], y[t]} /.
NDSolve[{x'[t] == fx[x[t], y[t]], y'[t] == fy[x[t], y[t]],
Thred[{x[time], y[time]} == point]}, {x, y}, {t, time,
time + 40}])[]}]

Out=
Function[point,
Function @@ {t, ({x[t], y[t]} /.
NDSolve[{Derivative[x][t] == fx[x[t], y[t]],
Derivative[y][t] == fy[x[t], y[t]],
Thred[{x[time], y[time]} == point]}, {x, y}, {t, time,
time + 40}])[]}] ## V. van der Pol Equation

splot = StreamPlot[{y, (1 - x^2) y - x}, {x, -4, 4}, {y, -3, 3}, StreamColorFunction -> "Rainbow"]; Manipulate[ Show[splot, ParametricPlot[ Evaluate[ First[{x[t], y[t]} /. NDSolve[{x'[t] == y[t], y'[t] == y[t] (1 - x[t]^2) - x[t] + 0.5 Cos[1.1 t], Thread[{x, y} == point]}, {x, y}, {t, 0, T}]]], {t, 0, T}, PlotStyle -> Red]], {{T, 20}, 1, 100}, {{point, {3, 0}}, Locator}, SaveDefinitions -> True]
Or just to show off (again a rip off from the documentation)
splot = LineIntegralConvolutionPlot[{{y, (1 - x^2) y - x}, {"noise", 1000, 1000}}, {x, -4, 4}, {y, -3, 3}, ColorFunction -> "BeachColors", LightingAngle -> 0, LineIntegralConvolutionScale -> 3, Frame -> False]; Manipulate[ Show[splot, ParametricPlot[ Evaluate[ First[{x[t], y[t]} /. NDSolve[{x'[t] == y[t], y'[t] == y[t] (1 - x[t]^2) - x[t]+0.5 Cos[1.1 t], Thread[{x, y} == point]}, {x, y}, {t, 0, T}]]], {t, 0, T}, PlotStyle -> White]], {{T, 20}, 1, 100}, {{point, {3, 0}}, Locator}, SaveDefinitions -> True]
You can solve the equation with (you might want to change the initial conditions) :
sol[t_] = NDSolve[{x''[t] - (1 - x[t]^2) x'[t] + x[t] == 0.5 Cos[1.1 t], x == 0, x' == 1}, x[t], {t, 0, 10}][[1, 1, 2]]
Now you can use the solution as any other function; in particular, you can plot it versus its derivative :
ParametricPlot[{sol[t], sol'[t]}, {t, 0, 10}]

1. D’Acunto M., Determination of limit cycles for a modified van der Pol oscillator. Mechanics Research Communications, 2006, 33:93–100. [doi:10.1016/j.mechrescom.2005.06.012]
2. Gottlieb, H.P.W., 2006. Harmonic balance approach to limit cycles for nonlinear jerk equations. Journal of Sound and Vibration, 297:243–250. [doi:10.1016/j.jsv.2006.03.047]
3. He, J.H., 2003. Determination of limit cycles for strongly nonlinear oscillators. Physical Review Letters, 90(17), Article No. 174301.
4. He, J.H., 2006a. Determination of limit cycles for strongly nonlinear oscillators. Physic Review Letter, 90:174–181.