Boundary conditions should always use solution object

[ErikQQY]

Fix #107
Fix #185

Supesed #208

Let's finish this! Now I have only changed MIRK methods when evaluating boundary conditions, the workaround is to change values during __mirk_loss but build a new one in __mirk_loss_bc where there would be dual numbers since we are using automatic differentiation to build sparse Jacobain.

using BoundaryValueDiffMIRK
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -9.81 * sin(θ)
end
function bc!(residual, u, p, t)
    residual[1] = u(pi/4)[1] + pi / 2
    residual[2] = u(pi/2)[1] - pi / 2
end
prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
sol = solve(prob, MIRK4(), dt = 0.05)

[github-actions]

Benchmark Results

	master	b5a7f5d...	master/b5a7f5d6640964...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	6.33 ± 0.17 ms	6.39 ± 0.18 ms	0.99
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	2.41 ± 0.1 ms	2.43 ± 0.11 ms	0.991
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	0.849 ± 0.048 ms	0.85 ± 0.047 ms	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	0.898 ± 0.063 ms	0.901 ± 0.074 ms	0.997
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	0.999 ± 0.11 ms	1 ± 0.095 ms	0.997
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.95 ± 0.57 ms	1.91 ± 0.59 ms	1.02
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.54 ± 0.89 ms	3.57 ± 0.35 ms	0.992
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0526 ± 0.0068 s	0.0535 ± 0.01 s	0.982
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0647 ± 0.0074 s	0.0694 ± 0.015 s	0.932
Simple Pendulum/IIP/Shooting(Tsit5())	0.237 ± 0.076 ms	0.254 ± 0.078 ms	0.932
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.0385 ± 0.0029 s	0.0406 ± 0.003 s	0.949
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	11.2 ± 0.8 ms	11.7 ± 0.76 ms	0.954
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	3.24 ± 0.21 ms	3.38 ± 0.2 ms	0.958
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	3.32 ± 0.25 ms	3.4 ± 0.25 ms	0.978
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	3.43 ± 0.27 ms	3.53 ± 0.25 ms	0.971
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.39 ± 2.5 ms	3.3 ± 2.6 ms	1.03
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	5.98 ± 4.4 ms	6.04 ± 4.6 ms	0.989
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0985 ± 0.0031 s	0.107 ± 0.009 s	0.917
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.119 ± 0.0072 s	0.12 ± 0.0055 s	0.995
Simple Pendulum/OOP/Shooting(Tsit5())	0.602 ± 0.032 ms	0.603 ± 0.034 ms	0.999
time_to_load	7.02 ± 0.12 s	6.95 ± 0.15 s	1.01

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ErikQQY]

OK, this PR only changed the boundary condition convention to solution-evaluation way in MIRK methods:

function bc!(residual, u, p, t)
	residual[1] = u(pi / 4)[1] + pi / 2
	residual[2] = u(pi / 2)[1] - pi / 2
end

If it's okay, I will replace other solvers as well. Funny that when I am benchmarking the changes in this PR, the solving is getting a bit faster🤯

This PR

julia> @benchmark sol = solve(prob, MIRK4(), dt = 0.05)
BenchmarkTools.Trial: 9976 samples with 1 evaluation.
 Range (min … max):  398.833 μs …  21.430 ms  ┊ GC (min … max):  0.00% … 97.59%
 Time  (median):     434.625 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   500.182 μs ± 563.331 μs  ┊ GC (mean ± σ):  12.65% ± 10.75%

  █▄                                                            ▁
  ██▇▅▄▃▁▁▁▃▁▁▁▃▁▃▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▃▆▇ █
  399 μs        Histogram: log(frequency) by time       4.41 ms <

 Memory estimate: 1.08 MiB, allocs estimate: 24342.

Master branch

julia> @benchmark sol = solve(prob, MIRK4(), dt = 0.05)
BenchmarkTools.Trial: 9520 samples with 1 evaluation.
 Range (min … max):  411.166 μs …  15.173 ms  ┊ GC (min … max):  0.00% … 96.72%
 Time  (median):     459.167 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   524.214 μs ± 598.530 μs  ┊ GC (mean ± σ):  11.43% ±  9.61%

  █▂                                                            ▁
  ██▄▃▃▁▃▄▄▁▃▁▁▁▄▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▁▃▁▁▃▁▁▁▃▁▃▁▁▃▁▁▁▃▆ █
  411 μs        Histogram: log(frequency) by time       5.38 ms <

 Memory estimate: 1.03 MiB, allocs estimate: 23762.