Systems of 1st order linear eq¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
In [2]:
matplotlib inline
Plot direction field, eigenvectors, and some solutions¶
Matrix [1 1; 4 1] --- saddle point¶
In [3]:
# ADJUST: Define the matrix right-hand side of the ODE system
P = np.array([[1,1],[4,1]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Find and plot the eigenvectors
r,v = np.linalg.eig(P);
v1 = v[:,0]; v2 = v[:,1];
# Scale to fit the axes
v1 = 2*v1 / np.max(np.abs(v1)); v2 = 2*v2 / np.max(np.abs(v2));
# Draw the eigenvectors
plt.plot([0, v1[0]], [0, v1[1]], 'k', linewidth=2);
plt.plot([0, -v1[0]], [0, -v1[1]], 'k', linewidth=2);
plt.plot([0, v2[0]], [0, v2[1]], 'k', linewidth=2);
plt.plot([0, -v2[0]], [0, -v2[1]], 'k', linewidth=2);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
epss = np.array([0,0.05])
# Basing good starting points on the eigenvectors
ics = [v1+epss, v1-epss, v2+epss, v2-epss, -v1+epss, -v1-epss, -v2+epss, -v2-epss, \
0.2*v1+epss, 0.2*v1-epss, 0.2*v2+epss, 0.2*v2-epss, -0.2*v1+epss, -0.2*v1-epss, -0.2*v2+epss, -0.2*v2-epss];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r')
Matrix [1 -2; 1 4] --- source (unstable node)¶
In [4]:
# ADJUST: Define the matrix right-hand side of the ODE system
P = np.array([[1,-2],[1,4]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Find and plot the eigenvectors
r,v = np.linalg.eig(P);
v1 = v[:,0]; v2 = v[:,1];
# Scale to fit the axes
v1 = 2*v1 / np.max(np.abs(v1)); v2 = 2*v2 / np.max(np.abs(v2));
# Draw the eigenvectors
plt.plot([0, v1[0]], [0, v1[1]], 'k', linewidth=2);
plt.plot([0, -v1[0]], [0, -v1[1]], 'k', linewidth=2);
plt.plot([0, v2[0]], [0, v2[1]], 'k', linewidth=2);
plt.plot([0, -v2[0]], [0, -v2[1]], 'k', linewidth=2);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
epss = np.array([0,0.05])
# Basing good starting points on the eigenvectors
ics = [[-0.4,0], [0.4,0], v1+epss, v1-epss, v2+epss, v2-epss, -v1+epss, -v1-epss, -v2+epss, -v2-epss, \
0.2*v1+epss, 0.2*v1-epss, 0.2*v2+epss, 0.2*v2-epss, -0.2*v1+epss, -0.2*v1-epss, -0.2*v2+epss, -0.2*v2-epss];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r')
Matrix [-2 -1; 3 -6] --- sink (stable node)¶
In [5]:
# ADJUST: Define the matrix right-hand side of the ODE system
P = np.array([[-2, -1], [3, -6]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Find and plot the eigenvectors
r,v = np.linalg.eig(P);
v1 = v[:,0]; v2 = v[:,1];
# Scale to fit the axes
v1 = 2*v1 / np.max(np.abs(v1)); v2 = 2*v2 / np.max(np.abs(v2));
# Draw the eigenvectors
plt.plot([0, v1[0]], [0, v1[1]], 'k', linewidth=2);
plt.plot([0, -v1[0]], [0, -v1[1]], 'k', linewidth=2);
plt.plot([0, v2[0]], [0, v2[1]], 'k', linewidth=2);
plt.plot([0, -v2[0]], [0, -v2[1]], 'k', linewidth=2);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
epss = np.array([0,0.05])
# Carefully chosen starting points
ics = [[-2,2], [-1,2], [0,2], [1,2], [2,-2], [-1,-2], [0,-2], [1,-2]];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r')
Matrix with a zero eigenvalue¶
In [6]:
# ADJUST: Define the matrix right-hand side of the ODE system
P = np.array([[4, -3], [8, -6]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Find and plot the eigenvectors
r,v = np.linalg.eig(P);
v1 = v[:,0]; v2 = v[:,1];
# Scale to fit the axes
v1 = 2*v1 / np.max(np.abs(v1)); v2 = 2*v2 / np.max(np.abs(v2));
# Draw the eigenvectors
plt.plot([0, v1[0]], [0, v1[1]], 'k', linewidth=3);
plt.plot([0, -v1[0]], [0, -v1[1]], 'k', linewidth=3);
plt.plot([0, v2[0]], [0, v2[1]], 'k', linewidth=3);
plt.plot([0, -v2[0]], [0, -v2[1]], 'k', linewidth=3);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
epss = np.array([0,0.05])
# Carefully chosen starting points
ics = [[-2,2], [-1,2], [0,2], [1,2], [2,-2], [-0.7,-2], [-1,-2], [0,-2], [1,-2]];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r')
Spiral point (no reason to plot eigenvectors)¶
In [7]:
# ADJUST: Define the matrix right-hand side of the ODE system
P = np.array([[-0.5, 1], [-1, -0.5]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
# Carefully chosen starting points
ics = [[-2,2], [-1,2], [0,2], [1,2], [2,-2], [-0.7,-2], [-1,-2], [0,-2], [1,-2]];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r')
A system depending on parameter, 7.6.11-15¶
In [8]:
# 7.6.12: changes at a = -sqrt(20), 0, sqrt(20)
# ADJUST: Change a
a = -1
P = np.array([[0, -5], [1, a]])
RHS = lambda Y, T: [P[0,0]*Y[0] + P[0,1]*Y[1], P[1,0]*Y[0] + P[1,1]*Y[1]];
# Plot the direction field
X, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-2,2,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.quiver(X, Y, DX/nn, DY/nn)
plt.grid(); plt.axis([-2,2,-2,2]);
# Find and plot the eigenvectors (if eigenvaluess are real)
r,v = np.linalg.eig(P);
if np.isreal(r[1]):
v1 = v[:,0]; v2 = v[:,1];
# Scale to fit the axes
v1 = 2*v1 / np.max(np.abs(v1)); v2 = 2*v2 / np.max(np.abs(v2));
# Draw the eigenvectors
plt.plot([0, v1[0]], [0, v1[1]], 'k', linewidth=3);
plt.plot([0, -v1[0]], [0, -v1[1]], 'k', linewidth=3);
plt.plot([0, v2[0]], [0, v2[1]], 'k', linewidth=3);
plt.plot([0, -v2[0]], [0, -v2[1]], 'k', linewidth=3);
# Compute solution numerically
tt = np.linspace(0, 20, 500)
epss = np.array([0,0.05])
# Carefully chosen starting points
ics = [[-1.5,1], [-1,1], [-0.5,1], [0,1], [0.5,1], [1,1], [1.5,1], \
[-1.5,-0.1], [-1,-0.1], [-0.5,-0.1], [0,-0.1], [0.5,-0.1], [1,-0.1], [1.5,-0.1]];
for ic in ics:
sol = scipy.integrate.odeint(RHS, ic, tt);
plt.plot(sol[:,0], sol[:,1], 'r', linewidth=1)
Larger system: mechanical vibration 7.6.25¶
In [9]:
m1=1; m2=4/3; k1=1; k2=3; k3 = 4/3;
RHS = lambda Y, T: [Y[2], Y[3], -(k1+k2)/m1*Y[0] + k2/m1*Y[1], k2/m2*Y[0] - (k2+k3)/m2*Y[1]];
tt = np.linspace(0, 20, 500)
sol = scipy.integrate.odeint(RHS, [2,1,0,0], tt);
plt.figure(); plt.grid();
plt.plot(tt, sol[:,0], 'r');
plt.plot(tt, sol[:,1], 'b');
In [10]:
# Decompose the solution into two fundamental modes
# Mode 1
sol1 = scipy.integrate.odeint(RHS, [10/7,10/7,0,0], tt);
plt.figure(); plt.grid();
plt.plot(tt, sol1[:,0], 'r');
plt.plot(tt, sol1[:,1], 'b');
# Mode 2
sol2 = scipy.integrate.odeint(RHS, [4/7,-3/7,0,0], tt);
plt.figure(); plt.grid();
plt.plot(tt, sol2[:,0], 'r');
plt.plot(tt, sol2[:,1], 'b');
# The sum of two modes: the original solution
plt.figure(); plt.grid();
plt.plot(tt, sol1[:,0] + sol2[:,0], 'r');
plt.plot(tt, sol1[:,1] + sol2[:,1], 'b');
In [ ]: