'...must return a column vector' errors when using ODE45

Hi all,

Here's my function code:

and here's my solution code with the stated error:

I've not included a large amount of the code in the function as I don't believe my errors lies in there but if it'd help you help me I can add it.

I'm aiming at a solution that provides a 16x300 matrix, i.e. the value of 16 variables contained in x over 300 time steps. But Matlab is saying it must return a column vector. Is this the nature of ode45 or can I adapt my code somehow? My function is rejecting the input arguments as shown. What should they be here? The independent variable is t which defines x and Q, while all the remaining terms are constant.

Thanks in advance.

EDIT: Full code posted below

[Star Strider]

That is difficult to read.

Your ‘xdot’ needs to be a column vector, that means ‘x’ also needs to be a column vector, not a (16x300) matrix. (I have no idea what ‘B*Q’ is.)

[Keelan Toal]

Sorry about the format. I need xdot to be a 16x300. Is there any way to do this? Perhaps ode45 is the wrong method? Trying to avoid having to solve 16 different ODEs. Q's components are defined in the 'if' loop shown so it's a 4x300 time-dependent matrix. Whereas A and B are 16x16 and 16X4 matrices of constants respectively.

[Walter Roberson]

Please post actual code. The image of the code is not clear enough for me to make out (I have vision problems.)

[Keelan Toal]

Function code:

function xdot = myode1(t,x)
h=0.05;lam=20;V=20; a=1.524; b=1.156; tau=(a+b)/V; w=(2*pi*V)/lam;
t=zeros(1,300);Q1=zeros(1,300);Q3=zeros(1,300);Q2=zeros(1,300);Q4=zeros(1,300); 
x=zeros(16,300);

for i=1:1:300
    t(i)=(i-1)/100;
    if t(i)<=2
        Q1(i)=(h/2)*(1-cos(w*(t(i))));
    else
        Q1(i)=0;
    end
    if (t(i)<=2.2 && t(i)>=0.2)
        Q3(i)=(h/2)*(1-cos(w*(t(i)-0.2)));
    else
        Q3(i)=0;
    end
    if (t(i)<=2+tau && t(i)>=tau)
        Q2(i)=(h/2)*(1-cos(w*(t(i)-tau)));
    else
        Q2(i)=0;
    end
    if (t(i)<=2+tau+0.2 && t(i)>=tau+0.2)
        Q4(i)=(h/2)*(1-cos(w*(t(i)-tau-0.2)));
    else
        Q4(i)=0;
    end
end

Kt=200000; Kp=0; Cp=0;K1=20000; K2=20000; K3=20000; K4=20000; C1=1000; C2=1000; C3=1000; C4=1000;
Mp=100; M=2160; M1=85; M3=85; M2=60; M4=60;
Ix=946; Iy=4140; W=1.45/2; Xp=0.234; Yp=0.375;

B1=[0,0,0,0,0,0,0,0,0,(Kt/M1),0,0,0,0,0,0]';
B2=[0,0,0,0,0,0,0,0,0,0,0,(Kt/M2),0,0,0,0]';
B3=[0,0,0,0,0,0,0,0,0,0,0,0,0,(Kt/M3),0,0]';
B4=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,(Kt/M4)]';
B=[B1,B2,B3,B4];
A1=[0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
A3=[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0];
A5=[0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0];
A7=[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0];
A9=[0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0];
A11=[0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0];
A13=[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0];
A15=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1];

A2=(-1/Mp)*[Kp,Cp,-Kp,-Cp,-Yp*Kp,-Yp*Cp,-Xp*Kp,-Xp*Cp,0,0,0,0,0,0,0,0];

A4=(-1/M)*[-Kp,-Cp,(K1+K2+K3+K4+Kp),(C1+C2+C3+C4+Cp),(W*(K1+K2-K3-K4)+Yp*Kp),(W*(C1+C2-C3-C4)+Yp*Cp),(-a*K1+b*K2-a*K3+b*K4+Xp*Kp),(-a*C1+b*C2-a*C3+b*C4+Xp*Cp),-K1,-C1,-K2,-C2,-K3,-C3,-K4,-C4];

A6=(-1/Ix)*[-Kp*Yp,-Cp*Yp,W*(K1+K2-K3-K4)+Yp*Kp,W*(C1+C2-C3-C4)+Yp*Cp,((W^2)*(K1+K2+K3+K4)-(Yp^2)*Kp),((W^2)*(C1+C2+C3+C4)-(Yp^2)*Cp),(-a*K1*W+b*K2*W-a*K3*W+b*K4*W-Xp*Yp*Kp),(-a*C1*W+b*C2*W-a*C3*W+b*C4*W-Xp*Yp*Cp),-K1*W,-C1*W,-K2*W,-C2*W,K3*W,C3*W,K4*W,C4*W];

