Calculate displacement of a mdof sytem with ODE45?

I am modelling a 2D half-car model in SimMechanics and MapleSim (to compare these two multibody software packages) for my project. My supervisor also want me to make a analytical model for comparing the outputs of both programs. So I have derived the mass matrix and stiffness matrix of this half-car model and I now want to calculate the displacement in vertical direction of the car. My supervisor told me I could do this with ODE45. I do not have any experience with MATLAB, so I have no idea what to do now.

This is my equation:

Mx"+Kx= 0 where M and K are a 4x4 matrix and they change every timestep(they are dependent of angle phi)

I have four generalised coordinates, so I have the initials conditions y0 which is a column of 8 scalars(displacement and velocity at the start)

My supervisor told me something like this: y = [x; xdot] so ydot=[xdot; xddot] where xddot is equal to: inv(M)-Kx, this is clear to me. This only has to be done for all the 4 generalised coordinates. so y = [x1; x1dot; phi; phidot; x2; x2dot; x3; x3dot] and for ydot the same. Am I right?

But what next? What do I have to do in MATLAB to generate the displacement using the ODE45 solver... Maybe I am a bit vague, so if you want to know something, just ask.

Answer 1

[Matt Tearle]

You're definitely on the right track. Write your equations in matrix form My' = -Ky where y is the 8-element vector you described. In this setting, M and K will be 8-by-8 matrices (which are functions of t and y).

Then in MATLAB

f = @(t,y) -[matrix K goes here]*y;
m = @(t,y) [matrix M goes here];

[t,y] = ode45(f,...,odeset('Mass',m));

You may also need to set other options related to the mass matrix. See doc odeset for details.

Answer 2

[Bas]

Matt, thank you for your answer, I have done what you wrote, but MATLAB gives the following error:

