DifferentialEquations.jl look at the ODEsystem

General Usage
1

I am solving a large system of ODEs with the DifferentialEquation.jl package. Since I am using a number of for loops to set up the RHS to the du, so I would like to see the RHS matrix before solving. Is there a way I can do that?

Lets say I have the following

tspan=[0.0,1.0]
function flux!(du,u,p,t)
for i in 1:n
     dummy=0.0
     for j in 1:n
         dummy=dummy+u[j]
     end
   du[i]= max(0,dummy)*u[i]+min(0,dummy)*u[i+1]
end

prob = ODEProblem(flux!, u0, tspan)
sol = solve(prob)

So basically I am trying to put a convolution like function on the right hand side, and trying to see if the correct ODE system is generated or not.

2

I’m a bit confused. Is this ModelingToolkit ODESystem or direct DifferentialEquations.jl?

3

This is DifferentialEquations.jl.

4

Then it’s just Julia code so you instrument it with @show. What did you try?

function flux!(du,u,p,t)
for i in 1:n
     dummy=0.0
     for j in 1:n
         dummy=dummy+u[j]
     end
   du[i]= max(0,dummy)*u[i]+min(0,dummy)*u[i+1]
end
@show du
end

you were missing an end btw

5

One debugging trick here is, because it’s linear, the matrix can be found by calculating the Jacobian. Just do ForwardDiff.jacobian on f and check the matrix.

6

Thanks for the answer. I now can print the du. Can I also, acces it as a variable(Matrix)?

7

du isn’t a matrix. But if you’re looking for the operator, in this case it’s equivalent to the Jacobian, so just calculate the Jacobian.

8

May be I am not being clear. So in the system

d_tu=A \cdot u

I want the matrix A, so that I can compute the coefficients by hand and see that they are indeed what they should be. I am really sorry, I am new to Julia and I don’t understand macros.

9

Yes, have you tried calculating the Jacobian?

10

Thanks for all the help, the issue was when I was writing the flux as a sum on a new line for each component, only the 1st component was being added somehow. So when I changed the coefficient for the 1st component to zero to see how the others affected the solution, the RHS essentially became zero, as the rest of the components were not being added. That’s why I wanted to look at the RHS.

It took me a while to realize that whatever I sum to the flux after putting a newline is not being added to the RHS.