Converting equations to first order system for ODE45

9 views (last 30 days)
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.

Accepted Answer

Steven Lord
Steven Lord on 23 Nov 2022
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.
  6 Comments
Steven Lord
Steven Lord on 26 Nov 2022
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."

Sign in to comment.

More Answers (1)

William Rose
William Rose on 23 Nov 2022
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.
  6 Comments
Torsten
Torsten on 23 Nov 2022
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
William Rose on 23 Nov 2022
@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.

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!