Search code examples
pythonmatplotlibnumericdifferential-equations

How do I generate a vector field plot for logistic equation K = 1 using matplotlib?


I'm going through Strogatz's Nonlinear Dynamics and Chaos and I've hit a snag in chapter 2 Exercise 2.8.1. (Educator flag: I've graduated so this isn't for a class, I'm just trying to get back into the numerical solving of differential equations) It's a pretty simple differential equation and I can plot individual solution curves given different initial conditions but I'm trying to use quiver or streamplot to superimpose individual solutions on top of the vector field.

My problem is in understanding how to translate the vector field plots for similar problems in the dy/dx form found here over to the dx/dt form that's primarily tackled in Strogatz's book.

Given that the x vector that's defined in the logistic function is only one dimensional I'm having a hard time reasoning out how express the u and v flows in quiver or streamplot since the problem only seems to have a u flow. It's probably super easy and is being over-thought but any guidance or assistance would be much appreciated!

So far I have the following:

# 2.8.1
# Plot the vector field and some trajectories for xdot = x(1-x) given 
# some different initial conditions for the logistic equation with carrying 
# capacity K = 1

# dx/dt = x(1-x)  

# Imports:
from __future__ import division
from scipy import *
import numpy as np
import pylab
import matplotlib as mp
from matplotlib import pyplot as plt  
import sys
import math as mt

def logistic(x,t):
    return np.array([x[0]*(1-x[0])])

def RK4(t0 = 0, x0 = np.array([1]), t1 = 5 , dt = 0.01, ng = None):  
    tsp = np.arange(t0, t1, dt)
    Nsize = np.size(tsp)
    X = np.empty((Nsize, np.size(x0)))
    X[0] = x0

    for i in range(1, Nsize):
        k1 = ng(X[i-1],tsp[i-1])
        k2 = ng(X[i-1] + dt/2*k1, tsp[i-1] + dt/2)
        k3 = ng(X[i-1] + dt/2*k2, tsp[i-1] + dt/2)
        k4 = ng(X[i-1] + dt*k3, tsp[i-1] + dt)
        X[i] = X[i-1] + dt/6*(k1 + 2*k2 + 2*k3 + k4)
    return X

def tplot():
    t0 = 0
    t1 = 10
    dt = 0.02
    tsp = np.arange(t0,t1,dt)
    X = RK4(x0 = np.array([2]), t1 = 10,dt = 0.02, ng = logistic)
    Y = RK4(x0 = np.array([0.01]), t1 = 10,dt = 0.02, ng = logistic)
    Z = RK4(x0 = np.array([0.5]), t1 = 10,dt = 0.02, ng = logistic)
    P = RK4(x0 = np.array([3]), t1 = 10,dt = 0.02, ng = logistic)
    Q = RK4(x0 = np.array([0.1]), t1 = 10,dt = 0.02, ng = logistic)
    R = RK4(x0 = np.array([1.5]), t1 = 10,dt = 0.02, ng = logistic)
    O = RK4(x0 = np.array([1]), t1 = 10,dt = 0.02, ng = logistic)
    pylab.figure()
    pylab.plot(tsp,X)
    pylab.plot(tsp,Y)
    pylab.plot(tsp,Z)
    pylab.plot(tsp,P)
    pylab.plot(tsp,Q)
    pylab.plot(tsp,R)
    pylab.plot(tsp,O)
    pylab.title('Logistic Equation - K=1')
    pylab.xlabel('Time')
    pylab.ylabel('Xdot')
    pylab.show()

print tplot()

image here


Solution

  • To graph a slope from a derivative (like, dx/dt), you can first find dx/dt, and then use a fixed dt to calculate dx. Then, at each (t, x) of interest, plot the little line segment from (t,x) to (t+dt, x+dx).

    Here's an example for your equation dx/dt = x(1-x). (The Strogatz picture doesn't have arrowheads so I removed them too.)

    import numpy as np
    import matplotlib.pyplot as plt
    
    times = np.linspace(0, 10,  20)
    x = np.linspace(0 ,2, 20)
    T, X = np.meshgrid(times, x)   # make a grid that roughly matches the Strogatz grid  
    
    dxdt = X*(1-X)            # the equation of interest
    dt = .5*np.ones(X.shape)  # a constant value (.5 is just so segments don't run into each other -- given spacing of times array
    dx = dxdt * dt            # given dt, now calc dx for the line segment
    
    plt.quiver(T, X, dt, dx, headwidth=0., angles='xy', scale=15.)
    plt.show()
    

    enter image description here