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

9 次查看（过去 30 天）
Paras
Paras 2024-7-12
编辑： Torsten 2024-7-14
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.

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2024-7-12
Ran in:
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')
  4 个评论
显示 2更早的评论
Paras
Paras 2024-7-14
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
Torsten 2024-7-14
编辑：Torsten 2024-7-14
"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.

请先登录，再进行评论。

更多回答（1 个）

Sam Chak
Sam Chak 2024-7-14
Ran in:
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)

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by