Joining two ODE solutions

Is there any way to join two ODE solutions so that I can use the final product as a function like an usual ODE solution?

My problem is that I have a function α(ϕ) defined in two regimes, so I have to solve two ODEs. The first starts with physical initial conditions at ϕ=0, the second starts when the first stops and continues until ϕ=2. The first step is then to join the two ODE solutions from ϕ=[0,2], the second to symmetrise this function so that, in the end, the solution is valid for ϕ=[-2,2].

Below is a MWE in which I solve the ODEs, then extract the solution and time in arrays and use append! to glue the solutions from the two regimes, and finally create a DiffEqArray to relate the complete solution with the values of ϕ. Plotting this array show the correct solution, but I cannot use it as a function, which is necessary for me to solve another ODE. Is there an alternative to what I have tried to do? Maybe reconstructing an ODESolution type from the two ODE solutions?

using DifferentialEquations, Plots, RecursiveArrayTools
const global lp=1.
s2 = 81.0
d = -0.0009
B = 1.0
β = 1.0
k = 1.0
α₀ = 0.083163
ϕ₀ = 2.0

function quantum!(dv,v,p,ϕ)
    s2,d,B,β,k=p
    α = v[1]
    #@show (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))
    dv[1] = (ϕ*s2*sin(2*d*α)+2*d*sinh(s2*α*ϕ))/(-α*s2*sin(2*d*α)+2*d*cos(2*d*α)+2*d*cosh(s2*α*ϕ))
end
# When dα/dϕ = 1, we reach the classical regime, and stop the integration
condition(v,ϕ,integrator) = (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))==1.0
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)
# Initial Condition at the bounce
a₀ = [α₀]
classspan = (0,ϕ₀)
probquant = ODEProblem(quantum!,a₀,classspan,(s2,d,B,β,k))
solquant = solve(probquant,Tsit5(),callback=cb)

function classic!(du,u,p,ϕ)
    αc = u[1]
    dαc = u[2]
    du[1] = dαc
    du[2] = 3*(-dαc^3/sqrt(2)+dαc^2+dαc/sqrt(2)-1)
end
# IC retrieved from end of quantum phase
init = [last(solquant);1.0]
classspan = (last(solquant.t),ϕ₀)
probclass = ODEProblem(classic!,init,classspan)
solclass = solve(probclass,Tsit5())

# I tried to join the two regime like this:
# α(ϕ) for ϕ>0
solu = append!(solquant[1,:],solclass[1,:]) # α
solt = append!(solquant.t,solclass.t) # ϕ

# α(ϕ) for ϕ<0
soloppu = reverse(solu)
soloppt = -reverse(solt)
pop!(soloppu) # Remove -0.0 (redundant)
pop!(soloppt)

# Join the two solutions
soltotu = append!(soloppu,solu)
soltott = append!(soloppt,solt)
soltot = DiffEqArray(soltotu,soltott)
plot(soltot) # Plot is fine, but how to use soltot as a function?

@Vaibhavdixit02 didn’t you have a solution for this?

Yeah it’s this code in DiffEqParameEstim.jl DiffEqParamEstim.jl/multiple_shooting_objective.jl at 227ececc2af8b24746e395c88973ad9f7f87e0af · SciML/DiffEqParamEstim.jl · GitHub

3 Likes

We should get something on this added to SciMLBase.

Sorry for the late reply!

So it seems that I should add a line similar to

sol_new = DiffEqBase.build_solution(prob, alg, sol_loss.t, sol_loss.u, retcode =:Success)

from Vaibhavdixit02’s answer. I’m a bit confused on the prob option though. Applying this to my case, I initially thought I could have used something like

soltot = DiffEqBase.build_solution(prob, Tsit5(), soltott, soltotu,retcode = :Success)

instead of the DiffEqArray at the end of my code. However, I don’t know what I should use instead of prob since I have both probquant and probclass as problems. This issue does not seem to appear in Vaibhavdixit02’s answer since prob is passed as an argument in function multiple_shooting_objective. Any idea?

The prob is used for the f and p in the interpolation of the lazy algorithms. For the most part, you should be okay merging like this.