打破一個大陣分成幾個有意義的變量在Matlab的ODE求解器

Question

我想提高我模擬用Matlab /倍頻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 

它之前幫我很多，但步驟2,3和6（以及在一定程度上4）隨着問題的複雜性增加（60多個州），變得難以管理和重複性管理。

更具體地說，每當我添加一個新的變量，我在3或4個地方傳播信息。我可以將步驟2,3和4捆綁在一起，代替每個變量，這是一個改進，因爲那時信息將在2個位置（開始和結束）。儘管如此，我的目標是將所有內容都放在1個位置，不必擔心索引

以下是我心中的一個僞代碼，我嘗試遵循更多面向對象的方法：

%%(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在你調用你的函數並將它們添加爲輸入之前。您可以根據具體問題調整。