'Sol' not defined upon opening Julia?

New to Julia
1

Hi. I have this code I am trying to make for graphing the trajectory of a thrown baton. The script involves a mass matrix (of the baton) being multiplied by dydt to show its trajectory.
My problem comes from the use of the name of the mass matrix, ‘M’ and the name of the solution, ‘sol’. Every time I close and re-open Julia it suddenly doesn’t recognise ‘M’ or ‘sol.t’ or ‘sol.u’ but once I figure it out it doesn’t have a problem identifying them.
Here is the code:


P = Dict(:m1 => 0.1, :m2 => 0.1, :L => 1, :g => 9.81)

params = (0.1, 0.1, 1.0, 9.81)

function dydt(du,u,p,t)
  m1, m2, L, g = p
  du[1] = u[2]
  du[2] = m2*L*u[6]^2*cos(u[5])
  du[3] = u[4]
  du[4] = m2*L*u[6]^2*sin(u[5])-(m1 + m2)*g
  du[5] = u[6]
  du[6] = -g*L*cos(u[5]) 
  nothing # This is to specify the default jacobian
end

function mass(t, u, P)
    M = [1.0 0 0 0 0 0;
         0 P[:m1] + P[:m2] 0 0 0 -P[:m2]*P[:L]*sin(sol.u[1][5]);
        0 0 1 0 0 0;
         0 0 0 P[:m1]+P[:m2] 0 P[:m2]*P[:L]*cos(sol.u[1][5]);
         0 0 0 0 1 0;
         0 -P[:L]*sin(u[1][5]) 0 P[:L]*cos(u[1][5]) 0 P[:L]^2]
   return M
end

 Smotion = ODEFunction(dydt,mass_matrix=M)
  prob_mm = ODEProblem(Smotion,[6.0,0.0,0.0,0.0,0.0,0.0],(0,35),(0.04,3e7,1))
  sol = solve(prob_mm,Rodas5(),reltol=1e-28,abstol=1e-8)
2

Can you clarify what you mean by this? Every time you open a new instance of Julia you’re in a fresh environment and won’t have access to old workspace results without re-calculating them.

1 Like
3

Yeah, so see at the bottom it says Smotion = ODEFunction(dydt,mass_matrix=M). And ‘M’ is defined above it. When I open this file in new successive times, it doesn’t recognise ‘M’ for some reason.

4

If you wrap the code in this Discourse post in backticks it will display as a code block, which will make it easier for us to read and understand.

It looks like you define M inside of a function. In that case you’re returning the contents of M, not the name itself.

Are you supposed to pass in a matrix or the function that returns the matrix?

5

Also, that seems like an incredibly tight relative tolerance setting. If I remember correctly, Float64 types are fundamentally limited to more like 1e-16 precision. Common values are more like 1e-8.

6

Thanks Mike. I’ve put the code into backticks, and changed the reltol to 16. I took the ‘M’ out of its function structure but now it is saying ‘u’ is not defined. Is this because ‘u’ is being defined at the end of the script inside prob_mm?

7

This is what it looks like now:


P = Dict(:m1 => 0.1, :m2 => 0.1, :L => 1, :g => 9.81)
params = (0.1, 0.1, 1.0, 9.81)

function dydt(du, u, p, t)
    m1, m2, L, g = params
    du[1] = u[2]
    du[2] = m2 * L * u[6]^2 * cos(u[5])
    du[3] = u[4]
    du[4] = m2 * L * u[6]^2 * sin(u[5]) - (m1 + m2) * g
    du[5] = u[6]
    du[6] = -g * L * cos(u[5])
    nothing
end

M = [1.0 0 0 0 0 0;
     0 P[:m1] + P[:m2] 0 0 0 -P[:m2] * P[:L] * sin(u[5]);
     0 0 1.0 0 0 0;
     0 0 0 P[:m1] + P[:m2] 0 P[:m2] * P[:L] * cos(u[5]);
     0 0 0 0 1.0 0;
     0 -P[:L] * sin(u[5]) 0 P[:L] * cos(u[5]) 0 P[:L]^2]
   

Smotion = ODEFunction(dydt, mass_matrix=M)
prob_mm = ODEProblem(Smotion, [0.0, 4.0, L, 20.0, -π/2, 2.0], (0.0, 35), (0.04, 3e7, 1))
sol = solve(prob_mm, Rodas5(), reltol=1e-16, abstol=1e-8)
8

