BoundaryValueDiffEq: How to restart to exploit initial guess?

Specific Domains
1

My boundary value problem has a parameter that renders the problem easy or hard.

Is there a way to restart the solve function for the hard problem by specifying the solution of the easy problem as initial guess?

Below a MWE (solves Burgers equation by harmonic balancing. If nonlinscaling=0 (see below), problem reduces to linear Helmholtz equation)

# set source function 
src(x) = 1. 

# set right-hand side 
function rhs!(du, u, p, x)
    c, omd = p
    C = u[1]; A = u[2]; B = u[3]
    dC = du[1]; dA = du[2]; dB = du[3]
    nonlinscaling = 0.2
    du[1] = nonlinscaling*(A*dA+B*dB+C*dC)
    du[2] = omd/c*B + nonlinscaling*(A*dC+C*dA) 
    du[3] = -omd/c*A - src(x)/c + nonlinscaling*(B*dC+C*dB) 
end

# set boundary conditions 
function bc(residual, u, p, x)
    Cleft = 0; Aleft = 0; Bleft=0; 
    residual[1] = u[1][1] - Cleft
    residual[2] = u[1][2] - Aleft 
    residual[3] = u[1][3] - Bleft 
end

# set constants 
c = 0.2
omd = 2*pi
p = [c, omd]

# set domain size 
xspan = (0.,1.)

# set initial guess 
start = [0., 1., 1.]

# define and solve the problem 
bvp = BVProblem(rhs!, bc, start, xspan,p)
sol = solve(bvp, MIRK4(), dt = 0.005)
2

BoundaryValueDifEq.jl does support using the solution from previous solving in the next solving, just pass the solution to the construction of BVProblem with BVProblem(the!, bc, sol.u, xspan, p) and call the solver on this new problem.

1 Like
3

I run into a method error.

Is there any example on this?

Thx!

4

Could you please share your error stacktrace?

5

Sure.

solrestart = sol;

nonlinscaling = 0.5 

# store parameter value in p
p = [c, omd, nonlinscaling]

# define and solve the problem 
bvp = BVProblem(rhs!, bc, solrestart.u, xspan,p)
sol = solve(bvp, MIRK4(), dt = 0.005)

results in

MethodError: no method matching zero(::Type{Vector{Float64}})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Dates/src/periods.jl:53
  zero(::G) where G<:Graphs.AbstractGraph at ~/.julia/packages/Graphs/FXxqo/src/interface.jl:346
  zero(::OffsetArrays.OffsetArray) at ~/.julia/packages/OffsetArrays/xdnHf/src/OffsetArrays.jl:407
  ...

Stacktrace:
 [1] zeros(#unused#::Type{Vector{Float64}}, dims::Tuple{Int64})
   @ Base ./array.jl:589
 [2] zeros(#unused#::Type{Vector{Float64}}, dims::Int64)
   @ Base ./array.jl:584
 [3] vector_alloc
   @ ~/.julia/packages/BoundaryValueDiffEq/xfc07/src/vector_auxiliary.jl:5 [inlined]
 [4] BoundaryValueDiffEq.BVPSystem(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(rhs!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(bc), SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, x::Vector{Float64}, order::Int64)
   @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/xfc07/src/collocation.jl:21
 [5] __solve(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(rhs!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(bc), SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::MIRK4; dt::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/xfc07/src/solve.jl:32
 [6] solve_call(_prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(rhs!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(bc), SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::MIRK4; merge_callbacks::Bool, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:221
 [7] solve_up(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(rhs!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(bc), SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Vector{Vector{Float64}}, p::Vector{Float64}, args::MIRK4; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:254
 [8] solve(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(rhs!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(bc), SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::MIRK4; sensealg::Nothing, u0::Nothing, p::Nothing, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/S7V8q/src/solve.jl:234
 [9] top-level scope
   @ In[13]:10
6

What version are you using? Looking at the stack trace, src/vector_auxiliary.jl and src/collocation.jl don’t seem to be in current BoundaryValueDiffEq.jl. But searching outside the repo in GitHub, they are in forks dating back about 5-7 years.

You seem to be on Julia 1.8, and BoundaryValueDiffEq.jl has changed a lot since then. Can report the Pkg status ]st?

2 Likes
7

You are using an old version of BoundaryValueDiffEq.jl, exploiting initial guess requires the latest version of BoundaryValueDiffEq.jl

2 Likes
8

Switching to Julia 1.10 solved the issue.

You are heros! Many thx!

1 Like