我正在研究Strogatz的《非线性动力学和混沌》,在第二章练习2.8.1中遇到了麻烦。 (教育工作者的旗帜:我毕业了,所以这不是一堂课,我只是想重新回到微分方程的数值求解中),这是一个非常简单的微分方程,在不同的初始条件下,我可以绘制各个解曲线但我正在尝试使用颤动或流图将单个解决方案叠加在向量字段之上。

我的问题在于理解如何将here中发现的dy / dx形式的矢量问题的矢量场图转换成Strogatz的书主要解决的dx / dt形式。

鉴于在logistic函数中定义的x向量仅是一维的,因此我很难理解为什么u和v如何在颤动或流图中表示,因为问题似乎只有u流动。这可能是超级容易的,并且正在考虑中,但是任何指导或帮助将不胜感激!

到目前为止,我有以下内容:

# 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

最佳答案

要从导数(例如dx / dt)绘制斜率,可以先找到dx / dt,然后使用固定的dt计算dx。然后,在每个感兴趣的(t,x)处绘制从(t,x)到(t + dt,x + dx)的小线段。

这是您的方程dx / dt = x(1-x)的示例。 (Strogatz图片中没有箭头,因此我也将其删除了。)

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()

关于python - 如何使用matplotlib为逻辑方程K = 1生成矢量场图?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/26834173/

10-11 22:48
查看更多