Was you earlier code based off of some prior/demo code? I’m not very experienced with DifferentialEquations.jl so it’s not super obvious to me what ODEFunction expects to be passed in. Your earlier code made sense if it’s expecting a function that can be computed based on the current conditions, but in that case you should be passing in the function name (mass).

9

Yes, I’m new to Julia so I was using the below code by Chris Rackauckas to base this similar problem off of (there is a multiplication of matrices in both). The main thing I want to see is the successful matrix multplication.

function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end
M = [1. 0  0
     0  1. 0
     0  0  0]
f = ODEFunction(rober,mass_matrix=M)
prob_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob_mm,Rodas5(),reltol=1e-8,abstol=1e-8)

plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))
10

It may be worth noting that I did get it to work using this structure a few times, it was when I restarted that I started getting error messages about variables it couldn’t recognise.

11

It seems like you’re trying to use a state-dependent mass matrix. In that case, you have to make it a SciMLOperator. This is mentioned in the ODEFunction documentation: ODE Problems · DifferentialEquations.jl (specifically the line “must be an AbstractArray or AbstractSciMLOperator”). See Home · SciMLOperators.jl to get started.

12

Does the code I based it off not use a state-dependent mass matrix?

13

No, as you can see, it just defines a constant M that has no dependence on u.

That said, SciMLOperators.jl has an unresolved issue where it conflates the vector it’s supposed to act on with the solver state it could depend on. I don’t know if it’s possible to accomplish what (I think) you’re attempting before that is fixed. It’s mentioned in several issues, such as Separate RHS and state vectors for `FunctionOperator` · Issue #57 · SciML/SciMLOperators.jl · GitHub, Overload operator application for `ComposedOperator` · Issue #159 · SciML/SciMLOperators.jl · GitHub, and How do I use this package to pass analytical jvp? · Issue #223 · SciML/SciMLOperators.jl · GitHub.

14

Actually, never mind, the bug only affects FunctionOperator and the documentation (but not implementation) of certain other operators. You can use MatrixOperator to accomplish what (I think) you want.

15

Here’s my best attempt to do what I think you’re trying to do. Note that I had to do some interpretation, such as making the mass matrix dependent on the current state u, not the initial state sol.u[1] of some previous solution.

using LinearAlgebra, SciMLOperators, OrdinaryDiffEq

function dydt(du,u,p,t)
  (; m1, m2, L, g) = p
  du[1] = u[2]
  du[2] = m2 * L * (u[6]^2) * cos(u[5])
  du[3] = u[4]
  du[4] = m2 * L * (u[6]^2) * sin(u[5]) - (m1 + m2) * g
  du[5] = u[6]
  du[6] = -g * L * cos(u[5]) 
  return nothing
end

# Function to update an existing mass matrix using the state u,
# assuming constant entries are correct
function updateM!(M, u, p, t)
    (; m1, m2, L, g) = p
    M[2, 6] = -m2 * L * sin(u[5])
    M[6, 2] = -L * sin(u[5])
    M[4, 6] = m2 * L * cos(u[5])
    M[6, 4] = L * cos(u[5])
    return nothing
end

# Function to create the mass matrix for the state u
function newM(_, u, p, t)
    (; m1, m2, L, g) = p
    # Create a 6x6 identity matrix and fill in the constant non-identity parts
    M = Matrix(1.0 * I(6))
    M[2, 2] = M[4, 4] = m1 + m2
    M[6, 6] = L^2
    # Use updateM! to fill in the state-dependent parts
    updateM!(M, u, p, t)
    return M
end

# Function to create a MatrixOperator representing the mass matrix
function mass_matrix(u, p)
    M = newM(nothing, u, p, 0.0)
    return MatrixOperator(M; update_func = newM, update_func! = updateM!)
end

params = (; m1=0.1, m2=0.1, L=1.0, g=9.81)
u0 = [0.0, 4.0, params.L, 20.0, -π / 2, 2.0]
smotion = ODEFunction(dydt; mass_matrix=mass_matrix(u0, params))
prob_mm = ODEProblem(smotion, u0, (0, 35), params)
sol = solve(prob_mm, Rodas5())

I’m not sure whether this accomplishes what you want, but it runs!

16

Thank you so much! That’s very helpful. So it was expecting a rigid mass matrix instead of a state-dependent one. Thanks Daniel.