运行DifferentialEquations.jl官方简谐振荡器示例时绘图报错

问题描述 #

运行DifferentialEquations.jl官方的简谐振荡器示例时，触发BoundsError: attempt to access 1-element Vector{Float64} at index [2]错误，原代码如下：

# Simple Harmonic Oscillator Problem
using OrdinaryDiffEq, Plots

#Parameters
ω = 1

#Initial Conditions
x₀ = [0.0]
dx₀ = [π / 2]
tspan = (0.0, 2π)

ϕ = atan((dx₀[1] / ω) / x₀[1])
A = √(x₀[1]^2 + dx₀[1]^2)

#Define the problem
function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end

#Pass to solvers
prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
sol = solve(prob, DPRKN6())

#Plot
plot(sol, vars = [2, 1], linewidth = 2, title = "Simple Harmonic Oscillator",
    xaxis = "Time", yaxis = "Elongation", label = ["x" "dx"])
plot!(t -> A * cos(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution x")
plot!(t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution dx")

核心错误信息：

ERROR: BoundsError: attempt to access 1-element Vector{Float64} at index [2]
...
[24] top-level scope
    @ c:\data\JuliaLanguage\Teste_DE.jl:25

错误原因分析 #

1. 绘图变量索引逻辑错误：SecondOrderODEProblem的解是ArrayPartition类型，结构为(du, u)，其中每个内部元素是长度为1的向量。原代码中vars = [2,1]试图访问向量的第2个元素，导致越界。
2. 初始条件计算隐患：x₀ = [0.0]会导致ϕ计算时出现除以0的情况，影响解析解的正确性。

修正后的代码 #

# Simple Harmonic Oscillator Problem
using OrdinaryDiffEq, Plots

#Parameters
ω = 1

#Initial Conditions（修正x₀避免除以0）
x₀ = [1.0]
dx₀ = [0.0]
tspan = (0.0, 2π)

#重新计算解析解参数
ϕ = atan((dx₀[1] / ω) / x₀[1])
A = √(x₀[1]^2 + (dx₀[1]/ω)^2)  #修正振幅公式，速度项需除以ω

#Define the problem
function harmonicoscillator(ddu, du, u, ω, t)
    ddu .= -ω^2 * u
end

#Pass to solvers
prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
sol = solve(prob, DPRKN6())

#Plot：修正vars参数，用(时间, 状态索引)指定
plot(sol, vars = [(0,1), (0,2)], linewidth = 2, title = "Simple Harmonic Oscillator",
    xaxis = "Time", yaxis = "Amplitude", label = ["x (位移)" "dx (速度)"])
#添加解析解
plot!(t -> A * cos(ω * t - ϕ), lw = 3, ls = :dash, label = "解析解x")
plot!(t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "解析解dx")

关键修改说明 #

1. 初始条件调整：将x₀设为[1.0]、dx₀设为[0.0]，避免除以0的问题，同时让波形更直观。
2. 绘图变量修正：vars = [(0,1), (0,2)]表示：
   • (0,1)：x轴为时间，y轴为解中的u（位移，对应SecondOrderODEProblem的第二个初始条件x₀）
   • (0,2)：x轴为时间，y轴为解中的du（速度，对应SecondOrderODEProblem的第一个初始条件dx₀）
3. 解析解参数修正：振幅A的正确公式为√(x₀² + (v₀/ω)²)，需将速度项除以角频率ω。

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