Matlab od45 question

Question

I am learning to use Matlab's ode45 to solve this equation and plot its time history:

$$I_gd\omega/dt + \omega \times (I_g\omega) = 0$$

It was suggested that I use:

omega_dot = fun(t,omega) ...

global I_G
I_G = I_g(:,:,m);
a = I_G*omega;
b = cross(omega,a);
omega_dot = -I_G/b;

I_g is a 3-by-3 inertia matrix predefined in another piece of code. My time vector is tspan = [0 100]. The initial state vector is omega = [1 0 0].

I understand that I need to set up my code in the form:

[t,y] = ode45(@vdpl,[0 20],[2 0]);
plot (t,y(:,1));

I need suggestions on how to make the information that I have work with the information that the examples for ode45 suggest. Any help would be greatly appreciated.

Answer 1

You're nearly there. First, don't use global variables. There's no need. And they're inefficient. The code below shows how to pass in parameters such as your I_G matrix (note that integrationfun now has a third input argument). Second, you didn't define m – it can be passed in just like I_G if it's constant (though then it'd make more sense to just pass in I_g(:,:,m)) or you can define it in your integration function. You can define one main function with a sub-function that describes the ODE:

function mainfun
tspan = [0 100];
omega0 = [1 0 0].';
I_G = ... % Define here, passed into integrationfun
[t,y] = ode45(@(t,y)integrationfun(t,y,I_G),tspan,omega0);
plot (t,y(:,1));

function omega_dot = integrationfun(t,omega,I_G)
m = ... % Define m
I_G = I_g(:,:,m);
a = I_G*omega;
b = cross(omega,a);
omega_dot = -I_G/b;

By the way, I'm not sure the last line of your integration function matches up with the differential equation you specified.