Skip to content

Commit e0a9c26

Browse files
authored
Merge pull request #279 from ErikQQY/qqy/initial_guess
2 parents a94a33e + 7689568 commit e0a9c26

File tree

9 files changed

+91
-12
lines changed

9 files changed

+91
-12
lines changed

docs/Project.toml

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
1515
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1616
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
1717
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
18+
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
1819
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
1920
SimpleBoundaryValueDiffEq = "be0294bd-f90f-4760-ac4e-3421ce2b2da0"
2021

@@ -34,5 +35,6 @@ LineSearch = "0.1.4"
3435
LinearAlgebra = "1.10"
3536
LinearSolve = "2.36.2"
3637
OrdinaryDiffEq = "6.90.1"
38+
Plots = "1"
3739
SciMLBase = "2.60.0"
3840
SimpleBoundaryValueDiffEq = "1.1.0"

docs/pages.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
pages = ["index.md",
44
"Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
5+
"Tutorials" => "tutorials/continuation.md",
56
"Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
67
"basics/solve.md", "basics/autodiff.md"],
78
"Solver Summaries and Recommendations" => Any[

docs/src/refs.bib

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,4 +16,14 @@ @article{ascher1994collocation
1616
pages={938--952},
1717
year={1994},
1818
publisher={SIAM}
19-
}
19+
}
20+
21+
@book{ascher1995numerical,
22+
author = {Ascher, Uri M. and Mattheij, Robert M. M. and Russell, Robert D.},
23+
title = {Numerical Solution of Boundary Value Problems for Ordinary Differential Equations},
24+
publisher = {Society for Industrial and Applied Mathematics},
25+
year = {1995},
26+
doi = {10.1137/1.9781611971231},
27+
URL = {https://epubs.siam.org/doi/abs/10.1137/1.9781611971231},
28+
eprint = {https://epubs.siam.org/doi/pdf/10.1137/1.9781611971231}
29+
}

docs/src/tutorials/continuation.md

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
# Solve BVP with Continuation
2+
3+
Continuation is a commonly used technique for solving numerically difficult boundary value problems, we exploit the priori knowledge of the solution as initial guess to accelerate the BVP solving by breaking up the difficult BVP into a sequence of simpler problems. For example, we use the problem from [ascher1995numerical](@Citet) in this tutorial:
4+
5+
```math
6+
\epsilon y'' + xy' = \epsilon \pi^2\cos(\pi x) - \pi x\sin(\pi x)
7+
```
8+
9+
for $\epsilon=10^{-4}$, on $t\in[-1,1]$ with two point boundary conditions $y(-1)=-2, y(1)=0$. With analitical solution of $y(x)=\cos(\pi x)+\text{erf}(\frac{x}{\sqrt{2\epsilon}})/\text{erf}(\frac{1}{\sqrt{2\epsilon}})$, this problem has a rapid transition layer at $x=0$, making it difficult to solve numerically. In this tutorial, we will showcase how to use continuation with BoundaryValueDiffEq.jl to solve this BVP.
10+
11+
We use the substitution to transform this problem into a first order BVP system:
12+
13+
```math
14+
y_1'=y_2\\
15+
y_2'= -\frac{x}{e}y_2 - \pi^2\cos(\pi x) - \frac{\pi x}{e}\sin(\pi x)
16+
```
17+
18+
Since this BVP would become difficult to solve when $0<\epsilon\ll 1$, we start the continuation with relatively bigger $\epsilon$ to first obtain a good initial guess for cases when $\epsilon$ are becoming extremely small. We can just use the previous solution from BVP solving as the initial guess `u0` when constructing a new `BVProblem`.
19+
20+
```@example continuation
21+
using BoundaryValueDiffEq, Plots
22+
function f!(du, u, p, t)
23+
du[1] = u[2]
24+
du[2] = -t / p * u[2] - pi^2 * cos(pi * t) - pi * t / p * sin(pi * t)
25+
end
26+
function bc!(res, u, p, t)
27+
res[1] = u[1][1] + 2
28+
res[2] = u[end][1]
29+
end
30+
tspan = (-1.0, 1.0)
31+
sol = [1.0, 0.0]
32+
e = 0.1
33+
for i in 2:4
34+
global e = e / 10
35+
prob = BVProblem(f!, bc!, sol, tspan, e)
36+
global sol = solve(prob, MIRK4(), dt = 0.01)
37+
end
38+
plot(sol, idxs = [1])
39+
```
40+
41+
In the iterative solving, the intermediate solutions are each used as the initial guess for the next problem solving.
42+
43+
## On providing initial guess
44+
45+
There are several ways of providing initial guess in `BVProblem`/`TwoPointBVProblem`:
46+
47+
1. Solution from BVP/ODE solving from SciML packages with `ODESolution` type.
48+
2. `VectorOfArray` from RecursiveArrayTools.jl.
49+
3. `DiffEqArray` from RecursiveArrayTools.jl.
50+
4. Function handle of the form `f(p, t)` for specifying initial guess on time span.
51+
5. `AbstractArray` represent only the possible initial condition.

lib/BoundaryValueDiffEqCore/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "BoundaryValueDiffEqCore"
22
uuid = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
33
authors = ["Qingyu Qu <[email protected]>"]
4-
version = "1.6.0"
4+
version = "1.7.0"
55

66
[deps]
77
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

lib/BoundaryValueDiffEqCore/src/utils.jl

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -204,6 +204,15 @@ function __extract_problem_details(
204204
return Val(true), eltype(u0), length(u0), Int(cld(t₁ - t₀, dt)), u0
205205
end
206206

207+
function __extract_problem_details(
208+
prob, u0::SciMLBase.ODESolution; dt = 0.0, check_positive_dt::Bool = false)
209+
# Problem passes in a initial guess function
210+
check_positive_dt && dt 0 && throw(ArgumentError("dt must be positive"))
211+
_u0 = first(u0.u)
212+
_t = u0.t
213+
return Val(true), eltype(_u0), length(_u0), (length(_t) - 1), _u0
214+
end
215+
207216
function __initial_guess(f::F, p::P, t::T) where {F, P, T}
208217
if hasmethod(f, Tuple{P, T})
209218
return f(p, t)
@@ -372,6 +381,7 @@ Takes the input initial guess and returns the value at the starting mesh point.
372381
@inline __extract_u0(u₀::DiffEqArray, p, t₀) = u₀.u[1]
373382
@inline __extract_u0(u₀::F, p, t₀) where {F <: Function} = __initial_guess(u₀, p, t₀)
374383
@inline __extract_u0(u₀::AbstractArray, p, t₀) = u₀
384+
@inline __extract_u0(u₀::SciMLBase.ODESolution, p, t₀) = u₀.u[1]
375385
@inline __extract_u0(u₀::T, p, t₀) where {T} = error("`prob.u0::$(T)` is not supported.")
376386

377387
"""
@@ -383,6 +393,8 @@ Takes the input initial guess and returns the mesh.
383393
@inline __extract_mesh(u₀, t₀, t₁, dt::Number) = collect(t₀:dt:t₁)
384394
@inline __extract_mesh(u₀::DiffEqArray, t₀, t₁, ::Int) = u₀.t
385395
@inline __extract_mesh(u₀::DiffEqArray, t₀, t₁, ::Number) = u₀.t
396+
@inline __extract_mesh(u₀::SciMLBase.ODESolution, t₀, t₁, ::Int) = u₀.t
397+
@inline __extract_mesh(u₀::SciMLBase.ODESolution, t₀, t₁, ::Number) = u₀.t
386398

387399
"""
388400
__has_initial_guess(u₀) -> Bool
@@ -392,6 +404,7 @@ Returns `true` if the input has an initial guess.
392404
@inline __has_initial_guess(u₀::AbstractVector{<:AbstractArray}) = true
393405
@inline __has_initial_guess(u₀::VectorOfArray) = true
394406
@inline __has_initial_guess(u₀::DiffEqArray) = true
407+
@inline __has_initial_guess(u₀::SciMLBase.ODESolution) = true
395408
@inline __has_initial_guess(u₀::F) where {F} = true
396409
@inline __has_initial_guess(u₀::AbstractArray) = false
397410

@@ -404,6 +417,7 @@ guess is supplied, it returns `-1`.
404417
@inline __initial_guess_length(u₀::AbstractVector{<:AbstractArray}) = length(u₀)
405418
@inline __initial_guess_length(u₀::VectorOfArray) = length(u₀)
406419
@inline __initial_guess_length(u₀::DiffEqArray) = length(u₀.t)
420+
@inline __initial_guess_length(u₀::SciMLBase.ODESolution) = length(u₀.t)
407421
@inline __initial_guess_length(u₀::F) where {F} = -1
408422
@inline __initial_guess_length(u₀::AbstractArray) = -1
409423

@@ -417,6 +431,7 @@ initial guess, it returns `vec(u₀)`.
417431
vec, hcat, u₀)
418432
@inline __flatten_initial_guess(u₀::VectorOfArray) = mapreduce(vec, hcat, u₀.u)
419433
@inline __flatten_initial_guess(u₀::DiffEqArray) = mapreduce(vec, hcat, u₀.u)
434+
@inline __flatten_initial_guess(u₀::SciMLBase.ODESolution) = mapreduce(vec, hcat, u₀.u)
420435
@inline __flatten_initial_guess(u₀::AbstractArray) = vec(u₀)
421436
@inline __flatten_initial_guess(u₀::F) where {F} = nothing
422437

@@ -435,6 +450,9 @@ end
435450
@inline function __initial_guess_on_mesh(u₀::DiffEqArray, mesh, p)
436451
return copy(u₀)
437452
end
453+
@inline function __initial_guess_on_mesh(u₀::SciMLBase.ODESolution, mesh, p)
454+
return copy(VectorOfArray(u₀.u))
455+
end
438456
@inline function __initial_guess_on_mesh(u₀::AbstractArray, mesh, p)
439457
return VectorOfArray([copy(vec(u₀)) for _ in mesh])
440458
end

lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -410,12 +410,11 @@ end
410410
-pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
411411
sol3 = solve(bvp3, RadauIIa5(), dt = 0.05)
412412

413-
# Needs a SciMLBase fix
414413
bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2),
415414
pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
416-
@test_broken solve(bvp4, RadauIIa5(), dt = 0.05) isa SciMLBase.ODESolution
415+
@test SciMLBase.successful_retcode(solve(bvp4, RadauIIa5(), dt = 0.05))
417416

418417
bvp5 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), DiffEqArray(sol3.u, sol3.t),
419418
(0, pi / 2), pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
420-
@test SciMLBase.successful_retcode(solve(bvp5, RadauIIa5(), dt = 0.05).retcode)
419+
@test SciMLBase.successful_retcode(solve(bvp5, RadauIIa5(), dt = 0.05))
421420
end

lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -435,14 +435,13 @@ end =#
435435
-pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
436436
sol3 = solve(bvp3, RadauIIa5(; nested_nlsolve = true), dt = 0.05)
437437

438-
# Needs a SciMLBase fix
439438
bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2),
440439
pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
441-
@test_broken solve(bvp4, RadauIIa5(; nested_nlsolve = true), dt = 0.05) isa
442-
SciMLBase.ODESolution
440+
@test_broken SciMLBase.successful_retcode(solve(
441+
bvp4, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
443442

444443
bvp5 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), DiffEqArray(sol3.u, sol3.t),
445444
(0, pi / 2), pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
446445
@test_broken SciMLBase.successful_retcode(solve(
447-
bvp5, RadauIIa5(; nested_nlsolve = true), dt = 0.05).retcode)
446+
bvp5, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
448447
end

lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -279,14 +279,13 @@ end
279279
-pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
280280
sol3 = solve(bvp3, MIRK4(), dt = 0.05)
281281

282-
# Needs a SciMLBase fix
283282
bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2),
284283
pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
285-
@test_broken solve(bvp4, MIRK4(), dt = 0.05) isa SciMLBase.ODESolution
284+
@test SciMLBase.successful_retcode(solve(bvp4, MIRK4(), dt = 0.05))
286285

287286
bvp5 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), DiffEqArray(sol3.u, sol3.t),
288287
(0, pi / 2), pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
289-
@test SciMLBase.successful_retcode(solve(bvp5, MIRK4(), dt = 0.05).retcode)
288+
@test SciMLBase.successful_retcode(solve(bvp5, MIRK4(), dt = 0.05))
290289
end
291290

292291
@testitem "Compatibility with StaticArrays" begin

0 commit comments

Comments
 (0)