1) Choose two natural interval extensions of the system under investigation.

These equations are mathematically equivalent, but they represent a different sequence of arithmetic operations, resulting in different results. This happens because some properties of conventional arithmetic are not valid in floating-point arithmetic.

F(X_n) = 2.6868X_n - 0.2462X_n³   

G(X_n) = 2.6868X_n - (0.2462X_n)X_n²

2) With exactly the same initial conditions, simulate the two natural interval extensions. In this case:

 x(1) = 0.1; y(1) = x(1);

for k=1:100

x(k+1,1) = 2.6868*x(k) - 0.2462x(k)^3;

y(k+1,1) = 2.6868*x(k) - (0.2463*y(k))*y(k)^2;

end

figure()

plot(1:60,x(1:60),'o-k')

hold on

plot(1:60,y(1:60),'x-r')

xlabel('n')

ylabel('F(X_n) ,G(X_n)')

3) The next step is to determine the lower bound error, which is described by the following equation:

where  are two pseudo-orbits derived from two interval natural extensions.

For a better visualization of the results, the error will be plotted in logarithm base 10, as shown in the algorithm.

for k=1:101

    lbe(k) = log10(abs(x(k)-y(k))/2);

end

figure()

plot(6:k,lbe(6:k),'-ok');

Evolution of relative error.