Uge 9b: Integration of Vector Fields: The Line Integral and Flux#
Demo by Christian Mikkelstrup, Hans Henrik Hermansen, Jakob Lemvig, Karl Johan Måstrup Kristensen, and Magnus Troen
from sympy import *
from dtumathtools import *
init_printing()
Brief Intro to Vector Fields in SymPy#
In general when working with vector fields in SymPy
it is often easiest to defined them as lambda/python functions, since we often are about to insert a vector function into \(\pmb{V}\).
x, y, z = symbols("x,y,z", real=True)
V = lambda x, y, z: Matrix([y * cos(x), y * sin(x), z])
V(x, y, z)
If we for instance have an expression of a parameter function
then we can easily insert \(\boldsymbol{r}_t\) into \(V\) by doing:
r1, r2, r3 = symbols("r1,r2,r3", cls=Function)
t = symbols("t")
r = Matrix([r1(t), r2(t), r3(t)])
# r inserted into V
V(*r)
If we want to use \(\pmb{V}\), we can analytically simply insert symbolic variables. In later/other courses or in the project period in this course one may get the need for computing vector field properties such as curl and divergence. This is done easily with dtumathtools
:
rotV = dtutools.rot(V(x, y, z), (x, y, z))
rotV
divV = dtutools.div(V(x, y, z), (x, y, z))
divV
![../_images/6ac6b1dd0b2a2d12ebb2d4e0599a6f3887d3d5ce08b486131444d8a3d0a60d5b.png](../_images/6ac6b1dd0b2a2d12ebb2d4e0599a6f3887d3d5ce08b486131444d8a3d0a60d5b.png)
If one is to carry out sufficiently many analytical computations on \(\pmb{V}\), it may though turn out to be more optimal to just define \(\pmb{V}\) as a SymPy object instead and then use SymPy’s \(\verb|.subs()|\) for inserting variables:
V = Matrix([y * cos(x), y * sin(x), z])
V
Line Integral of Vector Fields (Tangential Line Integral)#
We consider the vector field
as well as two curves \(K_1\) and \(K_2\) with the parametric representations
where \(u \in [0,4\pi]\).
x, y, z, u = symbols("x y z u", real=True)
r1 = Matrix([cos(u), sin(u), u / 2])
r2 = Matrix([1, 0, u / 2])
V = Matrix([-y, x, 2 * z])
u_range = (u, 0, 4 * pi)
K1 = dtuplot.plot3d_parametric_line(
*r1, u_range, show=False, rendering_kw={"color": "red"}, colorbar=False
)
K2 = dtuplot.plot3d_parametric_line(
*r2, u_range, show=False, rendering_kw={"color": "blue"}, colorbar=False
)
vectorfield_V = dtuplot.plot_vector(
V,
(x, -1, 1),
(y, -1, 1),
(z, 0, 6),
n=5,
quiver_kw={"alpha": 0.5, "length": 0.1, "color": "black"},
colorbar=False,
show=False,
)
combined = K1 + K2 + vectorfield_V
combined.legend = False
combined.show()
![../_images/79d95e2345dd7332666470a1e7a6af73d75abbd215e8bdf2f2f4f6de0eef09a7.png](../_images/79d95e2345dd7332666470a1e7a6af73d75abbd215e8bdf2f2f4f6de0eef09a7.png)
We want to compute the tangential line integral of \(\pmb{V}\) along each of the two curves from the point \(a = (1,0,0)\) to \(b = (1,0,2\pi)\). This corresponds for both of them to \(u \in [0,4\pi]\). First, we will find the tangent vectors:
r1d = r1.diff(u)
r2d = r2.diff(u)
r1d, r2d
The functions that are to be integrated are according to the theorem:
integrand1 = V.subs({x: r1[0], y: r1[1], z: r1[2]}).dot(r1d)
integrand2 = V.subs({x: r2[0], y: r2[1], z: r2[2]}).dot(r2d)
integrand1.simplify(), integrand2.simplify()
![../_images/9e9ee152c6ee988377abefdc2488b140c0741bd42af19bc5c40ee1b33dd63974.png](../_images/9e9ee152c6ee988377abefdc2488b140c0741bd42af19bc5c40ee1b33dd63974.png)
The tangential line integral of \(f\) along the curve \(K_1\) is then
integrate(integrand1, (u, 0, 4 * pi))
![../_images/eb2e323f28743cbb92d556527d655ff334140e357f0449e1588c0b1471797e15.png](../_images/eb2e323f28743cbb92d556527d655ff334140e357f0449e1588c0b1471797e15.png)
The tangential line integral of \(f\) along the curve \(K_2\) is then
integrate(integrand2, (u, 0, 4 * pi))
![../_images/e26e785982680e961d65484c538555ed7e9063e716b2dd40f2ad6ab24d0d8ce6.png](../_images/e26e785982680e961d65484c538555ed7e9063e716b2dd40f2ad6ab24d0d8ce6.png)
Perspective:#
If we interpret \(\pmb{V}\) as a force vector field, then the work (energy) with which \(\pmb{V}\) constributes when the particles is moved along the vertical line is \(4\pi\) less than if the particle was moved along the spiral curve. The magnitude of the work is hence in this example dependent on the path. \(\pmb{V}\) can hence not be a gradient vector field.
Seaching for Anti-Derivative via Stair Lines. Theory#
An anti-derivative to a smooth vector field \(\pmb{V}\) can, if an anti-derivative exists, be determined via the tangential line integral of \(\pmb{V}\) along an arbitrary curve from the origin \(\pmb{x}_0 = (0,0,0)\) to an arbitrary fixed point \(\pmb{x} = (x, y, z)\). One achives the easiest computations by using the stair line from the origin to \(P\) as the curve.
u_range = (u, 0, 1)
r1 = dtuplot.plot3d_parametric_line(
u, 0, 0, u_range, show=False, rendering_kw={"color": "red"}, colorbar=False
)
r2 = dtuplot.plot3d_parametric_line(
1, u, 0, u_range, show=False, rendering_kw={"color": "red"}, colorbar=False
)
r3 = dtuplot.plot3d_parametric_line(
1, 1, u, u_range, show=False, rendering_kw={"color": "red"}, colorbar=False
)
xyz = dtuplot.scatter(Matrix([1, 1, 1]), show=False, rendering_kw={"color": "black"})
combined = r1 + r2 + r3 + xyz
combined.legend = False
combined.camera = {"azim": 37, "elev": 16}
combined.show()
![../_images/78e98466b428b476a2a2dad96cd55aa9f99dd78fe9dd895b227d644e39c7cef0.png](../_images/78e98466b428b476a2a2dad96cd55aa9f99dd78fe9dd895b227d644e39c7cef0.png)
The plot above shows the three stair lines which the methods utilizes. We will start with the following parametrization of the lines:
r1 = Matrix([u, 0, 0])
r2 = Matrix([x, u, 0])
r3 = Matrix([x, y, u])
Here we have, respectively, \(u \in [0,x]\), \(u \in [0,y]\), and \(u \in [0,z]\). An advantage with this method is that the three tangent vectors are very simple.
r1d = r1.diff(u)
r2d = r2.diff(u)
r3d = r3.diff(u)
r1d, r2d, r3d
Now consider an arbitrary smoot vector field
The tangential line integral of \(\pmb{V}\) along the stair line consists of the sum of the tangential line integral of \(\pmb{V}\) along each of the three line segments of which the stair consists.
We consider the three integrands that are involved in the three integrals:
Hence we have the following formula for the tangential line integral of \(\pmb{V}\) along a stair line:
Searching for Anti-Derivative via Stair Lines. Example#
We are to investigate whether the following vector field \(\pmb{V}\) is a gradient vector field:
We compute the tangential line integral of \(\pmb{V}\) along the stair line by use of the formula:
We find the three integrands:
V = Matrix([y**2 + z, 2 * y * z**2 + 2 * y * x, 2 * y**2 * z + x])
integrand1 = V[0].subs({x: u, y: 0, z: 0})
integrand2 = V[1].subs({y: u, z: 0})
integrand3 = V[2].subs({z: u})
integrand1, integrand2, integrand3
![../_images/31bde68c89e9326a644578a6a2fbdd1e7e0ee473c27f83ba481ed1f5fd18c17c.png](../_images/31bde68c89e9326a644578a6a2fbdd1e7e0ee473c27f83ba481ed1f5fd18c17c.png)
The tangential line integral of \(\pmb{V}\) along the stair line is then:
F = (
integrate(integrand1, (u, 0, x))
+ integrate(integrand2, (u, 0, y))
+ integrate(integrand3, (u, 0, z))
)
F
![../_images/f10e8f5dfd7f81d745c6d093322bd8fad250b3d28cfa9da97bbe073967992552.png](../_images/f10e8f5dfd7f81d745c6d093322bd8fad250b3d28cfa9da97bbe073967992552.png)
If we consider the line integral as a function \(F\) of \((x,y,z)\), we can now test whether \(F\) really is an anti-derivative to \(\pmb{V}\).
V_F = dtutools.gradient(F)
V_F, Eq(V_F, V)
Since \(\pmb{V}\) is identical to the gradient of \(F\), \(F\) is an anti-derivative to \(\pmb{V}\)! This means that the tangential line integral of \(\pmb{V}\) along an arbitrary closed curve always will result in \(0\). Let us check this claim on the following knot:
t = symbols("t")
knot = (
Matrix(
[
-10 * cos(t) - 2 * cos(5 * t) + 15 * sin(2 * t),
-15 * cos(2 * t) + 10 * sin(t) - 2 * sin(5 * t),
10 * cos(3 * t),
]
)
* S(1)
/ 10
)
knot
dtuplot.plot3d_parametric_line(
*knot, (t, -pi, pi), rendering_kw={"color": "blue"}, legend=False
)
![../_images/016cacb5af93fdaa01bd881feeb1656b9608bf20a4e13f7eb3ce903ed9654916.png](../_images/016cacb5af93fdaa01bd881feeb1656b9608bf20a4e13f7eb3ce903ed9654916.png)
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f3ce59cfbe0>
## Note -- This integral takes a long time for SymPy to compute. It might easily take more than a minute, so be patient
#integrate(V.subs({x: knot[0], y: knot[1], z: knot[2]}).dot(knot.diff(t)), (t, -pi, pi))
The work which the force field \(\pmb{V}\) contributes with as the particle moves a trip around the know is thus \(0\). Just as expected, since \(\pmb{V}\) is a gradient vector field.
Surface Integral of Vector Fields (Flux through Surfaces)#
We are given a vector field by:
x, y, z = symbols("x y z", real=True)
V = lambda x, y, z: Matrix([x, y, 2])
V(x, y, z)
We consider the part of the standard paraboloid \(F_{parab}\) that is located below the plane \(z=1\). This can be parametrized by:
u, v = symbols("u v", real=True)
r = Matrix([u * cos(v), u * sin(v), u**2])
r
where \(u\in [0,1]\) and \(v\in [-\pi,\pi]\). The vector field and the surface are illustrated:
p_parab = dtuplot.plot3d_parametric_surface(
*r,
(u, 0, 1),
(v, -pi, pi),
use_cm=False,
camera={"elev": 25, "azim": -55},
show=False
)
p_field = dtuplot.plot_vector(
V(x, y, z),
(x, -1, 1),
(y, -1, 1),
(z, 0, 1.2),
n=4,
use_cm=False,
quiver_kw={"alpha": 0.5, "length": 0.1, "color": "red"},
show=False,
)
(p_parab + p_field).show()
![../_images/0d6bb0091db0d2be2909506af06da8926e1c9f5910cbd4d2361bcf91a115e288.png](../_images/0d6bb0091db0d2be2909506af06da8926e1c9f5910cbd4d2361bcf91a115e288.png)
We can do with just plotting the vector field on the surface itself:
# Vectors arrive from all corners. There are too many
# vectors when "p_parab" is used, so here we reduce n
dummy_surface = dtuplot.plot3d_parametric_surface(
*r, (u, 0, 1), (v, -pi, pi), n=8, show=False
)
p_field2 = dtuplot.plot_vector(
V(x, y, z),
(x, -1, 1),
(y, -1, 1),
(z, 0, 1),
slice=dummy_surface[0],
use_cm=False,
quiver_kw={"alpha": 0.5, "length": 0.1, "color": "red"},
show=False,
)
# Better camera angle for looking at the vectors with
p_parab.camera = {"elev": 60, "azim": -50}
(p_parab + p_field2).show()
![../_images/de7030678b9f90007689f9e7aabc5d15fbc6a530cc736c51f0ae8049882c0160.png](../_images/de7030678b9f90007689f9e7aabc5d15fbc6a530cc736c51f0ae8049882c0160.png)
Computation of Flux “down” through the Paraboloid \(F_{parab}\)#
This clearly looks as though there is a (positive) flux up through the surface. We would thus expect to see that the flux down through the surface is negative. We will start by finding the normal vector to the parabloid given by our chosen parametrization:
n_parab = r.diff(u).cross(r.diff(v))
simplify(n_parab)
Notice here that the normal vector \(n\) is pointing in the “wrong” direction (inwards towards the \(z\) axis and thus with a postive \(z\) coordinate); “wrong” in the sense that we actually wish to find the flux down through the surface. We thus simply have to remember to change the sign of the normal vector \(n_{parab}\):
n_parab = -simplify(n_parab)
n_parab
The length of the normal vector is only needed for arguing that the parametrization is regular:
n_parab.norm().simplify()
![../_images/cbbce2c44bc9837012f874cc999d38980fa418067ff1ab1e314435a004119db7.png](../_images/cbbce2c44bc9837012f874cc999d38980fa418067ff1ab1e314435a004119db7.png)
This expression is positive on the interior of \(\Gamma\), so it is indeed regular. We will now find the integrand:
integrand = n_parab.dot(V(*r))
simplify(integrand)
![../_images/2ef979e230cb8c3f1f42e774e4c0e76d9ce1fbc2a2fdb72e0805cabace2c6dc1.png](../_images/2ef979e230cb8c3f1f42e774e4c0e76d9ce1fbc2a2fdb72e0805cabace2c6dc1.png)
with which the flux is computed by doing
integrate(integrand, (u, 0, 1), (v, -pi, pi))
![../_images/498164797fbb562656ed470ccea29fa5cf189cf5b0f949dfa2ddb3bd2c21e64a.png](../_images/498164797fbb562656ed470ccea29fa5cf189cf5b0f949dfa2ddb3bd2c21e64a.png)
The wanted flux is thus \(-\pi\)!
Computation of the Flux through the “Lid” \(F_{lid}\) of the Paraboloid in the (Positive) Direction of the \(z\) axis#
r2 = Matrix([u * cos(v), u * sin(v), 1])
n2 = r2.diff(u).cross(r2.diff(v))
simplify(n2)
The direction of the normal vector fits with the wanted direction through the surface \(F_{lid}\).
integrand2 = n2.dot(V(*r2))
integrate(integrand2, (u, 0, 1), (v, -pi, pi))
![../_images/014a8f4145d73d4e78fd0f4ce71fbb60187bfd1bb1e5dc407404758e04c54a74.png](../_images/014a8f4145d73d4e78fd0f4ce71fbb60187bfd1bb1e5dc407404758e04c54a74.png)
The flux through the lid is \(\mathrm{Flux}(\pmb{V},F_{lid}) = 2\pi\).
The Flux through the closed Surface with Paraboloid and Lid \(F_{closed} = F_{parab} \cup F_{lid}\)#
p_lid = dtuplot.plot3d_parametric_surface(
*r2,
(u, 0, 1),
(v, -pi, pi),
use_cm=False,
n1=4,
n2=16,
camera={"elev": 25, "azim": -55},
show=False
)
p_fieldlid = dtuplot.plot_vector(
V(x, y, z),
(x, -1, 1),
(y, -1, 1),
(z, 0, 1),
slice=p_lid[0],
use_cm=False,
quiver_kw={"alpha": 0.5, "length": 0.1, "color": "red"},
show=False,
)
p_parab.camera = {"elev": 15, "azim": -60}
(p_parab + p_field2 + p_lid + p_fieldlid).show()
![../_images/084b1f754d804542d3a0c49f9f0dc83bddf2659f78e13c03eab86b63537d0e73.png](../_images/084b1f754d804542d3a0c49f9f0dc83bddf2659f78e13c03eab86b63537d0e73.png)
We do not need to continue with any further calculations! We here have the flux through all surfaces, from which the total flux through the closed surface is \(Flux(\pmb{V},F_{closed})=-\pi+2\pi=\pi\).
Flux through a Piece of a Spherical Surface#
We are given a vector field,
x, y, z = symbols("x y z")
V = Matrix([-y + x, x, z])
V
and a section of a spherical surface,
u, v, w = symbols("u v w", real=True)
radius = 2
r = radius * Matrix([sin(u) * cos(v), sin(u) * sin(v), cos(u)])
r
with \(u\in[\pi/6,\pi/2]\) and \(v\in[0,\pi]\),
p_surface = dtuplot.plot3d_spherical(
radius,
(u, pi / 6, pi / 2),
(v, 0, pi),
aspect="equal",
camera={"elev": 25, "azim": 55},
show=False,
)
dummy_surface = dtuplot.plot3d_spherical(
radius, (u, pi / 6, pi / 2), (v, 0, pi), show=False, n=8
)
p_field = dtuplot.plot_vector(
V,
(x, -1, 1),
(y, -1, 1),
(z, 0, 1),
slice=dummy_surface[0],
use_cm=False,
quiver_kw={"alpha": 0.5, "length": 0.2, "color": "red"},
show=False,
)
(p_surface + p_field).show()
![../_images/ffcacac5444dfd6297c0bd661a7de30b7de64cf5faf8eebadfe46592070047de.png](../_images/ffcacac5444dfd6297c0bd661a7de30b7de64cf5faf8eebadfe46592070047de.png)
The flux is computed to be
n = r.diff(u).cross(r.diff(v))
integrand = n.dot(V.subs({x: r[0], y: r[1], z: r[2]}))
integrate(integrand, (u, pi / 6, pi / 2), (v, 0, pi))
![../_images/4f9489d8500adb5cfd71971179239868ecd08d6d1bfaa7b0b75b1f7807263f99.png](../_images/4f9489d8500adb5cfd71971179239868ecd08d6d1bfaa7b0b75b1f7807263f99.png)
Flux through a Section of a spherical surface using Gauss’ Divergence Theorem (not Syllabus)#
We will compare the result from the flux computation with a result from using Gauss’ Theorem through a semi-sphere with a radius \(a\) (\(u\in[0,a]\), \(v\in[0,\pi/2]\), and \(w\in[-\pi,\pi]\)), with the flux through its two surfaces. This should give the same result!
V = Matrix([8 * x, 8, 4 * z**3])
r = u * Matrix([sin(v) * cos(w), sin(v) * sin(w), cos(v)])
V, r
First, using Gauss’ Theorem,
a = symbols("a", real=True, positive=True)
M = Matrix.hstack(r.diff(u), r.diff(v), r.diff(w))
# the determinatn is positive. The only term that is not
# squared is sin(v), which is positive within v's limits
jacobian = M.det()
divV = dtutools.div(V, [x, y, z])
divV_r = divV.subs({x: r[0], y: r[1], z: r[2]})
integrate(divV_r * jacobian, (u, 0, a), (v, 0, pi / 2), (w, -pi, pi))
![../_images/026406098406e3a8c891856880537f6a050f214e3b182044271283c049e4af45.png](../_images/026406098406e3a8c891856880537f6a050f214e3b182044271283c049e4af45.png)
and then through the two surfaces,
# The first surface corresponds to fixing u as the radius
r1 = r.subs(u, a)
n1 = r1.diff(v).cross(r1.diff(w)) # we have checked that it points outwards
integrand1 = n1.dot(V.subs({x: r1[0], y: r1[1], z: r1[2]}))
flux1 = integrate(integrand1, (v, 0, pi / 2), (w, -pi, pi))
# Second surface corresponds to fixing v to pi/2
# The remaining becomes the disc with radius a
r2 = r.subs(v, pi / 2)
n2 = -r2.diff(u).cross(r2.diff(w)) # must be negative in order to point outwards
integrand2 = n2.dot(V.subs({x: r2[0], y: r2[1], z: r2[2]}))
flux2 = integrate(integrand2, (u, 0, a), (w, -pi, pi))
flux1 + flux2
![../_images/026406098406e3a8c891856880537f6a050f214e3b182044271283c049e4af45.png](../_images/026406098406e3a8c891856880537f6a050f214e3b182044271283c049e4af45.png)
As expected!
Flow Curves / Integral Curves#
Vector fields can also “determine” a particle’s motion in the plane or in 3D space (a flow curve). We image that we at time \(t=0\) throw a particle into a force field at the point \(\pmb{x}_0\), and we wish to model the particle’s trajectory (curve) \(\pmb{r}(t)\) as a function of time. These curves are called integral curves (and sometimes flow curves). They are solutions to the differential equation system:
where \(\pmb{x}_0\) is the initial point.
We consider the vector field
and two points that indicate the starting points \(s_1\) and \(s_2\) of two particles within a vector field \(V\)
s1, s2 = Matrix([5, 0]), Matrix([-3, S(1)/2])
Can we now find the flow curves of these two particles? We denote them, respectively, \(r_1\) and \(r_2\).
In this case the differential equation system is linear, so we know from the solution model from Mathematics 1a. We will first find the system matrix:
and we can use the eigenvalues and corresponding eigenvectors of \(A\) to determine the flow curves
A = Matrix(2, 2, [S("-1/4"), S("1/2"), S("1/2"), S("1/4")])
A.eigenvects(), A.diagonalize()
Both \(\verb|A.eigenvects()|\) and \(\verb|A.diagonalize()|\) provide the wanted info, and which one uses is just a question of preference.
With \(\verb|A.eigenvects()|\):#
k = symbols("k1:3")
ev = A.eigenvects()
ev
The general solution is set up as we know it from system of first-order differential equations (with two different eigenvalues):
r_general = (
k[0] * exp(ev[0][0] * u) * ev[0][-1][0] + k[1] * exp(ev[1][0] * u) * ev[1][-1][0]
)
r_general
The given starting points are our initial conditions, and the equations
are solved in order to find our specific flow curves
ks_r1 = solve(r_general.subs(u, 0) - s1)
ks_r2 = solve(r_general.subs(u, 0) - s2)
ks_r1, ks_r2
![../_images/b536e2116856544b59c54ea0218a8f47cff5889857f8b65d339d0b0b59d85a36.png](../_images/b536e2116856544b59c54ea0218a8f47cff5889857f8b65d339d0b0b59d85a36.png)
These can be substituted directly into \(\boldsymbol r(u)\)
r1 = r_general.subs(ks_r1)
r2 = r_general.subs(ks_r2)
and we arrive at the flow curves
r1, r2
v_field = dtuplot.plot_vector(
A * Matrix([x, y]),
n=10,
scalar=False,
colorbar=False,
quiver_kw={"color": "black"},
show=False,
)
r1_plot = dtuplot.plot_parametric(
*r1, [u, 0, 3], rendering_kw={"color": "red"}, colorbar=False, show=False
)
r2_plot = dtuplot.plot_parametric(
*r2, (u, 0, 4), rendering_kw={"color": "blue"}, colorbar=False, show=False
)
(v_field + r1_plot + r2_plot).show()
![../_images/65848e059fb52273711efeb2e5b3de94340ba880fa1a32b1fad2ee59e41b500f.png](../_images/65848e059fb52273711efeb2e5b3de94340ba880fa1a32b1fad2ee59e41b500f.png)
Solution via \(\verb|A.diagonalize()|\):#
k = symbols("k1:3")
matrix_of_evs, LAMBDA = A.diagonalize()
matrix_of_evs, LAMBDA
r_general = k[0] * exp(LAMBDA[0, 0] * u) * matrix_of_evs.col(0) + k[1] * exp(
LAMBDA[1, 1] * u
) * matrix_of_evs.col(1)
r_general
eq_k_r1 = Eq
ks_r1 = solve(r_general.subs(u, 0) - s1)
ks_r2 = solve(r_general.subs(u, 0) - s2)
ks_r1, ks_r2
![../_images/b536e2116856544b59c54ea0218a8f47cff5889857f8b65d339d0b0b59d85a36.png](../_images/b536e2116856544b59c54ea0218a8f47cff5889857f8b65d339d0b0b59d85a36.png)
and we have the same flow curves as above
r1 = r_general.subs(ks_r1)
r2 = r_general.subs(ks_r2)
r1, r2
An Example with the Gauss Theorem (not Syllabus)#
We are given the vector field
V = Matrix([-x + cos(z), -y * x, 3 * z + exp(y)])
V
as well as the spatial set \(\Omega\) given by
Here the set is illustrated:
a = symbols("a")
# Since one cannot plot a 3D volumen (not easily, that is), we will plot its surfaces,
p = dtuplot.plot3d_parametric_surface(
x,
y,
y**2,
(x, 0, 3),
(y, 0, 2),
{"color": "royalblue"},
use_cm=False,
aspect="equal",
show=False,
)
p.extend(
dtuplot.plot3d_parametric_surface(
0,
y,
a * y**2,
(a, 0, 1),
(y, 0, 2),
{"color": "royalblue", "alpha": 0.5},
use_cm=False,
aspect="equal",
show=False,
)
)
p.extend(
dtuplot.plot3d_parametric_surface(
3,
y,
a * y**2,
(a, 0, 1),
(y, 0, 2),
{"color": "royalblue", "alpha": 0.5},
use_cm=False,
aspect="equal",
show=False,
)
)
p.extend(
dtuplot.plot3d_parametric_surface(
x,
2,
a * 4,
(x, 0, 3),
(a, 0, 1),
{"color": "royalblue", "alpha": 0.5},
use_cm=False,
aspect="equal",
show=False,
)
)
p.extend(
dtuplot.plot3d_parametric_surface(
x,
y,
0,
(x, 0, 3),
(y, 0, 2),
{"color": "royalblue", "alpha": 0.5},
use_cm=False,
aspect="equal",
show=False,
)
)
p_field = dtuplot.plot_vector(
V,
(x, 0, 3),
(y, 0, 2),
(z, 0, 4),
n=4,
use_cm=False,
quiver_kw={"alpha": 0.5, "length": 0.05, "color": "red"},
show=False,
)
(p + p_field).show()
![../_images/0dfcfbc42e6d720945dc2bc7dafc53c9d30c51d7f08adc3bb091b5e5a3a7b95e.png](../_images/0dfcfbc42e6d720945dc2bc7dafc53c9d30c51d7f08adc3bb091b5e5a3a7b95e.png)
We will now determine the flux of the vector field \(\pmb{V}\) out through the closed surface of \(\Omega\) by use of Gauss’ Theorem!
We will first determine a parametric representation of the solid volume,
u, v, w = symbols("u v w")
r = Matrix([u, v, w * v**2])
r
M = Matrix.hstack(r.diff(u), r.diff(v), r.diff(w))
M, det(M)
We can compute the determinant since
We must, though, not forget the absolute value (see Definition 25.19). It might, though, often be possible to achieve a nicer-looking expression by not computing the absolute value if one can garantee that the determinant is always positive. That is the case here, since \(v>0\), so
jacobian = M.det()
jacobian
![../_images/991c555945927c0ff64df2492db151ba10a33d599e2b8ab6e7210bcb6d7b4312.png](../_images/991c555945927c0ff64df2492db151ba10a33d599e2b8ab6e7210bcb6d7b4312.png)
The divergence now becomes
divV = dtutools.div(V, (x, y, z))
divV
![../_images/3d979c74bba76ec3142337dedd3a40b43d9ae1fd679a088983cc3c6317c66082.png](../_images/3d979c74bba76ec3142337dedd3a40b43d9ae1fd679a088983cc3c6317c66082.png)
The divergence boudned to the parametrized region is thus found to be,
divV_r = divV.subs({x: r[0], y: r[1], z: r[2]})
divV_r
![../_images/0622d75a4b87d2722a3e749f60a24c567ab631e2208ce0281f1d397b3d7e18de.png](../_images/0622d75a4b87d2722a3e749f60a24c567ab631e2208ce0281f1d397b3d7e18de.png)
Since
integrate(divV_r * jacobian, (u, 0, 3), (v, 0, 2), (w, 0, 1))
![../_images/6e262f67411aeb20db2e81d76081079096b3452f9bd38dc0ff032aaf1e7f2e5f.png](../_images/6e262f67411aeb20db2e81d76081079096b3452f9bd38dc0ff032aaf1e7f2e5f.png)
we get from Gauss’ Theorem that:
In this example it was particularly obvious to use Gauss’ Theorem, both since we would avoid computing the flux through \(5\) surfaces, after which we would add up those partial results, and since the divergence is mathematically simple in comparison to its vector field. This must, though, be considered every time a new exercise is given!