如何用Plots.jl绘制符合DynamicalSystems.jl签名的向量场

问题根源 #

原plot_flow!函数的核心问题在于：

• Plots.jl的quiver!要求坐标点与对应向量为一维扁平化数组，但原代码中us/vs是二维数组，与一维的xspace/yspace维度不匹配
• 嵌套循环生成向量场的方式冗余且效率低

修正后的完整代码 #

using DynamicalSystems
using Plots
using UnPack
using LaTeXStrings

"""
绘制双稳态系统的向量场（箭头图）
"""
function plot_flow!(myplot, α, β, xspace = 0:0.2:2, yspace = 0:0.2:2)
    # 生成网格坐标点的一维扁平化数组
    xs, ys = repeat(xspace, inner=length(yspace)), repeat(yspace, outer=length(xspace))
    # 用广播批量计算每个点的向量
    flows = cellular_switch_rule.([[x, y] for x in xspace, y in yspace], Ref(Dict(:α => α, :β => β)), nothing)
    # 将二维向量数组转为一维，匹配坐标维度
    us = vec([flow[1] for flow in flows])
    vs = vec([flow[2] for flow in flows])
    # 绘制向量场，缩放箭头避免重叠
    quiver!(myplot, xs, ys, quiver=(us, vs), scale=0.5, label="向量场")
end 

"""
    cellular_switch_rule(y, p, t)

大肠杆菌基因表达的生化速率方程，符合DynamicalSystems.jl函数签名要求
"""
function cellular_switch_rule(y, p, t)
    @unpack α, β = p
    u, v = y
    udot = α/(1 + v^β) - u
    vdot = α/(1 + u^β) - v
    return SVector(udot, vdot)
end 

cellular_switch(u0, p) = CoupledODEs(cellular_switch_rule, u0, p)

u_nullcline(v, α, β) = α/(1 + v^β)
v_nullcline(u, α, β) = α/(1 + u^β)

function plot_isoclines!(myplot, α, β, uspace = 0:0.1:2, vspace = 0:0.1:2)
    n_u = [u_nullcline(vi, α, β) for vi in vspace]
    n_v = [v_nullcline(ui, α, β) for ui in uspace]
    plot!(myplot, uspace, n_v, 
        linestyle = :dash,
        label = L"\frac{dv}{dt} = 0 \Leftrightarrow n_v(u)")
    plot!(myplot, n_u, vspace, 
        linestyle = :dash,
        label = L"\frac{du}{dt} = 0 \Leftrightarrow n_u(v)")
end 

# 绘制零倾线+向量场
α = 1
β = 2

myplot = plot(
    xlabel = "u", ylabel = "v", 
    xlim = (0, 2), ylim = (0, 2),
    legend = :right,
    size=(600,600))

plot_isoclines!(myplot, α, β)
plot_flow!(myplot, α, β)
display(myplot)

关键修改说明 #

1. 坐标扁平化：用repeat生成一维坐标数组xs/ys，确保与向量数组维度完全匹配
2. 广播优化：用cellular_switch_rule.()批量计算所有点的向量，替代冗余嵌套循环
3. 向量扁平化：用vec()将二维向量数组转为一维，解决维度不匹配问题
4. 显示优化：增加scale=0.5控制箭头大小，调整采样步长为0.2减少箭头密度，提升图面可读性
5. 强制显示：添加display(myplot)确保脚本运行时绘图窗口正常弹出

流线图替代方案 #

若需绘制流线图而非箭头图，可替换plot_flow!为以下实现：

function plot_flow!(myplot, α, β, xspace = 0:0.05:2, yspace = 0:0.05:2)
    # 定义适配streamplot的匿名函数
    f(u, v) = (α/(1 + v^β) - u, α/(1 + u^β) - v)
    streamplot!(myplot, f, xspace, yspace, label="流线图", linewidth=1)
end

内容的提问来源于stack exchange，提问作者Jared