# Variable rate jumps with DifferentialEquations.jl and matrices

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 + u
affect!(integrator) = (integrator.u = 0.25; integrator.u = 0.5;
integrator.u = 0.75; integrator.u = 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 = u; du = u` for the components of the matrix, but I would prefer to have just one equation with matrices.

It’s a bug which will be fixed by https://github.com/SciML/DiffEqJump.jl/pull/117 . I’ll release the fix tomorrow.

Thank you!

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 + u
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 = complex(1.0); integrator.u = complex(2.0);
integrator.u = complex(3.0); integrator.u = 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?

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

Yeah, I think that’s the issue, Maybe norm(u + u) 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

My bad! The rate should be real. But changing to ` norm(u + u) ` 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 + u)
affect!(integrator) = (integrator.u = complex(1.0); integrator.u = complex(2.0);
integrator.u = complex(3.0); integrator.u = 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 + u
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()`.