using DifferentialEquations

# 最简单的 ODE: dy/dt = -y, y(0) = 1
f(u, p, t) = -u
u0 = 1.0
tspan = (0.0, 5.0)

prob = ODEProblem(f, u0, tspan)
sol = solve(prob)

# 查看结果
sol.t        # 时间点
sol.u        # 解的值
sol(0.5)     # 在 t=0.5 处插值
sol(0.5:0.1:2.0)  # 批量插值

# 绘图
using Plots
plot(sol, xlabel="t", ylabel="y(t)", label="y = e^(-t)", linewidth=2)

# Lotka-Volterra 捕食者-猎物模型
function lotka_volterra!(du, u, p, t)
    α, β, γ, δ = p
    du[1] = α * u[1] - β * u[1] * u[2]    # 猎物
    du[2] = γ * u[1] * u[2] - δ * u[2]    # 捕食者
end

u0 = [1.0, 1.0]        # 初始条件 [猎物, 捕食者]
p = [1.5, 1.0, 3.0, 1.0]   # 参数 [α, β, γ, δ]
tspan = (0.0, 10.0)

prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob, Tsit5())

plot(sol, xlabel="t", ylabel="种群数量",
    label=["猎物" "捕食者"], linewidth=2)

# 默认求解器（自动选择）
sol = solve(prob)

# 指定求解器
sol = solve(prob, Tsit5())         # Tsitouras 5/4 阶，非刚性首选
sol = solve(prob, BS3())           # Bogacki-Shampine 3/2 阶
sol = solve(prob, Vern7())         # Verner 7/6 阶，高精度
sol = solve(prob, DP5())           # Dormand-Prince 5/4 阶
sol = solve(prob, DP8())           # Dormand-Prince 8/7 阶

# 设置精度
sol = solve(prob, Tsit5();
    abstol = 1e-10,     # 绝对误差容限
    reltol = 1e-8,      # 相对误差容限
    dt = 0.01,          # 初始步长
    adaptive = true,    # 自适应步长
    maxiters = 10000    # 最大迭代数
)

using DifferentialEquations, Plots

# Van der Pol 振荡器（刚性示例）
# 当 μ 很大时成为刚性方程
function vanderpol!(du, u, p, t)
    μ = p[1]
    du[1] = u[2]
    du[2] = μ * (1 - u[1]^2) * u[2] - u[1]
end

# 非刚性（μ = 1）
prob1 = ODEProblem(vanderpol!, [2.0, 0.0], (0.0, 10.0), [1.0])
sol1 = solve(prob1, Tsit5())

# 刚性（μ = 1000）
prob2 = ODEProblem(vanderpol!, [2.0, 0.0], (0.0, 3000.0), [1000.0])
sol2 = solve(prob2, Rodas5())     # 必须用刚性求解器

# 尝试用非刚性求解器求解刚性问题（会很慢或失败）
# sol_bad = solve(prob2, Tsit5())   # 不推荐

plot(sol2, xlabel="t", ylabel="u",
    title="Van der Pol 振荡器 (μ=1000)",
    label=["u₁" "u₂"])

using DifferentialEquations

# 示例：球弹跳模型
# 当球触地（y=0）时速度反向
function ball!(du, u, p, t)
    du[1] = u[2]      # 位置
    du[2] = -9.81      # 重力加速度
end

# 事件条件：u[1] == 0（触地）
condition(u, t, integrator) = u[1]

# 事件影响：速度反向并衰减
function affect!(integrator)
    integrator.u[2] = -0.8 * integrator.u[2]   # 恢复系数 0.8
end

cb = ContinuousCallback(condition, affect!)

u0 = [10.0, 0.0]    # 初始高度 10m，初速 0
prob = ODEProblem(ball!, u0, (0.0, 10.0))
sol = solve(prob, Tsit5(); callback=cb)

plot(sol, idxs=1, xlabel="t", ylabel="高度(m)",
    title="弹跳球", label="高度")

