How to solve 2nd order coupled differential equations?

30 次查看（过去 30 天）
Gloria
Gloria 2023-9-23
评论： Sam Chak 2023-9-28
I'd like to solve the 2 equation system below. I usually reduce the problem to a system of first order differencial equations and then use, for instance, ode113. However, here I don't know what to do because in each equation there are both the second derivatives of the variables.
Could anyone please help?
thanks a lot in advance

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

回答（2 个）

Torsten
Torsten 2023-9-23
M*x'' + A*x' + B*x = v
->
x'' = M^(-1) * (v - A*x' - B*x)
->
u1' = u2;
u2' = M^(-1) * (v - A*u2 - B*u1)
  2 个评论
Gloria
Gloria 2023-9-23
Hi Torsten, thank you for your answer. Unfortunately I don't think this solution works because if
u1=x
u2=x'
u3=y
u4=y'
in the equation for u2' there would be also a (u4)'
u2' = M^(-1) * (v - A*u2 - B*u1-Cu4-Du4')
and for this reason I don't know how I could do
Torsten
Torsten 2023-9-23
编辑：Torsten 2023-9-23
u1 = [theta;phi];
u2 = [dtheta/dt;dphi/dt];
u = [u1;u2]
du/dt = [du1/dt;du2/dt] = [u2;M^(-1)*(v-A*u2-B*u1)]
4 differential equations in the unknown vector of functions u = [theta;phi;dtheta/dt;dphi/dt]

请先登录，再进行评论。

Sam Chak
Sam Chak 2023-9-23
编辑：Sam Chak 2023-9-24
Ran in:
Hi @Gloria,
If there are time-dependent elements in the mass matrix, the idea is to define the generalized 'Mass' matrix in odeset(). The system
can be rearranged to
,
where and can be defined as a state vector :
.
In the 1st-order ODE form, the system can be expressed in terms of state variables as
.
and rewritten compactly in the state-space representation:
,
where is the identity matrix and ℳ is the generalized mass matrix of the state-space that contains time-dependent elements.
Here is an example since we don't know the parameters of the system.
% Parameters
lr = 1;
h = 1;
md = 1;
rd = 1;
ms = 1;
rs = 1;
w = 1;
% Generalized Mass matrix
M = @(t) [eye(2), zeros(2); zeros(2), [lr+2*md*h^2+(2*md*rd^2+ms*rs^2)*(cos(w*t))^2, 1/2*(2*md*rd^2+ms*rs^2)*(sin(w*t))^2; 1/2*(2*md*rd^2+ms*rs^2)*(sin(w*t))^2, lr+2*md*h^2+(2*md*rd^2+ms*rs^2)*(sin(w*t))^2]];
opt = odeset('Mass', M, 'MStateDependence', 'none');
% ode113 solver
tspan = linspace(0, 30, 3001);
x0 = [4 -3 -2 1];
[t, x] = ode113(@odefcn, tspan, x0, opt);
% Solution Plot
plot(t, x, 'linewidth', 1.5),
grid on, xlabel('t'), ylabel('\bf{x}')
legend({'$\theta$', '$\phi$', '$\dot{\theta}$', '$\dot{\phi}$'}, 'interpreter', 'latex', 'location', 'best', 'fontsize', 14)
function xdot = odefcn(t, x) % right-hand side of state-space
xdot = zeros(4, 1);
% Parameters
c = 1;
dc = 1;
w = 1;
md = 1;
rd = 1;
ms = 1;
rs = 1;
Iz = 1;
k = 1;
dk = 1;
h = 1;
% Matrices
C11 = c*dc^2 - w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
C12 = Iz*w + 2*w*(2*md*rd^2 + ms*rs^2)*(cos(w*t))^2;
C21 = - Iz*w - 2*w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
C22 = c*dc^2 + w*(2*md*rd^2 + ms*rs^2)*(sin(w*t))^2;
C = [C11, C12; C21, C22]; % damping matrix
K = diag([k*dk^2, k*dk^2]); % stiffness matrix
F = 2*md*rd*h*w^2*[cos(w*t); sin(w*t)]; % force matrix
% System
xdot(1:2) = x(3:4);
xdot(3:4) = F - C*x(3:4) - K*x(1:2);
end
  6 个评论
显示 4更早的评论
Torsten
Torsten 2023-9-27
编辑：Torsten 2023-9-27
x(1) = theta
x(2) = phi
x(3) = dtheta/dt
x(4) = dphi/dt
Thus
xdot(3:4) = d/dt ([dtheta/dt;dphi/dt]) = [d^2theta/dt^2;d^2phi/dt^2]
x(3:4) = [dtheta/dt;dphi/dt]
x(1:2) = [theta;phi]
so
[d^2theta/dt^2;d^2phi/dt^2] = M^(-1) * (F - C*[dtheta/dt;dphi/dt] - K*[theta;phi])
or
M * [d^2theta/dt^2;d^2phi/dt^2] + C*[dtheta/dt;dphi/dt] + K*[theta;phi] = F
Sam Chak
Sam Chak 2023-9-28
Hi @Gloria
If matrix remains on the left-hand side, then yes, and are still coupled.
However, when matrix is moved to the right-hand side of Eq. (1), it is essentially the same as solving a system of two algebraic equations:
,
,
where
.
In other words, and are considered to be fully decoupled in this form
,
which can be expressed in a relatively lengthy decoupled format:
So, the expression xdot(3:4) = M\(F - C*x(3:4) - K*x(1:2)) serves as a concise representation of the extensive decoupled form in MATLAB code, and it corresponds to Equation (2) mathematically.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by