ChrisRackauckas

That said, Symbolics was released in 2022. The thread you posted was from 2022. In 2026, it does support limits (https://github.com/SciML/SymbolicLimits.jl), integrals, (https://docs.sciml.ai/SymbolicIntegration/dev/) etc., and in fact that symbolic integration is much more advanced than SymPy, except on symbolic integration of ODEs (surprisingly, though that has been pretty deliberate since that domain almost never has an analytical solution so it's not the most useful to real practitioners and is more of a teaching tool)

But all of this is more about maturing the ecosystem to be more amenable to static compilation and analysis. The whole SciML stack had an initiative starting at the beginning of this summer to add JET.jl testing for type inference everywhere and enforcing this to pass as part of standard unit tests, and using AllocCheck.jl for static allocation-free testing of inner solver loops. With this we have been growing the surface of tools that have static and real-time guarantees. Not done yet, some had to be marked as `@test_broken` for having some branch that can allocate if condition number hits a numerical fallback and such, but generally it's getting locked down. Think of it as "prototype in Julia, deploy in Rust", except instead of re-writing into Rust we're just locking down the behavior with incrementally enforcing the package internals to satisfy more and more static guarantees.

With Julia v1.12's small binary generation, we plan to release forms via binaries with C ABIs over the next year as well.

Are those sufficient or should we consider supporting other deployments?

If people want a tl;dr, the main idea is you can construct an ODE where the forward pass is trivial, i.e. the ODE solution going forwards is exact, but its derivative is "hard". An easy way to do this is to make it so you have for example `x' = x - y, y' = y - x`, with initial conditions x(0)=y(0). If you start with things being the same value, then the solution to the ODE is constant since `x' = y' = 0`. But the derivative of the ODE solution with respect to its initial condition is very non-zero: a small change away from equality and boom the solution explodes. You write out the expression for dy(t)/dy(0) and what you get is a non-trivial ODE that has to be solved.

What happens in this case though is that automatic differentiation "runs the same code as the primal case". I.e., automatic differentiation has the property that it walks through your code and differentiates the steps of your code, and so the control flow always matches the control flow of the non-differentiated code, slapping the derivative parts on each step. But here the primal case is trivial, so the ODE solver goes "this is easy, let's make dt as big as possible and step through this easily". But this means that `dt` is not error controlled in the derivative (the derivative, being a non-trivial ODE, needs to have a smaller dt and take small steps in order to get an accurate answer), so the derivative then gets error due to this large dt (or small number of steps). Via this construction you can make automatic differentiation give as much error as you want just by tweaking the parameters around.

Thus by this construction, automatic differentiation has no bound to the error it can give, and no this isn't floating point errors this is strictly building a function that does not converge to the correct derivative. It just goes to show that automatic differentiation of a function that computes a nice error controlled answer does not necessarily give a derivative that also has any sense of error control, that is a property that has to be proved and is not true in general. This of course is a bit of a contrived example to show the point that the error is unbounded, but then it points to real issues that can show up in user code (in fact, this example was found because a user opened an issue with a related model).

Then one thing that's noted in here too is that the Julia differential equation solvers hook into the AD system to explicitly "not do forward-mode AD correctly", incorporating the derivative terms into the time stepping adaptivity calculation, so that it is error controlled. The property that you get is that for these solvers you get more steps to the ODE when running it in AD mode than outside of AD mode, and that is a requirement if you want to ensure that the user's tolerance is respected. But that's explicitly "not correct" in terms of what forward-mode AD is, so calling forward-mode AD on the solver doesn't quite do forward mode AD of the solver's code "correctly" in order to give a more precise solution. That of course is a choice, you could instead choose to follow standard AD rules in such a code. The trade-off is between accuracy and complexity.