ODE求解器在一起運行時緩慢運行

Question

我有一個類的一些ODE求解器的自定義實現。我遇到問題，如果我的時間步dt小於.2，程序將停止。但是，如果我註釋掉一個龍格 - 庫塔求解器，它會很快執行，我可以切換哪一個被註釋掉，這樣我就可以從兩個求解器中獲得解決方案。我想知道如何解決這個問題。我一直在試圖找到某種解決方法可能會干擾另一方的方法，但我不明白這是怎麼發生的。

實現：

global dt; 
dt = .5; % going below ~.25 makes the program take a very long time to exit 

g = 1; 
c_d = 2; 
m = 3; 

tf = 15; 

dudt = @(t, u) g - (c_d/m) * u.^2 

[t_euler, u_euler] = euler(dudt, [0, tf], 0); 
[t_rk4, u_rk4] = rk4(dudt, [0, tf], 0); % either of this one or rk2 
             % can be commented out to make the program 
             % run quickly, but rk2 and rk4 cannot be 
             % run at the same time 
[t_rk2, u_rk2] = rk2(dudt, [0, tf], 0); 


%% rk4.m %% 

function [t, u] = rk4(odefun, tspan, u0) 

t0 = tspan(1); 
t = [ t0 ]; 
t_new = t0; 
global dt; 
tf = tspan(2); 
u_new = u0; 

u = [ u0 ]; 

while (t_new < tf) 
    if (t_new + dt > tf) 
     dt = tf - t_new; 
    end 

    k1 = dt * odefun(t_new, u_new); 
    k2 = dt * odefun(t_new + dt/2, u_new + k1/2); 
    k3 = dt * odefun(t_new + dt/2, u_new + k2/2); 
    k4 = dt * odefun(t_new + dt, u_new + k3); 

    u_new = u_new + 1/6 * (k1 + 2*k2 + 2*k3 + k4); 
    % rk2.m is the same as rk4.m, except u_new = u_new + k2 
    % euler.m is also the same, except u_new = u_new + k1 
    u = [ u, u_new ]; 

    t_new = t_new + dt; 
    t = [ t, t_new ]; 
end 
end 

Answer 1

當計算的時間間隔長度的最後一個步驟，以滿足tf準確，可以減少全局變量dt，而在另一種方法是使用dt現在較小的值。有了您的數據，您可能會有新的dt小到1e-16甚至0。

您想將全局dt視爲一個常量，也許使用局部變量h而不是dt。或者通過dt作爲功能的參數。