Why Runge-Kutta methods cannot find the solution of Lorenz system? The solution of the following Lorenz system
s=10; r=28; b=8/3;
f = @(t,y) [-s*y(1)+s*y(2); -y(1)*y(3)+r*y(1)-y(2); y(1)*y(2)-b*y(3)];

in the interval $[0,8]$ with initial values $y_1(0)=1$,$y_2(0)=0$,$y_3(0)=0$ using MATLAB ode45() function are
>> [t y] = ode45(f,[0 8],[1 0 0]);
>> [t(end) y(end,:)]
ans =
    8.0000   -7.3560   -5.6838   27.8372

But the results are very poor with RK4 even with one million sub-interval (i.e. h=(8-0)/1e6)
ans = 
    8.0000   -7.4199   -5.6639   28.0052

My questions are:


*

*Why the results are different?

*Is ode45() the best function for Lorenz system?

*How can I improve the RK4 accuracy?

*Is Lorenz system an stiff equation?

 A: The Lorenz system
is notoriously unstable.
Try perturbing the initial conditions
by a very small amount
(say $10^{-6}$)
and see how much
the result changes.
You might try other solvers,
including a predictor-corrector
of some type.
Also,
make sure that
the computations are done
in double precision, at least.
A: *

*ode45 uses dynamical step size adaptation and is an order 5 method, thus it is more precise for the same effort than RK4 with constant step size (even if the step size control uses the "wrong" error estimate for the step size control). It is advisable to explicitly set the error tolerances, as the default tolerances might be too low. I get from

opts = odeset('AbsTol',1e-12,'RelTol',1e-10);
[t y] = ode45(f,[0 8],[1 0 0],opts);
[t(end) y(end,:)]

the result
    8.0000   -7.4186   -5.6635   28.0027

with relative and absolute error tolerances of 1e-10 and 1e-12. This is closer to your RK4 result. The dynamically selected step sizes range from 2e-4 to 6e-3.


*There is no best numerical method. But implicit methods with step size control might be better.


*Choose a smaller step size. But due to the increasing number of sub-steps the accumulated numerical noise will impose a limit in increasing precision.


*Yes. In general, non-linear ODE that are not sub-linear (i.e., there is no global Lipschitz constant) have regions where they may be stiff or produce dynamical blow-up.
