>>> from env_helper import info; info()
页面更新时间: 2023-06-24 12:31:34
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

2.6. 解常微分方程组

scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。下面让我们来看看如何用odeint 计算洛仑兹吸引子的轨迹。洛仑兹吸引子由三个微分方程定义:

洛仑兹吸引子三个微分方程

图 2.1 洛仑兹吸引子三个微分方程

洛仑兹吸引子的详细介绍: http://bzhang.lamost.org/website/archives/lorenz_attactor

这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分,就 可以计算出无质量点在此空间中的运动轨迹。其中 σ, ρ, β 为三个常数,不同的参数可以计算出不同的 运动轨迹: x(t), y(t), z(t)。 当参数为某些值时,轨迹出现馄饨现象:即微小的初值差别也会显著地影 响运动轨迹。下面是洛仑兹吸引子的轨迹计算和绘制程序:

>>> %matplotlib inline
>>> from scipy.integrate import odeint
>>> import numpy as np
>>> def lorenz(w, t, p, r, b):
>>>     # 给出位置矢量w,和三个参数p, r, b计算出
>>>     # dx/dt, dy/dt, dz/dt的值
>>>     x, y, z = w
>>>     # 直接与lorenz的计算公式对应
>>>
>>>     return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
>>>
>>> t = np.arange(0, 30, 0.01) # 创建时间点
>>> # 调用ode对lorenz进行求解, 用两个不同的初始值
>>> track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
>>> track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
>>>
>>> # 绘图
>>> from mpl_toolkits.mplot3d import Axes3D
>>> import matplotlib.pyplot as plt
>>>
>>> fig = plt.figure()
>>> ax = Axes3D(fig)
>>> ax.plot(track1[:,0], track1[:,1], track1[:,2])
>>> ax.plot(track2[:,0], track2[:,1], track2[:,2])
>>> plt.show()
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7fa457db5f50>]
<Figure size 640x480 with 0 Axes>

用odeint函数对洛仑兹吸引子微分方程进行数值求解所得到的运动轨迹

我们看到即使初始值只相差0.01,两条运动轨迹也是完全不同的。 在程序中先定义一个lorenz函数,它的任务是计算出某个位置的各个方向的微分值,这个计算直接根 据洛仑兹吸引子的公式得出。然后调用odeint,对微分方程求解,odeint有许多参数,这里用到的四 个参数分别为: 1. lorenz, 它是计算某个位移上的各个方向的速度(位移的微分) 2. (0.0, 1.0, 0.0),位移初始值。计算常微分方程所需的各个变量的初始值 3. t, 表示时间的数组,odeint对于此数组中的每个时间点进行求解,得出所有时间点的位置 4. args, 这些参数直接传递给lorenz函数,因此它们都是常量