# using ODE45 to solve simultaneously independent ODE's vs. calling the function twice - different results...

73 views (last 30 days)

Show older comments

Hi,

i am trying to solve for the displacement of a spring & mass with rotating eccentric mass in X and Y directions.

If i solve it using the following code:

[t,y] = ode45(@odefun, tspan, ic);

function dydt = odefun(t,y)

%SDOF (Uncoupled) Forced (Unbalance), Undamped

dydt = zeros(4,1);

dydt(1) = y(2);

dydt(2) = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);

dydt(3) = y(4);

dydt(4) = -(k/M)*y(3)+(m/M)*e*w^2*sin(w*t);

end

I get one result (which seems incorrect for the second equation).

If i solve it using the following code (which effectively separates it into two independent ODEs:

[t,y] = ode45(@odefun, tspan, ic);

function dydt = odefun(t,y)

%SDOF (Uncoupled) Forced (Unbalance), Undamped

dydt = zeros(2,1);

dydt(1) = y(2);

dydt(2) = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);

end

[t2,y2] = ode45(@odefun2, tspan, ic2);

function dydt2 = odefun2(t,y2)

dydt2 = zeros(2,1);

dydt2(1) = y2(2);

dydt2(2) = -(k/M)*y2(1)+(m/M)*e*w^2*sin(w*t);

end

I get a different result, which seems more in line with what i expect.

Why would i get different results for these two methods?

See attached images..

Any help would be greatly appreciated!

### Accepted Answer

Jan
on 2 Jun 2018

Edited: Jan
on 2 Jun 2018

This is a problem of the accuracy only. Decrease the tolerance to get equivalent trajectories for the simultaneous and separate evaluation:

function YourFcn

for method = 1:2

if method == 1

opt = odeset('RelTol', 1e-3, 'AbsTol', 1e-6);

else

opt = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);

end

% Solve simultaneously:

tspan = [0, 10];

ic = [0,0,0,0];

[t,y] = ode45(@odefun, tspan, ic, opt);

figure;

for k = 1:4

subplot(2,4,k);

plot(t, y(:,k));

end

% Solve separately:

ic = [0,0];

[t,y] = ode45(@odefun1, tspan, ic, opt);

subplot(2,4,5);

plot(t, y(:,1));

subplot(2,4,6);

plot(t, y(:,2));

[t,y] = ode45(@odefun2, tspan, ic, opt);

subplot(2,4,7);

plot(t, y(:,1));

subplot(2,4,8);

plot(t, y(:,2));

end

end

function dydt = odefun(t,y)

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dydt = [y(2); ...

-(k/M)*y(1)+(m/M)*e*w^2*cos(w*t); ...

y(4); ...

-(k/M)*y(3)+(m/M)*e*w^2*sin(w*t)];

end

function dydt = odefun1(t,y)

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dydt = [y(2); ...

-(k/M)*y(1)+(m/M)*e*w^2*cos(w*t)];

end

function dydt = odefun2(t,y)

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dydt = [y(2); ...

-(k/M)*y(1)+(m/M)*e*w^2*sin(w*t)];

end

High tolerances - different results: First row is the simultaneous solution, second row the separated:

Lower tolerances, which means a more accurate integration - equivalent results:

With the lower tolerances the trajectory shows the expected periodic behavior also, while for the less accurate integration the local accumulated discreatization error dominates the solutions. It is simply too coarse.

Note that the default of the absolute tolerance of 1e-6 is in the magnitude of your values. Then this is not a meaningful limit to determine a proper stepsize.

##### 4 Comments

Jan
on 5 Nov 2020

For interdependencies you have to solve the equations together. Then there is the special case, that the two systems have a very different frequency behaviour, e.g. two coupled pendulums and the eigenfrequency of one is 1e8 times higher than the other. In this case, the step size control of ODE45 considers the high frequency part of the system only and the accumulated rounding errors can dominate the solution for the low frequency part. There is no smart solution to handle this problem, because the Runge-Kutta-one step solvers are not designed for this case. You need a stiff solver for a reliable solution.

This is a problem for the simulation of a car: If you simulate the body, the wheels and consider the engine as a black box producing force, you will not catch the high-frequency vibrations which cause discomfort and noise. A complete model including all small parts is very expensive and you have to choose the integrator carefully. Even then the parameters of the integrator can influence the results massively.

