Problem with order of substitution with odeToVectorfield- Solving second order ode with ode45
조회 수: 6 (최근 30일)
이전 댓글 표시
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?
댓글 수: 0
답변 (4개)
Star Strider
2016년 6월 2일
편집: Star Strider
2016년 6월 2일
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]);
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.
댓글 수: 2
Star Strider
2016년 6월 2일
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
2017년 4월 7일
Hello!
I experience just the same issue. Hope support will reply to this topic. Because right now it feels just wrong.
댓글 수: 0
Victor Quezada
2019년 10월 29일
편집: Victor Quezada
2019년 10월 29일
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]')
댓글 수: 0
Matthew Ediz Beadman
2024년 3월 16일
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
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Equation Solving에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!