# Numerical Inaccuracies for Ranges of Parameters When Using ode45 for Two Identical (But Transformed) Systems

1 view (last 30 days)

Show older comments

Hi,

I have two systems of ODEs as defined below. A 4D system:

dxdt = [

w + (k/n)*(sin(x(2)-x(1)+a)-r*sin(2*(x(2)-x(1))) + sin(x(3)-x(1)+a)-r*sin(2*(x(3)-x(1))) + sin(x(4)-x(1)+a)-r*sin(2*(x(4)-x(1))));

w + (k/n)*(sin(x(1)-x(2)+a)-r*sin(2*(x(1)-x(2))) + sin(x(3)-x(2)+a)-r*sin(2*(x(3)-x(2))) + sin(x(4)-x(2)+a)-r*sin(2*(x(4)-x(2))));

w + (k/n)*(sin(x(1)-x(3)+a)-r*sin(2*(x(1)-x(3))) + sin(x(2)-x(3)+a)-r*sin(2*(x(2)-x(3))) + sin(x(4)-x(3)+a)-r*sin(2*(x(4)-x(3))));

w + (k/n)*(sin(x(1)-x(4)+a)-r*sin(2*(x(1)-x(4))) + sin(x(2)-x(4)+a)-r*sin(2*(x(2)-x(4))) + sin(x(3)-x(4)+a)-r*sin(2*(x(3)-x(4))));

];

And a 3D system:

