Error using ODE45. Function returns vector of length 1001 and initial conditions vector is 2.

I'm trying to solve an ODE but I'm getting the following error...

Error using odearguments
@(T,Q)[Q(2);(((21.*RHO.*A.*G)-(504.*E'.*I.*Q(1))-(504.*ETA'.*I.*Q(2)))./(RHO.*A.*L^4))] returns a
vector of length 1001, but the length of initial conditions vector is 2. The vector returned by
@(T,Q)[Q(2);(((21.*RHO.*A.*G)-(504.*E'.*I.*Q(1))-(504.*ETA'.*I.*Q(2)))./(RHO.*A.*L^4))] and the
initial conditions vector must have the same number of elements.

This is my code...

%ODE parameters
%A = A; already set above
rho = 1050;
g = 9.81;
L = 50e-03;
D = 2e-03;
I = 0.25*pi*(D/2)^4;
t_span = linspace(0,600,1000);
q0 = [0;0];

%Solving ODE

odefun = @(t,q) [q(2);(((21.*rho.*A.*g)-(504.*E'.*I.*q(1))-(504.*eta'.*I.*q(2)))./(rho.*A.*L^4))];

[t,q] = ode45(odefun,t_span,q0);

figure (4)
    plot(t,q(:,1),'-o')

Rather than creating a function I've created an anonymous function. The ODE is second order so I've created a system of linear equations. E and I are previously created 1x1000 vectors so I think I want my output q to be 1000 in size as well. Not sure if it's possible or not to solve an ODE with variables that are vectors this way or not. Displacement and velocity initial conditions are both 0.

Help?

[Jan]

"E and I are previously created 1x1000 vectors so I think I want my output q to be 1000 in size as well. Not sure if it's possible or not to solve an ODE with variables that are vectors this way or not."

The error message shows, that it is not possible. The output of the function to be integrated must have the same number of elements as the initial values, of course. I cannot imagine, what you expect as result for your code. The output of the integration is a matrix with the size tspan rows and the number of initial values as columns. A "size of 1000" might mean the number of time points, but why are the constants not scalars?

What do you want to calculate?

Answer 1

[Torsten]

You have to interpolate the values of E, I and eta to the time t supplied by the integrator and insert these interpolated values in your second-order ODE instead of the 1x1000 vectors. Use a function instead of a function handle to do this.

[Joshua Bowes-Smith]

Ok thanks that's helpful. As a result I've made the following function...

function dqdt = MyODE(t,q,E,eta,time,A,rho,g,L,I)
    E_intp = interp1(time,E,t);
    eta_intp = interp1(time,eta,t);
    dq1dt = q(2);
    dq2dt = (((21*rho*A*g)-(504.*E_intp'*I.*q(1))-(504.*eta_intp'*I.*q(2)))./(rho*A*L^4));
    dqdt = [dq1dt; dq2dt];
end

E and eta are the time-dependent 1x1000 vectors both calculated at the following time vector called 'time'...

tmin = 0;
tmax = 600;
N = 1000;
dt = (tmax - tmin)/(N-1);
time = tmin:dt:tmax;

A, rho, g, L, and I are all constant values in the dq2dt line. Apprently ode45 needs to output a vertical vectors/matrix so I tranposed E_intp and eta_intp. In my main code I used the function as follows...

t_span = [0 600];
q0 = [0;0];

[t,q] = ode45(@(t,q) MyODE(t,q,E,eta,time,A,rho,g,L,I),t_span,q0);

Have I created the function correctly? When I run it it gets stuck solving indefinitely. Can you help? Sorry, I'm pretty new to Matlab.

[Torsten]

You forgot I_intp.

[Joshua Bowes-Smith]

Sorry no I miswrote that originally. I is constant. It's just E and eta that need interpolating.

Any other insights?

[Torsten]

Any other insights?

No. Without knowing your arrays eta and E and the values of A,rho,g,L and I, there is nothing more I could say.

Maybe you could try ode15s instead of ode45.

And check your arrays for NaN or Inf values.

And take care that your "time" array covers the range from 0 to 600. If you interpolate outside the range of the "time" array, MATLAB returns NaN.

[Joshua Bowes-Smith]

I've attached the complete files. Would really appreciate if you could take a look.