Skip to content

Commit ca01c8c

Browse files
ErikQQYErikQQY
authored
Merge pull request #426 from SciML/qqy/singular
Collocation solvers handle singular BVP
2 parents 61772fb + e1c0a5d commit ca01c8c

9 files changed

Lines changed: 209 additions & 28 deletions

File tree

‎lib/BoundaryValueDiffEqCore/src/utils.jl‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -693,3 +693,20 @@ end
693693
@inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs, optimize_kwargs) = (;
694694
nlsolve_kwargs...,
695695
)
696+
697+
"""
698+
__add_singular_term!(K, singular_term, y, t)
699+
700+
Helper function to add the singular term contribution S * y / t to K for t > 0.
701+
Used in collocation residual computation for singular BVPs of the form y' = S*y/t + f(t,y).
702+
"""
703+
@inline function __add_singular_term!(K, singular_term::Nothing, y, t)
704+
return nothing
705+
end
706+
707+
@inline function __add_singular_term!(K, singular_term::AbstractMatrix, y, t)
708+
if t > 0
709+
mul!(K, singular_term, y, one(t) / t, one(t))
710+
end
711+
return nothing
712+
end

‎lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
2525
__build_solution, __Fix3, __split_kwargs, _sparse_like,
2626
get_dense_ad, __internal_optimization_problem,
2727
__internal_solve, __default_sparsity_detector, __build_cost,
28-
__tunable_part
28+
__tunable_part, __add_singular_term!
2929

3030
using ConcreteStructs: @concrete
3131
using DiffEqBase: DiffEqBase

‎lib/BoundaryValueDiffEqFIRK/src/collocation.jl‎

Lines changed: 19 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
function Φ!(residual, cache::FIRKCacheExpand, y, u, trait, constraint)
22
return Φ!(
33
residual, cache.fᵢ_cache, cache.k_discrete, cache.f, cache.TU, y, u, cache.p,
4-
cache.mesh, cache.mesh_dt, cache.stage, cache.f_prototype, trait, constraint
4+
cache.mesh, cache.mesh_dt, cache.stage, cache.f_prototype, cache.singular_term, trait, constraint
55
)
66
end
77

@@ -14,7 +14,7 @@ end
1414

1515
@views function Φ!(
1616
residual, fᵢ_cache, k_discrete, f!, TU::FIRKTableau{false}, y, u, p,
17-
mesh, mesh_dt, stage::Int, f_prototype, ::DiffCacheNeeded, ::Val{true}
17+
mesh, mesh_dt, stage::Int, f_prototype, singular_term, ::DiffCacheNeeded, ::Val{true}
1818
)
1919
(; c, a, b) = TU
2020
L_f_prototype = length(f_prototype)
@@ -59,7 +59,7 @@ end
5959

6060
@views function Φ!(
6161
residual, fᵢ_cache, k_discrete, f!, TU::FIRKTableau{false}, y, u, p,
62-
mesh, mesh_dt, stage::Int, f_prototype, ::DiffCacheNeeded, ::Val{false}
62+
mesh, mesh_dt, stage::Int, _, singular_term, ::DiffCacheNeeded, ::Val{false}
6363
)
6464
(; c, a, b) = TU
6565
tmp1 = get_tmp(fᵢ_cache, u)
@@ -81,7 +81,9 @@ end
8181
for r in 1:stage
8282
@. tmp1 = yᵢ
8383
__maybe_matmul!(tmp1, K, a[:, r], h, T(1))
84-
f!(residual[ctr + r], tmp1, p, mesh[i] + c[r] * h)
84+
t = mesh[i] + c[r] * h
85+
f!(residual[ctr + r], tmp1, p, t)
86+
__add_singular_term!(residual[ctr + r], singular_term, tmp1, t)
8587
residual[ctr + r] .-= K[:, r]
8688
end
8789

@@ -95,7 +97,7 @@ end
9597

