MIRK4 adaptive time steps

[sebastiantk]

Hi,

I am defining a TwoPointBVProblem and use

solve(prob, MIRK4(), dt = 1e-1, adaptive =true)

to solve it. However, the step size is not adjusted and remains at 0.01.

Am I specifying this incorrectly? From what I understood #89 added this functionality.

Thank you!

[ErikQQY]

The adaptivity in MIRK solvers can automatically adapt the mesh according to the quality of numerical solution, could you please share your MWE and ]st?

[sebastiantk]

I am not sure yet where exactly the problem lies in my case, so I will take the official example.

using BoundaryValueDiffEq
using Plots
const g = 9.81
L = 1.0
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g / L) * sin(θ)
end

function bc2a!(resid_a, u_a, p) # u_a is at the beginning of the time span
    resid_a[1] = u_a[1] + pi / 2 # the solution at the beginning of the time span should be -pi/2
end
function bc2b!(resid_b, u_b, p) # u_b is at the ending of the time span
    resid_b[1] = u_b[1] - pi / 2 # the solution at the end of the time span should be pi/2
end
bvp2 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), [pi / 2, pi / 2], tspan;
    bcresid_prototype=(zeros(1), zeros(1)))
sol2 = solve(bvp2, MIRK4(), dt=0.05,adaptive=true)
sol3 = solve(bvp2, MIRK4(), dt=0.5, adaptive=true)

plot(sol2.t, reduce(hcat, sol2.u)')
plot!(sol3.t, reduce(hcat, sol3.u)')

Using the larger 0.5 timestep the solver still returns Success but the solutions are completely different. The timestep is adapted to approximately 0.1.

What I was hoping for and what I think might be a problem in my use case is that the grid is adapted uniformly instead of making the grid denser in areas where the change in the function is greater.

My st

(@v1.10) pkg> st
Status `~/.julia/environments/v1.10/Project.toml`
  [764a87c0] BoundaryValueDiffEq v5.6.3
  [13f3f980] CairoMakie v0.11.9
  [0c46a032] DifferentialEquations v7.13.0
  [31c24e10] Distributions v0.25.107
  [634d3b9d] DrWatson v2.14.1
  [a3315474] EconPDEs v1.0.3
  [26cc04aa] FiniteDifferences v0.12.31
  [587475ba] Flux v0.14.13
  [f6369f11] ForwardDiff v0.10.36
  [713c75ef] Franklin v0.10.95
  [d8418881] Intervals v1.10.0
  [c8e1da08] IterTools v1.10.0
  [98e50ef6] JuliaFormatter v1.0.53
  [b964fa9f] LaTeXStrings v1.3.1
  [b2108857] Lux v0.5.23
  [dde8697e] MakiePublication v0.3.4
  [442fdcdd] Measures v0.3.2
  [961ee093] ModelingToolkit v9.4.0
  [54ca160b] ODEInterface v0.5.0
  [429524aa] Optim v1.9.2
  [3bd65402] Optimisers v0.3.2
  [7f7a1694] Optimization v3.23.0
  [3e6eede4] OptimizationBBO v0.2.1
  [e4316d97] OptimizationMultistartOptimization v0.2.0
  [4e6fcdb7] OptimizationNLopt v0.2.0
  [36348300] OptimizationOptimJL v0.2.2
  [500b13db] OptimizationPolyalgorithms v0.2.0
  [1dea7af3] OrdinaryDiffEq v6.74.0
  [f0f68f2c] PlotlyJS v0.18.13
  [91a5bcdd] Plots v1.40.2
  [c3e4b0f8] Pluto v0.19.40
  [49802e3a] ProgressBars v1.5.1
  [1fd47b50] QuadGK v2.9.4
  [fcd29c91] QuantEcon v0.16.6
  [f2b01f46] Roots v2.1.2
  [1ed8b502] SciMLSensitivity v7.56.1
  [90137ffa] StaticArrays v1.9.3
  [f3b207a7] StatsPlots v0.15.7
  [c3572dad] Sundials v4.24.0
  [24249f21] SymPy v2.0.1
  [0c5d862f] Symbolics v5.23.0
  [770da0de] UpdateJulia v0.4.4
  [e88e6eb3] Zygote v0.6.69
  [c771fb93] ODEInterface_jll v0.0.1+0
  [37e2e46d] LinearAlgebra
  [44cfe95a] Pkg v1.10.0

[ErikQQY]

Both the solutions you mentioned are the true solutions to the example problem.

IIRC, the mesh refinement routine in MIRK solvers indeed returns uniform mesh. The MIRK solvers don't support grid coarsing for now.

[sebastiantk]

Not sure I understand - I get two different solutions:

Both are correct?

[ErikQQY]

When we are constructing a boundary value problem using BVProblem or TwoPointBVProblem, u0=[pi/2, pi/2] is just the initial guess, the side condition of this problem is only the two boundary condition specified in $u=0$ and $u=\pi/2$, which don’t guarantee a unique solution. The returned numerical solutions you mentioned both satisfy the differential equations and boundary conditions, so the solver would return a successful return code and terminate the solving.

[sebastiantk]

Ah, I have erroneously assumed that the uniqueness of an initial value problem translates to a boundary value problem - thank you for pointing this out!

Thank you for answering my question!

[ChrisRackauckas]

Yeah this is just not unique. To match the IVP you'd just need to add another constraint equation to the initial point.