I stumbled upon a problem while reproducing some numerical solution of MATLAB (ode45) in Python (using solve_ivp) using the default Abstol (1e-6) and Reltol (1e-3) for both cases. And to my surprise for a particular case the solution computed in MATLAB and Python does not match. Each function uses (RK4(5)), all other parameters, initial condition are same, and both solutions look a valid solution, as they converge to a solution of periodic orbit but of different kind. As I plotted a variable for the whole-time length in Python, I found the trajectory is making many rounds near two regions, and I suspect one them is the solution computed by MATLAB. What exactly is going on here? Are there two co-existing solutions? How can I can correctly tell which one is correct? My system is a coupled non-linear set of ordinary differential equations. Any comment, suggestion in this regard will be appreciated.
I am attaching a few figures that I obtained in MATLAB as well as Python herewith.
Hello, I think there are two possible causes:
1. the conventional solvers not always work for Stiff Equation, and you can choose a solver with postfix "s" (in MatLab);
2. You system has a "chaotic nature" (like Three-body problem) and small changes in initial conditions and parameters cause large changes in behavior.
P.S. thus you can't trust to numerical solution sometimes
In the case of solving nonlinear, and specifically, chaotic systems, explicit Runge-Kutta methods cannot be called reliable. I recommend using geometric integrators, especially for simulating conservative systems. For stiff systems, however, implicit methods should be applied.
Nevertheless, the simulations performed with the same method and data type should be equal (or nearly equal) in various programming languages. I believe, the key reason for observed discrepancies is the different step control algorithms in two considered ODE solver implementations.
That's why I always recommend constructing finite-difference schemes yourself and not using built-in ODE solvers. This insures your experiments are reliable and reproducible. Second, one should remember that we investigate a discrete map of the continuous system in PC simulation, not the system itself. Some properties may be different, and some of them may not even exist in the original continuous system.
Please, feel free to contact me if you need any further help with this simulation.
Hello.
It is interesting question. Earlier I ran into the system for which all trajectories exhibit drifting behavior and very sensitive to choise of method and its options.
I can suppose that Matlab and Python implement the algorithm of error control in different ways. Are you try other solvers? For example from gsl: rkf45, implicit Gaussian of 2nd and 4th orders?
I think that there is an approach to check quality of numerical solution.
1. You can try to find partial case of given system whish is admits the first integral.
2. The next, it is convinient to consider the behavior of the initial system in the vicinity of integrable case. Also you can try the asymptotic methods for deriving the semi-analitycal solution in the vicinity of integrable case during the time interval of the interest.
Las respuestas deben ser iguales en ambos softwares. Deben dar resultados aceptables y los resultados se comprueban resolviendo el sistema de ecuaciones diferenciales manualmente, usando el algoritmo de Runge-Kutta de cuarto orden.
Chaotic solutins can be a problem. But I think this is not the case as the solutions are periodic.
Is supect that the problem is the tolerance. The defaults are not very accurate. Sometimes I use even 1e-10 or 1e-12. So the most easy solution would be to use more stringent tolerances. Maybe gradually, and check how the two solutions start to become similar.
- Badges
- Science topic
- Similar topics
- Engineering
- Control Theory