# Python vs. Julia ODE Solver

**URL:** <https://discourse.julialang.org/t/python-vs-julia-ode-solver/56123>\
**Category:** Performance\
**Tags:** python\
**Created:** [February 26, 2021, 11:03pm UTC](https://discourse.julialang.org/t/python-vs-julia-ode-solver/56123 "2021-02-26T23:03:06Z")\
**Posts on this page:** 1\
**Showing post:** 22

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 27, 2021, 10:26am UTC](https://discourse.julialang.org/t/python-vs-julia-ode-solver/56123/22 "2021-02-27T10:26:50Z")

</div>

> [@qtoloza](#):
>
> That’s odd, right?

It is. But you are not running exactly the “my” version, because I have much less allocations than you.

> [@baggepinnen](#):
>
> ```julia
> function bar(x, y)
> x[2:5] .= y .+ y
> nothing
> end
> 
> ```

Just to be clear, it was not possible to do this in the first version because of the `rand(Float64,(N,1))` definition I mentioned above:

```julia
julia> v = rand(Float64,(5,1));

julia> y = rand(5);

julia> z = zeros(5);

julia> z[1:5] .= v .+ y
ERROR: DimensionMismatch("cannot broadcast array to have fewer dimensions")

```

The problem is that `v` here is a two-dimension array (even though the second dimension is 1):

```julia
julia> v
5×1 Array{Float64,2}:
 0.23460063971299383
 0.977599317555971
 0.11413962191900962
 0.9259610441668678
 0.0253857101609829

```

That is solved if `v` is a vector:

```julia
julia> v = rand(Float64,5);

julia> z[1:5] .= v .+ y
5-element view(::Array{Float64,1}, 1:5) with eltype Float64:
 1.7707818802990465
 1.5662031379445822
 1.1979284388950557
 0.7127751475367243
 1.1919315206952428

```

It might be a parallel discussion why/if the first broadcasting should work.

The final code I am running is here, It seems to be faster than what you are running there still.

> **Code**
>
> ```julia
> using OrdinaryDiffEq, BenchmarkTools, LSODA, Plots, Profile
> using LinearAlgebra
> 
> function forward_RHS(du, u, p, t)
> v_s = view(u, 1:p[2])
> u_s = view(u, p[2]+1:2*p[2])
> w_s = view(u, 2*p[2]+1:3*p[2])
> s = view(u, 3*p[2]+1:4*p[2])
> 
> p[26] .= t .* view(p[23],:,1)
> p[26] .= sin.(p[26])
> mul!(p[24],p[20],p[26])
> 
> du[1:p[2]] .= ((p[3].*(v_s .- p[4]).^2 .+ p[5]).*(v_s .- 1.0)
> .- p[12].*u_s.*(v_s .- p[14]) .- p[13].*w_s.*(v_s .- p[15])
> .+ p[22]
> .+ p[24]
> .+ mul!(p[25],p[21],s))./p[16]
> du[p[2]+1:2*p[2]] .= (exp.(p[6].*(v_s .- p[7])) .+ p[8] .- u_s)./p[17]
> du[2*p[2]+1:3*p[2]] .= (p[9].*(v_s .- p[10]).^2 .+ p[11] .- w_s)./p[18]
> du[3*p[2]+1:4*p[2]] .= (-1 .* s .+ (v_s .> 0.25).*(v_s .< 0.3).*(view(du,1:p[2]) .> 0.0).*view(du,1:p[2])./0.05)./p[19]
> 
> nothing
> 
> end
> 
> function runprob(N)
> 
> parameters = (2, N,
> -500, 0.2, 0.0, 100.0, 0.7, 0.0, 15.0, 0.24, 0.02, 1700.0, 8.0, 0.19, -0.5,
> 25.0, 6.0, 12.0, 10.0,
> 2.0*1.4 .* (rand(Float64,N,2) .- 0.5), #20 
> 0.5 .* (rand(Float64,N,N) .- 0.8), #21
> 0.14 .* (rand(Float64,N) .- 0.5), #22
> rand(Float64,2,4), #23
> zeros(N), #24
> zeros(N), #25 
> zeros(2)) #26
> 
> prob = ODEProblem(forward_RHS, zeros(4*N), (0.0, 500.0), parameters)
> 
> solve(prob, lsoda(), saveat=0.1, reltol=1e-8, abstol=1e-8)
> 
> end
> 
> #@time solution = runprob(N)
> 
> ```

I get:

```julia
julia> @time runprob(20);
  0.132728 seconds (195.93 k allocations: 18.890 MiB)

julia> @time runprob(90);
  1.365390 seconds (512.74 k allocations: 54.033 MiB)

julia> @time runprob(100);
  5.902604 seconds (563.84 k allocations: 59.490 MiB, 0.33% gc time)

julia> @time runprob(120);
  6.846018 seconds (559.88 k allocations: 62.312 MiB, 0.26% gc time)

julia> @time runprob(200);
 16.374404 seconds (837.08 k allocations: 96.083 MiB, 0.10% gc time)

```

Perhaps there are indeed some “jumps” in the performance?

---

_[View the full topic](https://discourse.julialang.org/t/python-vs-julia-ode-solver/56123)._