Use finite-difference AD for BVP default nonlinear solvers

[ChrisRackauckas-Claude]

Ignore until reviewed by @ChrisRackauckas.

Summary

This changes the internal BVP-compatible default nonlinear solver polyalgorithms to use AutoFiniteDiff() when no autodiff backend is specified.

The master downgrade failures were coming from the downgraded NonlinearSolve stack selecting ForwardDiff through these internal defaults, which propagated dual numbers into user boundary-condition residual arrays and failed with Float64(::ForwardDiff.Dual...).

Validation

• git diff --check
• Runic: julia +1 --startup-file=no --project=/home/crackauc/sandbox/tmp_20260708_121627_10236/runic_env -e 'using Runic; exit(Runic.main(["--check", "--verbose", "."]))' checked 110 files successfully.
• Downgraded Shooting Core on Julia 1.10 after rebasing onto e2b50b00: BoundaryValueDiffEqShooting tests passed with Shooting Basic Problems Tests | 82 pass, Shooting NLLS Tests | 56 pass, and Shooting Orbital Tests | 40 pass.
• Downgraded MIRK Core on Julia 1.10 after rebasing onto e2b50b00: BoundaryValueDiffEqMIRK tests passed with MIRK Basic Tests | 200 pass, MIRK NLLS Tests | 72 pass, MIRK Ensemble Tests | 10 pass, MIRK Singular BVP Tests | 8 pass, MIRK VectorOfVector Initials Tests | 2 pass, and MIRK Dynamic Optimization Tests | 6 pass.
• After PR creation, master advanced to 4dfc3ef2 via release-prep version-only Project.toml bumps. I rebased over that and reran git diff --check plus Runic successfully.

Master-failure investigation: reproduced on clean master and bisected the Shooting downgrade failure to fe630c867aae5706a2e2a6edd1967519ec3c91f8, where the SciMLBase compat floor moved the downgraded stack to SciMLBase 3.30 / NonlinearSolveBase 2.33 / NonlinearSolveFirstOrder 2.1+.

[github-actions]

Benchmark Results (Julia v1.11)

Time benchmarks

	master	5ccf7f7...	master / 5ccf7f7...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	0.667 ± 0.012 s	0.646 ± 0.013 s	1.03 ± 0.028
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	13.3 ± 0.35 ms	13.2 ± 0.19 ms	1.01 ± 0.03
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	2.81 ± 0.073 ms	2.81 ± 0.075 ms	1 ± 0.037
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	3.38 ± 0.085 ms	3.39 ± 0.082 ms	0.997 ± 0.035
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	1.61 ± 0.046 ms	1.61 ± 0.053 ms	1 ± 0.044
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.55 ± 0.069 ms	1.55 ± 0.064 ms	1 ± 0.061
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.15 ± 0.12 ms	3.15 ± 0.13 ms	1 ± 0.055
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0413 ± 0.0031 s	0.0414 ± 0.0033 s	0.999 ± 0.11
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0631 ± 0.00067 s	0.0628 ± 0.00087 s	1 ± 0.018
Simple Pendulum/IIP/Shooting(Tsit5())	0.181 ± 0.083 ms	0.178 ± 0.083 ms	1.01 ± 0.66
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.796 ± 0.0031 s	0.786 ± 0.0093 s	1.01 ± 0.013
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	16.8 ± 5.4 ms	16.6 ± 5.4 ms	1.01 ± 0.46
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	3.41 ± 0.13 ms	3.44 ± 0.21 ms	0.99 ± 0.072
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	4.11 ± 0.2 ms	4.12 ± 0.19 ms	0.996 ± 0.067
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	1.96 ± 0.12 ms	1.95 ± 0.13 ms	1 ± 0.088
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.55 ± 0.62 ms	3.52 ± 0.9 ms	1.01 ± 0.31
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	6.73 ± 6.3 ms	6.84 ± 6.3 ms	0.984 ± 1.3
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0933 ± 0.0033 s	0.0927 ± 0.0037 s	1.01 ± 0.053
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.144 ± 0.0034 s	0.143 ± 0.0052 s	1 ± 0.043
Simple Pendulum/OOP/Shooting(Tsit5())	0.674 ± 0.065 ms	0.635 ± 0.046 ms	1.06 ± 0.13
time_to_load	7.45 ± 0.014 s	7.46 ± 0.023 s	0.999 ± 0.0036

Memory benchmarks

	master	5ccf7f7...	master / 5ccf7f7...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	0.389 M allocs: 0.0443 GB	0.389 M allocs: 0.0443 GB	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	0.0436 M allocs: 4.88 MB	0.0436 M allocs: 4.88 MB	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	15.7 k allocs: 1.65 MB	15.7 k allocs: 1.65 MB	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	22.1 k allocs: 2.03 MB	22.1 k allocs: 2.03 MB	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	12.8 k allocs: 1.05 MB	12.8 k allocs: 1.05 MB	1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	25.5 k allocs: 1.82 MB	25.5 k allocs: 1.82 MB	1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	0.049 M allocs: 3.39 MB	0.049 M allocs: 3.39 MB	1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.553 M allocs: 0.0535 GB	0.553 M allocs: 0.0535 GB	1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.833 M allocs: 0.0778 GB	0.833 M allocs: 0.0778 GB	1
