How to solve ODE with measured time series data?

I have data from measurement with time step dt = 1/6600, which are stored in psi_dot and psi_ddot, both are 1651x1 double. My time span is 0.25 second. I want to use the time dependent data to solve the below ODE.

function [z_dot] = MYODE(t,z,c,pd,pdd) % z(1) = theta % z(2) = thetadot

% Constants:

c1 = c(1);

c2 = c(2);

c3 = c(3);

c4 = c(4);

c5 = c(5);

% thetadot

z_dot(1,1) = z(2);

% thetadot_dot

z_dot(2,1) = c1.*pd.*pd*sin(2*z(1))...

            +c2.*pdd*cos(z(1))...

            +c3.*sign(pd).*pd.^2*(3.48*cos(z(1))).*(0.82*abs(sign(pd)*z(1)+pi/2)/pi+0.05)...

            +c4*abs(z(2))*z(2)...

            +c5*z(1);

end

I have my solver step up like this:

step = 1650; % N/A

t0 = 0; % Start time (sec)

t1 = 0.25; % Final time (sec)

t = linspace(t0 , t1 , (step+1));

[t,z1] = ode45(@(t,z) MYODE(t,z,c0,psi_dot,psi_ddot),t,z0);

where c0 are the constants passed into the function, psi_dot and psi_dot are the time series data.

It keeps throwing error at me saying Subscripted assignment dimension mismatch at z_dot(2,1) = c1.*pd.*pd*sin(2*z(1))...

Any idea where went wrong? I am not sure how the time series data passed into the function is indexed.

Answer 1

[Brian B]

The subscripted assignment dimension error is caused when you try to assign a vector of length 1651x1 to a 1x1 subarray of z_dot. The problem is that ode45 does not know that psi_dot and psi_ddot are to be interpreted as a time series, so it computes

   c1.*pd.*pd*sin(2*z(1))...
      +c2.*pdd*cos(z(1))...
      +c3.*sign(pd).*pd.^2*(3.48*cos(z(1))).*(0.82*abs(sign(pd)*z(1)+pi/2)/pi+0.05)...
      +c4*abs(z(2))*z(2)...
      +c5*z(1);

where pd and pdd are the full 1651x1 vectors psi_dot and psi_ddot.

To correct this, you need to compute pd and pdd from the data and from t. One way is as follows:

First, edit the definition of MYODE to accept an additional input tpsi which is a vector of times corresponding to the elements in psi_dot and psi_ddot. Then interpolate (I have used linear interpolation, but you could use nearest-neighbor interpolation or something else as appropriate) to get pd and pdd:

function [z_dot] = MYODE(t,z,c,tpsi,psi_dot,psi_ddot)
    % z(1) = theta % z(2) = thetadot

      % time-varying data
      pd  = interp1(tpsi, psi_dot, t, 'linear');
      pdd = interp1(tpsi, psi_ddot, t, 'linear');

      % Constants:
      c1 = c(1);
      c2 = c(2);
      c3 = c(3);
      c4 = c(4);
      c5 = c(5);

      % thetadot
      z_dot(1,1) = z(2);

      % thetadot_dot
      z_dot(2,1) = c1.*pd.*pd*sin(2*z(1))...
          +c2.*pdd*cos(z(1))...
          +c3.*sign(pd).*pd.^2*(3.48*cos(z(1))).*(0.82*abs(sign(pd)*z(1)+pi/2)/pi+0.05)...
          +c4*abs(z(2))*z(2)...
          +c5*z(1);
  end

Then, in your main script or function, you just need to pass in tpsi:

    tpsi = (0:1/6600:0.25)'; 
    step = 1650; % N/A
    t0 = 0; % Start time (sec)
    t1 = 0.25; % Final time (sec)
    t = linspace(t0 , t1 , (step+1));
    [t,z1] = ode45(@(t,z) MYODE(t,z,c0,tpsi, psi_dot, psi_ddot),t,z0);

[Fan]

This worked.

But I still do not understand, the pd constructed form the interpolation is exact the same as the original data series. Where does the tpsi come into play when the solver calculates the zdot(2,1)?

[Brian B]

No, pd is not the same as psi_dot. The latter is a vector, while the former is a scalar. The interpolation looks at tpsi and psi_dot as a time series and interpolates at the scalar time t to find pd.

To see this, just remove the semicolon to print out the value of pd at every (t,z) pair where ode45 evaluates MYODE:

   pd  = interp1(tpsi, psi_dot, t, 'linear')

You will see that every time, pd is just a singe number.

[Fan]

I think I understand know. Thanks!!