How to solve 2nd order coupled differential equations?

38 Ansichten (letzte 30 Tage)
Gloria
Gloria am 23 Sep. 2023
Kommentiert: Sam Chak am 28 Sep. 2023
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

Antworten (2)

Torsten
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)
  2 Kommentare
Gloria
Gloria am 23 Sep. 2023
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 am 23 Sep. 2023
Bearbeitet: Torsten am 23 Sep. 2023
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]

Melden Sie sich an, um zu kommentieren.


Sam Chak
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
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
Sam Chak am 28 Sep. 2023
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.

Melden Sie sich an, um zu kommentieren.

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!

Translated by