System of differential equations in Julia (BVP)

General Usage
,
1

Hello, i’m trying to solve a differential equation using DifferentialEquations.jl.
In the end, i’ll try to solve the equation

$$ \ddot{\lambda} + A\lambda = f ,\quad \dot{\lambda}(0) = a, \ \dot{\lambda}(1) = b $$

where \lambda is a matrix-valued function defined on [0,1], A is a linear operator, and different values of f that I only have numerically. I don’t mind vectorizing lambda to make it a simple vector, but I can’t understand the doc to make it work. Here’s what I wrote :


using LinearAlgebra
using DifferentialEquations

# y"+ y' = 0  , y 3x3 matrices
# y'(0) = 0, y'(1)=Id


function g2!(du, u, p, t)
    du[1] = u[2,]
    du[2] = -u[1]
end


myI = Matrix(Float64.(I(3)))

function bc2!(residual, u, p, t)
    residual[1] = u(0)[2]          # the derivative at the beginning of the time span should be 0
    residual[2] = u(1)[2] - myI # the derivative at the end of the time span should be Id
end

tspan = (0.0, 1.0)
u0= zeros(3,3) # initial condition
bvp2 = BVProblem(g2!, bc2!, u0 , tspan)

sol = solve(bvp2, MIRK4(), dt = 0.5)

stacktrace :

ERROR: MethodError: no method matching -(::ForwardDiff.Dual{ForwardDiff.Tag{BoundaryValueDiffEqMIRK.var"#42#48"{…}, Float64}, Float64, 9}, ::Matrix{Float64})
For element-wise subtraction, use broadcasting with dot syntax: scalar .- array
The function `-` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  -(::ChainRulesCore.NotImplemented, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/U6wNx/src/tangent_arithmetic.jl:49
  -(::Any, ::ChainRulesCore.NotImplemented)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/U6wNx/src/tangent_arithmetic.jl:50
  -(::ChainRulesCore.ZeroTangent, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/U6wNx/src/tangent_arithmetic.jl:101
  ...

Stacktrace:
  [1] bc2!(residual::Matrix{…}, u::BoundaryValueDiffEqCore.EvalSol{…}, p::SciMLBase.NullParameters, t::Vector{…})
    @ Main ~/Documents/Codes recherche/perso_augmented_lagrangian_quantum/simpleODE.jl:125
  [2] __vec_bc!
    @ ~/.julia/packages/BoundaryValueDiffEqCore/YqNKE/src/utils.jl:287 [inlined]
  [3] #20
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:83 [inlined]
  [4] eval_bc_residual!
    @ ~/.julia/packages/BoundaryValueDiffEqCore/YqNKE/src/utils.jl:123 [inlined]
  [5] __mirk_loss_bc!(resid::Vector{…}, u::Vector{…}, p::SciMLBase.NullParameters, pt::SciMLBase.StandardBVProblem, bc!::BoundaryValueDiffEqMIRK.var"#20#33"{…}, y::Vector{…}, mesh::Vector{…}, cache::BoundaryValueDiffEqMIRK.MIRKCache{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:289
  [6] #42
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:216 [inlined]
  [7] FixTail
    @ ~/.julia/packages/DifferentiationInterface/Vap8g/src/utils/context.jl:169 [inlined]
  [8] chunk_mode_jacobian(f!::DifferentiationInterface.FixTail{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/jacobian.jl:187
  [9] jacobian
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/jacobian.jl:42 [inlined]
 [10] jacobian(f!::BoundaryValueDiffEqMIRK.var"#42#48"{…}, y::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{…}, backend::AutoForwardDiff{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/Vap8g/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:470
 [11] __construct_nlproblem(cache::BoundaryValueDiffEqMIRK.MIRKCache{…}, y::Vector{…}, loss_bc::BoundaryValueDiffEqMIRK.var"#42#48"{…}, loss_collocation::BoundaryValueDiffEqMIRK.var"#44#50"{…}, loss::BoundaryValueDiffEqMIRK.var"#46#52"{…}, ::SciMLBase.StandardBVProblem)
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:357
 [12] __construct_nlproblem
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:238 [inlined]
 [13] __perform_mirk_iteration(cache::BoundaryValueDiffEqMIRK.MIRKCache{…}, abstol::Float64, adaptive::Bool; nlsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:161
 [14] __perform_mirk_iteration
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:159 [inlined]
 [15] solve!(cache::BoundaryValueDiffEqMIRK.MIRKCache{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/9OHWB/src/mirk.jl:140
 [16] #__solve#33
    @ ~/.julia/packages/BoundaryValueDiffEqCore/YqNKE/src/BoundaryValueDiffEqCore.jl:36 [inlined]
 [17] __solve
    @ ~/.julia/packages/BoundaryValueDiffEqCore/YqNKE/src/BoundaryValueDiffEqCore.jl:33 [inlined]
 [18] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/DG18l/src/solve.jl:635 [inlined]
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/DG18l/src/solve.jl:592 [inlined]
 [20] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Matrix{…}, p::SciMLBase.NullParameters, args::MIRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/DG18l/src/solve.jl:1142
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/DG18l/src/solve.jl:1120 [inlined]
 [22] solve(prob::BVProblem{…}, args::MIRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/DG18l/src/solve.jl:1057
 [23] macro expansion
    @ ./timing.jl:581 [inlined]
 [24] top-level scope
    @ ~/Documents/Codes recherche/perso_augmented_lagrangian_quantum/simpleODE.jl:132
Some type information was truncated. Use `show(err)` to see complete types.

Is there a way to make it work ?

2

I don’t know how to edit, sorry, the equation is
\ddot{\lambda} + A\lambda = f ,\quad \dot{\lambda}(0) = a, \dot{\lambda}(1) = b

3

Your issue seems to be that you’re trying to subtract a Matrix from a scalar without using broadcasting here: residual[2] = u(1)[2] - myI. Like the error message says, you might want to try using broadcasting, like residual[2] = u(1)[2] .- myI.