I wouldn’t use the phrase “evaluate a symbolics variable”. At this stage, the solutions are available and they are numeric. But anyway, maybe your phrase is ok.
It is important to understand the structure of the solution which I have called sol. The solution contains two fields t and u – field t holds the time instances where the solution has been found, while field u holds the computed values of the unknowns in these time instances.
In the case of Backstress:
> sol.t' # I use transpose `'` symbol to save vertical space
1×32 adjoint(::Vector{Float64}) with eltype Float64:
0.0 0.762987 8.39286 84.6916 … 44162.9 45314.9 47666.0 50000.0
> sol.u'
1×32 adjoint(::Vector{Vector{Float64}}) with eltype LinearAlgebra.Adjoint{Float64, Vector{Float64}}:
[1.0;;] [1.00003;;] [1.00033;;] [1.00329;;] … [70.7544;;] [70.751;;]
Observe that we have found the solution in 32 time instances. Also, observe that sol.u is a vector of vectors. Each element in the vector holds the solution of all unknowns in the corresponding time instance. So as an example, vector element sol.u[3] holds the solution of unknowns corresponding to time sol.t[3]. Two specific points:
- Why is there only one unknown with solution in
sol.u?
- The shape of
sol.u is not in a form that can be plotted.
sol.u and unknowns.
Field u of the solution only holds the solution to variables formally defined as unknowns. To find out which these are, and their order of storage in sol.u, you need to query:
> unknowns(sys)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
σ_T(t)
What about other “unknowns”? These have been stripped off the model and made into observed variables; the observed variables are found by:
> observed(sys) |> println
Equation
[ϵ°(t) ~ A°*Exp*
sinh((-σ_T(t) + σ - σ_d) / σ_ref)*(σ_T(t)^2)]
I’ll get back to how to compute this observed variable – which is what you ask about, I guess.
- The shape of
sol.u and how to make it plottable…
Perhaps the most efficient way to make sol.u plottable is:
> reduce(vcat,sol.u)' # Again: transpose to save space
1×32 adjoint(::Vector{Float64}) with eltype Float64:
1.0 1.00003 1.00033 1.00329 … 70.7375 70.7376 70.7544 70.751
Now you have the solution in a vector of elements, and not in a vector of vectors. The same should work if sol.u has multiple unknowns.
Plotting the result can be done by:
plot(sol.t, reduce(vcat,sol.u))
which produces:
Notice how jagged the solution is compared to doing:
plot(sol)
which produces
The reason for the jaggedness in the first method, plot(sol.t, reduce(vcat,sol.u)) , is that this draws straight lines between the datapoints of sol.t[j] and sol.u[j].
One of the fancy things with Julia solvers is that sol not only holds the solution in datapoints, i.e., sol.t, etc. The solution also includes interpolation functions. Interpolation in a certain time point t0 is given by syntax sol(t0), e.g.:
> sol(4e4)
1-element Vector{Float64}:
70.74021605315112
This is the solution σ_T(t=4e4).
Final question: how can you specify which variable you want to know the value of? And how can you find the value of observed variables? Simple:
> sol(4e4, idxs=σ_T)
70.74021605315112
> sol(4e4, idxs=ϵ°)
3.6413814395589957e-7
> sol(4e4, idxs=[σ_T, ϵ°])
2-element Vector{Float64}:
70.74021605315112
3.6413814395589957e-7
What about plotting observed variables (and unknowns) with interpolation?
plot(sol, idxs=σ_T)
produces.
It is also possible to scale time variable and results:
plot(sol, idxs=(t/60,[ϵ°*1e8,σ_T]), label=["epsilon*1e8" "sigma"], xlabel="time [min]", ylabel="")
produces:
If you want a parametric plot of the two dependent variables with time as parameter:
plot(sol, idxs=(σ_T, ϵ°*1e8))
produces:
