runaway 'unable to simplify to float approximation' in iterated numerical integral

[kcrisman]

See ​this ask.sagemath question. Apparently sometimes iterated numerical integrals don't recognize that they are just of one variable on the inside - or something. The following example is probably nowhere near minimal but is from the reporter:

var('T1, T2, T3, T4, T5, T6, T7, T8, r, r_0, r_2, r_a, r_s,theta_a, phi_a, theta_c, theta, phi, A');

r_a=2;

theta_a=pi/2;

phi_a=0;

r_s=vector([r, theta, phi]);

r_0=vector([r_a, theta_a, phi_a]);

T1=vector ([sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)]);

T2=vector ([cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta)]);

T3=vector ([-sin(phi), cos(phi), 0]);

T5=vector ([sin(theta_a)*cos(phi_a), sin(theta_a)*sin(phi_a), cos(theta_a)]);

T6=vector ([cos(theta_a)*cos(phi_a), cos(theta_a)*sin(phi_a), -sin(theta_a)]);

T7=vector ([-sin(phi_a), cos(phi_a), 0]);

var('L_s, SA_0, SA_1, SA_2, theta_h, I_0, I, rho, I_1, Phi_0, Phi_1, Phi_2, E_v') ;

L_s=20;

I_0=2.4; 

theta_h= 2*pi/45; 

Phi_0=I_0*SA_0; 

Phi_0;

T4=matrix ([T1, T2, T3]);

T8=matrix ([T5, T6, T7]);

(r_s*T4)*(r_0*T8);

A=((r_s*T4)*(r_0*T8));

theta_c=arccos(A/(r*r_a));

r_2=sqrt(r^2+r_a^2-2*r*r_a*cos(theta_c)); 

Phi_0 ;

Phi_0*(1-exp(-r/L_s));

f=diff(Phi_0*(1-exp(-r/L_s)), r); f

rho=1/(SA_0*r^2)*f; rho; SA_0; f;

y=rho/(r_2)^2*r^2*sin(theta)^2;

numerical_integral(lambda phi: numerical_integral(lambda theta: numerical_integral(lambda r: y, 0,5)[0], 0, pi)[1], 0, pi)

Note that numerical_integral(lambda r: y.subs(theta=1,phi=1),0,5) yields this problem while z = y.subs(theta=pi/2,phi=pi/2); numerical_integral(z,0,5) does not, and yet numerical_integral(lambda r: z,0,5) does.