How to solve 2nd order coupled differential equations?

38 ビュー (過去 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 日
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 件のコメント
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 日
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 ExchangeNumerical Integration and Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by