Error when I run ODE45 function

Could I know where I might be wrong in ode.m as I am getting following message as error when I run the ODE file. F1 and F2 are function of time.

Error using odearguments (line 93)

OUTPUT must return a column vector.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Ode (line 11)

[t, x] = ode45(@Output,tspan,IC);

State space model: (Output.m)

function [sdot] = Output(t,x)

Ze = zeros(3);
I = eye(3);

M = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
C = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
MA = M+Ad;
Q = C+Cm;


Fwt = 2.793977550000000e+04;    %Wind force on turbine tower

Fg = -79086000;                 %Gravitational force
Fbuoy = 7.844814740000000e+07;  %Buoyancy force
Fp = 2.712318560000001e+06;     %Heave force
Fmx = -334731.8545;             
Fmz = -3517000;                 

t = 0:0.1:10

for i = 1:length(t)
    
    Fh = hydro(t(i))
    F1(1,i) = Fmx+27939.6+6.5*10^5+Fh(1,i);
    F2(1,i) = Fmz+Fg+Fbuoy+Fp;
    F3(1,i) = -112510430.2+3.44*10^6+266.5*10^5+(Fh(1,i)*18);

 A = [Ze, I; -inv(MA)*Q, Ze];
 B = [Ze; inv(MA)];
 F = [F1;F2;F3];
  
 Svec = [x(1);x(2);x(3);x(4);x(5);x(6)];
% 
 sdot = A*Svec + B*F;
end

end 

ODE.m:

%Initial Condition 
x(10) = 0; x(20) = 0; x(30) = 0;
x(40) = 0; x(50) = 0; x(60) = 0;

IC = [x(10);x(20);x(30);x(40);x(50);x(60)];

%time span
t0 = 0; tf = 10;
tspan = t0:0.1:tf;

[t, x] = ode45(@Output,tspan,IC);
IC = [IC;x(end)];
%x(1) = x(:,1);
%x(2) = x(:,2);
%x(3) = x(:,3);
%x(4) = x(:,4);
%x(5) = x(:,5);
%x(6) = x(:,6);

%Plot Results 
figure(1), clf 
plot(t,x(:,1)), xlabel('time(s)'), ylabel('surge(m)')
title ('Surge vs Time')

[Star Strider]

I cannot run your code:

Unrecognized function or variable 'hydro'.

Fh = hydro(t(i))

so I cannot check to see if ‘A’, ‘B’ and ‘F’ are matrices. If they are, ‘sdot’ will likely not be a column vector.

[Satish Jawalageri]

Thanks for your reply. Here is the code:

hydro.m:

function [Fh] = hydro(t)

Cd = 0.6;
Ca = 0.9699;
R = 4.7;

t=0:0.1:10;

for i=1:length(t)
    vel = compute_wavevelocity(t(i));
    acc = compute_waveacceleration(t(i));
    fun = sym((0.5*997*Cd*2*R*abs(vel)*vel)+(Ca*997*pi*(R^2)*acc)+(pi*(R^2)*997*acc));
    Fh(1,i) = double(int(fun,0,-120));
    
end
plot(t,Fh), xlabel('time(s)'), ylabel('Fh')
title ('Fh vs Time')

end

compute_wavevelocity.m:

function velocity = compute_wavevelocity(t)
       T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
       vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega*t));
       
       velocity = mean(vel,1);
       
end

compute_waveacceleration.m:

function acceleration = compute_waveacceleration(t)
        T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
        acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega*t));
        
        acceleration = mean(acc,1);
               
    end

[Star Strider]

Before we go further with this, check to see if ‘A’, ‘B’ and ‘F’ are matrices. If they are, ‘sdot’ will likely not be a column vector.

The easiest way to do that is ether to add size calls without a closing semicolon so that they print to the Command Window, or just leave off the closing semicolon for ‘sdot’.

[Satish Jawalageri]

Yes, 'A', 'B' and 'F' are matrices. And I am getting sdot in order of 6x101 matrix.

[Star Strider]

Those are not column vectors.

You will have to decide how best to solve that, since I have no idea what you are doing.

[Satish Jawalageri]

Actually, sdot is column vector for each value of time (ie 0:0.1:10). I am trying to run the ODE45 code in such a way that the @output in ODE45 should take the column vector of sdot corresponding to the particular time.

[Walter Roberson]

Why are you looping over time there? why not just process for the one time that was passed into the function? The processing for any one t in the loop does not appear to depend on previous iterations.

[Satish Jawalageri]

Till output.m file I have used looping for getting values corresponding the time. So, do I need to loop for ODE to process according to time?If so, how can I do?

Thanks.

Answer 1

[Gifari Zulkarnaen]

You might see my example here to solve equation of motion using ODE:

https://www.mathworks.com/matlabcentral/fileexchange/74348-linear-dynamic-analysis-of-mdof-system-using-ode

[Satish Jawalageri]

I cant see any matlab code in zip file.

Thanks.

[Gifari Zulkarnaen]

I have revised the file

Thank you for your notice