Sol.interp and sol.k DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

Hello, please help me understand what it is sol.interp, sol.k and is there any way not to save them? I found that they take up quite a lot of RAM, which is why ubuntu 24.04 freezes.

Thank you for your helps!

Code:

using StaticArrays, DifferentialEquations

function FHN2_try3(u, p ,t)
    x1, y1, x2, y2, z= u
    ϵ, a, g, k, σ, α, k1, k2 = p

    I(ϕ_i) = g * (1.0/(1.0 + exp(k*(cos(σ/2) - cos(ϕ_i - α - σ/2)))))
    ρz = k1 + k2 * z ^ 2

    ϕ2 = atan(y2, x2)
    ϕ1 = atan(y1, x1)

    dx1dt = (x1 - x1 ^ 3 / 3 - y1 + I(ϕ2) + ρz * (x2 - x1) ) / ϵ
    dy1dt = x1 - a
    dx2dt = (x2 - x2 ^ 3 / 3 - y2 + I(ϕ1) + ρz * (x1 - x2) ) / ϵ
    dy2dt = x2 - a
    dzdt = x1 - x2
    return SVector(dx1dt, dy1dt, dx2dt, dy2dt, dzdt)
end

function FHN2_try3_params()
    ϵ = 0.01; a = -1.01;
    g = 0.1; k = 50.0; σ = 50.0 * pi / 180; α = 160.0 * pi / 180;
    k1 = 0.0; k2 = 0.0
    return [ ϵ, a, g, k, σ, α, k1, k2]
end

function jac_FHN(u, p, t)
    x1, y1, x2, y2, z = u
    ϵ, a, g, k, σ, α, k1, k2 = p

    pz = k1 + k2 * z^2

    exp_ϕ2 = exp( k * ( cos(σ/2) - cos( α + σ/2 - atan(y2, x2) ) ) )
    sin_ϕ2 = sin( α + σ/2 - atan(y2, x2) )
    exp_ϕ1 = exp( k * ( cos(σ/2) - cos( α + σ/2 - atan(y1, x1) ) ) )
    sin_ϕ1 = sin( α + σ/2 - atan(y1, x1) )
    
    x1x1 = (1.0 - x1^2 - pz) / ϵ
    x1y1 = -1.0/ϵ
    x1x2 = (-g * k * y2 * exp_ϕ2 * sin_ϕ2) / ( x2^2 * (1 + y2^2 / x2^2 ) * ( exp_ϕ2 + 1.0 )^2 ) + k1 + k2 * z^2
    x1x2 =  x1x2 / ϵ
    x1y2 = ( g * k * x2 * exp_ϕ2 * sin_ϕ2 ) / ( (x2^2 + y2^2) * ( 1.0 + exp_ϕ2 )^2 )
    x1y2 = x1y2 / ϵ
    x1z = 2.0 * k2 * z * ( -x1 + x2 ) / ϵ
    #-----------------------------------------------
    y1x1 = 1.0;
    y1y1 = 0.0; y1x2 = 0.0; y1y2 = 0.0; y1z = 0.0;
    #----------------------------------------------
    x2x2 = (1.0 - x2^2 - pz) / ϵ
    x2y2 = -1.0/ϵ
    x2x1 = (-g * k * y1 * exp_ϕ1 * sin_ϕ1) / ( x1^2 * (1 + y1^2 / x1^2 ) * ( exp_ϕ1 + 1.0 )^2 ) + k1 + k2 * z^2
    x2x1 =  x2x1 / ϵ
    x2y1 = ( g * k * x1 * exp_ϕ1 * sin_ϕ1 ) / ( (x1^2 + y1^2) * ( 1.0 + exp_ϕ1 )^2 )
    x2y1 = x2y1 / ϵ
    x2z = 2.0 * k2 * z * ( -x2 + x1 ) / ϵ
    #----------------------------------------------
    y2x2 = 1.0;
    y2y2 = 0.0; y2x1 = 0.0; y2y1 = 0.0; y2z = 0.0;
    #----------------------------------------------
    zx1 = 1.0; zy1 = 0.0; zx2 = -1.0; zy2 = 0.0; zz = 0.0;
    return SMatrix{5,5}(x1x1, y1x1, x2x1, y2x1, zx1,
                        x1y1, y1y1, x2y1, y2y1, zy1,
                        x1x2, y1x2, x2x2, y2x2, zx2,
                        x1y2, y1y2, x2y2, y2y2, zy2,
                        x1z,  y1z,  x2z,  y2z,  zz)
end

parameters = FHN2_try3_params()
tspan = (0.0, 10000.0)
parameters[7] = 0.09
parameters[8] = 75.74
u0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525] 
u0 = SVector{5}(u0)
prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
alg = Vern9();


sol = solve(prob, alg, adaptive = false, dt = 0.001, maxiters = 1e7) #true, abstol = abs_tolerance, reltol = rel_tolerance, maxiters = max_iters);

