在 Matlab 中为 ODE 求解器将一个大数组分解为几个有意义的变量

Question

我正在尝试改进使用 Matlab/Octave ode 求解器模拟动态系统的方式，更具体地说，是将主要变量数组拆分为更小、更有意义/可读变量的方式。

我实现ode45调用的odefun的标准方式如下（一个简单的例子）：

function d_states=model_as_I_usually_do(t,states)

% 1- parameter definition
m=1;
k=2;


% 2- index definition
pos_i=1:3; %position
vel_i=pos_i(end)+(1:3); %velocity

% 3- extracting data from states vector
pos=states(pos_i);
vel=states(vel_i);

% 4- even more specific renaming for convenince
x=pos(1);y=pos(2);z=pos(3);
u=vel(1);v=vel(2);w=vel(3);

% 5- dynamics
F=[1,1,1]';
a=F/m-k/m*[x y^2 z^3]'; %some equations

% 6- filling output (derivative) vector
d_states(pos_i)=vel;
d_states(vel_i)=a;

end

它以前对我有很大帮助，但是随着问题的复杂性增加（60 多个状态），步骤 2、3 和 6（以及在某种程度上 4）变得难以管理和重复。

更具体地说，每当我添加一个新变量时，我都会在 3 或 4 个位置传播信息。我可以将每个变量的步骤 2,3 和 4 捆绑在一起，这是一个改进，因为这样信息将位于 2 个位置（开头和结尾）。不过，我的目标是将所有内容集中在一个地方，而不必明确担心索引

以下是我想到的伪代码，我尝试遵循更面向对象的方法：

%%( some place before the call to the ode solver)
% new position representation
newstate.name="pos"
newstate.size=3;
newstate.alt_names={"x","y","z"}; % or alternatively, nothing
newstate.derivative="vel";
state_organizer.add(newstate);

% new velocity representation
newstate.name="vel";
newstate.size=3;
newstate.alt_names={"u","v","w"};
newstate.derivative="a";
state_organizer.add(newstate);


function d_states=model_I_would_prefer(t,states,state_organizer)
% 1- parameter definition (the same as before)

% does 2,3 and (possibly)4 without having to worry about indexes
state_organizer.extract(states);

% 5- dynamics ( the same as before)

% does the same as 6
state_organizer.update_derivatives(d_states);
end

在带有指针的语言中，这样做很简单。可能看起来像这样：

% encoding the velocity state
newstate.variable_ptr=&vel;
newstate.derivative_ptr=&a
...
% extracting from the main state array ( states)
for each state in state_organizer
  *(state.variable_ptr)=states(state.index)
end
...
% transferring to array d_states
for each state in state_organizer
  d_states(state.index)=*(state.derivative_ptr)
end

Matlab 没有指针（it has handles，但是每个变量都需要在对象内）。一种替代方法是使用函数“eval”not a good one，但它看起来像这样：

% adding a state ( vel)
organizer.state(n).index=organizer.last_index+(1:size);
organizer.state(n).derivative='a';
organizer.state(n).name='vel';
organizer.last_index=organizer.last_index+size;
...
% step 6 ( similar in essence to step 2)
for each state in organizer 
   eval(['d_states(' organizer.state(k).index ')=' organizer.state(k).derivative ';'])
end

考虑到这一点和我的目标（每个州的信息在一个地方，自动索引），我的主要问题是：

1. 有没有办法在 Matlab/Octave 中实现类似的东西？ （不考虑 Simulink）。
2. 是否有更好的替代方法（更有条理或更可行）？例如，使用字典或自动生成 M 文件。
3. 与 1 和 2 相同，但考虑其他语言或资源，例如 python（我已经考虑切换到它）或 C++。

Answer 1

如果我正确理解您的需求，那么问题出在第 2 步和第 4 步中的硬编码数字 1、2 和 3。将它们作为函数中的输入提供将使函数具有足够的通用性，可以回收用于更大的类别楷模。

说：

% 2- index definition
pos_i=1:n_pos; %position, or pos_i=1:nx+ny+nz
vel_i=pos_i(end)+(1:n_pos); %velocity

% 4- even more specific renaming for convenince
x=pos(1:nx);y=pos(1+nx:nx+ny);z=pos(1+nx+ny:nx+ny+nz);
u=vel(1:nu);v=vel(1+nu:nu+nv);w=vel(1+nu+nv:nu+nv+nw);

您也可以将nx, ny, nz 包装在一个数组中（速度也一样），然后说pos_ind

n_pos = numel(pos_ind);
nx = pos_ind(1);
ny = pos_ind(2);
nz = pos_ind(3);

那么您需要做的就是在调用函数并将它们添加为输入之前定义pos_ind 和vel_ind。您可以根据具体问题进行适当调整。