How to get acceleration data while using ODE solver for some governing equation of motion?

I am using ODE solver for a system of equations governing the motion of a model. In fact I am integrating the acceleration expressions by using ODE solver after converting them into state space form. and the Output I can get is position and velocities. How can I get the corresponding acceleration data? As the acceleration is alraedy definded in the form of expression in dydt vector, can't I store the value of thethis vector which can give the acceleration? Please suggest me any simpler way to get acceleration without using any numerical approach to calculate acceleration from velocity data or finding the 3rd order diff. equation and then integrate to get acceleration.

Answer 1

[Sam Chak]

Hi @Paras

Use the deval() command to acceleration data.

sol     = ode45(@(t, x) [x(2); -x(1)], [0 10], [1; 0]);
t       = linspace(0, 10, 1001);
[x, xp] = deval(sol, t);
plot(t, x, t, xp(2,:)), grid on % dotx = dx/dt, % ddotx = d²x/dt²
xlabel('t'), ylabel('\bf{x}')
legend('x', 'dotx', 'ddotx')

[Paras]

Thank you so much dear Sam. This is exactly what I was looking for and you caught it rightly. It has simplified things for me which otherwise be tidious. Thanks for your valuable conribution :)

[Torsten]

Isn't the acceleration -x(1) in the example above ? y'' usually cannot be retrieved from a call to "deval". You will have to call the function handle that you pass to "ode45" itself.

[Paras]

Dear Torsten, -x(1) is acceleration expression to govern the system. sol will not give us that acceleration data after using ode solver as it intergrate that acceleration expression to give position and velocity. So x is storing position x(1,:) and velocity x(2,:) but we can also get xp using deval which will give us first devivative of x that is velocity xp(1,:) and acceleration xp(2,:). Hope it address your query.

[Torsten]

"deval" is less accurate than the "correct" acceleration that you define via your ODE function. Thus after integration has finished, you should call your ODE function with the solution values for position and velocity to get acceleration.

Answer 2

[Sam Chak]

Hi @Paras

I believe I understand what @Torsten is trying to imply. He is an expert in solving dynamical systems. It is true that the 'deval()' function cannot be applied in all cases to obtain the 2nd time-derivative data of the solution in a 2nd-order system. The 'deval()' command can only be used if the dynamical system is in canonical form. The following example demonstrates the difference in the results between a non-canonical system and a canonical system, despite both systems being dynamically equivalent.

Fortunately, most uncomplicated mechanical motion systems are described in canonical form.

Non-canonical system

Canonical system

%% Non-canonical form (original system)
function dx = ode1(t, x)
    dx  = [- x(1) + x(2)
           - x(1) - x(2)];
end

%% Canonical form (after coordinate transformation)
function dz = ode2(t, z)
    dz  = [             z(2)    % dz1 = dx1 = - x(1) + x(2)
           - 2*z(1) - 2*z(2)];
end

%% Solving Non-canonical system
x10     = 1;            % initial value of x1
x20     = 0;            % initial value of x2
x0      = [x10; x20];
sol1    = ode45(@ode1, [0 10], x0);
t1      = linspace(0, 10, 1001);
[x, xp] = deval(sol1, t1);
subplot(211)
plot(t1, x, t1, xp(2,:)), grid on
xlabel('t'), ylabel('\bf{x}')
legend({'$x_{1}$', '$x_{2}$', '$\dot{x}_{2} \neq \ddot{x}_{1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 12)
title({'Non-Canonical: Incorrect usage of deval() to get $\ddot{x}_{1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 14)

%% Solving Canonical system
z10     = x10;          % initial value of z1, equivalent to x1(0)
z20     = - x10 + x20;  % initial value of z2(0) = dz1(0) = - x1(0) + x2(0)
z0      = [z10; z20];
sol2    = ode45(@ode2, [0 10], z0);
t2      = linspace(0, 10, 1001);
[z, zp] = deval(sol2, t2);
subplot(212)
plot(t2, z, t2, zp(2,:)), grid on
xlabel('t'), ylabel('\bf{x}')
legend({'$z_{1} = x_{1}$', '$z_{2} = \dot{x}_{1}$', '$\dot{z}_{2} = \ddot{x}_{1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 12)
title({'Canonical: Correct usage of deval() to get $\ddot{x}_{1}$'}, 'Interpreter', 'LaTeX', 'fontsize', 14)