sol_alg = sol.alg
sol_alg_choice = sol.alg_choice
sol_dense = sol.dense
sol_errors = sol.errors
sol_interp = sol.interp
sol_k = sol.k
sol_prob = sol.prob
sol_retcode = sol.retcode
sol_stats = sol.stats
sol_t = sol.t
sol_tslocation = sol.tslocation
sol_u = sol.u
sol_u_analytic = sol.u_analytic

varinfo()
1 Like
2

Do you need dense output? If you don’t, then just saveat?

1 Like
3

Thanks, I’ll try to use saveat. I don’t understand what dense output is used for, so I can’t tell if I need it.

4

If you don’t know what points you will need until after a solution, saying doing some analysis and playing around in the REPL, or generating new interpolations on-demand (rootfinding and such), then you’d want the dense continuous output. If you know what values you want, just use saveat and directly compute those. That will be cheaper since then it doesn’t have to save the whole continuous output.

5

After the solution has been received, I analyze the time series. I define spikes in the time series, the extreme events and calculate the time intervals between extreme events.
I tried using saveat and it didn’t seem to affect the result.
Thank you for your help

6

That kind of post analysis should probably be done on the continuous time dense version unless you have a very dense saveat save (which would take more memory than the continuous one). Though if you don’t have strict accuracy constraints then having a less dense saveat would be a way to cut memory.

7

one approach to look into here would be putting the analysis in a callback. writing the code to work well in a callback is a little tricky, but doing so would likely allow you to skip saving entirely

8

I need to calculate the trajectory for quite a long time, for example time = 5000000 , so that later, build a probability density function for time intervals between exteme events. Therefore, I try to use the accuracy as high as possible. I am using the Verne9() and abstol, reltol = 1e-14

9

In my case, this will not work. I need a complete time series.

10

I tried to calculate the solution by parts. Then I noticed that starting from a certain time in the time diagram, the spikes look too discrete. For example part of time series

image
image1920×1080 84.8 KB

Code for calculate solution by parts:

using StaticArrays, DifferentialEquations, Statistics, JLD2 

function clear_work_space(sol, len_sol, ttr)
    sol = nothing
    len_sol = nothing
    ttr = nothing
    return sol, len_sol, ttr
end

t_truncate(t) = floor(Int64, t / 2)

function save_data(iteration, data, path_to_save)
    file_name = "$(iteration)_sol_x1.jld"
    full_path = path_to_save * file_name
    jldsave(full_path; data)
end

function first_iteration(u0, parameters, tspan, integrator_setting, path_to_save)

    prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
    sol = solve(prob, integrator_setting.alg, adaptive = true,
    abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
    maxiters = integrator_setting.max_iters);
    last_point = sol[end]
    println("last point: $(last_point)"); flush(stdout)
    len_sol = length(sol.t)

    ttr = t_truncate(len_sol); tend = len_sol
    data_x1 = [sol[1, ttr:tend], sol.t[ttr:tend]]
    sol, len_sol, ttr = clear_work_space(sol, len_sol, ttr)
    GC.gc()

    save_data(1, data_x1, path_to_save)
    return last_point
end

function calculate_timeseris(u0_start, parameters, integrator_setting,
    t_point, count_iteration, path_to_save)

    global tspan = (0.0, t_point)
    global u0 = u0_start

    u0 = first_iteration(u0, parameters, tspan, integrator_setting, path_to_save)
    println("----------------------------------------------"); flush(stdout)
    println("");flush(stdout)

    for iteration in range(2, count_iteration)

        println("ITERATION: $(iteration)")

        tspan = (t_point*(iteration-1), t_point*iteration)
        prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
        println("u0: $(u0)"); flush(stdout)
        sol = solve(prob, integrator_setting.alg, adaptive = true,
        abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
        maxiters = integrator_setting.max_iters);
        u0 = sol[end]

        data_x1 = [sol[1, :], sol.t]
        sol = nothing; GC.gc();
        save_data(iteration, data_x1, path_to_save)

        println("last point: $(u0)"); flush(stdout)
        println("----------------------------------------------"); flush(stdout)
        println("");flush(stdout)
    end
end

alg = Vern9()
max_iters = 1e8
abs_tol = 1e-14; rel_tol = 1e-14
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)

path_to_save = "/home/sergey/MEGA/dynamical-systems/FHN_Korotkov/data/timeseries_k2_75_72/"
parameters = FHN2_try3_params()
parameters[7] = 0.09
parameters[8] = 75.72

u0_start = [1.7, 0.7, -1.4, 0.35, 0.7 - 0.35]
u0_start = SVector{5}(u0_start)

t_point = 50000.0
count_iteration = 30

calculate_timeseris(u0_start, parameters, integrator_setting,
    t_point, count_iteration, path_to_save)
11

If you need to decrease the memory, putting the analysis into a callback as Oscar describes would be a good way to go.

12

Thank you
I’ll try to save only local minimums and maximums using an callback