求解一个含时的微分方程并对它的解相对于它的参数用隐函数的方法求数值积分

  1. 我要求解一阶(相对于时间t)微分方程组,它同时依赖于一些参数,比如p,我想对它的解的绝对值相对于参数p求数值定积分,我写了一个程序,但是始终报错,不知道错误在何处,请高人指点一二。下面是元代码:
using DifferentialEquations
using QuadGK
using ForwardDiff

# Constants
g = 0.0045
Δ0 = 0
Ω = 5
ϕ = 0

# Define B3 function
function B3(t, q, α, g, Δ0)
    return -q - g * t - (Δ0 - α * t)
end

# Define the system of ODEs
function odefunc!(du, u, p, t)
    fp, f3, fm = u
    q, α, Ω, ϕ = p
    du[1] = Ω/2 * exp(1im * ϕ) - 1im * B3(t, q, α, g, Δ0) * fp + Ω/2 * exp(-1im * ϕ) * fp^2
    du[2] = B3(t, q, α, g, Δ0) + 2 * 1im * fp * Ω / 2 * exp(-1im * ϕ)
    du[3] = exp(-1im * f3) * Ω / 2 * exp(-1im * ϕ)
end

# Initial conditions and time span
u0 = [0.0, 0.0, 0.0]  # Initial conditions for fp, f3, fm
tspan = (0.0, 10.0)

# Function to solve the ODE for a given parameter p
function solve_ode(p)
    prob = ODEProblem(odefunc!, u0, tspan, p)
    sol = solve(prob, Tsit5())
    return sol
end

# Function to compute the integral of the solution with respect to p
function integrate_solution(p_values, tf, α, Ω, ϕ)
    integrand(p) = begin
        sol = solve_ode([p, α, Ω, ϕ])
        abs(sol(tf)[1])  # Assuming you want to integrate the first component
    end
    result, _ = quadgk(integrand, p_values[1], p_values[end])
    return result
end

# Example usage
p_values = [0.0, 3.0]  # Range of p values
tf = 10.0  # Final time
result = integrate_solution(p_values, tf, g, Ω, ϕ)
println("Result: ", result)

在vscode中运行上述代码,出现下列错误:ERROR: LoadError: InexactError: Float64(2.500000000015625 + 3.7499999999999997e-6im)

这里你的初值似乎应该设置为一个复数向量:u0 = Complex.([0.0, 0.0, 0.0])

2 个赞

谢谢啦,问题就出在这里!解决了。