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

9 views (last 30 days)
Tom Maar on 2 Jun 2016
Edited: Victor Quezada on 29 Oct 2019
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?

Star Strider on 2 Jun 2016
Edited: Star Strider on 2 Jun 2016
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.

Tom Maar on 2 Jun 2016
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 on 2 Jun 2016
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.

Uni Vang on 7 Apr 2017
Hello!
I experience just the same issue. Hope support will reply to this topic. Because right now it feels just wrong.

Victor Quezada on 29 Oct 2019
Edited: Victor Quezada on 29 Oct 2019
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]')