From 9094c05b59d86c890fa763425954faaea488df57 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 3 Apr 2025 00:21:59 +0800
Subject: [PATCH] Support derivative boundary conditions

---
 .../src/interpolation.jl                      | 12 +++++
 .../test/mirk_basic_tests.jl                  | 47 +++++++++++++++++++
 2 files changed, 59 insertions(+)

diff --git a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
index 0585de771..0b421fad8 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
@@ -132,3 +132,15 @@ function (s::EvalSol{C})(tval::Number) where {C <: MIRKCache}
     z .= z .* dt .+ u[ii]
     return z
 end
+
+function (s::EvalSol{C})(tval::Number, T::Type{Val{1}}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    z′ = zero(last(u))
+    ii = interval(t, tval)
+    dt = mesh_dt[ii]
+    τ = (tval - t[ii]) / dt
+    _, w′ = interp_weights(τ, alg)
+    __maybe_matmul!(z′, k_discrete[ii].du[:, 1:stage], w′[1:stage])
+    return z′
+end
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 2c86c3a9f..923596a53 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -381,3 +381,50 @@ end
             prob, MIRK4(), dt = 0.05, abstol = 1e-5, controller = error_control)
     end
 end
+
+@testitem "Derivative boundary conditions" begin
+    prob_bvp_nonlinear_14_f(y1, y2, y3, y4) = y2 * y3 - y1 * y4
+    function prob_bvp_nonlinear_14_f!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = u[3]
+        du[3] = u[4]
+        du[4] = prob_bvp_nonlinear_14_f(u[1], u[2], u[3], u[4])
+    end
+    function prob_bvp_nonlinear_14_bc!(res, u, p, t)
+        res[1] = u(0.0)[1]
+        res[2] = u(0.0)[2]
+        res[3] = u(1.0)[1] - 1
+        res[4] = u(1.0)[2]
+    end
+    prob_bvp_nonlinear_14_function = BVPFunction(
+        prob_bvp_nonlinear_14_f!, prob_bvp_nonlinear_14_bc!)
+    prob_bvp_nonlinear_14_tspan = (0.0, 1.0)
+
+    prob_bvp_nonlinear_normal = BVProblem(
+        prob_bvp_nonlinear_14_function, [1.0, 0.0, 0.0, 0.0], prob_bvp_nonlinear_14_tspan)
+
+    @test_nowarn sol1 = solve(
+        prob_bvp_nonlinear_normal, MIRK4(), dt = 0.01, adaptive = false)
+
+    prob_bvp_nonlinear_14_f(y1, y2, y3, y4) = y2 * y3 - y1 * y4
+    function prob_bvp_nonlinear_14_f!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = u[3]
+        du[3] = u[4]
+        du[4] = prob_bvp_nonlinear_14_f(u[1], u[2], u[3], u[4])
+    end
+    function prob_bvp_nonlinear_bc!(res, u, p, t)
+        res[1] = u(0.0)[1]
+        res[2] = u(0.0, Val{1})[1]
+        res[3] = u(1.0)[1] - 1
+        res[4] = u(1.0, Val{1})[1]
+    end
+    prob_bvp_nonlinear_function = BVPFunction(
+        prob_bvp_nonlinear_14_f!, prob_bvp_nonlinear_bc!)
+    prob_bvp_nonlinear_14_tspan = (0.0, 1.0)
+
+    prob_bvp_nonlinear_14 = BVProblem(
+        prob_bvp_nonlinear_function, [1.0, 0.0, 0.0, 0.0], prob_bvp_nonlinear_14_tspan)
+
+    @test_nowarn sol2 = solve(prob_bvp_nonlinear_14, MIRK4(), dt = 0.01, adaptive = false)
+end