Using ODE45 to solve a state space system.

Hello there!

I stucked in how to use ode45. My problem is the following:

I have a ODE that i want to solve, the only difference is that my initial conditions are vector 3x1.

    *function xdot = double_int2(t, y) 

        xd1 = y(2,:) % xdot = v
        xd2 = temp2 + Td % temp2 and Td are vector 1x3

    end*

and the other function

    *
     function [T,Y] = call_double_int2()
         x01 = [10 0 0 ];
         v01 = [1 0 0];
         t_span = [0 5];

         [T,Y]= ode45(@double_int2, t_span, [x01 v01])

      end*

So, I don't know how to implent in wat that MatLab understand it`s a row vector, I tried to declare the funcion as double_int2(t,y(2,3)), but it always take it as a element.

[Francisco]

Let me give more details about the problem.

I have 4 odes, they are the following:

xdot = v; vdot = g * e3 + Td;

Rdot = R x w; J * wdot = tau + w + J*w;

Where: 1) e3, Td. w and w are vectors 3x1; 2) J - inertia matrix 3x3;

x0 ( position ), v0 ( linear velocity ) and w0 (angular velocity) are expressed in cartesian coordinates thats why they are vectors 3x1.

Now that i gave some more details let me answer you.

1) All the examples I saw they use the initial condition as a element and that is the reason I am really confused, since I was oriented to use ode45 to solve it.

2) Thanks for idea of declaring xdot.

3) I used y(2,:) because as I said the problem initial conditions are in cartesian coordinates.

Thanks for you help.

[Sara]

From your description, you don't have four ode's but 12 (3 components * 4 variables), which is totally doable in matlab. Let's take the system you have posted originally for simplicity. It will have to be:

x = y(1:3); %not used but added for clarity
v = y(4:6);
xd1 = v;
xd2 = temp2 + Td;
xdot = [xd1,xd2];

since your v is in y(4:6). Is this clearer? If you can't still do it, post ALL the inputs plus whatever code you have.

[Francisco]

% Position

x01 = [10 0 0 ]; x02 = [0 10 0 ]; x03 = [0 0 0 ]; x04 = [0 -10 0 ]; x05 = [-10 0 0 ]; x0 = [ x01; x02; x03; x04; x05 ];

% Linear Velocity

v01 = [1 0 0]; v02 = [0 1 0]; v03 = [0 0 0]; v04 = [0 -1 0]; v05 = [-1 0 0]; v0 = [ v01; v02; v03; v04; v05 ];

%% Parameters global alpha n kq Kw e3 b m A J g omega

n = 1; % Number of Vehicles

alpha = 3*n ; %%%% Verify

kq=50; Kw=10*eye(3); %dimension of matrix eye

e3 = [ 0 0 1 ] ;

b = 2; %%%% Verify

g = 9.806;

A = [0 1 0 0 0; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 1; 0 0 0 1 0]; % Adjacency Matrix

m = 2; % mass

J = 1.2416*eye(n,n); omega = ones(1,3); % Vector that makes Ti no null

%% Thrust Vector

aux = zeros(1,3);

sum = zeros(1,3);

for i=1:n for j=1:n

        if i == j

            continue
        end

        xij = x0(j,:)-x0(i,:); %%x0 -> current x
        vij = v0(j,:)-v0(i,:); %%v0 -> current v
        aux = A(i,j) + ((xij)*tanh(norm(xij))/norm(xij) +  (vij*tanh(norm(vij))/(b*norm(vij)) ) );
        sum = aux + sum;
    end
    Ti = aux+alpha*omega;
    Td(i,:) = Ti;  % Recording Ti in a matrix [5,3]
    %[t_temp,y_temp] = ode45(@double_int2, [0 10], [x0(i,:) v0(i,:)]);

    aux = zeros(1,3); % clean aux for next iteration
end

%% t_span = [0, 5];

[t,y]= ode45(@double_int2, t_span, [x01 v01]);

%% Now the function %% %%

function xdot = double_int2(t, y) % x is the vector [ position , velocity ]

    global Td m g e3

    x=y(1:3);
    v=y(4:6);

    xd1 = v; % xdot = v

    xd2 = m*g*e3 + Td;  % m * vdot = m * g * e3 + Td

    xdot = [xd1, xd2 ];

[Francisco]

My xdot is that cconcatenated matrix, but ode45 expect it to return a 1D vector as I think that is the only problem now. Thanks for you attention.

Answer 1

[Sara]

First, the size of the array of initial conditions has to be the same of the number of equations you want to solve. So, in your case, you either have 6 equations ([x01 v01] is an array 1 by 6) or you want only 1 element from x01 and y01. Given your ODE function, I think it's the second, but you tell us.

Then:

function xdot = double_int2(t, y)

expects that somewhere you define xdot, maybe as xdot = [xd1,xd2]. In addition, y is a 1D array, so replace y(2,:) with y(2).

You will also need to pass temp2 and Td to double_int2. You can do:

[T,Y]= ode45(@(x)double_int2(x,temp2,Td), t_span, [x01 v01])
function xdot = double_int2(t, y,temp2,Td)