Converting equations to first order system for ODE45

Hi,

I'm trying to convert a set of equations of motion for a simple vehicle model into a set of first order differential equations for use with ODE45. This model is more complicated than previous ones I’ve used, with some derivatives appearing multiple times in the equations to be solved. This seems to give a problem with the order in which variables are presented in the code.

If I simplify the equations to make this discussion easier by removing all the constants, then they are:

1. v' + r + phi'' = alpha1 + alpha2
2. r' + phi'' = alpha1 + alpha2
3. phi'' + v' + r + r' + phi = phi'
4. alpha1' = v + r + phi + alpha1 + input
5. alpha2' = v + r + phi + alpha2

Since there are 6 derivatives, I think there are 6 states that need initial conditions:

1. v
2. r
3. phi
4. phi'
5. alpha1
6. alpha2

The 6 equations in the ODE are:

1. dx(1) = phi’
2. dx(2) = r + phi'' + alpha1 + alpha2 % v’
3. dx(3) = phi'' + alpha1 + alpha2 % r’
4. dx(4) = v' + r + r' + phi + phi' % phi’’
5. dx(5) = v + r + phi + alpha1 + input % alpha1’
6. dx(6) = v + r + phi + alpha2 % alpha2’

Running this returns the error “Index in position 1 exceeds array bounds (must not exceed 1).” because dx(2) contains a reference to dx(4) (phi’’) which isn’t defined as yet. Which ever way I arrange equations dx(2), dx(3) and dx(4) I will have a problem with precedence.

How should I arrange the equations to avoid such problems?

Thanks.

Answer 1

[Steven Lord]

Let's look at your equations and your definition for your state vector.

1. v' + r + phi'' = alpha1 + alpha2
2. r' + phi'' = alpha1 + alpha2
3. phi'' + v' + r + r' + phi = phi'
4. alpha1' = v + r + phi + alpha1 + input
5. alpha2' = v + r + phi + alpha2

So using the standard ODE notation of M*y' = f(t, y) if we define y as [v; r; phi; phi'; alpha1; alpha2] we have that y' is [v'; r'; phi'; phi''; alpha1'; alpha2'].

Written using a vector-vector product your first equation is [1 0 0 1 0 0]*y' = -y(2) + y(5)+y(6). So [1 0 0 1 0 0] is the first row in the Mass matrix for this problem and the first element of the vector returned from your ODE function f(t, y) is -y(2)+y(5)+y(6). Continue in this way building up the whole 6-by-6 Mass matrix and the function f(t, y).

Then if you're not sure how to specify a Mass matrix in a call to an ODE solver, look at one of the examples in the Summary of ODE Examples and Files section on this documentation page that says it uses the Mass option and use it as a guide for implementing the function to solve your system.

[Simon Aldworth]

Thank you to @William Rose, @Torsten and @Steven Lord for the huge amount of help. I have gone with the mass matrix approach as there is a neatness and brevity to it that my brain can deal with. The model now runs but I'd like to ask a few questions just to check my approach:

1. The example given here of the baton thrown into the air says that the function for the right side of the equation (f(t,q)) takes three inputs, like the mass matrix function (two required, one optional). However, my model has an input and I need to pass that to this function. As per my previous models, I've added two additional inputs; one is the input vector and the other is the integration time vector. The function then interpolates the input time history for the current integration time. It seems to run ok, so I assume what I've done is ok?
2. I've vectorised the function for the right side of the equations (i.e. used y(2,:) instead of y(2)) so that I can call the function again and return the derivative time histories. Again, it seems to work so I guess that's ok?
3. When filling the vector that is to be solved in f(t,q) the baton example has each element on a new line. I assume the spaces at the beginning of each line are telling matlab that each line is a new element? I assume, therefore, that I can't use spaces to separate the parts of each element as that will confuse matlab? If I use a comma between elements, can I then use spaces between parts of elements? For example, currently it looks like:

ydot = [a+b

c+d

e+f]

So I assume I cannot do this:

ydot = [a + b

c + d

e + f]

But perhaps I can do this:

ydot = [a + b,

c + d,

e + f]

Thanks again for everyone's help.

Regards.

[William Rose]

@Simon Aldworth,

1. As I understand your q.1, you have something like

function dydt = f(t,y,P1,P2)

That should work fine if you do

[t,y] = ode45(@(t,y) f(t,y,P1,P2),tspan,y0,opts);

assuming you have passed the mass matrix in opts.

2. You say you use y(2,:) instead of y(2). Of course y function f() is a local variable name. You could pass x or y or q to f(t,y,P) when you call it. I'm not sure if you are refering to y inside the function or outside. You said it works, so that's what matters. When I do simulations, I spend a lot of time writing test scripts for various special cases, to gain confidence that the code is doing what I expect.

3. The following all produce the same result:

dydt=[a+b
      c+d
      e+f]
