Variable rate jumps with DifferentialEquations.jl and matrices

General Usage
1

I’ve been going through some examples on differential equations with variable jumps and I did this example for myself and got it working.

function f(du, u, p, t)
    du .= u
end

prob = ODEProblem(f, [1.0 2.0; 3.0 4.0], (0.0, 1.0))
rate(u,p,t) = u[1] + u[2]
affect!(integrator) = (integrator.u[1] = 0.25; integrator.u[2] = 0.5;
                       integrator.u[3] = 0.75; integrator.u[4] = 1)
jump = ConstantRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob)
plot(sol)

Notice here that the rate is not constant and therefore I should use VariableRateJump() function instead of ConstantRateJump(). But if I replace it in the code it doesn’t work (Dimension mismatch). So what is the correct way to do this for a variable jump rate, where u is not a float but a matrix instead? I understand I could do du[1] = u[1]; du[2] = u[2] for the components of the matrix, but I would prefer to have just one equation with matrices.

2

It’s a bug which will be fixed by fix variable rate jump broadcast sizing by ChrisRackauckas · Pull Request #117 · SciML/JumpProcesses.jl · GitHub . I’ll release the fix tomorrow.

3

Thank you!

4

I updated the package and now that example works! But I ran into more problems. This example seems to run just fine:

function f(du, u, p, t)
    du .= u
end

prob = ODEProblem(f, Complex{Float64}[1.0 2.0; 3.0 4.0], (0.0, 1.0))
rate(u,p,t) = u[1] + u[2]
affect!(integrator) = (integrator.u = Complex{Float64}[1.0 2.0; 3.0 4.0])
jump = ConstantRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob)

but changing from ConstantRateJump() to VariableRateJump() breaks it. There are couple things going on here. I did some testing and seems like the VariableRateJump() problem breaks if the affect is defined directly as a matrix:

affect!(integrator) = (integrator.u = [1.0 2.0; 3.0 4.0])

And also it breaks if the affect is complex:

affect!(integrator) = (integrator.u[1] = complex(1.0); integrator.u[2] = complex(2.0);
                       integrator.u[3] = complex(3.0); integrator.u[4] = complex(4.0))

I hope you could look into these once you have time, thanks for the help! Also, if I run into new bugs what is the best method for reporting them?

5

What does it mean to have a complex valued rate in a jump process? I’m not sure any of the jump methods we have will work with complex rates…

1 Like
6

Yeah, I think that’s the issue, Maybe norm(u[1] + u[2]) or something would work, but I am not sure what it means to have a complex value as the rate for a jump process.

1 Like
7

My bad! The rate should be real. But changing to norm(u[1] + u[2]) did not seem to fix the issues. So to be clear both of these do not work

function f(du, u, p, t)
    du .= u
end

prob = ODEProblem(f, Complex{Float64}[1.0 2.0; 3.0 4.0], (0.0, 1.0))
rate(u,p,t) = norm(u[1] + u[2])
affect!(integrator) = (integrator.u[1] = complex(1.0); integrator.u[2] = complex(2.0);
                       integrator.u[3] = complex(3.0); integrator.u[4] = complex(4.0))
jump = VariableRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob)

and

function f(du, u, p, t)
    du .= u
end

prob = ODEProblem(f, [1.0 2.0; 3.0 4.0], (0.0, 1.0))
rate(u,p,t) = u[1] + u[2]
affect!(integrator) = (integrator.u = [1.0 2.0; 3.0 4.0])
jump = VariableRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob)

but they do work for ConstantRateJump().