How to save values from a differential equation not part of initial conditions (u) and then graph them?

New to Julia
1

Hello,

I have a food web that is driven by differential equations, that includes changes in plants, animals, and boat ‘population’. There are variables such as amount of fish harvested (catch_F) and fish price (fish_price) which influence the population of harvested fish species and ‘population’ of boats, but are not saved at each time step of the solution (like the populations).

In the end I would like to be able to plot fish price vs harvested fish population, and harvested fish vs boat ‘population’. But I can’t figure it out. I assume I need to save the value of fish harvested and fish price at each time step.

dBdT = function(dB, bioS, p, t)
     D, S, K, v, r, X, num_nutrient, num_plant, num_animal, num_vessel_types, bm, foodWeb, b, q, ω, c, h, e, target_species, catch_max, b0V, μ, elasticity, clearance_V_type = p

    results =  fill(NaN,2) # I need to pre-allocate an array for catch_F and  fish_price

        # rate of change in clearance rate
    if target_species != 0  # if there is some harvest
        clearance_V = create_clearance_V(clearance_max,bioS,β_V,b0V,clearance_min, target_species,clearance_V_type)
        catch_F = (catch_max*bioS[target_species]/(catch_max/clearance_V+bioS[target_species]))*bioS[(num_nutrient+num_plant+num_animal+num_vessel_types)]
        dB[num_nutrient+num_plant+num_animal+num_vessel_types+1] = 0.0
        fish_price = (catch_F/γ_fish)^(1/elasticity)
        revenue = fish_price*catch_F
        cost_V =  maintenance+clearance_cost_scale*clearance_V
        cost_F = cost_V*bioS[(num_nutrient+num_plant+num_animal+num_vessel_types)]
        results = [catch_F fish_price] 
        results2 =  vcat(results)
    else  # in absence of harvest
        dB[num_nutrient+num_plant+num_animal+num_vessel_types+1] = 0.0
        catch_F = 0.0
    end

... #rest of differential equations

    return results2 # I need to make results exist outside of the function
end

But it says results2 not defined?

And then how to graph it? I assume it will be

plot(sol, vars = ((num_nutrient+num_plant+num_animal+num_vessel_types),results2[1,]), title = "Clearance 1 - Vessels vs Catch", xlabel = "Vessel Coverage (1/m²)" ,ylabel = "Harvested Fish (g)", lw=1, legend=false)
2

The easiest thing to do is to just recompute those extra values using the post solution interpolation. I am working on a better interface for this though.

3

What is the post solution interpolation? Could you give me an example?

I managed to get results2 to save values, but it is very crude and ugly (time steps make large jumps).

