python - 使用 matplotlib 绘制 Duffing 振荡器的吸引力盆地
问题描述
我有一个常微分方程系统,它有两个吸引子,一个在 (1, 0),另一个在 (-1, 0)。我想在笛卡尔坐标中绘制一个吸引盆,其中有两种颜色,显示随着时间趋于正无穷大,每个坐标点上的一个点最终会成为哪个吸引子。但是,我不知道如何用 matplotlib 绘制这样的图。这是我现在所做的:
from scipy.integrate import ode
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
"""
The system of ODE:
x' = y
y' = x - x**3 - gamma*y
"""
# The system of equation
def f(t, r, arg):
return [r[1], r[0] - r[0] ** 3 - arg * r[1]]
# The Jacobian matrix
def jac(t, r, arg):
return [[0, 1], [1 - 3 * r[0] ** 2, -arg]]
# r is the vector (x,y)
# Initial condition, length of time evolution, time step, parameter gamma
def solver(r0, t0, t1, dt, gamma):
solution = ode(f, jac).set_integrator('dopri5')
# Set the value of gamma
solution.set_initial_value(r0, t0).set_f_params(gamma).set_jac_params(gamma)
return solution
# The function to find the fixed point each starting point ends at
def find_fp(r0, t0, t1, dt, gamma):
solution = solver(r0, t0, t1, dt, gamma)
error = 0.01
while solution.successful():
if norm(np.array(solution.integrate(solution.t+dt)) - np.array([1, 0])) < error:
return 1
elif norm(np.array(solution.integrate(solution.t+dt)) - np.array([-1, 0])) < error:
return -1
def fp(i, j, gamma):
t0, t1, dt = 0, 10, 0.1
return find_fp([i, j], t0, t1, dt, gamma)
我已经定义了几个函数。f
是定义方程组的函数,即系统jac
的雅可比矩阵,用作使用(Kutta-Runge 方法)dopri5
方法求解 ODE 的参数。scipy.integrate.ode
该find_fp
函数被定义为返回相空间中一个点将结束的吸引子,返回值1
表示该点将结束于(1, 0),并结束-1
于(-1, 0)。到目前为止,这些功能似乎运行良好。但是,我不知道如何使用我对matplotlib
模块所做的事情来绘制一个吸引力盆地。有什么好的想法吗?
解决方案
Quick'n'dirty:选择靠近静止点的初始点并向后计算一段时间的解。绘制解决方案并根据吸引子对它们进行着色。
gamma = 1.2
def solution(x,y):
return odeint(lambda u,t: -np.array(f(t,u,gamma)), [x,y], np.arange(0,15,0.01))
for i in range(-10,11):
for j in range(-10,11):
sol = solution(-1+i*1e-4, 0+j*1e-4);
x,y = sol.T; plt.plot(x,y,'.b',ms=0.5);
sol = solution(+1+i*1e-4, 0+j*1e-4);
x,y = sol.T; plt.plot(x,y,'.r',ms=0.5);
plt.grid(); plt.show();
这给出了图像
其他值gamma
或更长的积分间隔需要小心处理,因为三次项会导致反向时间积分中的超指数爆炸。
推荐阅读
- list - 放入ZStack时swiftui列表不起作用
- html - 尽管我的文件在同一个工作区中,但我在 vs 代码上不断收到“无法打开文件...找不到文件”消息
- django - 如何设置灵活的 django 模型
- c# - C# 将标题和新列添加到 CSV 文件
- kotlin - net.corda.core.CordaRuntimeException: kotlin.KotlinNullPointerException on startTrackedFlow returnValue
- list - 如何使用 katalon 工具和 groovy 脚本在 Web 服务请求的正文中附加字符串列表
- c - 如何允许用户使用C同时输入字符串和int
- java - 我的 Vector 丢失数据的任何原因?
- java - 如何使用可变参数作为 newInstance() 的参数来调用构造函数并创建实例?
- node.js - channel.send 不是函数吗?(不和谐.JS)