I've seen (too) many scientific calculations, in which a complex system of differential equation was "solved" by a standard integrator like ODE45. Without an analysis of the stiffness (e.g. caused be extremely different eigenfrequencies) and the sensitivity (small variation of initial values and parameters), the result is not reliable or valid. Numerical maths is a profession and it takes some years to learn how to apply and control the different methods. So if you have to integrate a conmplicated system, ask a scientist working in the field of numerical maths.

Btw. it was a development team of VW, who ran a simulation of a detailed car model using an ODE45 integrator written by their own. They failed tremendously due to bugs in the stepsize control and the stiffness of the system. Fortunately this was in the 1990'th and today experienced mathematicians are members of the team.

### More Answers (1)

Abraham Boayue
on 2 Jun 2018

Edited: Abraham Boayue
on 2 Jun 2018

Answers have been inserted into your codes (I just modified them a bit).

function Plot_odes

T = 10;

N = 1000;

t = 0 :T/(N-1):T;

ic = [0 0];

[t,x]= ode45(@odefun, t,ic );

[t,y]= ode45(@odefun2, t,ic );

figure

plot(t ,x(:,1),'linewidth',2,'color','r')

hold on

plot(t ,y(:,1),'linewidth',2,'color','b')

legend('x','y');

a = title('x(t) y(t)');

set(a,'fontsize',14);

a = ylabel('y');

set(a,'Fontsize',14);

a = xlabel('t [0 T]');

set(a,'Fontsize',14);

xlim([0 T])

grid

% Your two ODEs are completely independent of each other and

% therefore, the solution of one does not interfere with the other. In other

% words, they are not simultineous equations; you actually tricked ode45

% into believing that it was solving simultaneous equations in the first

% case. However, the two functions (odefun and odefun2) below are valid solutions of

% the odes because the vectors dxdt = [x' x''], dydt = [y' y''] and their

% solution vectors, X = [ x1 x2] and Y = [y1 y2] (where x1 = x(t) as the

% required solution for the displacement in the x direction and x2 its

% derivative), are uniform vectors with each element dependent on the other. On the otherhand, in your first method,

% you had dydt = [x' x'' y' y''] and Y = [x1 x2 y1 y2] which can not be right

% because x and y are independent variables.

function dxdt = odefun(t,y)

%SDOF (Uncoupled) Forced (Unbalance), Undamped

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dxdt_1 = y(2);

dxdt_2 = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);

dxdt =[dxdt_1; dxdt_2];

function dydt = odefun2(t,y2)

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dydt_1 = y2(2);

dydt_2 = -(k/M)*y2(1)+(m/M)*e*w^2*sin(w*t);

dydt =[dydt_1; dydt_2];

Here is your first method.

function Plot_odes

T = 10;

N = 1000;

t = 0 :T/(N-1):T;

ic = [0 0 0 0]; % Initial conditions

[t,z]= ode45(@odefun, t,ic );

figure

plot(t ,z(:,1),'linewidth',2,'color','r')

hold on

plot(t ,z(:,3),'linewidth',2,'color','b')

legend('x','y');

a = title('x(t) y(t)');

set(a,'fontsize',14);

a = ylabel('y');

set(a,'Fontsize',14);

a = xlabel('t [0 7]');

set(a,'Fontsize',14);

xlim([0 T])

grid

% Even though you were able to get some results, however,

% they are misleading because ode45 sees

% dydt = [y' y'' y''' y''''] and gives Y = [y1 y2 y3 y4].

% See above for further explanations

function dydt = odefun(t,y)

%SDOF (Uncoupled) Forced (Unbalance), Undamped

k = 104000;

M = 200;

m = 1/16;

w = 1;

e = 1;

dydt_1 = y(2);

dydt_2 = -(k/M)*y(1)+(m/M)*e*w^2*cos(w*t);

dydt_3 = y(4);

dydt_4 = -(k/M)*y(3)+(m/M)*e*w^2*sin(w*t);

dydt =[dydt_1; dydt_2;dydt_3; dydt_4];

##### 3 Comments

Chenko DUSHIMIMANA
on 4 Nov 2020

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!