# 离散回调：在特定时间点触发
function periodic_affect!(integrator)
    println("t = $(integrator.t): u = $(integrator.u)")
end
cb_periodic = PeriodicCallback(periodic_affect!, 1.0)   # 每 1 秒

# 组合多个回调
cb_combined = CallbackSet(cb, cb_periodic)

using DifferentialEquations, ForwardDiffSensitivity, SciMLSensitivity

# 定义 ODE 问题（带参数）
function f!(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = p[3] * u[1] * u[2] - p[4] * u[2]
end

u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
prob = ODEProblem(f!, u0, (0.0, 5.0), p)

# 解对参数的灵敏度
# Forward 模式（参数少时快）
sol = solve(prob, Tsit5(); sensealg=ForwardSensitivity())

# Adjoint 模式（参数多时快）
sol_adj = solve(prob, Tsit5(); sensealg=InterpolatingAdjoint())

# 计算目标函数对参数的梯度
function loss(p)
    prob_local = remake(prob; p=p)
    sol_local = solve(prob_local, Tsit5(); abstol=1e-6, reltol=1e-6)
    return sum(sol_local[end])    # 目标：最终状态之和
end

dp = ForwardDiff.gradient(loss, p)
println("梯度: $dp")

using DifferentialEquations, Plots

# 几何布朗运动：dS = μS dt + σS dW
function gbm!(du, u, p, t)
    μ, σ = p
    du[1] = μ * u[1]       # 漂移项
end

function gbm_noise!(du, u, p, t)
    σ = p[2]
    du[1] = σ * u[1]       # 扩散项
end

u0 = [100.0]           # 初始股价
p = [0.05, 0.2]        # μ=5%, σ=20%
tspan = (0.0, 1.0)

prob = SDEProblem(gbm!, gbm_noise!, u0, tspan, p)
sol = solve(prob, EM(); dt=0.001)    # Euler-Maruyama 方法

# Monte Carlo 模拟
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob, EM(); dt=0.001, trajectories=1000)

# 绘制 1000 条路径
plot(sim, xlabel="t", ylabel="S(t)",
    title="几何布朗运动 Monte Carlo",
    legend=false, alpha=0.1, linewidth=0.5)

using DifferentialEquations

# 简单的摆锤 DAE（含约束）
function pendulum_dae!(du, u, p, t)
    x, y, vx, vy, λ = u
    g, L = p
    # 微分方程
    du[1] = vx
    du[2] = vy
    du[3] = -λ * x        # 加速度 x 分量
    du[4] = -λ * y - g    # 加速度 y 分量
    # 代数约束：x² + y² = L²
    du[5] = x^2 + y^2 - L^2
end

u0 = [1.0, 0.0, 0.0, 0.0, 0.0]    # [x, y, vx, vy, λ]
p = [9.81, 1.0]    # g, L
differential_vars = [true, true, true, true, false]  # λ 是代数变量

prob = DAEProblem(pendulum_dae!, du0, u0, tspan, p;
                  differential_vars=differential_vars)
sol = solve(prob, DFBDF())

using DifferentialEquations

# 带延迟的种群模型
# du/dt = r * u(t) * (1 - u(t-τ)/K)
function delay_model!(du, u, h, p, t)
    r, K, τ = p
    u_delay = h(p, t - τ)      # 延迟项
    du[1] = r * u[1] * (1 - u_delay[1] / K)
end

# 历史函数（t < 0 时的值）
h(p, t) = [0.1]

u0 = [0.1]
p = [1.0, 1.0, 0.5]   # r, K, τ
tspan = (0.0, 30.0)

prob = DDEProblem(delay_model!, u0, h, tspan, p)
sol = solve(prob, MethodOfSteps(Tsit5()))

plot(sol, xlabel="t", ylabel="u(t)",
    title="延迟微分方程", label="种群")

