Autonomous nonlinear systems

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

matplotlib notebook

Plot direction field and some solutions

Nonlinear pendulum

# ADJUST: specify window for the plot
wndw = [-8, 8, -7, 7];
# ADJUST: specify initial points for solutions to plot
ics = [[1,0], [3,0], [-3,4]];
w2 = 10; gamma = 0.25;
RHS = lambda Y, T: [Y[1], -w2*np.sin(Y[0]) - gamma*Y[1]];

# Plot the direction field
X, Y = np.meshgrid(np.linspace(wndw[0],wndw[1],num=25), np.linspace(wndw[2],wndw[3],num=25))
DXY = RHS([X, Y], 0);
DX = DXY[0];
DY = DXY[1];
# ADJUST: power between 0 (no scaling) and 0.5 (all arrows equal)
nn = (DX**2 + DY**2)**(0.4) + 1e-7
# Plot the result
plt.figure();
plt.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis(wndw);

# Compute solution numerically
tt = np.linspace(0, 5, 500)
for ic in ics:
    sol = scipy.integrate.odeint(RHS, ic, tt);
    plt.plot(sol[:,0], sol[:,1], 'r')

Exercise 9.3.5

# ADJUST: specify window for the plot
wndw = [-2, 2, -1, 3];
# ADJUST or REMOVE: list of critical points
CPs = [[0,0], [1,0], [0,3/2], [-1,2]];
# ADJUST: specify initial points for solutions to plot
ics = [[0,0.01], [-0.3,0.1], [-0.28,0.1], [2,4], [0.2, 0.0001], [2,0.001]];
# x is Y[0], y is Y[1]
RHS = lambda Y, T: [Y[0] - Y[0]**2 - Y[0]*Y[1], 3*Y[1]-Y[0]*Y[1]-2*Y[1]**2];

# Plot the direction field
X, Y = np.meshgrid(np.linspace(wndw[0],wndw[1],num=25), np.linspace(wndw[2],wndw[3],num=25))
DXY = RHS([X, Y], 0);
DX = DXY[0];
DY = DXY[1];
# ADJUST: power between 0 (no scaling) and 0.5 (all arrows equal)
nn = (DX**2 + DY**2)**(0.45) + 1e-7
# Plot the result
plt.figure();
plt.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis(wndw);

# Compute solution numerically
tt = np.linspace(0, 5, 500)
for ic in ics:
    sol = scipy.integrate.odeint(RHS, ic, tt);
    plt.plot(sol[:,0], sol[:,1], 'r')
    
for cp in CPs:
    plt.plot(cp[0], cp[1], 'k*')

# ADJUST: specify window for the plot
wndw = [-1, 5, -3, 2];
# ADJUST or REMOVE: list of critical points
CPs = [[0,0], [2,1], [4,-2], [2,-2]];
# ADJUST: specify initial points for solutions to plot
ics = [[5, -0.7], [5, 0.2], [0,0.01],[-1,0.35], [-1,0.25], [4,1], [4, 0.71]];
# x is Y[0], y is Y[1]
RHS = lambda T, Y: [(2+Y[1])*(2*Y[1]-Y[0]), (2-Y[0])*(2*Y[1]+Y[0])];

# Plot the direction field
X, Y = np.meshgrid(np.linspace(wndw[0],wndw[1],num=25), np.linspace(wndw[2],wndw[3],num=25))
DXY = RHS(0, [X, Y]);
DX = DXY[0];
DY = DXY[1];
# ADJUST: power between 0 (no scaling) and 0.5 (all arrows equal)
nn = (DX**2 + DY**2)**(0.45) + 1e-7
# Plot the result
plt.figure();
plt.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis(wndw);

# Compute solution numerically
tt = np.linspace(0, 5, 500)
for ic in ics:
    sol = scipy.integrate.solve_ivp(RHS, [0,5], ic, max_step=0.02);
    plt.plot(sol.y[0,:], sol.y[1,:], 'r')
    
for cp in CPs:
    plt.plot(cp[0], cp[1], 'k*')