Multiple variable assignment and ODE solution extracting

Specific Domains Modelling & Simulations
1

Hi there!

I am quite new to Julia and have read manuals, tutorials, FAQs, and almost all question threads regarding the following, but have not yet found a solution.

I have a system of 17 ordinary differential equations that I can successfully solve with solve. However, I need the solutions of each function (the us) to be assinged each to a variable name (say v1, v2, v3, …, v17). So I am expecting them to be arrays of Nx1 (or 1xN).

Do you know a way to do this?

Thank you in advance!

2

Simple example:

function pendulum!(du,u,p,t)
    g = 9.81
    L = 2
    alpha = u[1]
    omega = u[2]
    du[1] = omega
    du[2] = -g/L*sin(alpha)
end
#
u0 = [pi/4,0.1]
tspan = (0.,5)
prob = ODEProblem(pendulum!,u0,tspan)
sol = solve(prob)

Here, sol is a fancy data structure that contains the time instances for the solution in field t, sol.t. It also contains the solution vector of variable 1 in sol[1,:], so I can plot the first variable as follows:

plot(sol.t,sol[1,:])

However, sol also contains an interpolation function, and by using function syntax, sol(0.5) contains the state vector at time = 0.5. A better way to do the plotting is the following:

plot(sol,vars=(0,1))

which (i) plots variable 1 as a function of variable 0 (= time), and (ii) uses interpolation so the solution looks smoother than if you use plot(sol.t,sol[1,:]).

You can extract the solutions as a matrix with variable 1 in row 1, variable 2 in row 2, etc.:

sol[1:end,:]

Here, column n corresponds to the time instance sol.t[n].

You can achieve the same by:

reduce(hcat,sol.(sol.t))

Here, sol.(sol.t)) computes the solution in the points given by time vector sol.t, and reduce(hcat,.) wraps the result into the same matrix as sol[1:end,:].

Finally, the solver decides at which points you want to compute the solution, i.e., the time instances in sol.t. If you, for some reason, want to compute the solution at other time instances, say at regular time intervals tsamp = tspan[1]:0.1:tspan[2] for my pendulum, you can create the solution matrix in those points by:

reduce(hcat,sol.(tsamp))
3

BLI,

Thank you for your reply. I liked the use you gave to reduce to build the matrix !

However, I would still need to know how to assign each one of the 17 solutions (say each row of the matrix) to a given variable name. So far the only solution I have is the brute force:

v1 = M[1,:];
v2 = M[2,:];
...
v17 = M[17,:]

with M being the Array that is the result of your reduce suggestion.

Any ideas?

Thank you once again.

4

Hm… can you use this? First, I create a dummy M matrix just for the illustration.

julia> M = rand(0:9,5,10)
5×10 Array{Int64,2}:
 4  3  5  5  4  0  1  0  5  6
 4  3  7  3  5  9  4  0  9  1
 1  9  6  3  6  4  2  6  1  8
 3  1  8  0  0  4  8  5  7  6
 2  6  4  1  0  7  2  4  3  5

[5\times 10 matrix of random values taken from the integers from 0, …, 9].

Next, I pick out vectors from the rows of M and given them names:

julia> for i in 1:5 
             variable_name = Symbol("v$i")
             vector = M[i,:]
             @eval $variable_name = $vector
        end
julia> v1
10-element Array{Int64,1}:
 4
 3
 5
 5
 4
 0
 1
 0
 5
 6

Is this what you want?

2 Likes
5

I can work with that! I will need to adapt it since my variables are not actually named vi… each one has a different name since they represent physical variables that are easier to track in named after what they represent.

But this is very good. Thank you!