ODEソルバーの実行速度が遅い場合

Question

クラスのODEソルバーのカスタム実装があります。私は時間のステップをdtより小さくして.2より小さくすると、プログラムが停止する問題があります。しかし、Runge-Kuttaソルバーの1つをコメントアウトすると、非常に迅速に実行され、どのソルバーからコメントを外すかを切り替えることができるので、両方のソルバーから解を得ることができます。私はこれを修正する方法が不思議です。私は、あるソルバが他のソルバを妨害しているかもしれない何らかの方法を見つけることを試みてきましたが、どうしてこのようなことが起こるのか分かりません。

実装：正確tfを満たすために、最後のステップの間隔の長さを計算する場合

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

、あなたはグローバル変数dtを削減し、他の方法でdtの今小さい方の値が使用されます。データでは、新しいdtは1e-16、さらには0と小さくてもかまいません。

グローバルdtを定数として扱いたい場合は、おそらくdtの代わりにローカル変数hを使用します。または、関数のパラメータとしてdtを渡します。