Two functions which are identical except for one very small summand. Why does the ODE solver take orders of magnitude longer to solve the one without the summand?

I have two functions, which are identical save for one very small summand. When I try to run ode113 (or ode45 for that matter) to solve ODEs given by the two functions, the one with the extra summand takes less than a minute, but the one without the extra summand takes two hours! I'm scratching my head to figure out why this must be the case. Here is the code (I've tried to be as minimal as I can) :

Some variables that have to be initialized:

    a=5; b=1; n_ellipse_points = 500;
    s=linspace(0,2*pi,n_ellipse_points);
    c = sqrt(a^2 - b^2);
    u_0 = acosh(a/c);
    rnd(1);
    u_pert = rand(size(s))*1e-5;
    Z = c*cosh(u_0+u_pert).*cos(s) + 1i*c*sinh(u_0+u_pert).*sin(s); Z=Z.'

And now the two functions:

function U = dZdt(Z,zeta)

    %%%Z needs to be a column vector

    Z = Z.';
    U = zeros(1,length(Z));

    tmp = [Z(end) Z Z(1)];
    dZds = (tmp(3:length(Z)+2) - tmp(1:length(Z)))/2;
    for k=1:length(Z)
    tmp2 = (Z-Z(k))./conj(Z-Z(k));
    tmp2 =      tmp2([1:k-1 k+1:length(Z)]);
    tmp3 = conj(dZds([1:k-1 k+1:length(Z)]));  
    U(k) = sum(tmp2.*tmp3);
    end

    U = U(:) * zeta / (4 * pi);

and

    function U = dZdt2(Z,zeta)

    %%%Z needs to be a column vector

    Z = Z.'; %Transposes but keeps +- signs the same
    U = zeros(1,length(Z));

    tmp = [Z(end) Z Z(1)];
    dZds = (tmp(3:length(Z)+2) - tmp(1:length(Z)))/2;

    for k=1:length(Z)
    tmp2 = (Z-Z(k))./conj(Z-Z(k));
    tmp3 = conj(dZds);  
    if k == 1
        tangent = (Z(2)-Z(end))/norm(Z(2)-Z(end));
    elseif k == length(Z)
        tangent = (Z(1)-Z(end-1))/norm(Z(1)-Z(end-1));
    else
        tangent = (Z(k+1)-Z(k-1))/norm(Z(k+1)-Z(k-1));
    end
    tmp2(k) = exp(2i*angle(tangent));
    U(k) = sum(tmp2.*tmp3); 
    end

    U = U(:) * zeta / (4 * pi);

Notice that the only difference between the two is the inclusion of the "exp(2*i*tangent)" summand in the second function. However, running

    [T,ZZ] = ode113(@(t,Z) dZdt2(Z,zeta), [0 20], Z)
takes under a minute, whereas
    [T,ZZ] = ode113(@(t,Z) dZdt(Z,zeta), [0 20], Z)
takes two hours.

If anyone can figure out what's going on, I'd be eternally grateful. Thanks!

Answer 1

[Steven Lord]

That modification most likely makes the second ODE a stiff ODE. As per the "Basic Solver Selection" table on this documentation page, ode113 is intended to be used on nonstiff problems. Try a stiffer solver like ode15s.

"You can identify a problem as stiff if nonstiff solvers (such as ode45) are unable to solve the problem or are extremely slow. If you observe that a nonstiff solver is very slow, try using a stiff solver such as ode15s instead."