How to solve ODE with measured time series data?

Question

Fan 2014-7-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/142001-how-to-solve-ode-with-measured-time-series-data

评论： Fan 2014-7-16

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Brian B 2014-7-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/142001-how-to-solve-ode-with-measured-time-series-data#answer_145292

编辑：Brian B 2014-7-16

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Brian B 2014-7-16

在 MATLAB Online 中打开

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 2014-7-16