using DifferentialEquations, BenchmarkTools

# 标准测试问题：Brusselator
function brusselator!(du, u, p, t)
    a, b = p
    du[1] = a + u[1]^2 * u[2] - (b + 1) * u[1]
    du[2] = b * u[1] - u[1]^2 * u[2]
end

u0 = [1.0, 1.0]
p = [1.0, 3.0]
prob = ODEProblem(brusselator!, u0, (0.0, 10.0), p)

# 性能比较
solvers = [Tsit5(), DP5(), Vern7(), BS3()]
names = ["Tsit5", "DP5", "Vern7", "BS3"]

for (solver, name) in zip(solvers, names)
    t = @elapsed for _ in 1:100
        solve(prob, solver; abstol=1e-8, reltol=1e-8)
    end
    println("$name: $(round(t, digits=3))s (100 次)")
end

using DifferentialEquations, Plots

function sir!(du, u, p, t)
    S, I, R = u
    β, γ = p
    du[1] = -β * S * I       # 易感者减少
    du[2] = β * S * I - γ * I感染者变化
    du[3] = γ * I             # 康复者增加
end

# 参数
β = 0.3       # 传染率
γ = 0.1       # 康复率
R0 = β / γ    # 基本再生数
println("基本再生数 R₀ = $R0")

u0 = [0.99, 0.01, 0.0]    # 初始：99%易感，1%感染
p = [β, γ]
tspan = (0.0, 100.0)

prob = ODEProblem(sir!, u0, tspan, p)
sol = solve(prob, Tsit5())

plot(sol,
    xlabel="天数", ylabel="人口比例",
    title="SIR 模型 (β=$β, γ=$γ, R₀=$R0)",
    label=["易感者 S" "感染者 I" "康复者 R"],
    linewidth=2
)

# 找到感染高峰
I_vals = [u[2] for u in sol.u]
peak_idx = argmax(I_vals)
println("感染高峰: 第 $(round(sol.t[peak_idx], digits=1)) 天, 比例 $(round(maximum(I_vals)*100, digits=1))%")

using DifferentialEquations

# 三分子反应（Brusselator 变体）
function chemical!(du, u, p, t)
    k1, k2, k3, k4 = p
    A, B, X, Y = u
    du[1] = k1 - k3 * A * X          # A 的消耗
    du[2] = k2 * X - k4 * B          # B 的变化
    du[3] = k3 * A * X - k2 * X + X^2 * Y - X  # X 的反应
    du[4] = k2 * X - X^2 * Y          # Y 的反应
end

u0 = [1.0, 3.0, 0.0, 0.0]
p = [1.0, 1.0, 1.0, 1.0]
prob = ODEProblem(chemical!, u0, (0.0, 50.0), p)
sol = solve(prob, Rodas5())   # 可能是刚性的

plot(sol, xlabel="时间", ylabel="浓度",
    label=["A" "B" "X" "Y"], linewidth=2)

using DifferentialEquations

# RLC 串联电路
# L * dI/dt + R * I + Q/C = V(t)
function rlc!(du, u, p, t)
    Q, I = u
    R, L, C, ω, V0 = p
    V = V0 * sin(ω * t)     # 外加交流电压
    du[1] = I                # dQ/dt = I
    du[2] = (V - R * I - Q / C) / L   # dI/dt
end

R = 10.0      # 电阻 (Ω)
L = 0.1       # 电感 (H)
C = 1e-4      # 电容 (F)
ω = 100.0     # 角频率 (rad/s)
V0 = 10.0     # 电压幅值 (V)

u0 = [0.0, 0.0]    # 初始：Q=0, I=0
p = [R, L, C, ω, V0]
prob = ODEProblem(rlc!, u0, (0.0, 0.5), p)
sol = solve(prob, Tsit5())

plot(sol, idxs=2, xlabel="t(s)", ylabel="I(A)",
    title="RLC 电路电流", label="电流", linewidth=1)