diff --git a/src/solve/mirk.jl b/src/solve/mirk.jl index c2ccc795..45857e73 100644 --- a/src/solve/mirk.jl +++ b/src/solve/mirk.jl @@ -303,13 +303,13 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo jac_prototype = vcat(init_jacobian(cache_bc), init_jacobian(cache_collocation)) jac = if iip - (J, u, p) -> __mirk_mpoint_jacobian!(J, u, p, jac_alg.bc_diffmode, + (J, u, p) -> __mirk_mpoint_jacobian!(J, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode, cache_bc, cache_collocation, loss_bcₚ, - loss_collocationₚ, resid_bc, resid_collocation, cache.M, L) + loss_collocationₚ, resid_bc, resid_collocation, L) else - (u, p) -> __mirk_mpoint_jacobian(u, p, jac_prototype, jac_alg.bc_diffmode, + (u, p) -> __mirk_mpoint_jacobian(jac_prototype, u, jac_alg.bc_diffmode, jac_alg.nonbc_diffmode, cache_bc, cache_collocation, loss_bcₚ, - loss_collocationₚ, cache.M, L) + loss_collocationₚ, L) end nlf = NonlinearFunction{iip}(loss; resid_prototype = vcat(resid_bc, resid_collocation), @@ -317,17 +317,17 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo return (L == cache.M ? NonlinearProblem : NonlinearLeastSquaresProblem)(nlf, y, cache.p) end -function __mirk_mpoint_jacobian!(J, x, p, bc_diffmode, nonbc_diffmode, bc_diffcache, +function __mirk_mpoint_jacobian!(J, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache, loss_bc::BC, loss_collocation::C, resid_bc, resid_collocation, - M::Int, L::Int) where {BC, C} + L::Int) where {BC, C} sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, resid_bc, x) sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache, loss_collocation, resid_collocation, x) return nothing end -function __mirk_mpoint_jacobian(x, p, J, bc_diffmode, nonbc_diffmode, bc_diffcache, - nonbc_diffcache, loss_bc::BC, loss_collocation::C, M::Int, L::Int) where {BC, C} +function __mirk_mpoint_jacobian(J, x, bc_diffmode, nonbc_diffmode, bc_diffcache, + nonbc_diffcache, loss_bc::BC, loss_collocation::C, L::Int) where {BC, C} sparse_jacobian!(@view(J[1:L, :]), bc_diffmode, bc_diffcache, loss_bc, x) sparse_jacobian!(@view(J[(L + 1):end, :]), nonbc_diffmode, nonbc_diffcache, loss_collocation, x) @@ -358,10 +358,10 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo jac_prototype = init_jacobian(diffcache) jac = if iip - (J, u, p) -> __mirk_2point_jacobian!(J, u, p, jac_alg.diffmode, diffcache, lossₚ, + (J, u, p) -> __mirk_2point_jacobian!(J, u, jac_alg.diffmode, diffcache, lossₚ, resid) else - (u, p) -> __mirk_2point_jacobian(u, p, jac_prototype, jac_alg.diffmode, diffcache, + (u, p) -> __mirk_2point_jacobian(u, jac_prototype, jac_alg.diffmode, diffcache, lossₚ) end @@ -370,12 +370,12 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo return (L == cache.M ? NonlinearProblem : NonlinearLeastSquaresProblem)(nlf, y, cache.p) end -function __mirk_2point_jacobian!(J, x, p, diffmode, diffcache, loss_fn::L, resid) where {L} +function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid) where {L} sparse_jacobian!(J, diffmode, diffcache, loss_fn, resid, x) return J end -function __mirk_2point_jacobian(x, p, J, diffmode, diffcache, loss_fn::L) where {L} +function __mirk_2point_jacobian(x, J, diffmode, diffcache, loss_fn::L) where {L} sparse_jacobian!(J, diffmode, diffcache, loss_fn, x) return J end diff --git a/test/runtests.jl b/test/runtests.jl index b7349b04..8b07648f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,6 +14,11 @@ const GROUP = uppercase(get(ENV, "GROUP", "ALL")) @time @safetestset "Orbital" begin include("shooting/orbital.jl") end + if VERSION ≥ v"1.10-" + @time @safetestset "Shooting NLLS Tests" begin + include("shooting/nonlinear_least_squares.jl") + end + end end end @@ -31,8 +36,10 @@ const GROUP = uppercase(get(ENV, "GROUP", "ALL")) @time @safetestset "Interpolation Tests" begin include("mirk/interpolation_test.jl") end - @time @safetestset "MIRK Nonlinear Least Squares Tests" begin - include("mirk/nonlinear_least_squares.jl") + if VERSION ≥ v"1.10-" + @time @safetestset "MIRK NLLS Tests" begin + include("mirk/nonlinear_least_squares.jl") + end end end end diff --git a/test/shooting/nonlinear_least_squares.jl b/test/shooting/nonlinear_least_squares.jl new file mode 100644 index 00000000..7108f7e5 --- /dev/null +++ b/test/shooting/nonlinear_least_squares.jl @@ -0,0 +1,96 @@ +using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test + +@testset "Overconstrained BVP" begin + SOLVERS = [ + Shooting(Tsit5(); nlsolve = LevenbergMarquardt()), + Shooting(Tsit5(); nlsolve = GaussNewton()), + MultipleShooting(10, Tsit5(); nlsolve = LevenbergMarquardt()), + MultipleShooting(10, Tsit5(); nlsolve = GaussNewton())] + + # OOP MP-BVP + f1(u, p, t) = [u[2], -u[1]] + + function bc1(sol, p, t) + t₁, t₂ = extrema(t) + solₜ₁ = sol(t₁) + solₜ₂ = sol(t₂) + solₜ₃ = sol((t₁ + t₂) / 2) + # We know that this overconstrained system has a solution + return [solₜ₁[1], solₜ₂[1] - 1, solₜ₃[1] - 0.51735, solₜ₃[2] + 1.92533] + end + + tspan = (0.0, 100.0) + u0 = [0.0, 1.0] + + bvp1 = BVProblem(BVPFunction{false}(f1, bc1; bcresid_prototype = zeros(4)), u0, tspan) + + for solver in SOLVERS + @time sol = solve(bvp1, solver; + nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), + verbose = false) + @test norm(bc1(sol, nothing, sol.t)) < 1e-4 + end + + # IIP MP-BVP + function f1!(du, u, p, t) + du[1] = u[2] + du[2] = -u[1] + return nothing + end + + function bc1!(resid, sol, p, t) + (t₁, t₂) = extrema(t) + solₜ₁ = sol(t₁) + solₜ₂ = sol(t₂) + solₜ₃ = sol((t₁ + t₂) / 2) + # We know that this overconstrained system has a solution + resid[1] = solₜ₁[1] + resid[2] = solₜ₂[1] - 1 + resid[3] = solₜ₃[1] - 0.51735 + resid[4] = solₜ₃[2] + 1.92533 + return nothing + end + + bvp2 = BVProblem(BVPFunction{true}(f1!, bc1!; bcresid_prototype = zeros(4)), u0, tspan) + + for solver in SOLVERS + @time sol = solve(bvp2, solver; + nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), + verbose = false) + resid_f = Array{Float64}(undef, 4) + bc1!(resid_f, sol, nothing, sol.t) + @test norm(resid_f) < 1e-4 + end + + # OOP TP-BVP + bc1a(ua, p) = [ua[1]] + bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] + + bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + + for solver in SOLVERS + @time sol = solve(bvp3, solver; + nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), + verbose = false) + @test norm(vcat(bc1a(sol(0.0), nothing), bc1b(sol(100.0), nothing))) < 1e-4 + end + + # IIP TP-BVP + bc1a!(resid, ua, p) = (resid[1] = ua[1]) + bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) + + bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + + for solver in SOLVERS + @time sol = solve(bvp4, solver; + nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), + verbose = false) + resida = Array{Float64}(undef, 1) + residb = Array{Float64}(undef, 2) + bc1a!(resida, sol(0.0), nothing) + bc1b!(residb, sol(100.0), nothing) + @test norm(vcat(resida, residb)) < 1e-4 + end +end diff --git a/test/shooting/shooting_tests.jl b/test/shooting/shooting_tests.jl index 1fa04d11..53012d06 100644 --- a/test/shooting/shooting_tests.jl +++ b/test/shooting/shooting_tests.jl @@ -80,101 +80,6 @@ using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test end end -@testset "Overconstrained BVP" begin - SOLVERS = [ - Shooting(Tsit5(); nlsolve = LevenbergMarquardt()), - Shooting(Tsit5(); nlsolve = GaussNewton()), - MultipleShooting(10, Tsit5(); nlsolve = LevenbergMarquardt()), - MultipleShooting(10, Tsit5(); nlsolve = GaussNewton())] - - # OOP MP-BVP - f1(u, p, t) = [u[2], -u[1]] - - function bc1(sol, p, t) - t₁, t₂ = extrema(t) - solₜ₁ = sol(t₁) - solₜ₂ = sol(t₂) - solₜ₃ = sol((t₁ + t₂) / 2) - # We know that this overconstrained system has a solution - return [solₜ₁[1], solₜ₂[1] - 1, solₜ₃[1] - 0.51735, solₜ₃[2] + 1.92533] - end - - tspan = (0.0, 100.0) - u0 = [0.0, 1.0] - - bvp1 = BVProblem(BVPFunction{false}(f1, bc1; bcresid_prototype = zeros(4)), u0, tspan) - - for solver in SOLVERS - @time sol = solve(bvp1, solver; - nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), - verbose = false) - @test norm(bc1(sol, nothing, sol.t)) < 1e-4 - end - - # IIP MP-BVP - function f1!(du, u, p, t) - du[1] = u[2] - du[2] = -u[1] - return nothing - end - - function bc1!(resid, sol, p, t) - (t₁, t₂) = extrema(t) - solₜ₁ = sol(t₁) - solₜ₂ = sol(t₂) - solₜ₃ = sol((t₁ + t₂) / 2) - # We know that this overconstrained system has a solution - resid[1] = solₜ₁[1] - resid[2] = solₜ₂[1] - 1 - resid[3] = solₜ₃[1] - 0.51735 - resid[4] = solₜ₃[2] + 1.92533 - return nothing - end - - bvp2 = BVProblem(BVPFunction{true}(f1!, bc1!; bcresid_prototype = zeros(4)), u0, tspan) - - for solver in SOLVERS - @time sol = solve(bvp2, solver; - nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), - verbose = false) - resid_f = Array{Float64}(undef, 4) - bc1!(resid_f, sol, nothing, sol.t) - @test norm(resid_f) < 1e-4 - end - - # OOP TP-BVP - bc1a(ua, p) = [ua[1]] - bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] - - bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) - - for solver in SOLVERS - @time sol = solve(bvp3, solver; - nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), - verbose = false) - @test norm(vcat(bc1a(sol(0.0), nothing), bc1b(sol(100.0), nothing))) < 1e-4 - end - - # IIP TP-BVP - bc1a!(resid, ua, p) = (resid[1] = ua[1]) - bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) - - bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) - - for solver in SOLVERS - @time sol = solve(bvp4, solver; - nlsolve_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiters = 1000), - verbose = false) - resida = Array{Float64}(undef, 1) - residb = Array{Float64}(undef, 2) - bc1a!(resida, sol(0.0), nothing) - bc1b!(residb, sol(100.0), nothing) - @test norm(vcat(resida, residb)) < 1e-4 - end -end - @testset "Shooting with Complex Values" begin # Test for complex values function f1!(du, u, p, t)