# 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
#@show (ϕ*s2*sin(2*d*v)+2*d*sinh(s2*v*ϕ))/(-v*s2*sin(2*d*v)+2*d*cos(2*d*v)+2*d*cosh(s2*v*ϕ))
dv = (ϕ*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)+2*d*sinh(s2*v*ϕ))/(-v*s2*sin(2*d*v)+2*d*cos(2*d*v)+2*d*cosh(s2*v*ϕ))==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
dαc = u
du = dαc
du = 3*(-dαc^3/sqrt(2)+dαc^2+dαc/sqrt(2)-1)
end
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.

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.