Python Forum

Full Version: Integrating for the volume of a torus in SciPy
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Hi,
I have been trying to perform a triple integral in cylindrical coordinates using tplquad in SciPy. The particular integral is for a torus with radius R = 4 and cross-sectional radius r = 1. This is my code:
import numpy as np
from scipy.integrate import tplquad

R = 4
r = 1

# Volume integral
# Order of integration is z, rho, theta
f1 = lambda z, rho, theta: 8 * 2 * rho
# set theta limits
a, b = 0, np.pi/2
# set rho limits
rho_low = lambda rho: R - r
rho_high = lambda rho: R + r
# set z limits
z_low = lambda rho, theta: 0
z_high = lambda rho, theta: np.sqrt(r**2 - (rho - R)**2)
# calculate integral
V, _ = tplquad(f1, a, b, rho_low, rho_high, z_low, z_high)
print(V)
I am receiving the following error: "RuntimeWarning: invalid value encountered in sqrt
z_high = lambda rho, theta: np.sqrt(r**2 - (rho - R)**2)"

I've seen in previous integrations that when an integration includes a square root of one of the variables then the integration has to be undertaken in the first octant only and then multiplied by 8. This has been my approach here which is why the upper theta limit is pi/2. I am aware that I could do this integration as a volume of revolution with only one variable to integrate over but I would like to understand why this method won't work and how it can be resolved. Grateful for any pointers, thanks!
My guess is that for z_low and z_high it should be lambda theta, rho:... because the documentation says qfun(x, y) and rfun(x, y)
(Jan-07-2020, 03:43 PM)Gribouillis Wrote: [ -> ]My guess is that for z_low and z_high it should be lambda theta, rho:... because the documentation says qfun(x, y) and rfun(x, y)

It turns out you are correct! Thanks!
I ran the code
import numpy as np
from scipy.integrate import tplquad

R = 4
r = 1

# Volume integral
# Order of integration is z, rho, theta
f1 = lambda z, rho, theta: 2 * rho
# set theta limits
a, b = 0, 2 * np.pi
# set rho limits
rho_low = lambda rho: R - r
rho_high = lambda rho: R + r
# set z limits
z_low = lambda theta, rho: 0
z_high = lambda theta, rho: np.sqrt(r**2 - (rho - R)**2)
# calculate integral
V, _ = tplquad(f1, a, b, rho_low, rho_high, z_low, z_high)
print(V)
which gives the correct result. I did not have to integrate in the first octant only and then multiply by 8 after all Think. How come the order 'lambda rho, theta' did not work and 'lambda theta, rho' does?