Converting equations to first order system for ODE45

조회 수: 1 (최근 30일)
Simon Aldworth
Simon Aldworth . 2022년 11월 23일
댓글: Simon Aldworth . 2022년 11월 27일
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?

채택된 답변

Steven Lord
Steven Lord 2022년 11월 23일
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
Simon Aldworth
Simon Aldworth 2022년 11월 27일
Thank you.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

William Rose
William Rose 2022년 11월 23일
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
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
William Rose
William Rose 2022년 11월 23일
@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.

댓글을 달려면 로그인하십시오.


Help CenterFile Exchange에서 Equation Solving에 대해 자세히 알아보기




Community Treasure Hunt

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

Start Hunting!

Translated by