When I run the simple MWE posted below containing BVProblem which specifies solver RadauIIa5 (notice the lower-case “a5”!) , it throws an error complaining that the two arrays provided by the solver to my boundary-condition function have incompatible lengths.
What specifically did I do wrong?
# This example using BVProblem with solver RadauIIA5 is based on
# an attmpt to reproduce the analytic solution...
# y(t) = exp(-t^2) + 10.
using BoundaryValueDiffEqFIRK
function f!(dy_dt, y, par, t)
# It seems that t is provided as a vector.
# (A vector might have been expected only for an "orthogonal collocation method".)
dy_dt = .-2.0*t*(y[1,:] .- 10.)
dy_dt
end
function initialGuess(t)
# And here, t seems to be a vector
# The arguments of this function appear to be undocumented.
# We use the correct solution but without the constant as the initial guess
[exp.(-t.*t)]
end
function bc!(resid, y, par, t) # Boundary conditions
# Here again t seems to be a vector, which seems inevitable,
# as the documentation says we can constrain the entire tspan.
# But here we would like to insist on the exact solution including the constant,
# but we get an error...
resid[1] = maximum(abs.(exp.(-(t.*t) ) .+ 10. .- y[1,:]) )
# Result: Error: Array t has size different from size of y[1,:]
# Using this version below, it runs but gives all- zero results
# resid[1] = abs(y[1,1] - 11.0) #
end
tspan = (0.0, 10.0)
bvp = BVProblem(f!, bc!, initialGuess , tspan, 0.0)
sol = solve(bvp, RadauIIa5(); dt = 0.1, abstol=1e-2, reltol=1e-2)
julia> include("../../../BVProblem,Example2.jl")
ERROR: LoadError: DimensionMismatch: arrays could not be broadcast to a common size: a has axes Base.OneTo(101) and b has axes Base.OneTo(601)
Stacktrace:
[1] _bcs1
@ ./broadcast.jl:523 [inlined]
[2] _bcs
@ ./broadcast.jl:517 [inlined]
[3] broadcast_shape
@ ./broadcast.jl:511 [inlined]
[4] combine_axes
@ ./broadcast.jl:492 [inlined]
[5] _axes
@ ./broadcast.jl:236 [inlined]
[6] axes
@ ./broadcast.jl:234 [inlined]
[7] combine_axes
@ ./broadcast.jl:493 [inlined]
[8] instantiate
@ ./broadcast.jl:302 [inlined]
[9] materialize
@ ./broadcast.jl:867 [inlined]
[10] bc!(resid::Vector{…}, y::BoundaryValueDiffEqCore.EvalSol{…}, par::Float64, t::Vector{…})
@ Main ~/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.jl:36
[11] eval_bc_residual!
@ ~/.julia/packages/BoundaryValueDiffEqCore/vnuZ6/src/utils.jl:123 [inlined]
[12] __firk_loss!(resid::Vector{…}, u::Vector{…}, p::Float64, y::Vector{…}, pt::SciMLBase.StandardBVProblem, bc!::typeof(bc!), residual::Vector{…}, mesh::Vector{…}, cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, eval_sol::BoundaryValueDiffEqCore.EvalSol{…})
@ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/ddek7/src/firk.jl:670
[13] #119
@ ~/.julia/packages/BoundaryValueDiffEqFIRK/ddek7/src/firk.jl:426 [inlined]
[14] NonlinearFunction
@ ~/.julia/packages/SciMLBase/tWwhl/src/scimlfunctions.jl:2469 [inlined]
[15] evaluate_f
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/utils.jl:174 [inlined]
[16] __init(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Float64, reltol::Float64, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{…})
@ NonlinearSolveFirstOrder ~/.julia/packages/NonlinearSolveFirstOrder/3kzAL/src/solve.jl:151
[17] __init
@ ~/.julia/packages/NonlinearSolveFirstOrder/3kzAL/src/solve.jl:121 [inlined]
[18] #__solve#146
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:5 [inlined]
[19] __solve
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:1 [inlined]
[20] macro expansion
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:173 [inlined]
[21] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, verbose::Bool, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{…})
@ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:130
[22] __generated_polysolve
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:130 [inlined]
[23] #__solve#154
@ ~/.julia/packages/NonlinearSolveBase/Kek5u/src/solve.jl:127 [inlined]
[24] __perform_firk_iteration(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…}, abstol::Float64, adaptive::Bool; nlsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{…})
@ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/ddek7/src/firk.jl:358
[25] __perform_firk_iteration
@ ~/.julia/packages/BoundaryValueDiffEqFIRK/ddek7/src/firk.jl:354 [inlined]
[26] solve!(cache::BoundaryValueDiffEqFIRK.FIRKCacheExpand{…})
@ BoundaryValueDiffEqFIRK ~/.julia/packages/BoundaryValueDiffEqFIRK/ddek7/src/firk.jl:309
[27] __solve(::BVProblem{…}, ::RadauIIa5{…}; kwargs::@Kwargs{…})
@ BoundaryValueDiffEqCore ~/.julia/packages/BoundaryValueDiffEqCore/vnuZ6/src/BoundaryValueDiffEqCore.jl:37
[28] __solve
@ ~/.julia/packages/BoundaryValueDiffEqCore/vnuZ6/src/BoundaryValueDiffEqCore.jl:34 [inlined]
[29] #solve_call#44
@ ~/.julia/packages/DiffEqBase/R2Vjs/src/solve.jl:634 [inlined]
[30] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Function, p::Float64, args::RadauIIa5{…}; kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/R2Vjs/src/solve.jl:1122
[31] solve_up
@ ~/.julia/packages/DiffEqBase/R2Vjs/src/solve.jl:1101 [inlined]
[32] solve(prob::BVProblem{…}, args::RadauIIa5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/R2Vjs/src/solve.jl:1038
[33] top-level scope
@ ~/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.jl:49
[34] include(fname::String)
@ Main ./sysimg.jl:38
[35] top-level scope
@ REPL[16]:1
in expression starting at /Users/andy 1/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.jl:49
Some type information was truncated. Use `show(err)` to see complete types.
julia>