Second order ODE with two time dependent data sets using ode45 and interp1.

function x = myODE(t,x)

    global d w 

    d = 60;
    w = 0.1;

            i = ['vdt' '3905' '.xls'] ;
            [num,raw]=xlsread(i) ;
            chestacc = num(1:1:610,4) ;     %first vector
            seatacc = num(1:1:610,7) ;      %second vector
            time = num(1:1:610,1);     %time corresponding to values in vectors

            fseat = interp1(time,seatacc,[1:1:610],'spline','extrap');
            fchest = interp1(time,chestacc,[1:1:610],'spline','extrap');

            x(1) = fchest*9.81/(w^2);

            dx = [x(2) ; fseat - 2*d*w*x(2) - (w^2)*fchest*9.81/(w^2)] ; %second order differential equation

end

Then using ode45 by inputting this into the command window:

[T,X] = ode45(@myODE,[-95:0.25:514],[0 0])

error: In an assignment A(I) = B, the number of elements in B and I must be the same. (for both of these lines: x(1) = fchest*9.81/(w^2); dx = [x(2) ; fseat - 2*d*w*x(2) - (w^2)*fchest*9.81/(w^2)] ;

Am I missing steps previous to this? Are there mistakes in the way I am using these functions? Thank you for your time and assistance.

Answer 1

[Jan]

How could you get two errors, when Matlab stops after the first one already?

In your code fchest and fseat are vectors with 610 elements. Then you cannot use them to define the scalars x(1) and x(2).

But it is rather strange, that you use interp1 to get the interpolation at the original time points, such that the output equals the input exactly. I guess you want to interpolate at the specific time point t instead.

It will slow down your program massively, when you import the XLS-file in each integration step. Importing this once is enough, when you store the output persistently:

function x = myODE(t,x)
d = 60;  % No need to make them global
w = 0.1;

persistent time chestacc seatacc
if isempty(time)  
  file = 'vdt3905.xls';
  num = xlsread(file) ;
  chestacc = num(1:610,4) ;     %first vector
  seatacc  = num(1:610,7) ;      %second vector
  time     = num(1:610,1);     %time corresponding to values in vectors
end

fseat = interp1(time,seatacc, t,'spline','extrap');
fchest = interp1(time,chestacc, t, 'spline','extrap');

% x(1) = fchest*9.81/(w^2);  % Not used?!
dx = [x(2) ; fseat - 2*d*w*x(2) - (w^2)*fchest*9.81/(w^2)] ;
end

The "t" in the interp1 call is a bold guess.