【混动动力学引论】暗线方程图像的绘制

求解超稳定周期轨道的mu值表达式的步骤

0. 声明迭代方程$f(\mu ,x)=1-\mu x^2$
1. 将$f(\mu ,x)$复合$p$次（$p$为周期），得$f^{(p)}(\mu ,x)$
2. 求解$\frac{\partial f^{(p)}(\mu ,x)}{\partial x}=0$方程的根，得出临界点$O$（多个）
3. 将临界点逐一代入$f^{(p)}(\mu ,x)$，求关于$\mu$的方程

绘制暗线方程代码

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# 声明变量
x = sp.symbols('x')
mu = sp.symbols('mu')
# 迭代方程
fx = 1 - mu * x**2

# 求给定超稳定周期轨道的mu值表达式
def transStableTrackMu(fx, x, period):
    # 复合n次
    compositeFx = fx
    for i in range(1, period):
        compositeFx = compositeFx.replace(x, fx)
    # 一阶偏导数
    dif = compositeFx.diff(x)
    # 求解临界点方程
    boundaryPoints = sp.solve(dif, x)
    return [compositeFx.replace(x, fp) for fp in boundaryPoints]

# 超稳定轨道表达式列表
muExprList = transStableTrackMu(fx, x, 6) # 这里是画周期1~6的，如果画更多周期的，则需要把6改成指定周期数
# 超稳定轨道函数
muFuncList = [sp.lambdify(mu, muExpr) for muExpr in muExprList]

# 画出暗线方程图像
for muFunc in muFuncList:
    muVar = np.linspace(0,2, 100)
    yVar = [muFunc(mu) for mu in muVar]
    plt.plot(muVar, yVar)

plt.show(!)

运行所得图像：