This is the plot made by Mathematica for the maximum time of 100

This is the plot made by Python for the maximum time of 100

You can compare the data set generated rather than he figures

@Asesh: The datas are different (I just used the figures to have a better sense of the difference)

To me it is not surprising that you have different results, since it is chaotic it is natural that small differences due to error in the numerical computation gives diverging results. However the differences in the results should not effect other parameter of the solution as fractal dimension of the chaotic attractor or its poincare map.

I usually use the solver "ode45" in matlab, to check the reliability of the solution I would decrease the allowed tolerance until I do not have significant changes with further decrease of it.

@Giuseppe: Thank you!

I will try ode45 in Matlab as well. I have a parameter in my equations which shows the adaptation and the values of the parameer demonstrate significant information about the dynamics. I want to report the results on the values of this parameter, but it is hard to decide when the results are reliable. So, although they are all chaotic, but it is also important for me how the state variables change.

If you are using Mathematica's NDSolve function, you should check out the Help section for NDSolve. In particular, examine the various options at your disposal when using NDSolve. For example, work with AccuracyGoal and PrecisionGoal. In NDSolve Help section there is an example of how these two options can be used to extend the accuracy of a solution trajectory from a ODE that exhibits sensitivity to initial conditions. Also look at WorkingPrecision and MaxSteps options.

@Brian: Thank you for the response. I have been playing with them as well and every change can lead to a new result.

Let me ask about my method of tackling the problem: As I want to integrate the ODEs for a long time, I run out of memory. So, to solve the issue, I integrate the equations for a short period of time, save the last values of all state variables, make those values the initial conditions of the next integration, integrate for another short time and so on. Do you think this approach might be problematic?

If your system of equations is sensitive to ICs, then you should be particularly careful. First of all I would make sure all the parameters in my model are integers. For example, if you have parameter value of say 0 .236, replace it with 236/1000. Then attempt to integrate with WorkingPrecision->20. When you stop your integration the values of your function should presumably have 20 digits of precision, and then you can restart the process using the stop conditions. But check to see what you get working with the other Options.

@Brian: Thank you for your response.

My system of equations is very sensitive to ICs (it is very similar to Lorenz Attractor).

In Mathematica, I changed all the parameters to integers (or rational numbers), as you suggested. I also increased the WorkingPrecision to 20 or even a higher number. I also set MaxSteps to Infinity and took the advantage of options such as AccuracyGoal. Despite all these, it seems that breaking the integration into two (or more) parts does not make the same result as doing the integration in one continuous chunk.

I also tried the same idea in Python and noticed a slight difference between doing it in a few chunks and one chunk. So, it seems that I have not been able to correctly implement the idea of breaking integration into sub-integrations.

@V. G. Kurbatov: Thank you!

Dear Roozbeh,

I would like to second Guiseppe's suggestion: if you have a chaotic dynamics that is sensitive to initial conditions, why do you decide to focus on the numerical values of their (long term) trajectories, i.e. on the most sensitive aspect of your system? Instead, I would focus on one of the robust quantities of their dynamics, such as fractional dimension, Lyapunov exponent, KS entropy, a.s.o. and analyze how any of those measures behaves under parameter changes, etc.

Find if your system is NOT chaotic for certain value of the parameters as it happens for the Lorenz equations: one way it to look for the Lie symmetries admitted by your equations.

@Harold: Thank you for your suggestions. I will consider them.

@Maria: In fact I know that the system starts with chaotic regime. I set the initial parameters so that it already generates chaos and then I want to show if the system avoids chaos or not.

I researched this problem in my article "Numerical and physical modeling of the dynamics of the Lorenz system". I used a modified method of power series. The error at the step can be decreased by varying the accuracy of the remainder term estimate in the series, which is not allowed in the Euler method or 4th order Runge-Kutta method. In the considered modification of the power series method, an advantage over the general scheme of the Taylor series method is that the expansion coefficients can be rapidly calculated by fast formulas in comparison to the procedure of symbol differentiation of the right-hand sides of the system equations. In addition, in the nonlinear case much memory is needed to store the symbol expressions in the calculation of the higher-order derivatives.

Article Numerical and physical modeling of the dynamics of the Lorenz system

compare graph of deterministic and stochastic i.e indeterministic approach...

mathematica commands of

Entropy of chaotic System with respect to changing parameter