dydt = [

(k/n)*(sin(-y(1)+a)-r*sin(-2*y(1)) + sin(y(2)-y(1)+a)-r*sin(2*(y(2)-y(1))) + sin(y(3)-y(1)+a)-r*sin(2*(y(3)-y(1))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

(k/n)*(sin(-y(2)+a)-r*sin(-2*y(2)) + sin(y(1)-y(2)+a)-r*sin(2*(y(1)-y(2))) + sin(y(3)-y(2)+a)-r*sin(2*(y(3)-y(2))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

(k/n)*(sin(-y(3)+a)-r*sin(-2*y(3)) + sin(y(1)-y(3)+a)-r*sin(2*(y(1)-y(3))) + sin(y(2)-y(3)+a)-r*sin(2*(y(2)-y(3))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

];

They are defined in such a way that the 3D system (dydt) is the reduced version of the 4D system (dxdt) via:

% y(1) = x(2) - x(1)

% y(2) = x(3) - x(1)

% y(3) = x(4) - x(1)

%%%%%

% -y(1) = x(1) - x(2)

% -y(2) = x(1) - x(3)

% -y(3) = x(1) - x(4)

%%%%%

% dy(1)/dt = dx(2)/dt - dx(1)/dt

% dy(2)/dt = dx(3)/dt - dx(1)/dt

% dy(3)/dt = dx(4)/dt - dx(1)/dt

That is, plotting the difference between x(2)-x(1) should recover the same result as y(1) (and so on). This holds true for most values of the parameters, however I have found that there are discrepancies for few ranges in parameter, a. Below is a script that I put together very quickly to demonstrate this. I have found that without incredibly low tolerances and manually set step sizes, the two systems will not agree. Even with these settings, the code below will not produce identical results.

figure % 3D

options = odeset('RelTol',3e-14,'AbsTol',1e-15,'Stats','on','InitialStep',1e-4,'MaxStep',1e-3);

[t,x]=ode45(@hmm,[0 15000],[0 0.1 0.2 0.3],options);

plot(t,x(:,4)-x(:,1))

title('4D')

figure % 4D

[t,y]=ode45(@hmmPD,[0 15000],[0.1 0.2 0.3],options);

plot(t,y(:,3))

title('3D')

% 4D system

function dxdt = hmm(t,x)

w=1;

k=1/5;

a=1.3;

r=1/4;

n=4;

dxdt = [

w + (k/n)*(sin(x(2)-x(1)+a)-r*sin(2*(x(2)-x(1))) + sin(x(3)-x(1)+a)-r*sin(2*(x(3)-x(1))) + sin(x(4)-x(1)+a)-r*sin(2*(x(4)-x(1))));

w + (k/n)*(sin(x(1)-x(2)+a)-r*sin(2*(x(1)-x(2))) + sin(x(3)-x(2)+a)-r*sin(2*(x(3)-x(2))) + sin(x(4)-x(2)+a)-r*sin(2*(x(4)-x(2))));

w + (k/n)*(sin(x(1)-x(3)+a)-r*sin(2*(x(1)-x(3))) + sin(x(2)-x(3)+a)-r*sin(2*(x(2)-x(3))) + sin(x(4)-x(3)+a)-r*sin(2*(x(4)-x(3))));

w + (k/n)*(sin(x(1)-x(4)+a)-r*sin(2*(x(1)-x(4))) + sin(x(2)-x(4)+a)-r*sin(2*(x(2)-x(4))) + sin(x(3)-x(4)+a)-r*sin(2*(x(3)-x(4))));

];

end

% 3D System

% y(1) = \psi_{21}

% y(2) = \psi_{31}

% y(3) = \psi_{41}

% THUS

% y(2) - y(1) = x(3) - x(2)

% y(3) - y(1) = x(4) - x(2)

% y(1) - y(2) = x(2) - x(3)

% y(1) - y(3) = x(2) - x(4)

function dydt = hmmPD(t,y)

k=1/5;

a=1.3;

r=1/4;

n=4;

dydt = [

(k/n)*(sin(-y(1)+a)-r*sin(-2*y(1)) + sin(y(2)-y(1)+a)-r*sin(2*(y(2)-y(1))) + sin(y(3)-y(1)+a)-r*sin(2*(y(3)-y(1))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

(k/n)*(sin(-y(2)+a)-r*sin(-2*y(2)) + sin(y(1)-y(2)+a)-r*sin(2*(y(1)-y(2))) + sin(y(3)-y(2)+a)-r*sin(2*(y(3)-y(2))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

(k/n)*(sin(-y(3)+a)-r*sin(-2*y(3)) + sin(y(1)-y(3)+a)-r*sin(2*(y(1)-y(3))) + sin(y(2)-y(3)+a)-r*sin(2*(y(2)-y(3))) - sin(y(1)+a)+r*sin(2*y(1))-sin(y(2)+a)+r*sin(2*y(2))-sin(y(3)+a)+r*sin(2*y(3)));

];

end

I have spent a fair bit of time now with the documentation and have tried a few of the other solvers, but they haven't appeared to be better suited to this task. ode45() with these settings seems to work for the majority of a, just not around 1.12<=a<=1.22 and 1.28<=a<=1.37. Is this simply a case of needing to continually decrease the step sizes until it is accurate or is there a better way of going about it? It is getting quite computationally expensive at this point. Thank you for your time.

##### 1 Comment

Torsten
on 26 Jul 2022

Edited: Torsten
on 26 Jul 2022

### Answers (1)

Jan
on 28 Jul 2022

Edited: Jan
on 28 Jul 2022

I confirm the observed instabilities. The result depends critically on the initial value and the parameters. The trajectory has some fix points - the diagrams looks, like there are 2 of them. The radius of convergence of these fix points seems to be very small, so it looks randomly, when they are entered.

This means, that your function is instable. Numerical effects can increase such instabilities. You can check this by using the final as initial value and integrate in the opposite direction: change the order of the time span. Then vary the initial value to find the readius of convergence. Then examine the rate of growing apart from the order of convergence. My guess is an exponential behaviour in a local area.

Compare your function with simulating pool billiard without friction. On a microscopic level the surface of the table is smooth including the holes. But if a ball gets too near to one of the holes, in stucks in there. Which hole it is, depends critically on the initial values and the parameters, physically. For the simulation the accuracy and stepsize of the integrator matters in addition: If a step is calculated over a tiny hole, its effects are not detected. Limiting the step size manually has a strong effect then.

This means, that your function is mathematically and/or numerically instable. This happens for simulating a pen standing on its tip or a Galton board. Also you get a final value, the output of the integration should not be called a "result" of a simulation, because the instabilities do not allow to create a proper simulation model of the physical object.

This is neither a problem of Matlab or its integrators, but the validation of the model fails. In consequence your code is an expensive pseudo-random number generator with 1 bit output.

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!