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.
采纳的回答
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 个评论
Fan
2014-7-16
Brian B
2014-7-16
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 MATLAB 的更多信息
Translated by 
