Help with solving a system of coupled equations of motions in time domain?

Hi,

I've posted a similar post here, and I learned much through that post and comments. Now, I'm building on it and the system of equations has become a bit more complex again. I will explain it here again. Following were the set of equations I had, as I mentioned in the earlier post (linked).

But now, suppose I have another time-varying function added to each equation above. For instance, each DoF noted above includes another function f_1(t) (in fact, my system includes several of these, but I'm using one here for explanation). For example, as below.

So, I understand that these additional "functions" affects the system of equations at each time step. For example, f_1(t) here is the excitation loads due to air, as the object moves. i.e., one has to model it separately, and add a function/script to calculate that at the respective time step. Also, f_1(t) depends on the state variables which I am trying to solve.

For example, following what I had posted earlier (above linked), the equation system now changed to,

function sdot = learnF(A, B, s, data) % Equations to solve

sdot = [s(2) + f_1(inputs); % Added arbitrarily only to explain!!
        data.m*s(4)*s(12) + Ctn*A(2,2,1)*s(4)*s(12) + A(2,4,1)*s(8)*s(12) + A(2,6,1)*s(12)^2;
        s(4);
        -data.m*s(2)*s(12) - A(1,1,1)*s(2)*s(12) - A(2,4,1)*s(8)
        s(6);
        -data.m*data.g;
        s(8);
        data.I(4,4)*s(10)*s(12);
        s(10);
        -data.I(4,4)*s(8)*s(12);
        s(12);
        (data.I(4,4)-data.I(6,6))*s(8)*s(10) + A(1,1,1)*s(2)*s(4) ...
        - A(2,2,1)*s(2)*s(4) - A(2,4,1)*s(2)*s(8) - A(2,6,1)*s(2)*s(8) - A(6,4,1)*s(8)];
end

I added f_1(inputs) just arbitrarily (without deriving things) and I do know the placing could be wrong. The point is, now I have time-varying function within the function, which I have to solve using an ode solver. So my questions are;

1. Did I understand the approach correctly here? If yes, how can I model this f_1(t) to use inputs from sdot? (I can feel that I am wrong here).
2. I can, instead, use a for loop (time loop) and solve the equations step by step, but then I cannot use, for example, ode45 because I cannot solve ODE for one time step using that. Is there a work around for this?
3. IF I can still use ode45 and solve this, can I get outputs of f_1(t) after the simulation?

My sincere apologies for the long and possibly confusing post, but I'm kind of helpless here. I'd appreicate any feedback!

Answer 1

[Torsten]

How did you come up with the sdot vector for the equations so far ?

If you use a mass matrix, you simply have to add the additional force terms to the sdot vector.

A force term f1(t) in equation 1, e.g., has to be reflected as

sdot(2) = data.m*s(4)*s(12) + Ctn*A(2,2,1)*s(4)*s(12) + A(2,4,1)*s(8)*s(12) + A(2,6,1)*s(12)^2 + f1(t)

A force term f2(t) in equation 2 has to be reflected as

sdot(4) = -data.m*s(2)*s(12) - A(1,1,1)*s(2)*s(12) - A(2,4,1)*s(8) + f2(t)

and so on.

[Jake]

Hi @Torsten, Thank you for the response!

The sdot vector is derived following what I have done in the previous post, and with the help on Matlab example, "Solve Equations of Motion for Baton Thrown into Air". I use a mass matrix and so far, my mass matrix is still the same (although, I understand that these have to be changed accordingly), as I have noted in the earlier post. i.e., My mass mtarix is,

function M = learnMass(A, B, data) % Mass matrix function
m = data.m;
I = data.I;
M = zeros(12,12);
M(1,1)  = 1; M(3,3) = 1; M(5,5) = 1; M(7,7) = 1; M(9,9) = 1; M(11,11) = 1;
M(2,2)  = m + A(1,1,1);
M(2,6)  = A(3,1,1);
M(2,10) = A(5,1,1);
M(4,4)  = m + A(2,2,1);
M(4,8)  = A(4,2,1);
M(4,12) = A(6,2,1);
M(6,2)  = A(1,3,1);
M(6,6)  = m + A(3,3,1);
M(6,10) = I(5,5) + A(5,3,1);
M(8,8)  = I(4,4) + A(4,4,1);
M(10,2) = A(1,5,1);
M(10,6) = A(3,5,1);
M(10,10)= A(5,5,1);
M(12,4) = A(2,6,1);
M(12,8) = A(4,6,1);
M(12,12)= I(6,6) + A(6,6,1);
end

I agree that the f1(t) and f2(t) terms have to be reflected in the respective sdot lines as you mentioned. I just added the f1(t) in sdot(1) just to explain, but I agree it was a poor explanation.

But my question is how to model the f1, f2 parts. What I meant by f1(t) is that f1 is a function that depends on time - which I also wanted to implicitly mean that it depends on the velocity and/or displacement (which are the state variables that I get as sdot). I'm not sure if I make much sense here, so I'll try to explain further.

Suppose the EOM I'm solving are a motion of a car, and the module f1 represents an autopilot system. Then, at a particular time-step (t = t), autopilot system needs the inputs, which are the solutions of the EOM (velocity, displacement etc.) at the previous time-step (t = t-1). Further, f1 is modelled by another equation, which I believe I have to write a separate script. Am I missing something?

[Torsten]

Use an if-statement to adapt the inputs fi to the events that happen at certain times t.

Or integrate up to the times when events happen, stop the integration and restart again from the last time step of the previous calculation and with the new settings for the fi.

[Jake]

Thanks. I don't quite follow the first alternative, but the second one should be fine to implement.

[Torsten]

For the first (quick-and-dirty) alternative:

if t < t0
    f1 = ...;
else if t >= t0 && t <= t1
    f1 = ...;
else if ...
    ...
else
    ...
end

But the second way is recommended.