# Boundary Value Problems in ModelingToolkit
**URL:** https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491
**Category:** Modelling & Simulations
**Tags:** question, diffeq, modelingtoolkit
**Created:** [October 28, 2025, 6:44pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491 "2025-10-28T18:44:15Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@bradyelster](https://discourse.julialang.org/u/bradyelster)
#### Post date: [October 28, 2025, 6:44pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/1 "2025-10-28T18:44:15Z")
Is it currently possible to symbolically setup and then numerically solve boundary-value problems using ModelingToolkit? This [issue](https://github.com/SciML/ModelingToolkit.jl/issues/924) tells me this is possible through the dynamic optimization interface but I’m not sure how to set it up this way. As a test case I’m trying this simple example from the BoundaryValueDiffEq docs:
```julia-auto
using BoundaryValueDiffEq
function f!(du, u, p, t)
du[1] = u[2]
du[2] = u[1]
end
function bc!(res, u, p, t)
res[1] = u(0.0)[1] - 1
res[2] = u(1.0)[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
prob = BVProblem(f!, bc!, u0, tspan)
sol = solve(prob, MIRK4(), dt = 0.01)
```
Here is my attempt at setting this up in ModelingToolkit:
```julia-auto
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using BoundaryValueDiffEq # for BVProblem
@variables y(..)
eqs = [D(D(y(t))) ~ y(t)]
(ti, tf) = (0.0, 1.0)
cons = [y(ti) ~ 1.0, y(tf) ~ 0.0]
@named bvpsys = System(eqs, t; constraints=cons)
bvpsys = mtkcompile(bvpsys)
u0 = [0.0, 0.0]
prob = BVProblem(bvpsys, u0, (ti, tf))
# or do I need JuMPDynamicOptProblem() ?
sol = solve(prob, MIRK4(), dt=0.01)
```
My attempt above leads to a seg fault and my REPL crashes so something is horribly wrong. I’m guessing `BVProblem` is not the one to use here and I need a dynamic optimization solver like in [these](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/dynamic_optimization/) examples, but I’m not sure how. What is the proper syntax for this?
---
### Author:  [@Sevi](https://discourse.julialang.org/u/Sevi)
#### Post date: [November 3, 2025, 9:42am UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/2 "2025-11-03T09:42:08Z")
Welcome to the forum 👋
Did you end up trying it with `JuMPDynamicOptProblem` ?
From the docs page you linked it sounds like that is the recommended way to go and also shows how to use the problem. It should look roughly like this I guess.
```julia
using InfiniteOpt, Ipopt, DiffEqDevTools
jprob = JuMPDynamicOptProblem(bvpsys, [u0map; pmap], (ti, tf); dt = 0.001)
jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()));
```
where `[u0map; pmap]` have to be adjusted accordingly.
---
### Author:  [@bradyelster](https://discourse.julialang.org/u/bradyelster)
#### Post date: [November 3, 2025, 7:08pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/3 "2025-11-03T19:08:01Z")
Thanks for the recommendation. I tried it this way by supplying a guess for the value of the derivative. However, the solution seems to rapidly course-correct after the initial timestep in order to satisfy the boundary conditions. This does not match the correct solution found with BoundaryValueDiffEq.jl. Not sure what I’m doing wrong.
```julia-auto
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables y(..)
eqs = [D(D(y(t))) ~ y(t)]
@parameters deriv_guess
(ti, tf) = (0.0, 1.0)
cons = [y(tf) ~ 0.0]
@named bvpsys = System(eqs, t; constraints=cons)
bvpsys = mtkcompile(bvpsys)
u0map = [y(t) => 1.0, D(y(t)) => deriv_guess]
parammap = [deriv_guess => 1.0]
using InfiniteOpt, Ipopt, DiffEqDevTools
jprob = JuMPDynamicOptProblem(bvpsys, [u0map; parammap], (ti, tf); dt=0.001)
jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()));
using Plots
Plots.plot(jsol.sol)
```
I see [here](https://docs.sciml.ai/ModelingToolkit/stable/API/problems/#SciMLBase.BVProblem) that the docs have `BVProblem` defined for the `System` type, so maybe there is a way to set up the problem without using dynamic optimization? If there is, I’d happily contribute documentation. It would be great to leverage the abilities of ModelingToolkit for these types of problems.
---
### Author:  [@Sevi](https://discourse.julialang.org/u/Sevi)
#### Post date: [November 3, 2025, 10:07pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/4 "2025-11-03T22:07:31Z")
Something seems strange indeed. I’m not an expert on ModelingToolkit.jl or BoundaryValueDiffEq.jl, but I can reproduce your issue on my machine.
> **versioninfo()**
>
> ```julia-auto
> Julia Version 1.12.1
> Commit ba1e628ee49 (2025-10-17 13:02 UTC)
> Build Info:
> Official https://julialang.org release
> Platform Info:
> OS: Linux (x86_64-linux-gnu)
> CPU: 64 × AMD Ryzen Threadripper PRO 3975WX 32-Cores
> WORD_SIZE: 64
> LLVM: libLLVM-18.1.7 (ORCJIT, znver2)
> GC: Built with stock GC
> Threads: 1 default, 1 interactive, 1 GC (on 64 virtual cores)
> Environment:
> JULIA_PKG_USE_CLI_GIT = true
> JULIA_MAX_NUM_PRECOMPILE_FILES = 50
>
> ```
> **Packages**
>
> ```julia-auto
> [764a87c0] BoundaryValueDiffEq v5.18.0
> [f3b72e0c] DiffEqDevTools v2.48.0
> [20393b10] InfiniteOpt v0.5.9
> [b6b21f68] Ipopt v1.12.1
> [961ee093] ModelingToolkit v10.26.1
> [91a5bcdd] Plots v1.41.1
>
> ```
## BVProblem
This version only using BoundaryValueDiffEq.jl (which I think matches your original set of equations?) seems to work fine and identifies the initial condition for `dy ≈ -1.31`.
```julia
using BoundaryValueDiffEq
function testproblem!(du, u, p, t)
du[1] = u[2]
du[2] = u[1]
end
function bc!(residual, u, p, t)
residual[1] = u(ti)[1] - 1
residual[2] = u(tf)[1]
end
prob = BVProblem(testproblem!, bc!, [1.0, 1.0], (ti, tf))
sol = solve(prob, MIRK4(), dt = 0.001)
Plots.plot(sol)
```
## JuMPDynamicOptProblem
The code you tried with the Jump problem runs, but does the wrong thing:
```julia
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables y(..)
eqs = [D(D(y(t))) ~ y(t)]
@parameters deriv_guess
(ti, tf) = (0.0, 1.0)
cons = [y(tf) ~ 0.0]
@named bvpsys = System(eqs, t; constraints=cons)
bvpsys = mtkcompile(bvpsys)
u0map = [y(t) => 1.0, D(y(t)) => deriv_guess]
parammap = [deriv_guess => -1.0]
using InfiniteOpt, Ipopt, DiffEqDevTools
jprob = JuMPDynamicOptProblem(bvpsys, [u0map; parammap], (ti, tf); dt=0.001)
jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructRadauIIA5()))
```
As you mentioned the corresponding solution is clearly incorrect:
But I did get this warning which seems highly relevant:
```plaintext
┌ Warning: The control problem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by op) exceeds the total number of states. The solvers will default to doing a nonlinear least-squares optimization.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/b28X4/src/systems/optimal_control_interface.jl:99
```
So (as the plot also suggests) the solution is actually just a least squares fit with three datapoints (`y(ti)`, `y(tf)`, and the specified initial `dy(ti)` which was not treated as guess but rather as a fixed initial condition. Since there is no solution of the equation that satisfies these three points, `y(t)` immediately jumps to what would be the correct initial condition for the other two points.
It reminds me of a recent thread here about this warning about overdetermined systems:
> [@Overdetermined system when using modelingtoolkitize](https://discourse.julialang.org/t/overdetermined-system-when-using-modelingtoolkitize/132917):
>
> Hi, I would like to convert a numeric differential-algebraic equation to a symbolic model using modelingtoolkitize. However, when then trying to convert it back to an ODEProblem, I get an overdetermined warning and when trying to solve it an InitialFailue. The overdetermined warning is expected since I have 2 initial conditions + 1 algebraic equation but only 2 unknowns. I would like to get rid of the initial condition on u[2] but am not sure how to do this. Here’s and MWE with some comments ex…
It’s a different problem type, but the solution there was to specify the initial condition to be determined as `nothing`. I tried that here, but it didn’t work:
```julia
u0map = [y(t) => 1.0, D(y(t)) => nothing]
parammap = []
jprob = JuMPDynamicOptProblem(bvpsys, [u0map; parammap], (ti, tf); dt=0.001)
```
```julia-auto
ERROR: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:
yˍt(t)
```
… but specifying the “initial guess” leads to the above warning.
## MTK into BVProblem
Your original approach actually “works” for me sometimes, in the sense that I could create the problem once, but also crashes with a segfault either at the `solve` step or the construction of the `BVProblem`. I couldn’t really find out more so far …
```julia
prob = BVProblem(bvpsys, [1.0, 1.0], (ti, tf))
sol = solve(prob, MIRK4(), dt=0.001)
```
* * *
Overall I feel like this should work. Perhaps there is really something wrong with how the conversion should happen (then the documentation should probably be improved). Or it’s a bug somewhere 🤷
---
### Author:  [@bradyelster](https://discourse.julialang.org/u/bradyelster)
#### Post date: [November 4, 2025, 6:00pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/5 "2025-11-04T18:00:49Z")
Okay, I think I got it now. My main issue was with how initial values were being passed to the solver. ModelingToolkit requires a vector of `Pair` objects mapping unknowns to their initial values. Guesses for unknown initial conditions or parameters are supplied via the `guesses` keyword argument. The following code produces the correct solution now! 🙂
```julia-auto
using ModelingToolkit, BoundaryValueDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables y(..)
eq = [D(D(y(t))) ~ y(t)]
# right boundary condition treated as constraint
constr = [y(1.0) ~ 0.0]
@named bvp = System(
eq, t;
constraints=constr
)
bvpsys = mtkcompile(bvp)
# pass uncertain parameters/initial conditions to solver via `guesses`
prob = BVProblem(bvpsys, [y(t) => 1], (0.0, 1.0), guesses=[D(y(t)) => 0.0])
sol = solve(prob, MIRK4(), dt=0.01)
plot(sol)
```
Thank you @Sevi for taking the time to help!
---
### Author:  [@Sevi](https://discourse.julialang.org/u/Sevi)
#### Post date: [November 4, 2025, 6:45pm UTC](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491/6 "2025-11-04T18:45:43Z")
That’s great, glad you found the solution!