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

73 views (last 30 days)
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!
  2 Comments

Sign in to comment.

Accepted Answer

Jan
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
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.

Sign in to comment.

More Answers (1)

Abraham Boayue
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
Chenko DUSHIMIMANA on 4 Nov 2020
Hello, Thanks the detail explanation about ode45. What if the differential equations to be solved by ode45 are interconnected(i.e. they depend on each other, that is, to find the solution of one system of equations(say, system 1) requires the knowledge of the solution from the other system of equations(say, system 2)). Where, the "tspan" for both system may differ or be similar. How Can we solve these two systems? simultaneously or separately?

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!