9698
@views function Φ!(
9799
residual, fᵢ_cache, k_discrete, f!, TU::FIRKTableau{false}, y, u, p, mesh,
98-
mesh_dt, stage::Int, f_prototype, ::NoDiffCacheNeeded, ::Val{false}
100+
mesh_dt, stage::Int, _, singular_term, ::NoDiffCacheNeeded, ::Val{false}
99101
)
100102
(; c, a, b) = TU
101103
tmp1 = similar(fᵢ_cache)
@@ -117,7 +119,9 @@ end
117119
for r in 1:stage
118120
@. tmp1 = yᵢ
119121
__maybe_matmul!(tmp1, K, a[:, r], h, T(1))
120-
f!(residual[ctr + r], tmp1, p, mesh[i] + c[r] * h)
122+
t = mesh[i] + c[r] * h
123+
f!(residual[ctr + r], tmp1, p, t)
124+
__add_singular_term!(residual[ctr + r], singular_term, tmp1, t)
121125
residual[ctr + r] .-= K[:, r]
122126
end
123127

@@ -266,7 +270,7 @@ end
266270
function Φ(cache::FIRKCacheExpand, y, u, trait)
267271
return Φ(
268272
cache.fᵢ_cache, cache.k_discrete, cache.f, cache.TU, y, u,
269-
cache.p, cache.mesh, cache.mesh_dt, cache.stage, trait
273+
cache.p, cache.mesh, cache.mesh_dt, cache.stage, cache.singular_term, trait
270274
)
271275
end
272276

@@ -279,7 +283,7 @@ end
279283

280284
@views function Φ(
281285
fᵢ_cache, k_discrete, f, TU::FIRKTableau{false}, y,
282-
u, p, mesh, mesh_dt, stage::Int, ::DiffCacheNeeded
286+
u, p, mesh, mesh_dt, stage::Int, singular_term, ::DiffCacheNeeded
283287
)
284288
(; c, a, b) = TU
285289
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
@@ -302,7 +306,9 @@ end
302306
for r in 1:stage
303307
@. tmp1 = yᵢ
304308
__maybe_matmul!(tmp1, K, a[:, r], h, T(1))
305-
residuals[ctr + r] = f(tmp1, p, mesh[i] + c[r] * h)
309+
t = mesh[i] + c[r] * h
310+
residuals[ctr + r] = f(tmp1, p, t)
311+
__add_singular_term!(residuals[ctr + r], singular_term, tmp1, t)
306312
residuals[ctr + r] .-= K[:, r]
307313
end
308314

@@ -317,7 +323,7 @@ end
317323

318324
@views function Φ(
319325
fᵢ_cache, k_discrete, f, TU::FIRKTableau{false}, y,
320-
u, p, mesh, mesh_dt, stage::Int, ::NoDiffCacheNeeded
326+
u, p, mesh, mesh_dt, stage::Int, singular_term, ::NoDiffCacheNeeded
321327
)
322328
(; c, a, b) = TU
323329
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
@@ -340,7 +346,9 @@ end
340346
for r in 1:stage
341347
@. tmp1 = yᵢ
342348
__maybe_matmul!(tmp1, K, a[:, r], h, T(1))
343-
residuals[ctr + r] = f(tmp1, p, mesh[i] + c[r] * h)
349+
t = mesh[i] + c[r] * h
350+
residuals[ctr + r] = f(tmp1, p, t)
351+
__add_singular_term!(residuals[ctr + r], singular_term, tmp1, t)
344352
residuals[ctr + r] .-= K[:, r]
345353
end
346354