Simple Pendulum/IIP/Shooting(Tsit5())	4.64 k allocs: 0.224 MB	4.64 k allocs: 0.224 MB	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.89 M allocs: 0.984 GB	0.89 M allocs: 0.984 GB	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	0.0932 M allocs: 24.8 MB	0.0932 M allocs: 24.8 MB	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	0.0324 M allocs: 3.96 MB	0.0324 M allocs: 3.96 MB	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	0.045 M allocs: 4.99 MB	0.045 M allocs: 4.99 MB	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	25.3 k allocs: 2.17 MB	25.3 k allocs: 2.17 MB	1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	0.142 M allocs: 10.2 MB	0.142 M allocs: 10.2 MB	1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	0.266 M allocs: 18.7 MB	0.266 M allocs: 18.7 MB	1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	2.52 M allocs: 0.279 GB	2.52 M allocs: 0.279 GB	1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	3.82 M allocs: 0.404 GB	3.82 M allocs: 0.404 GB	1
Simple Pendulum/OOP/Shooting(Tsit5())	0.0373 M allocs: 1.69 MB	0.0373 M allocs: 1.69 MB	1
time_to_load	0.159 k allocs: 11.2 kB	0.159 k allocs: 11.2 kB	1

[ChrisRackauckas-Claude]

Verified against the clean current ModelingToolkit master downstream failure from SciML/ModelingToolkit.jl#4731.

Baseline environment:

• Julia 1.12.6
• ModelingToolkitBase 1.51.0 from SciML/ModelingToolkit.jl@86a4146139
• BoundaryValueDiffEqMIRK 1.17.2
• BoundaryValueDiffEqCore 2.7.3
• BoundaryValueDiffEqAscher 1.15.2
• clean depot and a disposable local environment (no test/Project.toml)

Baseline focused command:

timeout 3600 env JULIA_DEPOT_PATH="$PWD/../mtkb-bvp-repro-depot" \
  julia +1.12 --startup-file=no --project=../mtkb-bvp-repro-env \
  lib/ModelingToolkitBase/test/bvproblem.jl

It reproduced the CI error in the first Lotka-Volterra MIRK solve after 2m00.2s:

MethodError: no method matching Float64(::ForwardDiff.Dual...)
update_eval_sol! @ BoundaryValueDiffEqMIRK/src/interpolation.jl:380
__mirk_loss! @ BoundaryValueDiffEqMIRK/src/mirk.jl:450
Test Summary:  | Error  Total     Time
Lotka-Volterra |     1      1  2m00.2s

I then developed only this PR's lib/BoundaryValueDiffEqCore at version 2.7.3, restored released MIRK 1.17.2, and reran the same complete focused test file. It exited 0 with:

Test Summary:  | Pass  Total     Time
Lotka-Volterra |    4      4  2m55.7s
Test Summary:                  | Pass  Total     Time
Boundary Condition Compilation |    2      2  1m30.8s
Test Summary:        | Pass  Broken  Total     Time
ODE with constraints |    6       3      9  1m45.5s
Test Summary:             | Pass  Broken  Total     Time
Cost function compilation |    6       2      8  1m04.5s
Test Summary:        | Pass  Total     Time
Parameter estimation |    2      2  2m35.8s

This covers both errors from the ModelingToolkitBase InterfaceII job: Lotka-Volterra at bvproblem.jl:14 and cost-function compilation at bvproblem.jl:308. The existing assertions were unchanged.

For comparison, MIRK PR #542's type-adaptive interpolation cache removed the conversion error and made Lotka pass, but both cost-function assertions subsequently evaluated as NaN ≈ NaN. This PR is the only tested candidate that fixes both downstream errors and preserves finite cost results.

[ChrisRackauckas-Claude]

CI audit for the red FIRK AD Julia 1.12 job: this is an exact clean-master dependency regression, not introduced by #559. Clean master 4dfc3ef also signal-11s in the residual-comprehension allocation frame with Enzyme 0.13.182. The hosted boundary is exact: Enzyme 0.13.181 passed FIRK AD 9/9; 0.13.182 crashed without a FIRK source change. The post-CI 0.13.183 release now resolves by default, and I ran the native FIRK AD group locally on clean master: 9/9 passed in 18m05.5s. Existing focused draft #547 already hardens that exact allocation frame; independent reproduction/boundary evidence is recorded at #547 (comment). No #559 code change is indicated.

[ChrisRackauckas-Claude]

Final audit for the red FIRK NESTED job (86740825722): I reproduced its first failing configuration exactly and ran both nonlinear solvers on both revisions.

Environment: Julia 1.12.6, LinearSolve 4.3.0, BoundaryValueDiffEqCore 2.7.3, BoundaryValueDiffEqFIRK 1.17.2, SciMLBase 3.35.2, NonlinearSolveBase 2.34.1, NonlinearSolveFirstOrder 2.2.0, ForwardDiff 1.4.1. Case: UnderconstrainedProbArr[1], LobattoIIIa4, nested_nlsolve=true, dt=0.1, abstol=reltol=0.1.

Results:

• Use finite-difference AD for BVP default nonlinear solvers #559 source, NewtonRaphson(): retcode=Success in 1228.423704 s.
• Clean master, NewtonRaphson(): retcode=Success in 1244.463348 s.
• Use finite-difference AD for BVP default nonlinear solvers #559 source, GaussNewton(): retcode=Success in 1246.962855 s.
• Clean master, GaussNewton(): retcode=Success in 1228.765616 s.

The hosted job had already passed the basic FIRK NESTED block (223 pass / 26 broken) before reporting SingularException(72) for this same underconstrained case, then was cancelled after four hours. Four exact isolated successes rule out a deterministic #559 regression in the first reported case. The remaining distinction is full-suite ordering/resource state on the hosted runner; this evidence does not justify changing or suppressing any assertion.