dydt=[a + b
      c + d
      e + f]
dydt=[a+b;c+d;e+f] 
dydt=[a + b; c + d; e + f]

[Simon Aldworth]

Thank you!

[Simon Aldworth]

Hi,

I have a few more questions to ask, partly due to my confusion and partly because some parts of my code weren't working as well as I thought.

In previous uses of functions that were then passsed to ODE45, it was possible to obtain the derivatives by passing the results from the solver back to the function. I thought this might be possible here (hence my previous Q2 referring to vectorising the function) but since I'm using the mass matrix method the function f(t,y) only contains the right hand side of the equations. This Stack Overflow example also suggests passing the solver results back to the function, but I can't see how the variables on the left hand side (i.e. the derivatives factored by the mass matrix) get included in the results returned from the function. For example, my first equation is [1 0 0 1 0 0]*y' = -y(2) + y(5)+y(6), and y = [v; r; phi; phi'; alpha1; alpha2] and y' = [v'; r'; phi'; phi''; alpha1'; alpha2']. If I pass the solver results back to the function then the first element of the return would be -y(2) + y(5)+y(6) and so would be missing the contributions from v' and phi''. Is that correct? If so, what's the best way for me to obtain the derivatives at each time step? Or is it simpler to just calculate the derivatives of y using a function like 'gradient'?

I constructed the mass matrix from my five equations, following the advice from @Steven Lord. For the sixth equation, I put a 1 in element [6,3] and put y(4) in the last row of f(t,y), so saying phi' = phi'. Is that correct? I want to check this because the model is behaving well in all degrees of freedom except the one associated to phi, phi', etc.

Thanks.

[Steven Lord]

If you call the ODE solver with one output to obtain a solutions struct, call deval to evaluate that solution at other times. You can call it with two outputs to also return "the first derivative of the numeric solution produced by the solver."

[Simon Aldworth]

Thank you.

Answer 2

[William Rose]

@Simon Aldworth,

I think you are corect tht a state vector with six elements, as you have defined above, will suffice. If we call the state vector x, then its six elements are:

