BVProblem with solver RadauIIa5: a very simple MWE yields all-zero results (?)

General Usage
1

When I run the very simple MWE posted below containing BVProblem which specifies solver RadauIIa5, it runs with a return code “Unstable”, and produces unexpected all-zero results.
In view of the simplicity of this problem, I wonder what I might have done wrong.
I would appreciate it if you could identify any mistake(s).

I would also appreciate the identification of any benchmark run in the Julia language with this solver, together with its numerical results. All the benchmarks in the online documentation seem to have used other languages.

Please also consider the question appearing as one of the comments in this MWE.
thanks.


#  This example using BVProblem with solver RadauIIa5 is 
# an attempt to reproduce the analytic solution... 
#   y(t) = 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

[dy_dt]     # These square brackets are not required, and don't seem to matter.

end

function initialGuess(t)

And here, t again 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

[t .* t]      # If these square brackets are absent ==> error msg!  Why the difference?

end

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

Here again t seems to be a vector.

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

end

tspan = (1.0, 10.0)

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

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

Result below…

julia> include(“…/…/…/BVProblem,Example3.jl”)
retcode: Unstable
Interpolation: FIRK Order 9 Interpolation
t: 2881-element Vector{Float64}:
1.0
1.003125
1.00625
1.009375
1.0125
1.015625
1.0187499999999998
1.0218749999999999
1.025
1.028125
1.03125
1.034375
1.0375
1.0406250000000001
1.0437500000000002
1.046875
1.05
1.053125
1.05625
1.059375
1.0625
1.065625
1.06875
1.0718750000000001
1.0750000000000002
⋮
9.928125000000001
9.93125
9.934375
9.9375
9.940625
9.94375
9.946874999999999
9.95
9.953125
9.956249999999999
9.959374999999998
9.962499999999999
9.965625
9.96875
9.971875
9.975
9.978124999999999
9.98125
9.984375
9.9875
9.990625000000001
9.99375
9.996875
10.0
u: 2881-element Vector{Vector{Float64}}:
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
⋮
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]
[0.0]

julia>

I get the same result with this modification:

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

Here again t seems to be a vector.

resid[1] = abs(sol[1,1] - 11.0)

end

And this gives a syntax error (?) ...

resid[2] = abs(sol1 - 110.0) # But this gives a syntax error (?)

2

I don’t know how the comments got exploded. It looks like I quoted the pasted text correctly.

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

dy_dt[1] = 2.0 .* t

end
4

Thank you! That worked.