ode45 keeps saying that my function must return a column vector

So I am working on code for a project at school, and I cant seem to figure out why ode45 keeps telling me that my function must return a column vector. I tested my function on it own outside of ode45 and the value it returned was indeed a column vector.

My main code is

clc
close all

Cars_in = [0 0 0 0 0 4 10 12 8 15 25 53 56 48 24 22 32 18 12 0 0 0 0];
Cars_out = [0 0 0 0 0 0 4 6 0 0 12 47 48 47 12 16 32 34 20 12 49 0 0];
Cars = Cars_in + Cars_out;
Force = Cars*1818*9.81; %N
hours=0:length(Cars)-1;
random_num = 3600;
Time=0:(length(Cars)-1)*random_num-1;


F_cont=interp1 (hours,Force,Time/random_num,'spline');

Qin_avg = 5.36e-7; %L/s
Qin = abs(randn(1,length(Cars))*Qin_avg);
Q_cont = interp1 (hours,Qin,Time/random_num,'spline');


for i = 1:length(F_cont)
    if F_cont(i)<0
        F_cont(i) = 0;
    end
end

x0 = [0 0 0 0 0 0 0];

[t,y]= ode45(@(t,y)dy(t,y,F_cont,Q_cont),Time,x0);

and my function being called is

function dy=spline_interp(t,x,F_cont,Q_cont)

  
%%inertial elements
rj = 1; %m
rm = 0.5; %m
H = 0.1; %m
density_steel = 8050; %kg/m3
mj = density_steel*H*pi*rj^2; %kg
mm = density_steel*H*pi*rm^2;
J1 = 0.5*mj*rj^2;
Jm = 0.5*mm*rm^2;

%%torsional spring
k1 = 1872897.912/1000; %N-m/deg
k2 = 2809346.869/1000; %N-m/deg

%%viscous friction
bt = 50; 

%%generators
ke1 = .26; %deg/(s*Volt)
km = 3.59; %Nm/Amp
ke2 = 0.6;
kt2 = 0.000327741; %m^3/rev

%%Circuit
R = 20000; %ohms
Rl = 380;
L = 4300e-6; %Henry
C = 10000e-6; %farad

%%fluid elements
alpha = 2;
density_water = 997; %kg/L
delta_x = 6.096; %m
A  = 0.1^2;
If = alpha*density_water*delta_x/A;
bulk = 2.2e9; %N/m^2
volume = 0.15^2*pi*0.3048;
Cf = volume/bulk;

%%linear system
A = [0 0 -1/J1 0 0 0 0;
     0 -bt/Jm 1/Jm -km 0 0 0;
     k1+k2 -(k1+k2)*x(2) 0 0 0 0 0;
     0 ke1/L 0 -R/L -1/L 0 ke2*kt2;
     0 0 0 1/C -1/(C*Rl) 0 0;
     0 0 0 0 0 0 -1/Cf;
     0 0 0 0 0 1/If 0];
 
 B = [rj/J1 0; 0 0; 0 0; 0 0; 0 0; 0 1/Cf; 0 0];
 
  U = [F_cont; Q_cont];

dy=A*x+B*U;
end

Any insight would be greatly appreciated, but as the project is due on Thursday, November 29, there is not much time left.

Thank You

[Torsten]

The name of your function is "spline_interp", not "dy".

[Julien Gibbons]

the function is called spline_interp, but the file is called dy

Answer 1

[Torsten]

I think you mean

function main
Cars_in = [0 0 0 0 0 4 10 12 8 15 25 53 56 48 24 22 32 18 12 0 0 0 0];
Cars_out = [0 0 0 0 0 0 4 6 0 0 12 47 48 47 12 16 32 34 20 12 49 0 0];
Cars = Cars_in + Cars_out;
Force = Cars*1818*9.81; %N
hours=0:length(Cars)-1;
random_num = 3600;
Time=0:(length(Cars)-1)*random_num-1;

F_cont=interp1(hours,Force,Time/random_num,'spline');
Qin_avg = 5.36e-7; %L/s
Qin = abs(randn(1,length(Cars))*Qin_avg);
Q_cont = interp1(hours,Qin,Time/random_num,'spline');

for i = 1:length(F_cont)
   if F_cont(i)<0
        F_cont(i) = 0;
   end
