# MATLAB is faster for my large system of ODEs?

**URL:** <https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610>\
**Category:** Numerics\
**Tags:** matlab, ode, speed-optimization, solver\
**Created:** [November 13, 2024, 11:06pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610 "2024-11-13T23:06:02Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![martru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martru/32/213400_2.png) [@martru](https://discourse.julialang.org/u/martru)\
**Post date:** [November 13, 2024, 11:06pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/1 "2024-11-13T23:06:02Z")

</div>

Excuse the deliberately provocative title 🙂  
I have a very large system of ODEs of the form:

dot{x} = Ax + Bu

A is a sparse matrix 5505x5505. B\*u is a 5505 element vector. Using ode45 in matlab, it solves in about 59 minutes. In Julia I have the problem set up like so:

```julia
function sys!(dx,x,p,t)
    dx.= A*x + Bu
    nothing
end

x0 = zeros(5505)
prob = ODEProblem(sys!, x0, (0.0,5.0))
@btime sol = solve(prob, CVODE_BDF(), save_everystep = false)

```

I’ve tried CVODE\_BDF and QNDF but I cut them both off after neither solved within an hour. I don’t have a lot of experience with solving large systems, so any tips for accelerating solving would be much appreciated.

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 13, 2024, 11:21pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/2 "2024-11-13T23:21:59Z")

</div>

The rough equivalent of ode45 would be `Tsit5` You are trying to use implicit solvers which are going to be a lot slower if your matrices aren’t stiff. Also, you probably want `bu` to be part of `p` rather than a global variable, i.e.

```julia
function sys!(dx,x,p,t)
    Bu = p[1]
    dx.= A*x + Bu
    nothing
end

x0 = zeros(5505)
prob = ODEProblem(sys!, x0, (0.0,5.0), (Bu,))
@btime sol = solve(prob, Tsit5(), save_everystep = false)
```

---

<div class="post-metadata">

**Author:** ![martru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martru/32/213400_2.png) [@martru](https://discourse.julialang.org/u/martru)\
**Post date:** [November 13, 2024, 11:28pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/3 "2024-11-13T23:28:28Z")

</div>

I’ll give that a shot, thanks! Generally speaking, does Tsit5 tend to be any faster than ode45? Are there any other solvers you’d recommend for this case?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 13, 2024, 11:36pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/4 "2024-11-13T23:36:20Z")

</div>

> [@martru](#):
>
> `dx.= A*x + Bu`

Also, the above isn’t really in-place, it allocates two temporary arrays and then copies the result into dx.

I am not that familiar with sparse linear algebra, but look at the `mul!` function for in-place, non-allocating operations.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [November 13, 2024, 11:52pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/5 "2024-11-13T23:52:33Z")

</div>

This post was temporarily hidden by the community for possibly being off-topic, unfocused, inappropriate, or spammy.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 13, 2024, 11:52pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/6 "2024-11-13T23:52:46Z")

</div>

In general, `Tsit5` is very similar to ode45 (also known as `DP5`), but is generally ~20% faster due to using more optimized coefficients. There are tons of methods that are potentially better, but knowing which to look into requires knowing more about the properties of `A`

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [November 14, 2024, 6:48am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/7 "2024-11-14T06:48:25Z")

</div>

Shouldn’t `A` also be passed within the parameters, rather than being global?

BTW, rather than wait an hour to test a simulation from t = 0 to 5.0, a shorter time span is probably sufficient to get an idea of relative speed.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 14, 2024, 10:13am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/8 "2024-11-14T10:13:52Z")

</div>

