1
$\begingroup$

I am solving two different ODE whose solutions need to be matched. I am currently doing this by hand, which works great, but I would like to automatise this process. The second ODE takes one of its initial conditions and its initial time from the last values of the first solution. I put a MWE below. Now, I have naively tried to put together the first and second ODEs, as I show in the second code below (giving a StackOverflowError). It seems this is not the way to go, though. Does someone know how to do that?

EDIT: the second code is now the solution

using DifferentialEquations, Plots

# First ODE, which depends on some parameters

function param_dadp!(s2,d,a0,ϕ₀)
    #ϕ₀ = abs(find_time!(s2,d))
    #ϕ₀=2.0
    function def_dadp!(dv,v,p,ϕ)
        s2,d=p
        α = 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
    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!)
    α₀ = [a0]
    tspan = (0,ϕ₀)
    probdadp = ODEProblem(def_dadp!,α₀,tspan,(s2,d))
    soldadp = solve(probdadp,Tsit5(),callback=cb)
end
plot(param_dadp!(81.0,-0.0009,0.083163,2.0))

# Second ODE, which needs initial α and ϕ from the previous solution
function test!(dv,v,p,ϕ)
    α = v[1]
    dα = v[2]
    dv[1] = dα
    dv[2] = 3*(-dα^9/sqrt(2)+dα^8+dα^7/sqrt(2)-dα^6)
end
init = [1.06,1.0] # Put IC by hand
testspan = (0.75,3.0) # Put initial ϕ by hand
prob = ODEProblem(test!,init,testspan)
sol = solve(prob,Tsit5())
plot!(sol,vars=(0,1),xlims=(0,2))
using DifferentialEquations, Plots

function param_dadp!(s2,d,a0,ϕ₀)
    function def_dadp!(dv,v,p,ϕ)
        s2,d=p
        α = 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
    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!)
    α₀ = [a0]
    tspan = (0,ϕ₀)
    probdadp = ODEProblem(def_dadp!,α₀,tspan,(s2,d))
    soldadp = solve(probdadp,Tsit5(),callback=cb)

    function classic!(du,u,p,ϕ)
        αc = u[1]
        dαc = u[2]
        du[1] = dαc
        du[2] = 3*(-dαc^9/sqrt(2)+dαc^8+dαc^7/sqrt(2)-dαc^6)
    end
    init = [last(soldadp);1.0]
    classspan = (last(soldadp.t),ϕ₀)
    probclass = ODEProblem(classic!,init,classspan)
    solclass = solve(probclass,Tsit5())

    solu = append!(soldadp[1,:],solclass[1,:])
    solt = append!(soldadp.t,solclass.t)
    sol = DiffEqArray(solu,solt)
end
plot(param_dadp!(81.0,-0.0009,0.083163,2.0))
```
$\endgroup$
3
$\begingroup$

last(soldadp) is an array, so [last(soldadp),1.0] is a Vector{Union{Vector{Float64},Float64}}, i.e. it's like [[2.0],1.0]. That's not what you wanted: you wanted [2.0,1.0] which means you wanted to use the concatenation operator ;. init = [last(soldadp);1.0] fixes your error.

Now, I don't know what you mean by sol = soldadp+solclass: those two time series are not aligned in time, and they hold different numbers of vectors. Did you want to append the plot of the second after the other? Then you can use plot! to mutate the plot. Or you can build an array that concatenates the two in the way you're interpreting, but the + operator doesn't seem to make much sense here. Either way, it's a 1 x N object and a 2 x M object, so it's going to need some user input on what you actually mean by "put these together".

$\endgroup$
7
  • $\begingroup$ Thanks for pointing out the error! Indeed, sol as such does not make sense. What I really want as output is to extend the function α(ϕ) computed on [0,last(soldadp)] with the α(ϕ) computed on [last(soldadp), ϕ₀] (I will need this function later). In the end, I'm interested in u[1] only in the second ODE, so I thought I could just add the first component of the two solutions using sol = append!(soldadp[1,:],solclass[1,:]), but the plot in output is wrong. I guess I'm not accessing the solutions correctly... $\endgroup$
    – Free_ion
    Jun 19 at 8:45
  • $\begingroup$ That should be fine. Just make sure to plot with t too. Also, the default plot uses the interpolation to generate more values, so either use saveat directly or interpolate to get a denser plot. $\endgroup$ Jun 19 at 12:08
  • $\begingroup$ How can I specify what t is? I see that the plot generated now goes from 0 to ~40 instead of 0 to 2. Did I lose time somehow when I used append? $\endgroup$
    – Free_ion
    Jun 19 at 14:38
  • $\begingroup$ How are you generating the data? soldadp(0.0:0.1:soldadp.t) etc.? $\endgroup$ Jun 19 at 18:54
  • 1
    $\begingroup$ You're on the right track. And something that may be helpful is the RecursiveArrayTools.jl DiffEqArray for having the t paired with the u for plotting. $\endgroup$ Jun 19 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.