??? Error: File: newtest.m Line: 26 Column: 5 Expression or statement is incorrect--possibly unbalanced (, {, or [.

My m-file:

x1=1
xdot=0
phi=0
phidot=0
x2=1
x2dot=0
x3=1
x3dot=0

y=[x1; x1dot; phi; phidot; x2; x2dot; x3; x3dot]

M=[5 2 1 0;
    2 5 1 1;
    1 1 2.5 0;
    0 1 0 2.5]

K=[50 5 1 1;
    5 50 1 1;
    1 1 25 0;
    1 1 0 25]

f = @(t,y) -[K]*y;
m = @(t,y) [M];

[t,y] = ode45(f,...,odeset('Mass',m))    (this is line 25, the error said something about line 26..)

In this test file, I have made the M and K matrices independent, so they do not change in time now, because I first want to run this model correctly and the step to make the matrices dependent of 'phi' will follow afterwards. Do you have any idea what I am doing wrong?

[Bas]

I have changed line 25 with this ';' on the end. So:

[t,y] = ode45(f,...,odeset('Mass',m));

and then it gives this error message:

??? Error: File: newtest.m Line: 26 Column: 1 This statement is incomplete.

Do you have any idea? Sorry, I think the answer will be very easy, but I have no idea what to do, this is my first MATLAB assignment at the university..

[Matt Tearle]

Sorry, I was skipping bits of the code (doc style). Check the documentation for the full syntax for ode45. There should be more stuff (a timespan vector and initial conditions) in place of the "...". MATLAB treats a literal "..." as "this line is incomplete, see the next line for the rest" (and ignores the rest of the line). That's what's causing the error.

[Bas]

That is what I already thought, I will change my code and will inform you about the results. Thanks again

Answer 3

[Bas]

It is not working at the moment, I added a timespan and the initial conditions to the ode45 function, but I do not think that is the real problem. For instead, look to the error messages:

??? Error using ==> mtimes Inner matrix dimensions must agree.

Error in ==> @(t,y)-K*y

Error in ==> odearguments at 98 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 172 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> newtest at 28 [t,y] = ode45(f,[0 10],y,odeset('Mass',m));

The matrix multiplication -K*y is not possible, because K is still a 4x4 matrix and y is a 8x1 column. Furthermore, I do not know what the other errors mean..

This is my M-file till now:

x1=1;
x1dot=0;
phi=0;
phidot=0;
x2=1;
x2dot=0;
x3=1;
x3dot=0;
y=[x1; x1dot; phi; phidot; x2; x2dot; x3; x3dot];
M=[5 2 1 0;
 2 5 1 1;
 1 1 2.5 0;
 0 1 0 2.5];
K=[50 5 1 1;
  5 50 1 1;
  1 1 25 0;
  1 1 0 25];
f = @(t,y) -K*y;  
m = @(t,y) M;   
[t,y] = ode45(f,[0 10],y,odeset('Mass',m));
figure
plot(t,y)

[Matt Tearle]

Don't worry about the intermediate error messages. That's just the full trace from your initial script to the final error. The error is, of course, happening when you call ode45, so you're seeing an error in a function inside a function inside a function...

The problem is what you've identified: K is 4-by-4 and y is 8-by-1. Note in my original answer I said K and M would have to be 8-by-8 because you need to write the equations as My' = -Ky where y is the 8-element vector [x1; x1dot; phi ...]

That means that y' = [x1dot; x1dotdot; phidot; phidotdot; ...] = [y2; 50*x1 + 5*phi + x2 + x3; ...] = [y2; 50*y1 + 5*y3 + y5 + y7;...]. So K should be something like

[0  1 0  0 0  0 0  0]
[50 0 5  0 1  0 1  0]
[0  0 0  1 0  0 0  0]
[5  0 50 0 1  0 1  0]
[0  0 0  0 0  1 0  0]
[1  0 1  0 25 0 0  0]
[0  0 0  0 0  0 0  1]
[1  0 1  0 0  0 25 0]

(I think, based on what you have for K currently). And similarly for M.

[Bas]

Hi Matt, I have changed my model a little bit, with some help of my supervisor. I have only one question left, how do I make my mass and stiffness matrix time depended? For example phi changes every time step. So the matrices have to be calculated again every step. Thanks for your time! Here is my m-file which I now have, it are two files, one normal m-file and one function file:

    % %initial conditions:
%q = [u_car ; phi ; u_1 ; u_2];
u_car_0 = -0.2;
u_1_0 = -0.1;
u_2_0 = -0.1;
phi_0 = 0.0
u_car_0dot = 0;
phi_0dot = 0;
u_1_0dot= 0;
u_2_0dot= 0;
y0=[u_car_0; phi_0; u_1_0; u_2_0; u_car_0dot; phi_0dot; u_1_0dot; u_2_0dot];
% solve diff vgl
[t, y]=ode45(@odefun,[0 10],y0);
figure
plot(t,y(:,1))

And the function file:

    function dydt=odefun(t,y)
M=[5 2 1 0;            
   2 5 1 1;
   1 1 2.5 0;
   0 1 0 2.5];
K=[50 5 1 1;
    5 50 1 1;
    1 1 25 0;
    1 1 0 25];
dydt(1:4,1)= -inv(M)*K*y(5:8,1);
dydt(5:8,1)= y(1:4,1);

The M and K matrix are numbers in this case, but in my real assignment they are k_spring1*cos(phi) etcetera. So what do I have to prescribe for 'phi' for example to make this a time depended variable?

[Matt Tearle]

Inside odefun you have t and y, so you can calculate anything based on those values at the current timestep. It looks like phi is y_2, so

K = [50*cos(y(2)) 5 1 1;...]

Also, don't use inv(M). Use the Mass option in ode45. But if you really don't want to do that, at least use -M\(K*y(5:8)) instead of inv(M).

I also think you may have your indices for y back-to-front somewhere. In the calling script, it seems like y_2 is phi and y_6 is phi'. If that's the case, then, in odefun, dydt(2) would be phi', so dydt(2) = y(6). And dydt(8) would be phi'' (calculated from M\K*y). In other words, in the script, y(1:4) are the primitive variables and y(5:8) are the derivatives; but in the function it's the other way around.

[Bas]

Hey Matt, sorry for the late reaction! I already solved the problems, it is working fine now, the issue about the indices is completely right! My supervisor was a bit in a hurry when he told me how to do this, and I already thought this has to be wrong! So I changed it myself after I posted this here.

Only one question left, can you explain why you should not use inv(M), my supervisor told me to definitely use inv(M) and not the options you give, because that was my first idea, if I am honest.... I thought okay fine, my supervisor should be right, and I did not asked why I should use inv(M).. But now I am very confused, because of your remark.

Thanks for the help!

[Matt Tearle]

Your supervisor may have been thinking about just getting it working, and using inv is simple, conceptually and practically. Furthermore, in your example, the inverse is probably pretty painless.

But inv is, in general, slow and numerically sensitive. There are almost always better ways to do it, in terms of accuracy and efficiency. If an alternative is offered (such as the Mass property), you can bet that it's better.

(But, again, "better" can be a subjective term. If your goal is just to make this work, the simplest implementation may be the best.)

Answer 4

[Bas]

M=[5*cos(phi) 2 1 0;            
 2 5 1*cos(phi) 1;
 1 1*cos(phi) 2.5 0;
 0 1 0 2.5];

I mean something like this but than much more elements with a time depedent variable. So what do you have to fill in for phi(t) to make this matrix time dependent?