# calculate the flux of a vector field through a surface

How to calculate the flux of a vector field through a surface in mathematica? I've this field:

F = (x, x^2 * y, y^2 * z)


and this surface:

S = { (x,y,z) ∈ R^3 | 2 * Sqrt[x^2+y^2] <= z <= 1 + x^2 + y^2}


So, I'm trying:

region = ImplicitRegion[2 * Sqrt[x^2+y^2] <= z <= 1 + x^2 + y^2, {x, y, z}];

Integrate[#, {x,y,z} ∈ region]& /@ ({x, x^2 * y, y^2 * z} . {x, y, z})


I expect Pi/30 as a result, but it comes out "Infinity"...

• Volume[region] is infinite. Have a look at Plot3D[{2*Sqrt[x^2 + y^2], 1 + x^2 + y^2}, {x, -5, 5}, {y, -5, 5}]. Your region is all the 3D space (z) in between those surfaces for all x,y. It is not a surface but an infinite volume. You need to put some bounds on the x,y and decide where you want to cut it off, and also to reformulate your region as a surface and not a volume. Sep 5, 2020 at 16:38
• By hand, I solve by transforming into spherical coordinates. How can I get mathematica to do it?
– Teo7
Sep 5, 2020 at 17:17
• You have provided a region that is not a surface and is infinitely large. Please address this issue first. Describe the surface in plain English and maybe I can suggest the correct representation. Sep 5, 2020 at 17:25
• I've the surface right in this form... specifically, I've to calculate the flux of the field F leaving the surface of the solid S. Where S is 2 * Sqrt[x^2+y^2] <= z <= 1 + x^2 + y^2
– Teo7
Sep 5, 2020 at 17:49
• Right, but that surface has an infinite surface area and extends in all directions x,y. Do you want a cylindrical-like 'slice' of that solid , i.e break it into three parts the top,edge,bottom, or do you want a box-like surface where the x,y ranges are capped? Sep 5, 2020 at 17:55

The problem is the region is infinite. We need to restrict the radius like so:

region = ImplicitRegion[2*Sqrt[x^2+y^2] < z < x^2+y^2+1 && x^2+y^2 <= 1, {x,y,z}]


This region describes the volume in between a height 2 inverted cone of base radius 1, and a curved cap surface:

Plot3D[{2 Sqrt[x^2 + y^2], x^2 + y^2 + 1}, {x, -1, 1}, {y, -1, 1},
PlotStyle -> {Opacity[.5]},
RegionFunction -> Function[{x, y, z}, x^2 + y^2 < 1]]


We can then apply the divergence theorem to calculate the flux:

f = {x, x^2 y, y^2 z};
Integrate[Div[f, {x, y, z}], {x, y, z} ∈ region]
(* result Pi/5 *)


... or we can calculate the surface integral directly:

f = {x, x^2 y, y^2 z};
coneSurface = {x, y, 2 Sqrt[x^2 + y^2]};
capSurface = {x, y, x^2 + y^2 + 1};
(* multiply by -1 because the normal is pointing into the cone but we need it pointing out *)
coneNormal = -1*ResourceFunction["UnitNormal"][coneSurface, {x, y}];
capNormal = ResourceFunction["UnitNormal"][capSurface, {x, y}];
coneRegion = ImplicitRegion[x^2 + y^2 <= 1 && z == 2 Sqrt[x^2 + y^2], {x,y,z}];
capRegion = ImplicitRegion[x^2 + y^2 <= 1 && z == x^2 + y^2 + 1, {x,y,z}];
Integrate[coneNormal.f, {x,y,z} ∈ coneRegion] + Integrate[capNormal.f, {x,y,z} ∈ capRegion]

(* result: Pi/5 *)


The book result of $$\pi/30$$ is not correct. The reason why is because $$\nabla\cdot F$$ is $$1+x^2+y^2$$ not $$x^2+y^2$$, as the $$F$$ has a non-constant $$x$$ component.

If we make $$F_x$$ constant, e.g f = {0, x^2 y, y^2 z};, then both approaches above give $$\pi/30$$.