How to solve 2nd order coupled differential equations?
38 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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

0 Kommentare
Antworten (2)
Torsten
am 23 Sep. 2023
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)
Sam Chak
am 23 Sep. 2023
Bearbeitet: Sam Chak
am 24 Sep. 2023
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 Kommentare
Torsten
am 27 Sep. 2023
Bearbeitet: Torsten
am 27 Sep. 2023
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
am 28 Sep. 2023
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.
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!