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.

*No products are associated with this question.*

Answer by Matt Tearle
on 6 Sep 2012

Accepted answer

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

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)

Show 2 older comments

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.)

Opportunities for recent engineering grads.

## 0 Comments