end

x0 = [0 0 0 0 0 0 0];

[t,y]= ode45(@(t,y)spline_interp(t,y,Time,F_cont,Q_cont),Time,x0);

end

function dy=spline_interp(t,x,Time,F_cont,Q_cont)

%%inertial elements

rj = 1; %m
rm = 0.5; %m
H = 0.1; %m
density_steel = 8050; %kg/m3
mj = density_steel*H*pi*rj^2; %kg
mm = density_steel*H*pi*rm^2;
J1 = 0.5*mj*rj^2;
Jm = 0.5*mm*rm^2;

%%torsional spring

k1 = 1872897.912/1000; %N-m/deg
k2 = 2809346.869/1000; %N-m/deg


%%viscous friction

bt = 50; 

%%generators

ke1 = .26; %deg/(s*Volt)
km = 3.59; %Nm/Amp
ke2 = 0.6;
kt2 = 0.000327741; %m^3/rev

%%Circuit

R = 20000; %ohms
Rl = 380;
L = 4300e-6; %Henry
C = 10000e-6; %farad

%%fluid elements

alpha = 2;
density_water = 997; %kg/L
delta_x = 6.096; %m
A  = 0.1^2;
If = alpha*density_water*delta_x/A;
bulk = 2.2e9; %N/m^2
volume = 0.15^2*pi*0.3048;
Cf = volume/bulk;

%%linear system

A = [0 0 -1/J1 0 0 0 0;
     0 -bt/Jm 1/Jm -km 0 0 0;
     k1+k2 -(k1+k2)*x(2) 0 0 0 0 0;
     0 ke1/L 0 -R/L -1/L 0 ke2*kt2;
     0 0 0 1/C -1/(C*Rl) 0 0;
     0 0 0 0 0 0 -1/Cf;
     0 0 0 0 0 1/If 0];

B = [rj/J1 0; 0 0; 0 0; 0 0; 0 0; 0 1/Cf; 0 0];

U = [interp1(Time,F_cont,t,'spline'); interp1(Time,Q_cont,t,'spline')];

dy=A*x+B*U;

end

[Julien Gibbons]

I appreciate the help Torsten, but I must ask, have you tested this code, and if so, approximately how long did it take to run?

[Julien Gibbons]

Unfortunatey this did not solve my issue, as the code is still taking forever to run. Which in itself is strange as I have a code that is very similar to my current one that works fine.

[Steven Lord]

You have constants in your function that vary by many orders of magnitude, which makes me suspect your ODE is stiff. Try one of the stiffer solvers listed on this documentation page and see if performance improves. See this post on Cleve's blog for more information about stiffness.

[Julien Gibbons]

Thank you so much steven, this solved my issues exactly. Torsten, Jan, thank you for all your help. Only issue now though is that the program keeps experiencing failure at a certain value of t

example

Warning: Failure at t=3.972685e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.411381e-14) at time t.

> In ode23s (line 401)

In odetest (line 22)

[Torsten]

Then you should inspect the results up to this time and see what is going wrong.

Answer 2

[Jan]

This is a job for the debugger. Type this in te command window:

dbstop if error

Run the code again until Matlab stops at the error. Now check the failing line. I guess Torsten hit the point: You are running the function dy, but adjust spline_interp. Then:

[t, y]= ode45(@(t,y) spline_interp(t, y, F_cont, Q_cont), Time, x0)

[Julien Gibbons]

Thank you for your help, unfortunately this did not help solve my problem. I have tried renaming my function and its file to be consistent to no avail. I have also tried making F_cont and Q_cont global variables. Everytime I run to code I get the same error

Error using odearguments (line 93)

@(T,Y)SPLINE_INTERP(T,Y,F_CONT,Q_CONT) must return a column vector.

Error in ode45 (line 115)

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

Error in calc_ode45 (line 28)

[t,y]= ode45(@(t,y)spline_interp(t,y,F_cont,Q_cont),Time,x0);

I believe the issue is that I am passing t into ode45, but not using it within my created function. If I try to preform the interpolation inside dy while calling ode45, the code runs without giving me an error, but takes forever. With only 22 points the code was still running 8hrs later, we need 79200 points. I have tried to use t inside ode45 as a means of determining the index of F_cont and Q_cont that are required but to no avail.