ODE solver for BVProblem never calls the provided function for boundary conditions

General Usage
1

My experience with BVProblem is becoming kaftaesque. I hope I’m just missing something.

My previous discourse 125374 on another aspect of the problem did not produce a solution.
So I investigated further, simplifying the example even more… Now it seems as if this
BVProblem code may never have worked. I note that many of its benchmarks on the web
have 404 links to the source code, and the “reference” links mention only 3 other languages,
not including Julia. So might BVProblem be just a place-holder during work-in progress?

I suspect my problem really needs the boundary conditions specified over various points
distributed over the tspan, and BVProblem with a FIRK solver seems like the only choice.
The only FIRK solver for Julia that I can find is RadauIIA5, which has performed well for others, but who knows if the version that was ported to Julia has been checked out?
Many other FIRK solvers are said to be included in the NSDERungeKutta package, which I installed only to see this:

Failed to precompile NSDERungeKutta [3978d399-9848-4322-9a48-ca26b732c7dc] to "/Users/andy 1/.julia/compiled/v1.11/NSDERungeKutta/jl_kXY7PR".
ERROR: LoadError: UndefVarError: `AbstractInitialValueSolver` not defined in `NSDERungeKutta`
~/Applications/Julia/NSDERungeKutta.jl-master/src/abstract.jl:1

but this AbstractInitialValueSolver appears to be totally unknown to google.

Any suggestions? <<<<<<<<

So here’s my simplified yet diversified example from Kafka:

using BoundaryValueDiffEq, StaticArrays, LinearAlgebra, ArgCheck   

function f!(dy_dt, y, par, t)     ! Calculate derivatives 
#  It seems that t is not a vector.
#  (A vector might be expected only for an "orthogonal collocation method".)
    y1, y2 = y
    par1, par2 = par 

# Sign ambiguity of sqrt() should force use of boundary conditions, but nope...  
    dy_dt[1] = 2*y2 + sqrt(t*par1) 
    dy_dt[2] = 3*y1 - sqrt(t/par2)

    dy_dt 
end

function initialGuess(t)
#  And here, t might be  expected to be a vector (if it could be called without complaints). 
# It seems necessary here to have access to parameters in par, as in bc! function. 
# But the arguments of this function appear to be undocumented.  
# A clue about work-in-progress?

    guess1 = t.*2 
    guess2 = t.*3 

    (guess1, guess2)
end

function bc!(resid, y, par, t)   #  Boundary conditions 
# And here, one might expect t to be a vector?  Why not?
# The documentation says we can constrain the entire tspan.
# But bc never gets called! At least its arguments are documented.  
# Someone wrote y[1,i] & y[2,i] can be expected to have been evaluated at t[i]

    @argcheck 0>1  #  Do we ever hit this tripwire?  Nope.  No matter what I do. 
    println("Argument t: type = $(typeof(t) ); length = $(length(t)) ; value = $(t)")

    par1, par2 = par 

    resid[1] = abs(y[1,1] - par1)
    resid[2] = abs(y[2,end] - par2) 
end

tspan = (0.0, 50.0)

par = (1.0, 2.0)

# bvp = BVProblem(f!, bc!, initialGuess, tspan, par)
#  Fails: can't "iterate" the function initialGuess.  See discourse 125374. 

# For this next work-around attempt with fixed numerical initGuess, I had to assume that
# the solver could accept the entire span to be covered by bc, but in my own uniform steps. 

t1, t2 = tspan
initGuess  = initialGuess([t for t = t1:0.1:t2])  
# But nope...
# bvp = BVProblem(f!, bc!, initGuess, tspan, par)
# Fails -- can't zero(wrongType) for bc=[2,:] matrix, 
# just like discourse 40395  -- with very nice discussion there, but I don't see 
# how static arrays can work if dynamic grid refinement may be needed... 

# ... or fails on bad state type as a 2-tuple of vectors
# ERROR: LoadError: Tuple type used as a state. 
# Since a tuple does not have vector properties, 
# it will not work as a state type in equation solvers.

# But if the problem is stripped down to useless...
bvp = BVProblem(f!, bc!, [1.0, 2.0], tspan, par)
#  Reports "success", producing some nonzero result, but UNinfluenced by bc 

# bvp = BVProblem(f!, bc!, nothing, tspan, par)
#  Reports "success", producing "nothing"s, only at ends of tspan, incompatible with bc   

sol = solve(bvp, RadauIIA5(); abstol=1e-2, reltol=1e-2)