A8=(-1/Iy)*[-Kp*Xp,-Cp*Xp,(-a*K1+b*K2-a*K3+b*K4-Xp*Kp),(-a*C1+b*C2-a*C3+b*C4-Xp*Cp),(-a*K1*W+b*K2*W+a*K3*W-b*W*K4-Xp*Yp*Kp),(-a*C1*W+b*C2*w+a*C3*W-b*W*C4-Xp*Yp*Cp),(a^2)*(K1+K2+K3+K4)-(Xp^2)*Kp,(a^2)*(C1+C2+C3+C4)-(Xp^2)*Cp,a*K1,a*C1,-b*K2,-b*C2,a*K3,a*C3,-b*K4,-b*C4];

A10=(-1/M1)*[0,0,-K1,-C1,-K1*W,C1*W,a*K1,a*C1,(Kt+K1),C1,0,0,0,0,0,0];

A12=(-1/M2)*[0,0,-K2,-C2,-K2*W,-C2*W,-b*K2,-b*C2,0,0,(Kt+K2),C2,0,0,0,0];

A14=(-1/M3)*[0,0,-K3,-C3,K3*W,C3*W,a*K3,a*C3,0,0,0,0,(Kt+K3),C3,0,0];

A16=(-1/M4)*[0,0,-K4,-C4,K4*W,C4*W,-b*K4,-b*C4,0,0,0,0,0,0,(Kt+K4),C4];

A=[A1;A2;A3;A4;A5;A6;A7;A8;A9;A10;A11;A12;A13;A14;A15;A16];
Q=[Q1;Q2;Q3;Q4];
xdot=A*x+B*Q;
end

So the values of Q look like this wrt time:

Solution code:

Tspan=0:0.01:2.99; IC=[zeros(1,16)]';
[T,X]=ode45(@(t,x) 'myode1(t,x)',Tspan,IC);

Resultant error:

Error using odearguments (line 90)
@(T,X)'MYODE1(T,X)' must return a column vector.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Optimisation (line 2)
[T,X]=ode45(@(t,x) 'myode1(t,x)',Tspan,IC);

And what I need is x and xdot to be 16x300 matrices, i.e. x1,x2...x16 for 300 instances of time.

[Walter Roberson]

[T,X]=ode45(@(t,x) myode1(t,x),Tspan,IC);

Notice no string!!

[Keelan Toal]

Thanks for the suggestion but it seems to produce the same result.

[Walter Roberson]

xdot = xdot(:);

right you compute xdot.

[Torsten]

A*x is 16x1 and B*Q is 16x300. Incompatible.

Best wishes

Torsten.

[Walter Roberson]

Ah. We could show how to expand the A*x to be 16x300 itself, but I suspect that the poster wants a result of length 16 not a result of length 16*300 ...

[Keelan Toal]

We have progress! But I get the error:

Error using odearguments (line 92)

Error using odearguments (line 92)
@(T,X)MYODE1(T,X) returns a vector of length 4800, but the length of initial conditions vector is 16. The vector returned
by @(T,X)MYODE1(T,X) and the initial conditions vector must have the same number of elements.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Optimisation (line 2)
[T,X]=ode45(@(t,x) myode1(t,x),Tspan,IC);

If I change the length of IC to 4800 (which is wrong), the could does run but x is a 16x300 of zeros.

[Torsten]

function xdot = myode1(t,x)

Here, t is 1x1 and x is 16x1.

So, you have to evaluate Q1, Q2, Q3 and Q4 at this special time t. If you do this, Q will be 4x1 and everything is fine.

Best wishes

Torsten.

Answer 1

[Titus Edelhofer]

Keelan,

I think there is a conceptual misunderstanding of what your myfun should do: it looks as if your 300 is related to 300 time steps that you want to do? The time stepping is done by the ode solver, inside myfun you should only work with 16x1 vectors, no loop on t (t is given by the ode solver as first input argument!).

Titus

[Keelan Toal]

That's exactly what I what I was trying to do. Thanks for pointing that out.

How would you suggest I input the constraints on Q1-4? I.e. when they are zero. Also I do not get a matrix for x, only xdot, is this correct?

[Titus Edelhofer]

You can leave the computation of Q1-4 more or less unchanged. Simply use the t that you get as parameter (without indexing (i)). Remove the for loop and just write

if t<=2
  Q1 = (h/2) * (1-cos(w*t));
else
  Q1 = 0;
end

Titus