‎lib/BoundaryValueDiffEqFIRK/src/firk.jl‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
defect
2828
nest_prob
2929
resid_size
30+
singular_term
3031
nlsolve_kwargs
3132
optimize_kwargs
3233
kwargs
@@ -62,6 +63,7 @@ Base.eltype(::FIRKCacheNested{iip, T}) where {iip, T} = T
6263
fᵢ₂_cache
6364
defect
6465
resid_size
66+
singular_term
6567
nlsolve_kwargs
6668
optimize_kwargs
6769
kwargs
@@ -263,7 +265,7 @@ function init_nested(
263265
return FIRKCacheNested{iip, T, typeof(diffcache), fit_parameters}(
264266
alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p,
265267
alg, TU, ITU, f_prototype, bcresid_prototype, mesh, mesh_dt, k_discrete,
266-
y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size,
268+
y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size, prob.singular_term,
267269
nlsolve_kwargs, optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)
268270
)
269271
end
@@ -414,7 +416,7 @@ function init_expanded(
414416
return FIRKCacheExpand{iip, T, typeof(diffcache), fit_parameters}(
415417
alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p,
416418
alg, TU, ITU, f_prototype, bcresid_prototype, mesh, mesh_dt, k_discrete,
417-
y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size, nlsolve_kwargs,
419+
y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size, prob.singular_term, nlsolve_kwargs,
418420
optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...)
419421
)
420422
end

‎lib/BoundaryValueDiffEqFIRK/test/expanded/singular_bvp_tests.jl‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
@testitem "Singular BVP" tags = [:singular] begin
2+
using BoundaryValueDiffEqFIRK
3+
using LinearAlgebra
4+
5+
nested = false
6+
7+
for stage in (2, 3, 4, 5)
8+
s = Symbol("LobattoIIIa$(stage)")
9+
@eval lobattoIIIa_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...)
10+
end
11+
12+
for stage in (2, 3, 4, 5)
13+
s = Symbol("LobattoIIIb$(stage)")
14+
@eval lobattoIIIb_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...)
15+
end
16+
17+
for stage in (2, 3, 4, 5)
18+
s = Symbol("LobattoIIIc$(stage)")
19+
@eval lobattoIIIc_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...)
20+
end
21+
22+
for stage in (2, 3, 5, 7)
23+
s = Symbol("RadauIIa$(stage)")
24+
@eval radau_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...)
25+
end
26+
27+
# Lane-Emden equation of index 1:
28+
# y'' + (2/t)*y' + y = 0, y(0) = 1, y'(0) = 0
29+
# The exact solution is y(t) = sin(t)/t (with limit y(0) = 1)
30+
#
31+
# In first-order form: y[1]' = y[2], y[2]' = -y[1] - (2/t)*y[2]
32+
# This can be written as y' = S*y/t + f(t,y) where:
33+
# S = [0 0; 0 -2] and f(t,y) = [y[2]; -y[1]]
34+
35+
function lane_emden!(du, u, p, t)
36+
du[1] = u[2]
37+
du[2] = -u[1] # The -2*u[2]/t term is handled by singular_term
38+
end
39+
function lane_emden(u, p, t)
40+
return [u[2]; -u[1]]
41+
end
42+
43+
function lane_emden_bc_a!(resid, ua, p)
44+
resid[1] = ua[1] - 1.0 # y(0) = 1
45+
end
46+
function lane_emden_bc_b!(resid, ub, p)
47+
resid[1] = ub[1] - sin(1.0) # y(1) = sin(1) ≈ 0.8415
48+
end
49+
50+
lane_emden_bc_a(ua, p) = ua[1] - 1.0
51+
lane_emden_bc_b(ub, p) = ub[1] - sin(1.0)
52+
53+
# The singular term matrix S
54+
S = [0.0 0.0; 0.0 -2.0]
55+
56+
tspan = (0.0, 1.0)
57+
u0 = [1.0, 0.0]
58+
bcresid_prototype = (zeros(1), zeros(1))
59+
60+
prob_iip = TwoPointBVProblem(
61+
lane_emden!, (lane_emden_bc_a!, lane_emden_bc_b!), u0, tspan; bcresid_prototype
62+
)
63+
prob_oop = TwoPointBVProblem(
64+
lane_emden, (lane_emden_bc_a, lane_emden_bc_b), u0, tspan; bcresid_prototype
65+
)
66+
67+
# Test with different FIRK stages
68+
# Note: tolerance varies by order - FIRK2 is lower order so needs larger tolerance
69+
@testset "FIRK solvers with singular term" for prob in (prob_iip, prob_oop)
70+
for stage in (2, 3, 4, 5)
71+
sol = solve(prob, lobattoIIIa_solver(Val(stage)), dt = 0.01)
72+
@test SciMLBase.successful_retcode(sol)
73+
end
74+
for stage in (3, 4, 5)
75+
sol = solve(prob, lobattoIIIb_solver(Val(stage)), dt = 0.01)
76+
@test SciMLBase.successful_retcode(sol)
77+
end
78+
for stage in (3, 4, 5)
79+
sol = solve(prob, lobattoIIIc_solver(Val(stage)), dt = 0.01)
80+
@test SciMLBase.successful_retcode(sol)
81+
end
82+
for stage in (2, 3, 5, 7)
83+
sol = solve(prob, radau_solver(Val(stage)), dt = 0.01)
84+
@test SciMLBase.successful_retcode(sol)
85+
end
86+
end
87+
end

