在Julia中求解受迫振子微分方程的函数设置疑问

要解决这个二阶非齐次微分方程，核心是把它转化为一阶微分方程组——这是Julia中DifferentialEquations.jl求解器的标准输入形式。下面是完整的实现步骤：

1. 方程转化 #

原方程：m*y''(t) + c*y'(t) + k*y(t) = -m*g

令：

• y₁(t) = y(t)（位移）
• y₂(t) = y'(t)（速度）

则一阶方程组为：

• y₁'(t) = y₂(t)
• y₂'(t) = (-c*y₂(t) - k*y₁(t) - m*g)/m

你之前把-mg移到左侧的操作没问题，但这依然是非齐次方程（因为存在常数项-mg），不能直接套用齐次方程的解法，必须在一阶系统中保留这个常数项。

2. 完整代码实现 #

首先确保安装了必要的包：

using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")

然后编写求解代码：

using DifferentialEquations, Plots

# 定义微分方程（in-place 形式，符合DifferentialEquations.jl要求）
function forced_oscillator!(du, u, p, t)
    m, c, k, g = p  # 提取参数：质量、阻尼系数、劲度系数、重力加速度
    y, y_prime = u   # 当前状态：位移y，速度y'
    du[1] = y_prime  # 位移的导数是速度
    du[2] = (-c*y_prime - k*y - m*g)/m  # 速度的导数（加速度）由原方程推导而来
end

# 设置参数、初始条件和时间区间
m = 1.0    # 质量
c = 0.5    # 阻尼系数
k = 2.0    # 劲度系数
g = 9.81   # 重力加速度
p = (m, c, k, g)

u0 = [0.0, 0.0]  # 初始状态：[初始位移, 初始速度]
tspan = (0.0, 10.0)  # 求解的时间范围：从t=0到t=10

# 构建ODE问题并求解
prob = ODEProblem(forced_oscillator!, u0, tspan, p)
sol = solve(prob)  # 默认用高效的自适应求解器

# 绘制结果
plot(sol, 
     label=["位移 y(t)" "速度 y'(t)"],
     xlabel="时间 t",
     ylabel="数值",
     title="受迫振子的位移与速度变化")

3. 关键说明 #

• DifferentialEquations.jl要求微分方程函数是in-place的（即修改du数组来存储导数，而非返回新数组），这样能提升性能。
• 参数p把所有常数打包传递，方便后续调整参数时不用修改函数本身。
• 求解器会自动选择合适的算法，如果你需要指定特定求解器（比如Tsit5()），可以写solve(prob, Tsit5())。

内容的提问来源于stack exchange，提问作者empty-void