MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# Calculate displacement of a mdof sytem with ODE45?

Asked by Bas on 6 Sep 2012

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.

## Products

No products are associated with this question.

Answer by Matt Tearle on 6 Sep 2012

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 by Bas on 6 Sep 2012
Edited by Sean de Wolski on 6 Sep 2012

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 on 6 Sep 2012

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 on 6 Sep 2012

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 on 6 Sep 2012

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

Answer by Bas on 7 Sep 2012

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 on 12 Sep 2012

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 on 17 Sep 2012

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 on 18 Sep 2012

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 by Bas on 10 Sep 2012
```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?