# 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:** 1
**Showing post:** 4

<div class="post-metadata">

### Author: ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@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")

</div>

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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/b/db6fc00ee746100158b2af72ce7d984608d29cef.png)

## 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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/4/145e736c81517902de5709ce92910f8d335886cd.png)

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 🤷

---

_[View the full topic](https://discourse.julialang.org/t/boundary-value-problems-in-modelingtoolkit/133491)._