How to apply ode solvers to ode systems written in matrix form skipping expansion of their right- hand-sides? Show a representative example.

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.

Answer 1

[Aykut Satici]

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