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 = 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 = 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));
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!
