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

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)
``````

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.

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.

`sol(t)`.

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

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`.

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
``````

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)
``````

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

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.

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

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?

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)
``````

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)
``````

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

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
``````

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

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)
``````

`map(fdemand_sol,0.0:1.0:tmax)`