(I needed to define results and results2 outside of the differential equation

result = fill(NaN,1,6)
result2 = deepcopy(result)

dBdT = function(dB, bioS, p, t)
        #things inside equation
        results = [catch_F fish_price] 
        global result2 = vcat(result2,result) # you want to stick the new results onto the collection of results
        # rest of equation

I am currently piecing a SavingCallback (link) together with help from @jonniedie on Gitter.

4

sol(t).

A SavingCallback is fine. What you’re doing here is scary though since it assumes time is monotonic, but it isn’t.

5

Could you elaborate on what is assuming time is monotonic? Is it the Callback, because it can go a step backwards and save data at that time point? Why could this be problematic?

When I use sol(t) it gives me the values of u (in this case my bioS) and time t but I cannot get values that are not part of u.

6

you’re pushing result to the global, but that step might get rejected so result isn’t necessarily a real result, and the time points aren’t necessarily unique either.

You can write a function of t that uses sol(t) to compute any values that can be computed from the solution.

function fsol(t) 
  u = sol(t)
  [u[2],-u[1],u[1]*u[3]
end
7

Thank you for your reply. How do I specify that I want a range of t values (from 0.0 to 50000.0? I can only specify one value for t like so

function fdemand_sol(t)
    bioS = demand_sol(t)
    [1.0/(clearance_max*bioS[target_species]^β_V/(b0V^β_V+bioS[target_species]^β_V)+clearance_min)]
end

clearance_demand = fdemand_sol(1000.0)

Edit: The calculated values do not change if I change t. It also does not change if I specify a different sol.
What am I doing incorrectly?

function fdemand_sol(t)
    bioS = demand_sol(t)
    [clearance_V, catch_F, fish_price, bioS_vessel, bioS_target]
end
clearance_demand = fdemand_sol(5000.0)

function fstatic_sol(t)
    bioS = static_sol(t) #this is a different differential equation
    [clearance_V, catch_F, fish_price, bioS_vessel, bioS_target]
end
clearance_static = fstatic_sol(5000.0)
8

You’re not using the output of sol(t)…

9

I’m afraid I don’t understand then.

My thought process is:
sol() is the differential equation that I want to get values from
t is the time interval I want to use
I need to run the differential equation (sol()), then create a function (fsol()) that will use the outputs of sol() to calculate the values I am interested in at t.

I don’t understand what I am doing incorrectly.

10

Yes, but you don’t use bioS in your calculations.

11

I use bioS here

function fdemand_sol(t)
    bioS = demand_sol(t)
    [1.0/(clearance_max*bioS[target_species]^β_V/(b0V^β_V+bioS[target_species]^β_V)+clearance_min)]
end

Is this incorrect?

12

That’s fine, but that’s not what you’re showing in:

function fdemand_sol(t)
    bioS = demand_sol(t)
    [clearance_V, catch_F, fish_price, bioS_vessel, bioS_target]
end
clearance_demand = fdemand_sol(5000.0)

function fstatic_sol(t)
    bioS = static_sol(t) #this is a different differential equation
    [clearance_V, catch_F, fish_price, bioS_vessel, bioS_target]
end
clearance_static = fstatic_sol(5000.0)
13

Okay, I have written it properly now. One more question - how to specify the time range I want (0.0 - 50000.0) ? I’ve tried two ways and neither works.

function fdemand_sol(tspan)
    bioS = demand_sol(tspan)
    clearance_V = bioS[num_nutrient+num_plant+num_animal+num_vessel_types+1]
    catch_F = (catch_max*bioS[target_species]/(catch_max/clearance_V+bioS[target_species]))*bioS[(num_nutrient+num_plant+num_animal+num_vessel_types)]
    [clearance_V, catch_F]
end

tmax = 50000.0
tspan = (0.0,tmax)

clearance_demand = fdemand_sol(tspan)
clearance_demand = fdemand_sol(0.0-50000.0)
14

0.0:dt:tmax and apply it at each point. There’s nothing special going on here.

15

It got confusing really fast from errors of dimension mismatch. Here is a work around for when t is a single or range value.

function fdemand_sol(t)
    if length(t) > 1
        bioSall = demand_sol(t) # this is if you give a range
    else
        bioSall = [demand_sol(t)] # this is if you give a single t value
    end
    
    clearance_V = Array{Float64}(undef, length(bioSall))
    catch_F = Array{Float64}(undef, length(bioSall))

    for i in 1:length(bioSall)
        bioS = bioSall[i]

        clearance_V[i] = bioS[(num_nutrient+num_plant+num_animal+num_vessel_types+1)]
        catch_F[i] = (catch_max*bioS[target_species]/(catch_max/clearance_V[i]+bioS[target_species]))*bioS[(num_nutrient+num_plant+num_animal+num_vessel_types)]
    end
    hcat(clearance_V, catch_F)
end

clearance_demand = fdemand_sol(0.0:1.0:tmax) # range
clearance_demand = fdemand_sol(100.0) # single value
16

It’s a function. Broadcast it or use a map.

17

How do I do it correctly? If I do this I get

Vector{Array{Float64,1}} with 2 elements
Vector{Float64} with 84 elements
Vector{Float64} with 84 elements

instead of 50001×2 Array{Float64,2}

function fdemand_sol(tspan)
    bioS = demand_sol(tspan)
    clearance_V = @. bioS[num_nutrient+num_plant+num_animal+num_vessel_types+1]
    catch_F = @. (catch_max*bioS[target_species]/(catch_max/clearance_V+bioS[target_species]))*bioS[(num_nutrient+num_plant+num_animal+num_vessel_types)]
    [clearance_V, catch_F]
end

tmax = 50000.0
clearance_demand = fdemand_sol(0.0:1.0:tmax)
18

map(fdemand_sol,0.0:1.0:tmax)