‎lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
2424
DiffCacheNeeded, NoDiffCacheNeeded, __split_kwargs,
2525
__concrete_kwargs, __FastShortcutNonlinearPolyalg,
2626
__construct_internal_problem, __internal_solve,
27-
__default_sparsity_detector, __build_cost
27+
__default_sparsity_detector, __build_cost, __add_singular_term!
2828

2929
using ConcreteStructs: @concrete
3030
using DiffEqBase: DiffEqBase

‎lib/BoundaryValueDiffEqMIRK/src/collocation.jl‎

Lines changed: 19 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
11
function Φ!(residual, cache::MIRKCache, y, u, trait, constraint)
22
return Φ!(
33
residual, cache.fᵢ_cache, cache.k_discrete, cache.f, cache.TU, y, u, cache.p,
4-
cache.mesh, cache.mesh_dt, cache.stage, cache.f_prototype, trait, constraint
4+
cache.mesh, cache.mesh_dt, cache.stage, cache.f_prototype, cache.singular_term, trait, constraint
55
)
66
end
77

88
@views function Φ!(
99
residual, fᵢ_cache, k_discrete, f!, TU::MIRKTableau, y, u, p, mesh,
10-
mesh_dt, stage::Int, f_prototype, ::DiffCacheNeeded, ::Val{true}
10+
mesh_dt, stage::Int, f_prototype, singular_term, ::DiffCacheNeeded, ::Val{true}
1111
)
1212
(; c, v, x, b) = TU
1313
L_f_prototype = length(f_prototype)
@@ -43,7 +43,7 @@ end
4343

4444
@views function Φ!(
4545
residual, fᵢ_cache, k_discrete, f!, TU::MIRKTableau, y, u, p, mesh,
46-
mesh_dt, stage::Int, _, ::DiffCacheNeeded, constraint::Val{false}
46+
mesh_dt, stage::Int, _, singular_term, ::DiffCacheNeeded, constraint::Val{false}
4747
)
4848
(; c, v, x, b) = TU
4949

@@ -60,7 +60,9 @@ end
6060
for r in 1:stage
6161
@. tmp = (1 - v[r]) * yᵢ + v[r] * yᵢ₊₁
6262
__maybe_matmul!(tmp, K[:, 1:(r - 1)], x[r, 1:(r - 1)], h, T(1))
63-
f!(K[:, r], tmp, p, mesh[i] + c[r] * h)
63+
t = mesh[i] + c[r] * h
64+
f!(K[:, r], tmp, p, t)
65+
__add_singular_term!(K[:, r], singular_term, tmp, t)
6466
end
6567

6668
# Update residual
@@ -71,7 +73,7 @@ end
7173

7274
@views function Φ!(
7375
residual, fᵢ_cache, k_discrete, f!, TU::MIRKTableau, y, u, p,
74-
mesh, mesh_dt, stage::Int, _, ::NoDiffCacheNeeded, ::Val{false}
76+
mesh, mesh_dt, stage::Int, _, singular_term, ::NoDiffCacheNeeded, ::Val{false}
7577
)
7678
(; c, v, x, b) = TU
7779

