Problem with order of substitution with odeToVectorfield- Solving second order ode with ode45

Hello,

I have a 12 by 12 system of second order odes and I know how to make a first order problem out of it using odeToVectorField and calculate a solution using ode45. However, odeToVectorField seems to change the order of appearance of my variables. I defined a vector U and matrices M, K, C and a vector P are given.

syms u1z(t) u2z(t) u1x(t) u2x(t) u1y(t) u2y(t) u1rx(t) u2rx(t) u1ry(t) u2ry(t) u1tor(t) u2tor(t) 
U=[u1z(t); u2z(t); u1x(t); u2x(t); u1y(t); u2y(t); u1rx(t); u2rx(t); u1ry(t); u2ry(t); u1tor(t); u2tor(t)];

[V, V_subs]=odeToVectorField(M*diff(U, 2) + C*diff(U, 1) + K*U == P);

and V_subs is

V_subs =

      u2z
     Du2z
      u1z
     Du1z
      u1x
     Du1x
      u2x
     Du2x
     u2ry
    Du2ry
      u1y
     Du1y
      u2y
     Du2y
     u2rx
    Du2rx
     u1rx
    Du1rx
     u1ry
    Du1ry
    u1tor
   Du1tor
    u2tor
   Du2tor

As you can see, the order of appearance of my variables is changed compared to vector U. As a consequence, the solution I get from ode45 is changed. I mean it is still correct, but the order of my columns in the solution Matrix is changed. Is there any way to prevent odeToVectorField from changing order while doing its substitutions?

I already changed the column order of my solution matrix with a rather simple for-expression later on, but is there a more 'elegant' way?

Answer 1

[Star Strider]

Use the matlabFunction function to create an anonymous function or a function file. You can specify the order of the variables, and it saves you the trouble of putting the odeToVectorField symbolic code into executable code for the numeric ODE solvers.

From the documentation Specify Input Arguments for Generated Function (for a function file but this works for anonymous functions as well):

matlabFunction(r,'File','myfile','Vars',[y z x]);

You can also specify the output variables and their order. See Specify Output Variables.

Remember to put your independent variable as the first argument to the function to use it in ode45, if matlabFunction doesn’t do that for you.

EDIT — Added links to the online documentation.

[Tom Maar]

Hello, thank you for your help, but I don't really see how this solves my problem. I still have to rewrite my system of odes to a system of first order odes. I can do this either manually -which would be insane considering that it is a 12 by 12 system - or using odeToVectorField.

In fact, I do already use matlabfunction, in order make executable code out of symbolic code, just as you said. My further code is actually like this (inital values y0 are given):

syms u1z(t) u2z(t) u1x(t) u2x(t) u1y(t) u2y(t) u1rx(t) u2rx(t) u1ry(t) u2ry(t) u1tor(t) u2tor(t) 
U=[u1z(t); u2z(t); u1x(t); u2x(t); u1y(t); u2y(t); u1rx(t); u2rx(t); u1ry(t); u2ry(t); u1tor(t); u2tor(t)];

[V, V_subs]=odeToVectorField(M*diff(U, 2) + C*diff(U, 1) + K*U == P);

Vs=matlabFunction(V, 'vars', {'t', 'Y'});
tspan=[0 50];

[t, U]= ode45(Vs, tspan, y0);

I tried to specify the 'Outputs', but it did not really work. I mean the order of rows of V is already determined.

Did you want to say that I can avoid using odeToVectorField at all?

[Star Strider]

Since my approach doesn’t seem to be what you want to do, I would keep track of the variables in your calling code, and just use one variable, for example simply ‘u(1)’ to ‘u(24)’ (I believe I counted them correctly), in your ODE to present to ode45. You can assign them to their appropriate variables later, where necessary in your code.

What you want to do is similar to the ‘dynamically-created variables’ problem, never a recommended approach.

I would definitely use both odeToVectorField and matlabFunction.

Answer 2

[Uni Vang]

Hello!

I experience just the same issue. Hope support will reply to this topic. Because right now it feels just wrong.

Answer 3

[Victor Quezada]

Define your equations upside down, it worked for me. For example:

clear
syms x1(t) x2(t) f1(t) f2(t)
m=[2 3];c=[3 1 1];k=[3 2 1];f1=cos(2*t);f2=2*cos(3*t);
q=[diff(x2,2)==(f2-(c(2)+c(3))*diff(x2)+c(2)*diff(x1)...
    -(k(2)+k(3))*x2+k(2)*x1)/m(2);...
    diff(x1,2)==(f1-(c(1)+c(2))*diff(x1)+c(2)*diff(x2)...
    -(k(1)+k(2))*x1+k(2)*x2)/m(1)];