The imbedded comments indicate the options tried (and failed or succeeded insignificantly).

So next I will try other FIRK solvers if they emerge.
Or if anyone know of other problem types allowing distributed boundary conditions,
please let me know.

thanks

Andy

2

BoundaryValueDiffEq.jl does support multi-points BVP, and if you remove the @argcheck inyour MWE, the code does work with MIRK or FIRK methods, I believe you should use the RadauIIa5 methods from BoundaryValueDiffEqFIRK:

using BoundaryValueDiffEq, StaticArrays, LinearAlgebra, ArgCheck   

function f!(dy_dt, y, par, t)
#  It seems that t is not a vector.
#  (A vector might be expected only for an "orthogonal collocation method".)
    y1, y2 = y
    par1, par2 = par 

# Sign ambiguity of sqrt() should force use of boundary conditions, but nope...  
    dy_dt[1] = 2*y2 + sqrt(t*par1) 
    dy_dt[2] = 3*y1 - sqrt(t/par2)
end

function initialGuess(t)

    guess1 = t.*2 
    guess2 = t.*3 

    [guess1, guess2]
end

function bc!(resid, y, par, t)   #  Boundary conditions 

    #@argcheck 0>1  #  Do we ever hit this tripwire?  Nope.  No matter what I do. 
    println("Argument t: type = $(typeof(t) ); length = $(length(t)) ; value = $(t)")

    par1, par2 = par 

    resid[1] = abs(y[1,1] - par1)
    resid[2] = abs(y[2,end] - par2) 
end

tspan = (0.0, 50.0)

par = (1.0, 2.0)
t1, t2 = tspan
initGuess  = initialGuess([t for t = t1:0.1:t2])  

# But if the problem is stripped down to useless...
bvp = BVProblem(f!, bc!, [1.0, 2.0], tspan, par)
#  Reports "success", producing some nonzero result, but UNinfluenced by bc 

# bvp = BVProblem(f!, bc!, nothing, tspan, par)
#  Reports "success", producing "nothing"s, only at ends of tspan, incompatible with bc   

sol = solve(bvp, RadauIIa5(); dt = 0.01, abstol=1e-2, reltol=1e-2)
3

Thanks for your prompt reply.

I tried your recommendeRadauIIa5from BoundaryValueDiffEqFIRK,
(which required a lot of compiling struggles), while leaving only
BoundaryValueDiffEqFIRK in my “using” statement.

But of course, I did NOT first remove the @argcheck tripwire
from bc!, because that would have removed the main source
of information about the main problem — failure to call bc!.

And this different RadauIIa5 that you suggested produced
exactly the same results, while never calling bc!.

So please tell me — have you ever observed reportable
evidence that BVProblem with RadauIIa5 actually produced
correct results, based on calling the argument function to
fetch boundary conditions?

If so, could you please tell me where I might see this evidence?

(I’m still hoping for information about finding this mysterious
AbstractInitialValueSolver
in order to try other solvers.)

thanks again

Andy

4

But then I realized that I should have first used pkg to remove the
package BoundaryValueDiffEq.jl that was previously supplying RadauIIa5,
to ensure that I am now using the one in BoundaryValueDiffEqFIRK.
So I did this, but the result did not change.
So now we have learned that RadauIIa5 is present in both packages.

But now I wonder — besides the RadauIIa5 in OrdinaryDiffEqFIRK,
(which announced itself to be incompatible with BVProblem), maybe
the one in BoundaryValueDiffEq.jl and the one in BoundaryValueDiffEqFIRK,
both being compatible with BVProblem, should also be different
versions, but by package configuration error, are the same version with
neither supporting FIRK by calling the BC! function.

I’m speculating that the version in BoundaryValueDiffEq.jl doesn’t
call a function, but just expects constants for boundary conditions.
And BoundaryValueDiffEqFIRK was accidentally loaded with that
same version. Or the FIRK version wasn’t available yet…

What do you think?

It seems a pity that we have multiple modules with the same name.

5

I see you might have some misunderstandings here.

  1. BoundaryValueDiffEqFIRK.jl is actually a subpackage of BoundaryValueDiffEq.jl, so basically BoundaryValueDiffEq.jl is the superset of BoundaryValueDiffEqFIRK, BoundaryValueDiffEqMIRK and BoundaryValueDiffEqFIRK etc, if you using BoundaryValueDiffEq, it will automatically load BoundaryValueDiffEqFIRK package.
  2. The boundary conditions are coupled with the collocation equations in the residual control during the nonlinear system solving, so your statement(RadauIIa5 not calling bc!) is incorrect.
  3. AbstractInitialValueSolver in NSDERungeKutta.jl has nothing to do with BoundaryValueDiffEq.jl, so I could not provide further assistance here
