# Is Julia reliable for solving ordinary and stochastic differential equations?
**URL:** https://discourse.julialang.org/t/is-julia-reliable-for-solving-ordinary-and-stochastic-differential-equations/98374
**Category:** General Usage
**Created:** [May 5, 2023, 8:30am UTC](https://discourse.julialang.org/t/is-julia-reliable-for-solving-ordinary-and-stochastic-differential-equations/98374 "2023-05-05T08:30:03Z")
**Posts on this page:** 1
**Showing post:** 12
### Author:  [@Palli](https://discourse.julialang.org/u/Palli)
#### Post date: [May 6, 2023, 3:50pm UTC](https://discourse.julialang.org/t/is-julia-reliable-for-solving-ordinary-and-stochastic-differential-equations/98374/12 "2023-05-06T15:50:47Z")
> [@soldin](#):
>
> I know C++, Matlab, and R. Matlab is too slow for solving ODEs
I see from last month:
2023: MIT Center for Computational Science & Engineering MathWorks Prize for Outstanding Master’s Research in Computational Science & Engineering. For Songchen Tan (MIT Julia Lab) for his work on TaylorDiff.jl for use NeuralPDE.jl physics-informed neural networks (PINNs)
So it’s a Prize (also) in the name of MathWorks, i.e. the maker of MATLAB. The SciML ecosystem if for sure doing something (well man things) right. The only language (of those) potentially competing on speed with Julia is C++. The others not unless you can wrap fast libraries from other fast languages, and that can be done for e.g. BLAS, but I understand more challenging for solving equations.
I believe all the issues Yuri brought up, are solvable, if not already, and the only really interesting intriguing issue he brought up that time was OffsetArrays.jl, that I very much doubt is an issue for solving your equations:
> [@Discussion on "Why I no longer recommend Julia" by Yuri Vishnevsky](https://discourse.julialang.org/t/discussion-on-why-i-no-longer-recommend-julia-by-yuri-vishnevsky/81151/17):
>
> I do not think the article is entirely fair, while it is of much greater quality than the average article criticizing Julia. The comparison is done to older and more mainstream libraries. While this is a legitimate viewpoint for choosing what to use for your work just now, well, it is not fair in the more general sense. The comparison would need to be between ecosystems/languages at the same level of maturity. The more interesting question to me is if Julia will have these problems when it rea…
> There is a lot of problems with OffsetArrays.jl but most other languages do not even something like OffsetArrays.jl or the expectation that most code written would automatically work with custom indexes.
I (still) have full confidence in Julia core developers and e.g. Chris, for math/technical computing, such as you bring up. Julia however isn’t bug-free (no software ever will be), but issues are handled well.
Yuri no longer has any open PR in Julia, all of his 11 have been merged or closed. But he also has open issues, including 2 from 2 weeks ago, rather interesting:
[Issues · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues?q=author%3Ayurivish+is%3Aopen) [after clicking the link the + needs to be changed to a space; I’m not sure how to put in a correct link for Discource, it breaks if I do it, so I’m not sure it’s possible…]
Note, I think all integer and floating-point division is correct in Julia for all types (handled by CPU instructions, i.e. assuming correctly), from his issue, it’s about floored division (e.g. `div`)
>
>
> Julia has a \[\`div\`\](https://docs.julialang.org/en/v1/base/math/#Base.div) functi…on that is used to implement \[floor division\](https://docs.julialang.org/en/v1/base/math/#Base.fld) (\`fld\`) and \[ceil division\](https://docs.julialang.org/en/v1/base/math/#Base.cld) (\`cld\`). I think I’ve found a bug that causes all of these division functions to return incorrect results.
>
> Floor division is documented to return the largest integer less than or equal to \`x / y\`. It should never return a number that is \*greater\* than \`x / y\`.
>
> But it does:
>
> \`\`\`julia
> julia\> 514 / Float16(0.75)
> Float16(685.5)
>
> julia\> div(514, Float16(0.75)) # should round down, but rounds up instead
> Float16(686.0)
>
> julia\> fld(514, Float16(0.75)) # likewise
> Float16(686.0)
>
> julia\> fld(514, Float16(0.75)) ≤ 514 / Float16(0.75)
> false
> \`\`\`
>
> Similarly, ceil division should never return a number that is \*smaller\* than regular division, but it does:
>
> \`\`\`julia
> julia\> 515 / Float16(0.75)
> Float16(686.5)
>
> julia\> cld(515, Float16(0.75)) # should round up, but rounds down instead
> Float16(686.0)
>
> julia\> cld(515, Float16(0.75)) ≥ 515 / Float16(0.75)
> false
> \`\`\`
>
> This behavior is not limited to 16-bit floats. Here’s a case where \`fld\` produces an incorrect result for \`Float32\` inputs:
>
> \`\`\`julia
> julia\> 4\_194\_307 / Float32(0.75) # = 5592409.5
> 5.5924095f6
>
> julia\> fld(4\_194\_307, Float32(0.75)) # = 5592410, incorrectly rounded up
> 5.59241f6
>
> julia\> fld(4\_194\_307, Float32(0.75)) ≤ 4\_194\_307 / Float32(0.75)
> false
> \`\`\`
>
> And here’s the same for \`cld\`:
>
> \`\`\`julia
> julia\> 4\_194\_308 / Float32(0.75) # = 5592410.5
>
> julia\> cld(4\_194\_308, Float32(0.75)) # = 5592410, incorrectly rounded down
> 5.59241f6
>
> julia\> cld(4\_194\_308, Float32(0.75)) ≥ 4\_194\_308 / Float32(0.75)
> false
> \`\`\`
>
> The equivalent operations in Python produce the correct results:
>
> \`\`\`python
> \# For 16-bit floats:
> \>\>\> np.float16(514) / np.float16(0.75) # regular division
> 685.5
> \>\>\> np.float16(514) // np.float16(0.75) # floor division
> 685.0
>
> \# For 32-bit floats:
> \>\>\> np.float32(4\_194\_307) / np.float32(0.75) # regular division
> 5592409.5
> \>\>\> np.float32(4\_194\_307) // np.float32(0.75) # floor division
> 5592409.0
> \`\`\`
>
> Examples of this incorrect behavior are not hard to find – for most floats, you can find a divisor that will make either \`fld\` or \`cld\` return the wrong answer.
>
> Here are some examples for \`Float16\` where either \`fld\` or \`cld\` is incorrect:
>
> \- \`cld(1, Float16(0.000999)) \< 1 / Float16(0.000999)\`
> \- \`cld(2, Float16(0.001999)) \< 2 / Float16(0.001999)\`
> \- \`cld(3, Float16(0.002934)) \< 3 / Float16(0.002934)\`
> \- \`cld(4, Float16(0.003998)) \< 4 / Float16(0.003998)\`
> \- \`fld(5, Float16(0.004925)) \> 5 / Float16(0.004925)\`
>
> And here are some for \`Float32\`:
>
> \- \`fld(5, Float32(6.556511e-7)) \> 5 / Float32(6.556511e-7)\`
> \- \`fld(10, Float32(1.3113022e-6)) \> 10 / Float32(1.3113022e-6)\`
> \- \`fld(11, Float32(1.4305115e-6)) \> 11 / Float32(1.4305115e-6)\`
> \- \`cld(16, Float32(2.8014183e-6)) \< 16 / Float32(2.8014183e-6)\`
> \- \`cld(17, Float32(2.2053719e-6)) \< 17 / Float32(2.2053719e-6)\`
>
> For simplicity I’ve presented examples where the first argument is an integer; this bug also occurs for non-integral inputs.
>
> A divisor producing the wrong result can be found for over 51% of all possible 16-bit floats. I have not evaluated how widespread this is for \`Float32\`, but the results above suggest that it is similarly easy to find failures there too.
>
> I’ve tracked the invocations down to \[this definition\](https://github.com/JuliaLang/julia/blob/72aec423c2ab9f80c249d63fdd68b35833cfd7ed/base/div.jl#L370) of \`div\`, which has existed at least as far back as \[Julia 1.6\](https://github.com/JuliaLang/julia/blob/f9720dc2ebd6cd9e3086365f281e62506444ef37/base/div.jl#L279):
>
> \`\`\`julia
> div(x::T, y::T, r::RoundingMode) where {T\<:AbstractFloat} =
> convert(T, round((x - rem(x, y, r)) / y))
> \`\`\`
> That error rate suggests that this bug manifests once out of every ~1,700 randomly chosen Float16 divisions
> […]
> Based on a billion samples, the bug manifests once out of every ~10,000 randomly chosen Float32 divisions and once out of every ~74,000 Float64 divisions.
Division (and e.g. square root) CAN be correct, in general, can only be correctly rounded, so you always have a small error. According to chaos theory, small errors can blow up. I don’t think this is too worrying, i.e. to have a tiny bit smaller error, than the correctly rounded result. If you’re sensitive to it, then most likely the small error in the correctly rounded error too?
I’m curious, for differential equations, you can have arbitrary operators, e.g. division common, but would you have floored division often (or ever), or `div`? If not, then you don’t need to worry about that issue at least.
Since Yuri was posting new issues recently, he at least cares about Julia, and likely still uses, or why would he even have know or posted those two weeks ago.
---
_[View the full topic](https://discourse.julialang.org/t/is-julia-reliable-for-solving-ordinary-and-stochastic-differential-equations/98374)._