Phase portraits in nonlinear systems

import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
from scipy.integrate import solve_ivp
from mayavi import mlab

• In the previous lectures we have looked at behaviour of linear system and their phase plots and dependance of stability on the eigen values

• We now look at nonlinear systems

• Let us analyze the following system: $$\dot{x} = x + e^{-y}$$ $$\dot{y} = -y$$

• We reuse the code from prev lecture and plot the flow and trajectories

x = np.linspace(-3, 3, 20); y = np.linspace(-3, 3, 20);
X, Y = np.meshgrid(x, y);
U = X + np.exp(-Y);
V =-Y;
plt.quiver(X, Y, U, V);
ax = plt.gca(); ax.set_aspect(1);

x1 = np.linspace(-3, 3, 100); y1 = np.linspace(-3, 3, 100);
X1, Y1 = np.meshgrid(x1, y1);
U1 = X1 + np.exp(-Y1);
V1 =-Y1;
#seedpoints = np.array([[0, 1.5], [2.5, 0]]);
plt.streamplot(X1, Y1, U1, V1, integration_direction='forward', density=1.5)
plt.ylim(-3, 3);
plt.contour(X1, Y1, U1, levels=[0.0]);
plt.contour(X1, Y1, V1, levels=[0.0],colors='g');

• Let us now consider the following system of non-linear equations: $$ \dot{x} = -x + x^3$$ $$ \dot{y} = -2y$$
• The fixed points are (0,0), (-1,0),(1,0)
• We can analyze the behaviour near the fixed points by linearizing the equations ,forming the Jacobian and finding the eigen values and eigen vectors

def show_traj(x0=0.9,y0=1.0):
    x = np.linspace(-3, 3, 20); y = np.linspace(-3, 3, 20);
    X, Y = np.meshgrid(x, y);
    U = -X+X**3;
    V =-2*Y;
    plt.quiver(X, Y, U, V);
    ax = plt.gca(); ax.set_aspect(1);

    x1 = np.linspace(-3, 3, 100); y1 = np.linspace(-3, 3, 100);
    X1, Y1 = np.meshgrid(x1, y1);
    U1 = -X1 + X1**3;
    V1 =-2*Y1;
    #seedpoints = np.array([[0, 1.5], [2.5, 0]]);
    plt.streamplot(X1, Y1, U1, V1, integration_direction='forward', density=1)
    plt.ylim(-3, 3);



    def mysys(t, x): # returns the RHS
        return [-x[0] + x[0]**3,  -2*x[1]];

    tspan = [0,0.7]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True);

    tout = np.linspace(0, np.max(tspan), 100);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    plt.plot(xout, yout)
w = interactive(show_traj, x0 = (-1.5, 1.5, 0.1), y0 = (-1.5, 1.5, 0.1))
w

• Let us consider the problem: $$ \dot{x} = -y + ax(x^2+y^2)$$ $$ \dot{y} = x + ay(x^2 + y^2)$$
• The fixed point of the system is (0,0). When a=0 the trajectory is a closed path.
• When a>0, the trajectory spirals out and when a<0, the trajectory spirals in towards the fixed point.
• But if we linearize this system, then the eigen vectors indicate that the fixed point should have closed orbits around it irrespective of the value of a.
• So this case shows us that we cannot always linearize the equation to get the full picture.

def show_traj(x0=0.9,y0=1.0, a = 0.0):
    x = np.linspace(-3, 3, 20); y = np.linspace(-3, 3, 20);
    X, Y = np.meshgrid(x, y);
    U = -Y+a*X*(X**2+Y**2);
    V = X + a*Y*(X**2 + Y**2);
    plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    ax = plt.gca(); ax.set_aspect(1);

    def mysys(t, x): # returns the RHS
        return [-x[1] + a*x[0]*(x[0]**2 + x[1]**2),  x[0] + a*x[1]*(x[0]**2 + x[1]**2)];

    tspan = [0,4]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True);

    tout = np.linspace(0, np.max(tspan), 100);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    plt.plot(xout, yout)
w = interactive(show_traj, x0 = (-1.5, 1.5, 0.1), y0 = (-1.5, 1.5, 0.1), a = (-0.4, 0.4, 0.05))
w

• We now look at a conservative nonlinear system.
• We take a case of a oscillator with a non linear forcing function $$ m\ddot{x} = F $$ $$ F = x-x^3$$
• Let $y = \dot{x}$, now we have the following system of equations: $$ \dot{x} = y$$ $$ \dot{y} = x-x^3$$
• For this system we have 3 fixed points (0,0),(-1,0) and (1,0)
• We can analyze the behaviour near fixed points by making use of Jacobians

def show_traj(x0=0.9,y0=1.0):
    x = np.linspace(-3, 3, 15); y = np.linspace(-3, 3, 15);
    X, Y = np.meshgrid(x, y);
    U = Y;
    V = X-X**3;
    plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    ax = plt.gca(); ax.set_aspect(1);

    def mysys(t, x): # returns the RHS
        return ([x[1], x[0]-x[0]**3]);

    tspan = [0,20]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True);

    tout = np.linspace(0, np.max(tspan), 100);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    
    
    plt.plot(xout, yout)
w = interactive(show_traj, x0 = (-1.5, 1.5, 0.1), y0 = (-1.5, 1.5, 0.1))
w

• The energy of the system remains constant

def mysys(t, x): # returns the RHS
    return ([x[1], x[0]-x[0]**3]);

tspan = [0,20]
x0 = -1.1; y0 = 2.5
ics = [x0, y0]
sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

tout = np.linspace(0, np.max(tspan), 200);
xout = sol.sol(tout)[0];
yout = sol.sol(tout)[1];
E = 1/2*yout**2 - 1/2*xout**2 + 1/4*xout**4;

#plt.plot(xout, yout)
plt.plot(tout, E); 
plt.ylim(1, 3)

(1.0, 3.0)

• We now vary the initial y coordinate and while keeping x coodinate equal to 1
• We plot in 3D by taking the z axis as the energy of the contour that the particle follows
  • We color code the orbits according to energy level

def mysys(t, x): # returns the RHS
    return ([x[1], x[0]-x[0]**3]);



tspan = [0,20]

for y0 in np.linspace(0.1, 2.5, 30):

    x0 = -1.1;
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

    tout = np.linspace(0, np.max(tspan), 200);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    E = 1/2*yout**2 - 1/2*xout**2 + 1/4*xout**4;
    print(np.mean(E))
    #plt.plot(xout, yout)
    mlab.plot3d(xout, yout, E, color=(1/3.3*(np.mean(E)+0.3),1/3.3*(np.mean(E)+0.3),1-1/3.3*(np.mean(E)+0.3)))
    
mlab.show()

-0.23397561489269308
-0.22227546486453623
-0.20372641162887078
-0.17832827181792993
-0.14608047337486074
-0.1069837956706559
-0.061038589614078015
-0.008244054583729031
0.05139969792260138
0.11789225271146407
0.19123364880943797
0.27142389908567005
0.35846287585896996
0.4523512207578056
0.5530878351597842
0.6606742630205524
0.7751093029198924
0.8963935242642245
1.0245263741712352
1.1595086701142945
1.3013397392966255
1.4500196354592356
1.6055485985385343
1.7679265531714963
1.9371536360452677
2.1132296873383094
2.296154744755877
2.485928679617646
2.6825514636532533
2.886023361384373