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

**URL:** <https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183>\
**Category:** New to Julia\
**Tags:** plotting\
**Created:** [May 26, 2020, 10:32am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183 "2020-05-26T10:32:44Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 26, 2020, 10:32am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/1 "2020-05-26T10:32:44Z")

</div>

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.

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 28, 2020, 2:00pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/2 "2020-05-28T14:00:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 28, 2020, 3:02pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/3 "2020-05-28T15:02:07Z")

</div>

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

```julia
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](https://docs.sciml.ai/stable/features/callback_library/#saving_callback-1)) together with help from @jonniedie on Gitter.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 28, 2020, 3:34pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/4 "2020-05-28T15:34:46Z")

</div>

> [@HelgavonLichtenstein](#):
>
> What is the post solution interpolation? Could you give me an example?

`sol(t)`.

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

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [May 29, 2020, 8:20am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/5 "2020-05-29T08:20:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 30, 2020, 5:05am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/6 "2020-05-30T05:05:21Z")

</div>

> [@HelgavonLichtenstein](#):
>
> Could you elaborate on what is assuming time is monotonic?

> [@HelgavonLichtenstein](#):
>
> ```julia
> results = [catch_F fish_price] 
> global result2 = vcat(result2,result)
> 
> ```

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.

> [@HelgavonLichtenstein](#):
>
> 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 can write a function of `t` that uses `sol(t)` to compute any values that can be computed from the solution.

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

```

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 2, 2020, 11:53am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/7 "2020-06-02T11:53:05Z")

</div>

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

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2020, 12:56pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/8 "2020-06-02T12:56:28Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 2, 2020, 1:42pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/9 "2020-06-02T13:42:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2020, 1:48pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/10 "2020-06-02T13:48:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 2, 2020, 1:56pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/11 "2020-06-02T13:56:55Z")

</div>

I use `bioS` here

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2020, 3:01pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/12 "2020-06-02T15:01:48Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 5, 2020, 11:32am UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/13 "2020-06-05T11:32:42Z")

</div>

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.

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 5, 2020, 12:25pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/14 "2020-06-05T12:25:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 5, 2020, 2:40pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/15 "2020-06-05T14:40:40Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 5, 2020, 2:49pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/16 "2020-06-05T14:49:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HelgavonLichtenstein](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgavonlichtenstein/32/5005_2.png) [@HelgavonLichtenstein](https://discourse.julialang.org/u/HelgavonLichtenstein)\
**Post date:** [June 5, 2020, 3:06pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/17 "2020-06-05T15:06:47Z")

</div>

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

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

```

instead of `50001×2 Array{Float64,2}`

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 5, 2020, 3:46pm UTC](https://discourse.julialang.org/t/how-to-save-values-from-a-differential-equation-not-part-of-initial-conditions-u-and-then-graph-them/40183/18 "2020-06-05T15:46:38Z")

</div>

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