Linear 2D flows contd

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

In general, for any matrix relating $\mathbf{x}$ to $\dot{\mathbf{x}}$, say $\mathbf{A}$, we can find out the eigenvalues and eigenvectors. Based on that we can analyze various trajectories in the phase space and understand the role of the eigenproperties of the matrix.

def show2dphase(a=1, b = 1, c=0.05, d=1):
    #a = 1; b = 1; c = 4; d = -2;
    A = np.array([[a, b], [c, d]]);
    lam, v = np.linalg.eig(A);
    print(lam);
    print(v[:,0]);
    print(v[:,1]);

    x = np.linspace(-3, 3, 20); y = np.linspace(-3, 3, 20);
    X, Y = np.meshgrid(x, y);
    U = a*X + b*Y;
    V = c*X + d*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 = a*X1 + b*Y1;
    V1 = c*X1 + d*Y1;
    plt.streamplot(X1, Y1, U1, V1, integration_direction='forward', density=1)
    plt.ylim(-3, 3);
    plt.plot(x1, x1*(v[1,0]/v[0,0]));
    plt.plot(x1, x1*(v[1,1]/v[0,1]));

w = interactive(show2dphase, a = (-3, 3), b = (-3, 3), c = (-3, 4,0.01), d = (-3, 3))
w

• We see that for complex eigen values, with positive real part, the values spiral out and the magnitude grows.
• For real part negative, the values spiral in.
• When the real part is zero, we have a closed orbit
• When both eigen values become equal, all the vectors move out from the origin. They either grow exponentially or decay depending upon the sign of eigen value.

• Let us now plot the trajectory for a=1, b=1, c=-2, d=4
• We make use of solve_ivp routine to solve the initial value problem
  • We provide the initial condition and the time span

from scipy.integrate import solve_ivp

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

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

    tspan = [0, 1]
    #x0 = -2.5; 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);
    plt.xlim(-3,3); plt.ylim(-3, 3);
    plt.plot(x, x, '--k');
    plt.plot(x, -4*x, '-k');
    plt.plot(x0, y0, 'ok');
w = interactive(showtraj, x0 = (-2.5, 2.5, 0.05), y0 = (-2.5, 2.5, 0.05));
w