在特定條件下用標誌改變ODE計算中的常數

Question

我有一個用於計算酸性變化的ODE。一切工作都很好，只要酸度達到臨界點時我就想改變常數。它應該是我想模擬的某種不可逆轉的效果。

我的常量來自一個結構文件（c）我在ODE函數中加載一次。

[Time,Results] = ode15s(@(x, c) f1(x, c),[0 c.length],x0,options); 

主要的問題我這裏是不會告訴Matlab來改變常數，但記得如果在模擬過程中已經發生過一次。所以Matlab應該採用不可逆轉的常數而不是我在開始時提供的常數。

我想寫一個在ODE運行時保存的標誌和一個if條件，「如果標誌存在，改變常量」。怎麼做？

更新我： 問題就迎刃而解了

這裏第一個自己的解決方案，它不是拋光和需要結構文件的方法。這意味着，事件應該突然改變的常量必須是結構文件，這些結構文件將ODE函數交給應該評估的函數（請參閱上面的ODE語法）。該函數接受的輸入是這樣的：

function [OUTPUT] = f1(t, x, c) 

% So here, the constants all start with c. followed by the variable name. (because they are structs instead of globals) 

%% write a flag when that event happens 

if c.someODEevent <= 999 && exist ('flag.txt') == 0 
    dlmwrite ('flag.txt',1); 
end 

%% next iteration will either write the flag again or not. more importantly, if it exists, the constant will change because of this. 

if exist ('flag.txt') == 2 
     c.changingconstant = c.changingconstant/2; 
end 



end 

請看Horchlers樣的回答，你必須要小心，這樣的步驟可能引入不準確，你必須要小心檢查，如果你的代碼做什麼是應該去做。

Answer 1

要準確地做到這一點，您應該在ODE求解器中使用event detection。我不能給你一個具體的答案，因爲你只提供ode15s在你的問題中調用它，但你需要編寫一個events函數，然後通過odeset指定它。事情是這樣的：

function acidity_main 
% load c data 
... 
x0 = ... 
options = odeset('Events',@events); % add any other options too 

% integrate until critical value and stop 
[Time1,Results1] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options); 

x0 = Results(end,:); % set new initial conditions 
% pass new parameters -it's not clear how you're passing parameters 
% if desired, change options to turn off events for faster integration 
[Time2,Results2] = ode15s(@(x,c)f1(x,c),[0 c.length],x0,options); 

% append outputs, last of 1 is same as first of 2 
Time = [Time1;Time2(2:end)]; 
Results = [Results1;Results2(2:end,:)]; 
... 

function y=f1(x,c) 
% your integration function 
... 

function [value,isterminal,direction] = events(x,c) 
value = ... % crossing condition, evaluates to zero at event condition 
isterminal = 1; % stop integration when event detected 
direction = ... % see documentation 

你要使用事件進行正確的整合到了「acidicity達到一個臨界點」，並停止整合點。然後再次用新值調用ode15s並繼續進行整合。這可能看起來很粗糙，但它是如何準確地完成這種事情的。

您可以看到一個基本事件檢測示例here。在您的命令窗口中鍵入ballode以查看此代碼。您可以在my answer here中看到這個演示的稍微更復雜的版本。這裏的an example使用事件在指定的時間（而不是您指定的狀態值的情況下）準確地更改ODE。

注意：我覺得很奇怪，你將所謂的「常量」c作爲第二個參數傳遞給ode15s。該函數具有嚴格的輸入參數要求：第一個是自變量（通常是時間），第二個是狀態變量數組（與您的初始條件向量相同）。另外，如果f1只有兩個參數，則@(x,c)f1(x,c)是多餘的 - 通過@f1就足夠了。