代码之家  ›  专栏  ›  技术社区  ›  natemcintosh

微分方程把“什么都没有”交给我的微分方程

  •  0
  • natemcintosh  · 技术社区  · 6 年前

    我想解两体微分方程。D级微分方程.jl以前很容易就解决了我的方程,但到了某个时候,它开始交上来了 nothing 对于微分方程的初始条件。

    这是否可能是最近频繁出现的范围界定问题?

    using Unitful, RecursiveArrayTools, DifferentialEquations
    
    function propagateTrajNewWay(df::DataFrame, idx::Int)
        x  = df.E[idx]u"m"
        y  = df.F[idx]u"m"
        z  = df.G[idx]u"m"
        r0 = [x, y, z]
        vx = df.dE[idx]u"m/s"
        vy = df.dF[idx]u"m/s"
        vz = df.dG[idx]u"m/s"
        v0 = [vx, vy, vz]
        rv0 = ArrayPartition(r0, v0)
        tspan = (df.time[idx]u"s", df.time[end]u"s")
        μ = 3.986004418e14u"m^3/s^2"
        # prob = ODEProblem(twoBody, u0, tspan)
        # println(prob)
        # sol = solve(prob, VCABM(), reltol = 1e-12, abstol = 1e-12)
        prob = ODEProblem((t, y, dy) -> twoBodyNew(t, y, dy, μ), rv0, tspan)
        sol = solve(prob, VCABM())
        return sol
    end
    
    function twoBodyNew(t, y, dy, μ)
        show(y)
        r = norm(y.x[1])
        dy.x[1] .= y.x[2]
        dy.x[2] .= -μ .* y.x[1] / r^3
    end
    
    
    function propWholeTraj(df)
        for ii = 1:size(df, 1)
            propagateTrajNewWay(df, ii)
        end
    end
    
    propWholeTraj(df)
    

    这里的代码基本上是从 this page 使用静态数组。我也有自己版本的双体函数和设置并调用它的函数,但这也有同样的错误。

    你有没有想过这里会出什么问题?

    1 回复  |  直到 4 年前
        1
  •  1
  •   Chris Rackauckas    6 年前

    函数签名应为 twoBodyNew(t, y, μ, dy) ((du),u,p,t) . 链接页面上的文档已修复。