6

I am wondering why are you using @argcheck 0>1 in the bc! definition? that is very weird and causing errors

7

Thanks for asking about argcheck.
I didn’t realize this was unclear.

Yes, I’m taking a weird measure to clarify the weird situation that BVProblem solver can present its result to me in REPL session at the same time that argcheck is sitting in the BC function waiting to kill the calculation whenever BC is called.
I did this simply to prove that BC is never called before a result appears from the solver.
Don’t you think that is what is really weird?

You wrote:
` The boundary conditions are coupled with the collocation equations in the residual control during the nonlinear system solving, so your statement(RadauIIa5 not calling bc!) is incorrect."

I guess you wrote this because I failed to post the latest REPL session result together with the version of my example code. I will do so now below, but first let me try to better understand what you wrote. Is it correct that the “coupling with the collocation equations in the residual control during the nonlinear system solving”, if any, has already occurred when the solution result appears in my REPL session? Do you perhaps believe that the solution that appeared did not require “nonlinear solving”? Is the misunderstanding perhaps because my example problem is too simple to cause “nonlinear solving”? Anyway, here is the code containing the “weird argcheck” dog that never barked…

# using  BoundaryValueDiffEqFIRK, LinearAlgebra, ArgCheck 

function f!(dy_dt, y, par, t)
#  It seems that t is not a vector.
#  (A vector might be expected only for an "orthogonal collocation method".)
    y1, y2 = y
    par1, par2 = par 

# Sign ambiguity of sqrt() should force use of boundary conditions, but nope...  
    dy_dt[1] = 2*y2 + sqrt(t*par1) 
    dy_dt[2] = 3*y1 - sqrt(t/par2)

    dy_dt 
end

function initialGuess(t)
#  And here, t might be  expected to be a vector (if it could be called without complaints). 
# It seems necessary here to have access to parameters in par, as in bc! function
# But the arguments of this function appear to be undocumented. 

    guess1 = t.*2 
    guess2 = t.*3 

    (guess1, guess2)
end

function bc!(resid, y, par, t)   #  Boundary conditions 
# And here, one might expect t to be a vector?  Why not?
# The documentation says we can constrain the entire tspan.
# But bc never gets called! At least its arguments are documented.  
# Someone wrote y[1,i] & y[2,i] should have been evaluated at t[i]

    @argcheck 0>1  #  Is this function ever called?  Nope. No matter what I do. 
    println("Argument t: type = $(typeof(t) ); length = $(length(t)) ; value = $(t)")

    par1, par2 = par 

    resid[1] = abs(y[1,1] - par1)
    resid[2] = abs(y[2,end] - par2) 
end

tspan = (0.0, 50.0)

par = (1.0, 2.0)

# bvp = BVProblem(f!, bc!, initialGuess, tspan, par)
#  Fails: can't "iterate" the function initialGuess -- see error in discourse 125374.

# For this next work-around attempt with fixed numerical initGuess, I have to assume that
# the solver could accept the entire span to be covered by bc in my uniform steps.

t1, t2 = tspan
initGuess  = initialGuess([t for t = t1:0.1:t2])  
# But nope...
# bvp = BVProblem(f!, bc!, initGuess, tspan, par)
# Fails -- can't zero(wrongType) for bc=[2,:] matrix, or on bad state type as a 2-tuple of vectors

bvp = BVProblem(f!, bc!, [1.0, 2.0], tspan, par)
#  Reports "success", produces some nonzero result, uninfluenced by bc 

# bvp = BVProblem(f!, bc!, nothing, tspan, par)
#  Reports "success", producing "nothing"s, only at ends of tspan, incompatible with bc   

sol = solve(bvp, RadauIIA5(); abstol=1e-2, reltol=1e-2)

and here is the solution presented on REPL, obviously without the argcheck being encountered, because the bc function wasn’t ever called while calculationg this result…

julia> include("../../../BVProblem,Example1.jl")
retcode: Success
Interpolation: 3rd order Hermite
t: 132-element Vector{Float64}:
  0.0
  0.10479183976232473
  0.25100294534130374
  0.6393469920463636
  1.0276910387514233
  1.4160350854564832
  1.804379132161543
  2.192723178866603
  2.5810672255716627
  2.9694112722767225
  3.3577553189817824
  3.746099365686842
  4.134443412391902
  4.522787459096962
  4.911131505802022
  5.299475552507081
  5.6878195992121405
  6.0761636459172
  6.464507692622259
  6.852851739327319
  7.241195786032378
  7.6295398327374375
  8.017883879442497
  8.406227926147556
  8.794571972852616
  ⋮
 41.415471896077634
 41.80381594278269
 42.19215998948775
 42.580504036192806
 42.968848082897864
 43.35719212960292
 43.74553617630798
 44.13388022301304
 44.522224269718095
 44.91056831642315
 45.29891236312821
 45.68725640983327
 46.075600456538325
 46.46394450324338
 46.85228854994844
 47.2406325966535
 47.628976643358556
 48.01732069006361
 48.40566473676867
 48.79400878347373
 49.182352830178786
 49.570696876883844
 49.9590409235889
 50.0
u: 132-element Vector{Vector{Float64}}:
 [1.0, 2.0]
 [1.4784761627787666, 2.370778510242625]
 [2.338949976872162, 3.1554108530522784]
 [6.540514649735073, 7.791398085818848]
 [17.043932491135138, 20.34090778009886]
 [44.02166684484558, 53.168099600509024]
 [113.73656711923445, 138.39003455203607]
 [294.14995378894423, 359.21862154270167]
 [761.2118925861841, 931.1362766113605]
 [1970.4929487492575, 2412.0947509020557]
 [5101.577520805165, 6246.780686905033]
 [13208.709723208689, 16175.861925882084]
 [34200.1096124359, 41884.88914844597]
 [88552.18661479125, 108452.23890176378]
 [229283.61665110654, 280812.26196291856]
 [593673.3513300174, 727096.6497726047]
 [1.5371718101042858e6, 1.8826414802019214e6]
 [3.980131306890526e6, 4.874643529280528e6]
 [1.0305580065949526e7, 1.2621704392795416e7]
 [2.668378944104268e7, 3.2680832266144663e7]
 [6.909117473102453e7, 8.46190598496398e7]
 [1.7889477332675815e8, 2.19100454032248e8]
 [4.6320445627814394e8, 5.673072800576137e8]
 [1.1993551549686558e9, 1.468904072795267e9]
 [3.1054381471821814e9, 3.803369441901954e9]
 ⋮
 [1.5806270187719904e44, 1.935864834823972e44]
 [4.0926488053551184e44, 5.012450634765595e44]
 [1.0596917580839402e45, 1.297852045969242e45]
 [2.743813787984171e45, 3.3604718648871426e45]
 [7.104437725121664e45, 8.701123667963687e45]
 [1.8395211661653428e46, 2.252944114077279e46]
 [4.762992163060103e46, 5.833450224186192e46]
 [1.2332608486730967e47, 1.5104298995004137e47]
 [3.1932286865085376e47, 3.910890456981831e47]
 [8.268088178840666e47, 1.0126298593248529e48]
 [2.1408201179549935e48, 2.621958460037315e48]
 [5.543132436794436e48, 6.788923023408346e48]
 [1.4352591772723885e49, 1.757826316482069e49]
 [3.716254174753007e49, 4.551463241316328e49]
 [9.622335331528868e49, 1.1784905848100067e50]
 [2.4914694441895247e50, 3.0514144239999736e50]
 [6.451053489053346e50, 7.900894675790903e50]
 [1.6703432271939853e51, 2.0457443009695093e51]
 [4.324947082468293e51, 5.296956758293058e51]
 [1.1198397408162554e52, 1.371517979345195e52]
 [2.899552344108043e52, 3.5512118627777856e52]
 [7.50768479612472e52, 9.194998450078359e52]
 [1.943932176720242e53, 2.3808209637712046e53]
 [2.149083892016969e53, 2.6320794749380586e53]

julia>

I hope this clarifies!

8

And if you would prefer an even simpler non-linear problem, I offer the one below, which unfortunately will not run because of an “iterator” error. I can’t solve this error, perhaps because I can’t find any documentation for the arguments of the initial guess function.
Perhaps you might help with this?


#  This example using BVProblem with solver RadauIIA5 is based on
# an attmpt to reproduce the exact analytic solution
#   y(t) = exp(-t^2) + 10.
# by solving the associated ODE.

using BoundaryValueDiffEqFIRK, ArgCheck 

function f!(dy_dt, y, par, t)

#  It seems that t is not provided as a vector.
#  (A vector might be expected only for an "orthogonal collocation method".)

    dy_dt = -2.0*t*(y - 10.)  

    dy_dt 
end

function initialGuess(t)
#  And here, t might be  expected 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^2) 

end

function bc!(resid, y, par, t)   #  Boundary conditions 
# Here, one might expect t to be a vector?  Why not?
# The documentation says we can constrain the entire tspan.

#     @argcheck 0>1  #  Option for a test: Is this function ever called?  

    resid[1] = abs(y[begin] - 11.0)

end

tspan = (0.0, 50.0)


bvp = BVProblem(f!, bc!, initialGuess, tspan, 0.0)

sol = solve(bvp, RadauIIA5(); abstol=1e-2, reltol=1e-2) 

ERROR: LoadError: MethodError: no method matching iterate(::typeof(initialGuess))
The function `iterate` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  iterate(::Expronicon.JLMatch)
   @ Expronicon ~/.julia/packages/Expronicon/Oa1rJ/src/match.jl:139
  iterate(::Expronicon.JLMatch, ::Any)
   @ Expronicon ~/.julia/packages/Expronicon/Oa1rJ/src/match.jl:139
  iterate(::DataStructures.EnumerateAll)
   @ DataStructures ~/.julia/packages/DataStructures/95DJa/src/multi_dict.jl:98
  ...

Stacktrace:
  [1] isempty(itr::Function)
    @ Base ./essentials.jl:1121
  [2] norm(itr::Function, p::Int64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.11.2+0.x64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:602
  [3] ODE_DEFAULT_NORM(u::Function, t::Float64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/slKc5/src/common_defaults.jl:83
  [4] __init(prob::BVProblem{…}, alg::RadauIIA5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Float64, reltol::Float64, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{…}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/hCoet/src/solve.jl:326
  [5] __solve(::BVProblem{…}, ::RadauIIA5{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/hCoet/src/solve.jl:6
  [6] solve_call(_prob::BVProblem{…}, args::RadauIIA5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/slKc5/src/solve.jl:612
  [7] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Function, p::Float64, args::RadauIIA5{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/slKc5/src/solve.jl:1080
  [8] solve(prob::BVProblem{…}, args::RadauIIA5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/slKc5/src/solve.jl:1003
  [9] top-level scope
    @ ~/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.v0.jl:43
 [10] include(fname::String)
    @ Main ./sysimg.jl:38
 [11] top-level scope
    @ REPL[443]:1
in expression starting at /Users/andy 1/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.v0.jl:43
Some type information was truncated. Use `show(err)` to see complete types.

julia>
9

OK, I see what’s going wrong, from your error stacktrace, the Radau solver you are trying to use is the RadauIIA5 solver from OrdinaryDiffEq.jl, RadauIIa5 in BoundaryValueDiffEq.jl should be the lowercase RadauIIa5.

10

Now what? Are you planning to correct this “misconfiguration”?
Or have you already done so?

11

In the MWE you provided in this post, actually it should error since you are doing the solve incorrectly, can you make sure that this example have the successful return code?

12

I provided two MWEs. Which one do you mean?

And what mistake do you see in how I’m doing the solve?

13

Oh, sorry I am not being specific, can you make sure this MWE actually work?

using  BoundaryValueDiffEqFIRK, LinearAlgebra, ArgCheck 

function f!(dy_dt, y, par, t)
#  It seems that t is not a vector.
#  (A vector might be expected only for an "orthogonal collocation method".)
    y1, y2 = y
    par1, par2 = par 

# Sign ambiguity of sqrt() should force use of boundary conditions, but nope...  
    dy_dt[1] = 2*y2 + sqrt(t*par1) 
    dy_dt[2] = 3*y1 - sqrt(t/par2)

    dy_dt 
end

function initialGuess(t)
#  And here, t might be  expected to be a vector (if it could be called without complaints). 
# It seems necessary here to have access to parameters in par, as in bc! function
# But the arguments of this function appear to be undocumented. 

    guess1 = t.*2 
    guess2 = t.*3 

    (guess1, guess2)
end

function bc!(resid, y, par, t)   #  Boundary conditions 
# And here, one might expect t to be a vector?  Why not?
# The documentation says we can constrain the entire tspan.
# But bc never gets called! At least its arguments are documented.  
# Someone wrote y[1,i] & y[2,i] should have been evaluated at t[i]

    @argcheck 0>1  #  Is this function ever called?  Nope. No matter what I do. 
    println("Argument t: type = $(typeof(t) ); length = $(length(t)) ; value = $(t)")

    par1, par2 = par 

    resid[1] = abs(y[1,1] - par1)
    resid[2] = abs(y[2,end] - par2) 
end

tspan = (0.0, 50.0)

par = (1.0, 2.0)

# bvp = BVProblem(f!, bc!, initialGuess, tspan, par)
#  Fails: can't "iterate" the function initialGuess -- see error in discourse 125374.

# For this next work-around attempt with fixed numerical initGuess, I have to assume that
# the solver could accept the entire span to be covered by bc in my uniform steps.

t1, t2 = tspan
initGuess  = initialGuess([t for t = t1:0.1:t2])  
# But nope...
# bvp = BVProblem(f!, bc!, initGuess, tspan, par)
# Fails -- can't zero(wrongType) for bc=[2,:] matrix, or on bad state type as a 2-tuple of vectors

bvp = BVProblem(f!, bc!, [1.0, 2.0], tspan, par)
#  Reports "success", produces some nonzero result, uninfluenced by bc 

# bvp = BVProblem(f!, bc!, nothing, tspan, par)
#  Reports "success", producing "nothing"s, only at ends of tspan, incompatible with bc   

sol = solve(bvp, RadauIIA5(); abstol=1e-2, reltol=1e-2)
14

How can I make sure if BoundaryValueDiffEqFIRK contains the wrong solver?

And why do you think that the reason for the failure is my MWE?

Why can’t it be just another problem with the BoundaryValueDiffEqFIRK package?

Have you corrected the problem about the wrong RadauIIA5 solver?

15

How can I make sure if BoundaryValueDiffEqFIRK contains the wrong solver?

The correct FIRK solver in BoundaryValueDiffEqFIRK is RaduaIIa5 instead of RadauIIA5.

And why do you think that the reason for the failure is my MWE?

I can’t reproduce the same result as you got, and I am pretty sure you have some mistakes in your MWE.

16

So now I’m guessing that you mean that I mis-spelled “RadauIIa5”.
If so, why don’t you tell me how to spell it correctly?
I have never seen any other spelling in the documentation.

17

If so, why don’t you tell me how to spell it correctly?

In ODE solver for BVProblem never calls the provided function for boundary conditions - #9 by ErikQQY :

the Radau solver you are trying to use is the RadauIIA5 solver from OrdinaryDiffEq.jl, RadauIIa5 in BoundaryValueDiffEq.jl should be the lowercase RadauIIa5 .

I have never seen any other spelling in the documentation.

Documentations for FIRK methods on solving BVP is

and

18

You still haven’t identified which MWE you’re talking about, or what “result” you got.
Of course, the result will depend on which packages you have installed.
I’m using only BoundaryValueDiffEqFIRK. Have you removed all other relevant packages?

And since you’re the expert, what mistake do you see in my very simple MWE?
How can you accuse it of being mistaken, unless you see the mistake?

19

You still haven’t identified which MWE you’re talking about, or what “result” you got.

The MWE in your first post: ODE solver for BVProblem never calls the provided function for boundary conditions

I’m using only BoundaryValueDiffEqFIRK. Have you removed all other relevant packages?

Latest version of BoundaryValueDiffEqFIRK.jl should working fine with BVP no matter if initial guesses are provided.

And since you’re the expert, what mistake do you see in my very simple MWE?
How can you accuse it of being mistaken, unless you see the mistake?

The mistakes are obvious, and if you execute the MWE you have provided(without the commented part), errors are expected:

  1. FIRK methods don’t have automatic step size choosing, so you must provide a dt in your solve
  2. To solve a BVP, you should use RadauIIa5 instead of RadauIIA5
  3. As for the commented part where you are trying to provide an initial guess, they should be provided as a function returning AbstractArray instead of a Tuple
20

When I do that, while also saying “using BoundaryValueDiffEqFIRK, ArgCheck”,
I get:

```ERROR: LoadError: UndefVarError: RadauIIa5 not defined in Main
Suggestion: check for spelling errors or missing imports.
Stacktrace:
[1] top-level scope
@ ~/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.jl:47
[2] include(fname::String)
@ Main ./sysimg.jl:38
[3] top-level scope
@ REPL[446]:1
in expression starting at /Users/andy 1/Applications/Julia/RealisticSeismology-main/BVProblem,Example2.jl:47

So which package contains "RadauIIa5"?