Can't understand numerical errors on julia

Hello! First of all, I want to thank you all. This is the first time that I ask something here, but I always read you, and the community is great. I hope that, with enough experience, I can start answering things aswell :slight_smile: .

The reason why I’m here is that I can’t understand the behavior I get on a physical simulation im doing.
Previous research shows that, if I have a big perturbation only along one dimension of my system, then small perturbations should start growing along the transverse directions, but they never appear in my simulations using julia. The method im using to solve the equations basically relies on fft (FFTW.jl) and operator splitting.
In recent papers, I have seen the same system of equations solved using “noisier” methods, and they actually see the transversal growth. Even more, I have seen the system solved using the same method and library but using FORTRAN, and they see transversal growth.

At this point I should state that the system of equations has the following property: If there is exactly zero perturbation along some direction, then a perturbation can not be developed. A perturbation may only change its size, for example from numerical noise to dominate the system, but it cant be created from nothing.

Due to all I said, I started thinking: “Maybe, for some reason that I can’t understand, my simulations are too clean”, and to prove it I added Gaussian noise on the order of 10 to the -15, and voila! The simulations matched all previous research: simulations and experimental measures.

So, I am here to ask: How does julia handle the numerical errors, so that my simulation NEVER develops a component along the second dimension? I mean, julia returns exactly 0.0 amplitude in that direction. Not even something like 1e-310.

Whats the difference between julia and FORTRAN, such that the latter generates a noise component, and julia doesnt?

Many thanks in advance!

4 Likes

Welcome to Julia’s discourse!
Sincerely, i have no idea. A float64 is a float64 in julia and fortran, do you use any other Julia packages?

It is very hard to say much about this in the abstract, without seeing any code. In principle Julia does not do anything special, just IEEE 754, unless you use something like @fastmath.

Investing the effort into making an MWE could get you much more specific responses though.

4 Likes