How to apply ode solvers to ode systems written in matrix form skipping expansion of their right- hand-sides? Show a representative example.
1 次查看(过去 30 天)
显示 更早的评论
Let DY=A(t)Y be a system of linear odes in matrix form, where Y is an n*1-matrix and A is n*n-matrix. How to apply ode-solvers directly to this system without expending it?
Assume now that a system is nonlinear:
DY=f(t,Y); Y is an n*1-matrix. The question is how to apply ode-solvers without expending this equation.
0 个评论
回答(1 个)
Aykut Satici
2014-8-19
编辑:Aykut Satici
2014-8-19
If you have the right hand side of your ODE as an n-by-1 vector, then you can just assign the output of the ode function to this vector.
As an example, below I have implemented the simulation of a rotating rigid body with a constant angular velocity about the body x-axis. Note that the equations of motion of a rotating rigid body, as expressed in the body frame, are
R'(t) = R(t)*Omega(t);
I*omega'(t) = I*omega(t) x omega(t)
where R is the rotation matrix, Omega is the 3-by-3 skew-symmetric angular velocity matrix, omega is the representation of the angular velocity as a 3-vector (a vector in R^3), and I is the moment of inertia of the body expressed in the body frame. The primes denote time differentiation and "x" denotes the cross product in R^3.
In particular, the function "odefun" is doing what you are asking. It creates the vector field "f" and assigns it as the output of "odefun".
function odeSystemSample()
% Simulation of a rotating rigid-body
%%Set-up the problem and solve
par.I = diag([1,2,3]); % Inertia of the body
q0 = [reshape(eye(3),9,[]); 1; 0; 0]; % Initial condition, body is
% rolling with constant angular
% velocity
ti = 0; tf = 2;
opt = odeset('AbsTol',1.0e-07,'RelTol',1.0e-07);
[t,q] = ode45( @odefun, [ti:0.001:tf], q0 ,opt, par);
end
function dq = odefun(t,q,par)
% This is how you assign an n-by-1 vector to the output of the function
% that defines the differential equations.
R = reshape(q(1:9),3,3);
w = q(10:12);
dq = zeros(12,1);
dq(1:9) = reshape(R*phi(w,0),9,[]); % This creates a 9-by-1 vector
dq(10:12) = par.I \ cross(par.I*w, w); % This is a 3-by-1 vector
end
function x = phi(y,flag)
% This function is problem specific; you do not need to pay attention to it.
% This is the isomorphism between the Lie algebra so(3) and R^3.
% If the flag is on, then the function goes from so(3) to R^3.
% If the flag is off, then the function goes from R^3 to so(3).
if flag
x = [-y(2,3); y(1,3); -y(1,2)];
else
x = [0, -y(3), y(2); y(3), 0, -y(1); -y(2), y(1), 0];
end
end
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!