@@ -88,7 +90,9 @@ end
8890
for r in 1:stage
8991
@. tmp = (1 - v[r]) * yᵢ + v[r] * yᵢ₊₁
9092
__maybe_matmul!(tmp, K[:, 1:(r - 1)], x[r, 1:(r - 1)], h, T(1))
91-
f!(K[:, r], tmp, p, mesh[i] + c[r] * h)
93+
t = mesh[i] + c[r] * h
94+
f!(K[:, r], tmp, p, t)
95+
__add_singular_term!(K[:, r], singular_term, tmp, t)
9296
end
9397

9498
# Update residual
@@ -100,13 +104,13 @@ end
100104
function Φ(cache::MIRKCache, y, u, trait)
101105
return Φ(
102106
cache.fᵢ_cache, cache.k_discrete, cache.f, cache.TU, y, u,
103-
cache.p, cache.mesh, cache.mesh_dt, cache.stage, trait
107+
cache.p, cache.mesh, cache.mesh_dt, cache.stage, cache.singular_term, trait
104108
)
105109
end
106110

107111
@views function Φ(
108112
fᵢ_cache, k_discrete, f, TU::MIRKTableau, y, u,
109-
p, mesh, mesh_dt, stage::Int, ::DiffCacheNeeded
113+
p, mesh, mesh_dt, stage::Int, singular_term, ::DiffCacheNeeded
110114
)
111115
(; c, v, x, b) = TU
112116
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
@@ -123,7 +127,9 @@ end
123127
for r in 1:stage
124128
@. tmp = (1 - v[r]) * yᵢ + v[r] * yᵢ₊₁
125129
__maybe_matmul!(tmp, K[:, 1:(r - 1)], x[r, 1:(r - 1)], h, T(1))
126-
K[:, r] .= f(tmp, p, mesh[i] + c[r] * h)
130+
t = mesh[i] + c[r] * h
131+
K[:, r] .= f(tmp, p, t)
132+
__add_singular_term!(K[:, r], singular_term, tmp, t)
127133
end
128134

129135
# Update residual
@@ -136,7 +142,7 @@ end
136142

137143
@views function Φ(
138144
fᵢ_cache, k_discrete, f, TU::MIRKTableau, y, u, p,
139-
mesh, mesh_dt, stage::Int, ::NoDiffCacheNeeded
145+
mesh, mesh_dt, stage::Int, singular_term, ::NoDiffCacheNeeded
140146
)
141147
(; c, v, x, b) = TU
142148
residuals = [safe_similar(yᵢ) for yᵢ in y[1:(end - 1)]]
@@ -153,7 +159,9 @@ end
153159
for r in 1:stage
154160
@. tmp = (1 - v[r]) * yᵢ + v[r] * yᵢ₊₁
155161
__maybe_matmul!(tmp, K[:, 1:(r - 1)], x[r, 1:(r - 1)], h, T(1))
156-
K[:, r] .= f(tmp, p, mesh[i] + c[r] * h)
162+
t = mesh[i] + c[r] * h
163+
K[:, r] .= f(tmp, p, t)
164+
__add_singular_term!(K[:, r], singular_term, tmp, t)
157165
end
158166

159167
# Update residual

‎lib/BoundaryValueDiffEqMIRK/src/mirk.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
errors
2929
new_stages
3030
resid_size
31+
singular_term
3132
nlsolve_kwargs
3233
optimize_kwargs
3334
kwargs
@@ -228,8 +229,8 @@ function SciMLBase.__init(
228229
return MIRKCache{iip, T, use_both, typeof(diffcache), fit_parameters}(
229230
alg_order(alg), stage, N, size(X), f, bc, prob_, prob.problem_type, prob.p, alg,
230231
TU, ITU, f_prototype, bcresid_prototype, mesh, mesh_dt, k_discrete, k_interp, y,
231-
y₀, residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size, nlsolve_kwargs,
232-
optimize_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...)
232+
y₀, residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size, prob.singular_term
233+
, nlsolve_kwargs, optimize_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...)
233234
)
234235
end
235236

0 commit comments

Comments
 (0)