MATLAB: How to solve ODE with measured time series data

ode45

I have data from measurement with time step dt = 1/6600, which are stored in psi_dot and psi_ddot, both are 1651×1 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.

Best Answer

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);
Related Question