x=[v, r, phi, phi', alpha1, alpha2]

v, r, phi are related by equations that are not obviously separateable into explicit equaitons for their derivatives. The derivatives are defined implicitly by a set of algebraic equaitons. Therefore you should read about ode15i() and the example for the Robertson system. I think you will be able to solve your system with this approach.

With the ode15i approach, we dfine a vector which is the derivative of the state vector:

xp=[v', r', phi', phi'', alpha1', alpha2']

You have explicit equations for the derivatives of aplha1, alpha2 and phi':

dx(3) = x(4)                       % v'
dx(5) = x(1)+x(2)+x(3)+x(5)+input  % alpha'=v+r+phi+alpha1+input  
dx(6) = x(1)+x(2)+x(3)+x(6)        % alpha2'=v+r+phi+alpha2 

The equaitons above can be rewritten as

0=xp(3)-x(4)                       % v'
0=xp(5)-(x(1)+x(2)+x(3)+x(5)+input)  % alpha'=v+r+phi+alpha1+input  
0=xp(6)-(x(1)+x(2)+x(3)+x(6))        % alpha2'=v+r+phi+alpha2

I recommend that you rearrange your equaitons 1,2,3 so that phi'' only appears in one of the equations. You had:

1. v' + r + phi'' = alpha1 + alpha2
2. r' + phi'' = alpha1 + alpha2
3. phi'' + v' + r + r' + phi = phi'

You could write these as

1a. v' + r - r' = 0

2a. (phi')' - alpha1 - alpha2 + r' = 0

3a. alpha1 + alpha2 + v' + r + phi - phi' = 0

These can be rewritten in terms of x() and xp(): (recall x=[v, r, phi, phi', alpha1, alpha2])

0=xp(1)+x(2)-xp(2);                         % v'+r-r'=0
0=xp(4)-x(5)-x(6)+xp(2);                    % phi''-alpha1-alpha2+r'=0
0=x(5)+x(6)-xp(2)+xp(1)+x(2)+x(3)-xp(3);    % alpha1+alpha2+v'+r+phi-phi'=0

Now you follow the example for the Robertson problem in the ode15i help: combine the equaitons into a function which returns a vector. Each row of the vector should be an equaiton that equals zero.

function res = simonidae(t,y,yp)
%SIMONIDAE: Simon's implicit differential algebraic equation
%Assume input=sin(t); adjust that part as desired.
res = [xp(3)-x(4);                          % v'
       xp(5)-(x(1)+x(2)+x(3)+x(5)+sin(t));  % alpha'=v+r+phi+alpha1+input  
       xp(6)-(x(1)+x(2)+x(3)+x(6));         % alpha2'=v+r+phi+alpha2
       xp(1)+x(2)-xp(2);                      % v'+r-r'=0
       xp(4)-x(5)-x(6)+xp(2);                 % phi''-alpha1-alpha2+r'=0
       x(5)+x(6)-xp(2)+xp(1)+x(2)+x(3)-xp(3)] % alpha1+alpha2+v'+r+phi-phi'=0
end

Then you will set options, find a consistent initial conditions x0 and xp0, set options, set the time span, and finally, call ode15i to solve the system:

[t,x] = ode15i(@simonidae,tspan,x0,xp0,options);

Good luck.

[William Rose]

@Simon Aldworth, check all my equations and code carefully. It is alway possible that I made a typo somewhere.

[William Rose]

Note that simonidae() has six independent equations, each of which equals zero. This equals the number of elements of the state vector. That is as it should be.

[Simon Aldworth]

@William Rose: thanks for this comprehensive response. I do have one concern though: I removed all the constants just for clarity. If I put them back in, then the substitutions you recommend approx. half way down can't be done. The constants are not the same for each variable in each equation. For example, subtracting the (alpha1 + alpha2 - phi'') in equation 2 from the (alpha1 + alpha2 - phi'') in equation 1 would not be zero. Sorry if I misled you.

[William Rose]

@Simon Aldworth,

[edit: I forgot to include a closing bracket in funciton simonidae. I have added it now.]

You're welcome.

You are referring to equations

1. v' + r + phi'' = alpha1 + alpha2
2. r' + phi'' = alpha1 + alpha2
3. phi'' + v' + r + r' + phi = phi'

which I rearranged. I rearranged them because I suspect the solver will work better with a simplified set of equations. But you do not have to rearrange them, and you can still apply the approach which I recommended earlier.

1. 0 = v' + r + phi'' - alpha1 - alpha2
2. 0 = r' + phi'' - alpha1 - alpha2
3. 0 = phi'' + v' + r + r' + phi - phi'

Recall x=[v, r, phi, phi', alpha1, alpha2], xp=[v', r', phi', phi'', alpha1', alpha2']

0 = xp(1)+x(2)+xp(4)-x(5)-x(6)
0 = xp(2)+xp(4)-x(5)-x(6)
0 = xp(4)+xp(1)+x(2)+xp(2)+x(3)-x(4)

I have a feeling that in the last equation, it is important to use "-x(4)" rather than "-xp(3)". Now you use the three equations above in simonidae(), instead of the ones in my previous posting:

function res = simonidae(t,y,yp)
%SIMONIDAE: Simon's implicit differential algebraic equation
%Assume input=sin(t); adjust that part as desired.
res = [xp(3)-x(4);                      % v'
    xp(5)-(x(1)+x(2)+x(3)+x(5)+sin(t)); % alpha'=v+r+phi+alpha1+input  
    xp(6)-(x(1)+x(2)+x(3)+x(6));        % alpha2'=v+r+phi+alpha2
    xp(1)+x(2)+xp(4)-x(5)-x(6);         % v'+r+phi''-alpha1-alpha2
    xp(2)+xp(4)-x(5)-x(6);              % r'+phi''-alpha1-alpha2
    xp(4)+xp(1)+x(2)+xp(2)+x(3)-x(4)]   % phi''+v'+r+r'+phi-phi'
end

Try something like that.

[Torsten]

Maybe of help to get the system explicit in the derivatives:

syms xp1 xp2 xp4 x5 x6 x2 x3 x4
eqn1 = xp1 + x2 + xp4 - x5 - x6 ==0;
eqn2 = xp2 + xp4 - x5 - x6 == 0;
eqn3 = xp4 + xp1 + x2 + xp2 + x3 - x4 == 0;
sol = solve([eqn1,eqn2,eqn3],[xp1,xp2,xp4])

sol = struct with fields:
    xp1: x4 - x3 - x2 - x5 - x6
    xp2: x4 - x3 - x5 - x6
    xp4: x3 - x4 + 2*x5 + 2*x6

[William Rose]

@Simon Aldworth,

@Torsten and @Steven Lord are giving you excellent suggestions. They are such good resources on this site.

The mass matrix approach of @Steven Lord is an alternative way of specifying and solving implicit differential algebraic equations. It works when the implicit equaitons are linear in the derivative of the state vector. And your equations are.

When the implicit equations are linear, then they are solvable, i.e. the explicit equaitons can be found. The solution may not be very practical, which is why the mass matrix approach can be useful. @Torsten shows that in your case, the set of three equations has a nice simple solution. Therefore, using @Torsten's result, you can write six explicit equations for the derivative of the state vector, and solve in the standard way, with ode45() for example.

See here for a nice explanation of implicit, explicit, mass matrix, and which solvers work with what.

You said you have given a simplified version of the equations, with some constants removed. Maybe the explicit solution of @Torsten will not be so simple when the true equations are used.

My approach and the approaches of @Steven Lord and @Torsten should give the same results in the end.