`ode45` 和 tspan 错误尝试访问

Question

我正在使用 ode45 求解二阶微分方程。时间跨度是根据txt文件中的数字多少来确定的，因此，时间跨度定义如下

i    = 1;
t(i) = 0;
dt   = 0.1;
numel(theta_d)
while ( i < numel(theta_d) )
   i = i + 1;
   t(i) = t(i-1) + dt;
end

现在时间元素不应超过txt的大小（即numel(theta_d)）。在 main.m 中，我有

x0 = [0; 0];
options= odeset('Reltol',dt,'Stats','on');
[t, x] = ode45('ODESolver', t, x0, options);

和 ODESolver.m header 是

function dx = ODESolver(t, x)

如果我运行代码，我会收到此错误

Attempted to access theta_d(56); index out of bounds because numel(theta_d)=55.

Error in ODESolver (line 29)
theta_dDot  = ( theta_d(i) - theta_dPrev ) / dt;

为什么 ode45 没有固定时间跨度？

编辑：这是完整的代码

main.m

clear all
clc

global error theta_d dt;
error = 0;


theta_d = load('trajectory.txt');




i    = 1;
t(i) = 0;
dt   = 0.1;
numel(theta_d)
while ( i < numel(theta_d) )
  i = i + 1;
  t(i) = t(i-1) + dt;
end

x0 = [pi/4; 0];
options= odeset('Reltol',dt,'Stats','on');
[t, x] = ode45(@ODESolver, t, x0, options);

 %e = x(:,1) - theta_d; % Error theta


plot(t, x(:,2), 'r', 'LineWidth', 2);
title('Tracking Problem','Interpreter','LaTex');
xlabel('time (sec)');
ylabel('$\dot{\theta}(t)$', 'Interpreter','LaTex');
grid on

和ODESolver.m

function dx = ODESolver(t, x)

persistent i theta_dPrev

if isempty(i)
    i = 1;
    theta_dPrev = 0;
end

global error theta_d dt ; 

dx = zeros(2,1);

%Parameters:
m = 0.5;       % mass (Kg)
d = 0.0023e-6; % viscous friction coefficient
L = 1;         % arm length (m)
I = 1/3*m*L^2; % inertia seen at the rotation axis. (Kg.m^2)
g = 9.81;      % acceleration due to gravity m/s^2


% PID tuning
Kp = 5;
Kd = 1.9;
Ki = 0.02;

% theta_d first derivative
theta_dDot  = ( theta_d(i) - theta_dPrev ) / dt;
theta_dPrev = theta_d(i);


% u: joint torque
u = Kp*(theta_d(i) - x(1)) + Kd*( theta_dDot - x(2)) + Ki*error;
error = error + (theta_dDot - x(1));

dx(1) = x(2);
dx(2) = 1/I*(u - d*x(2) - m*g*L*sin(x(1)));

i = i + 1;
end

这是错误

Attempted to access theta_d(56); index out of bounds because numel(theta_d)=55.

Error in ODESolver (line 28)
theta_dDot  = ( theta_d(i) - theta_dPrev ) / dt;

Error in ode45 (line 261)
    f(:,2) = feval(odeFcn,t+hA(1),y+f*hB(:,1),odeArgs{:});

Error in main (line 21)
[t, x] = ode45(@ODESolver, t, x0, options);

Answer 1

这里的问题是因为你有离散时间点的数据，但是ode45需要能够在你的时间范围内的任何时间点计算导数。一旦解决了问题，它就会将结果插回到您想要的时间点。所以它会计算导数比你指定的时间点多很多次，因此你的 i 计数器根本不起作用。

由于您有离散数据，因此继续 ode45 的唯一方法是将 theta_d 插值到任何时间 t。您有一个值列表 theta_d 对应于时间 0:dt:(dt*(numel(theta_d)-1))，因此要插入到特定时间 t，请使用 interp1(0:dt:(dt*(numel(theta_d)-1)),theta_d,t)，我将其转换为匿名函数以提供theta_p 在给定时间 t

的内插值

那么你的导数函数看起来像

function dx = ODESolver(t, x,thetaI)

dx = zeros(2,1);

%Parameters:
m = 0.5;       % mass (Kg)
d = 0.0023e-6; % viscous friction coefficient
L = 1;         % arm length (m)
I = 1/3*m*L^2; % inertia seen at the rotation axis. (Kg.m^2)
g = 9.81;      % acceleration due to gravity m/s^2

% PID tuning
Kp = 5;
Kd = 1.9;
Ki = 0.02;

% theta_d first derivative
dt=1e-4;
theta_dDot  = (thetaI(t) - theta(I-dt)) / dt;
%// Note thetaI(t) is the interpolated theta_d values at time t

% u: joint torque
u = Kp*(thetaI(t) - x(1)) + Kd*( theta_dDot - x(2)) + Ki*error;
error = error + (theta_dDot - x(1));

dx=[x(2);  1/I*(u - d*x(2) - m*g*L*sin(x(1)))];    
end

并且您必须在使用 [t, x] = ode45(@(t,x) ODESolver(t,x,thetaI, t, x0, options);.

调用 ode45 之前定义 thetaI=@(t) interp1(0:dt:(dt*(numel(theta_d)-1)),theta_d,t);

我从 ODESolver 中删除了一些内容并更改了导数的计算方式。

请注意，我无法对此进行测试，但它应该可以帮助您。

`ode45` 和 tspan 错误尝试访问

`ode45` and tspan error Attempted to access

matlab

ode