[V,Vsubs] = odeToVectorField(q)
M = matlabFunction(V,'vars', {'t','Y'});
[t,y] = ode23(M,[0 20],[0.2 1 0 0]);
figure % graficamos las soluciones de las ecs. difs. originales
subplot(2,1,1),plot(t,y(:,1)),title('solución 1'),ylabel('x(t) [m]')
subplot(2,1,2),plot(t,y(:,3)),title('solución 2'),xlabel('t [s]'),ylabel('x(t) [m]')

Answer 4

[Matthew Ediz Beadman]

Hi,

I stumbled on this question as I was having the same issue. The crude solution I have found is to define a wrapper function for your state derivative function which you get from the vector field, i.e. the function Vs in the following:

[V, V_subs]=odeToVectorField(M*diff(U, 2) + C*diff(U, 1) + K*U == P);
Vs=matlabFunction(V, 'vars', {'t', 'Y'});

What you want to do is define a function that:

• extracts the state variables from a correctly ordered state variable X
• reorders them into Y (which is the state variables ordered such that they can be passed to the function Vs)
• Pass Y into Vs to obtain Y_dot
• extract the derivative of the state variables from Y_dot
• reorder them into the ordering you want and define this vector as X_dot (which is the output of the wrapper function)

Then you can pass this wrapper function to your ODE solver, and it will accept the state variables in the order that you define it to, and not the order defined by odeToVectorField.

I have applied this to my project which has 8 state variables, and the solution looks like so:

% x, x_d, phi, phi_d, theta, theta_d, psi_, psi_d
X0 = [0,0,0,0,-0.1,0,pi+0.1,0];

[t,X] = ode45(@(t,X) dXdt(t,X,ode_fcn_euler), [0 2000], X0);

plots_4state(t,X,'euler ode');

function X_dot = dXdt(t,X,ode_fcn_euler)
x = X(1);
x_d = X(2);
phi = X(3);
phi_d = X(4);
theta = X(5);
theta_d = X(6);
psi_ = X(7);
psi_d = X(8);

% Y_state_order = [7,8,1,2,3,4,5,6];
%   psi_ Dpsi_ x Dx phi Dphi theta Dtheta
Y = [psi_, psi_d, x, x_d, phi, phi_d, theta, theta_d];
% use euler function
Y_dot = ode_fcn_euler(t,Y);

psi_d = Y_dot(1);
psi_dd = Y_dot(2);
x_d = Y_dot(3);
x_dd = Y_dot(4);
phi_d = Y_dot(5);
phi_dd = Y_dot(6);
theta_d = Y_dot(7);
theta_dd = Y_dot(8);

% inv_Y_state_order = [3,4,5,6,7,8,1,2]
X_dot = [x_d; x_dd; phi_d; phi_dd; theta_d; theta_dd; psi_d; psi_dd];
end

The above is quite ugly code, but has been written to be explicit. If you want the wrapper function to be more concise you can define a list that grabs the index of the state variables as they should be when reordered to Y, and then define a function that can invert that list to grab the correctly order state derivative X_dot from Y_dot.

Like so:

% x, x_d, phi, phi_d, theta, theta_d, psi_, psi_d
X0 = [0,0,0,0,-0.1,0,pi+0.1,0];

Y_state_order = [7,8,1,2,3,4,5,6];
inv_Y_state_order = invertStateOrder(Y_state_order); % [3,4,5,6,7,8,1,2]

[t,X] = ode45(@(t,X) dXdt(t,X,ode_fcn_euler,Y_state_order,inv_Y_state_order), [0 2000], X0);

function X_dot = dXdt(t,X,ode_fcn_euler, Y_state_order, inv_Y_state_order)
Y = X(Y_state_order);
Y_dot = ode_fcn_euler(t,Y);

X_dot = Y_dot(inv_Y_state_order);
end

function inv_order = invertStateOrder(state_order)
find = 1;
inv_order = zeros([1, numel(state_order)]);
i = 0;
while find <= numel(state_order)
    index = mod(i,numel(state_order))+1;
    if state_order(index) == find
            inv_order(find) = index;
            find = find+1;
    end
    i = i+1;
end
end