State Dependent Mass Matrix

General Usage
1

Hello,

I’ve looked at the documentation for Handling Mass Matrices for DifferentialEquations.jl. In the example was a singular matrix. I would like to know how to implement a state dependent mass matrix, where the matrix M in equation
Mu’ =f(u,p,t)
has the dependent variables in the matrix itself.

The image below is the equation I’m trying to replicate:

Capture
Capture811×85 7.77 KB

I’m trying to integrate from Matlab to Julia for this type of problem. I thought replacing the mass_matrix option with a function would work but it produces the method error MethodError: no method matching issingular(::typeof(function goes here))

An example or helpful tip would be appreciated. The link relating to mass matrices is https://docs.sciml.ai/latest/tutorials/advanced_ode_example/#Handling-Mass-Matrices-1 Thank you!

2

I think this is not implemented right now ( but will be, I am sure ).

I planned on doing this after a discussion with @ChrisRackauckas and @YingboMa on Slack, but never found the time.

If I encounter such problems right now, I tend to use something along the lines of

function dudt(u, p, t)
    m = M(u,p,t)
    f = F(u,p,t)
    inv(m)*f
end

or solve it manually, if possible.

3

It’s implemented. Examples:

We just need to clean up our documentation of DiffEqOperators to make it more clear how to define dependent DiffEqOperators, since that’s all it is (the mass matrix can be any DiffEqOperator)

4

Ah, I did not know that :+1:

5

Thank you for that. I’m still new to Julia, could you link me to the equation your code snippet is solving to help visualize it for me? I can try to figure out the syntax from there.

6

I am really really slow with making matrices in LaTeX, so does this simplification help?

function update_func(A,u,p,t)
    A[1,1] = cos(t)
    A[2,1] = sin(t)*u[1]
    A[3,1] = t^2
    A[1,2] = cos(t)*sin(t)
    A[2,2] = cos(t)^2 + u[3]
    A[3,2] = sin(t)*u[2]
    A[1,3] = sin(t)
    A[2,3] = t^2
    A[3,3] = t*cos(t) + 1
end
M = DiffEqArrayOperator(ones(3,3),update_func=update_func)
u0 = ones(3)
tspan = (0.0, 1.0)
f(du,u,p,t) = (@. du = u + t; nothing)
prob = ODEProblem(ODEFunction(f, mass_matrix=M), u0, tspan)
sol = solve(prob,RadauIIA5())

The mass matrix is just the A defined in the update_func. I’ll trade you for LaTeX and together we can get this in the docs :slight_smile:

2 Likes
7

Has the solution to this problem updated? DifferentialEquations.jl documentation says that state dependent mass matrices are “not directly supported” and there is no documentation in DiffEqOperators.jl either. Would the current implementation support Complex variables? @ChrisRackauckas @Julius_Martensen

8

State-dependent mass matrices are inherently pretty unstable, so I would recommend a DAE formulation instead.

1 Like
9

Sorry for making revive this thread.

I have a question of curiosity. What are the reasons that may prevent from adding a mass matrix as a parameter of a solver ? Is it for a technical reason ?

Currently, I include the mass matrix in the definition of the problem using du = M\u. I am happy with that, but I notice that Matlab allows for the definition of a mass matrix via odeset for a couple of its solvers. I guess that in this case some numerical strategies are implemented to solve the problem efficiently (but I may be wrong).

Thank you for helping me to fill my ingnorance on this topic.

10

Your setup does not work when M is singular (i.e. a DAE).

Also you can specialize the nonlinear solver so the M part is effectively free.

1 Like