> [@martru](#):
>
> ```julia
> function sys!(dx,x,p,t)
> dx.= A*x + Bu
> nothing
> end
> 
> ```

```julia
function sys!(dx,x,p,t)
    mul!(dx, A, x)
    mul!(dx, B, u, 1, 1)
    nothing
end

```

this is the non-allocating version (depending a bit on where `u` is coming from).

If `u` is constant or piecewise constant, you can solve this differential equation by computing a single matrix exponential of `exp([A B; 0 0])` and iterate the corresponding system

\begin{align} x(k+1) = A\_d x(k) + B\_d u(k) \end{align}

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [November 14, 2024, 12:41pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/9 "2024-11-14T12:41:03Z")

</div>

@martru, for convenience of other readers of this forum, you may want to consider reflecting the essence of your problem better in the title. Perhaps changing it to something like “Matlab faster in solving a huge linear ODE?”

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [November 14, 2024, 1:00pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/10 "2024-11-14T13:00:59Z")

</div>

> [@martru](#):
>
> Excuse the deliberately provocative title

Best way to get help 😅 but I edited to be a bit more clear what it’s about to future searchers

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [November 14, 2024, 1:02pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/11 "2024-11-14T13:02:01Z")

</div>

Oh, didn’t see this suggestion before I acted 😉. I just got this power to edit titles and I might be drunk with it

---

<div class="post-metadata">

**Author:** ![martru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martru/32/213400_2.png) [@martru](https://discourse.julialang.org/u/martru)\
**Post date:** [November 14, 2024, 6:43pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/12 "2024-11-14T18:43:43Z")

</div>

u is indeed constant. I’ll try implementing this!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 14, 2024, 7:41pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/13 "2024-11-14T19:41:20Z")

</div>

It’s implemented here

> <https://github.com/JuliaControl/ControlSystems.jl/blob/46086894dc85dcb736776c2c1202954e7e90be70/lib/ControlSystemsBase/src/discrete.jl#L51>

You could either use this functionality directly, or just find the definition of `Ad, Bd` there

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [November 14, 2024, 8:39pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/14 "2024-11-14T20:39:11Z")

</div>

> [@baggepinnen](#):
>
> `exp!`

I was just wondering where does this `exp!` function come from and what does it do?

---

<div class="post-metadata">

**Author:** ![martru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martru/32/213400_2.png) [@martru](https://discourse.julialang.org/u/martru)\
**Post date:** [November 14, 2024, 10:02pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/15 "2024-11-14T22:02:28Z")

</div>

[This](https://web.mit.edu/2.14/www/Handouts/StateSpaceResponse.pdf) might clear it up.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [November 14, 2024, 10:08pm UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/16 "2024-11-14T22:08:15Z")

</div>

Clear up what?

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [November 15, 2024, 12:00am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/17 "2024-11-15T00:00:02Z")

</div>

> [@zdenek\_hurak](#):
>
> Clear up what?

I believe OP was answering your question about `exp!`, with a reference about using matrix exponentials to solve state-space systems. Perhaps your question was more about the `!` than the `exp`. If so, `exp!` is in-place matrix exponential, specific to LinearAlgebra:

> <https://github.com/JuliaLang/julia/blob/100e305b9fd0fa6f05795bdad7f46838eff36f97/stdlib/LinearAlgebra/src/dense.jl#L713>

EDIT: In the example above `exp!` is used in a way I never dared. A matrix is composed within the in-place function call and then assigned, e.g. `M = exp!([a b; c d])`. My mindset was always that you have to allocate an array before calling in place, otherwise how do you know what the result is? But by assigning to LHS you know the result. I know this is trivial, but this pattern wasn’t obvious to me, and it’s great!

` M = exp!([A*Ts B*Ts; zeros(nu, nx + nu)])`

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [November 15, 2024, 12:42am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/18 "2024-11-15T00:42:37Z")

</div>

Indeed, the `!` symbol in `exp!` was what made me wondering. I am aware of the general habit of naming functions that modify their argument in this way, but it is just that I was unable to find any description of the function `exp!` – neither in the [online documentation](https://docs.julialang.org/en/v1/), nor through the `?` help in REPL.

---

<div class="post-metadata">

**Author:** ![martru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martru/32/213400_2.png) [@martru](https://discourse.julialang.org/u/martru)\
**Post date:** [November 15, 2024, 1:31am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/19 "2024-11-15T01:31:13Z")

</div>

This essentially sped things up by a factor of 170! Thanks!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 15, 2024, 8:06am UTC](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610/20 "2024-11-15T08:06:28Z")

</div>

That’s very cool 🙂

If you find yourself in need of speeding it up even further, you may be able to perform all matrix operations on the GPU. Fig 7 in [the ControlSystems.jl paper](https://github.com/JuliaControl/CDC2021/blob/master/CDC2021.pdf) indicates this being worthwhile above state dimensions of about 100. You may also use the trick employed in `ControlSystems.lsim` that `B*u` does not need to be iterated and can be computed in a single operation and added to the resulting state by broadcasting.

[Next page](https://discourse.julialang.org/t/matlab-is-faster-for-my-large-system-of-odes/122610.md?page=2)