Demo 11: Integration af vektorfelter gennem flader (flux)

Demo af Christian Mikkelstrup, Hans Henrik Hermansen, Jakob Lemvig, Karl Johan Måstrup Kristensen og Magnus Troen. Revideret marts 2026 af shsp.

from sympy import *
from dtumathtools import *

init_printing()

Fluxen gennem en flade \(F\) er jf. lærebogen et ortogonalt fladeintegral, og formlen forenkles til følgende:

(5)\[\begin{equation} \int_F\boldsymbol V\cdot \mathrm d\boldsymbol S = \int_a^b\int_c^d \pmb{V}(\boldsymbol r(u,v))\cdot \boldsymbol N\,\mathrm{d}u\,\mathrm{d}v. \end{equation}\]

Denne formel bruger den kanoniske normalvektor \(\boldsymbol N\), udregnet som et krydsprodukt af de partielt afledte af fladeparametriseringen, \(\boldsymbol N=\frac{\partial \boldsymbol r}{\partial u}\times \frac{\partial \boldsymbol r}{\partial v}\). Normalvektorens retning afgør fluxens fortegn, da vektorfeltvektorer langs med \(\boldsymbol N\) giver positive bidrag, da prikproduktet da bliver positivt, osv. Normalvektoren er derfor en indikator for vores valg af positiv orientering, så vi kan fortolke fortegnet på den resulterende flux korrekt.

Flux gennem lukket paraboloide

Vi får givet vektorfeltet:

x, y, z = symbols("x y z", real=True)
V = lambda x, y, z: Matrix([x, y, 2])
V(x, y, z)

\[\begin{split}\displaystyle \left[\begin{matrix}x\\y\\2\end{matrix}\right]\end{split}\]

Lad os betragte den paraboloide-formede container med låg, som vi i demo 8 parametriserede som:

u, v = symbols("u v", real=True)
r_vaeg = Matrix([u * cos(v), u * sin(v), u**2])
r_laag = Matrix([u * cos(v), u * sin(v), 1])
r_vaeg, r_laag

\[\begin{split}\displaystyle \left( \left[\begin{matrix}u \cos{\left(v \right)}\\u \sin{\left(v \right)}\\u^{2}\end{matrix}\right], \ \left[\begin{matrix}u \cos{\left(v \right)}\\u \sin{\left(v \right)}\\1\end{matrix}\right]\right)\end{split}\]

begge defineret for \(u\in [0,1],v\in [-\pi,\pi]\). Vektorfeltet og den lukkede overflade i samme plot:

p_vaeg = dtuplot.plot3d_parametric_surface(*r_vaeg,(u, 0, 1),(v, -pi, pi),
    use_cm=False,
    camera={"elev": 25, "azim": -55},
    show=False
)
p_laag = dtuplot.plot3d_parametric_surface(*r_laag,(u, 0, 1),(v, -pi, pi),
    show=False
)
p_felt = 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_vaeg + p_laag + p_felt).show()

Normalvektor

Lad os beregne fluxen (det ortogonale fladeintegral) af vektorfeltet ud igennem denne lukkede paraboloideformede container. Altså, vi vil have normalvektorerne til at pege udad fra ethvert punkt på begge dele af overfladen. Første skridt er derfor at tjekke normalvektorretningen for begge dele af overfladen.

Normalvektoren tilhørende den paraboloiodeformede væg:

n_vaeg = r_vaeg.diff(u).cross(r_vaeg.diff(v))
simplify(n_vaeg)

\[\begin{split}\displaystyle \left[\begin{matrix}- 2 u^{2} \cos{\left(v \right)}\\- 2 u^{2} \sin{\left(v \right)}\\u\end{matrix}\right]\end{split}\]

peger indad i containeren, så lad os bytte rundt på rækkefølgen af vektorerne:

n_vaeg = r_vaeg.diff(v).cross(r_vaeg.diff(u))
simplify(n_vaeg)

\[\begin{split}\displaystyle \left[\begin{matrix}2 u^{2} \cos{\left(v \right)}\\2 u^{2} \sin{\left(v \right)}\\- u\end{matrix}\right]\end{split}\]

Sådan. Alternativt kunne man blot ændre fortegnet på selve normalvektoren. Eller måske nemmere: husk indtil vi er færdige, at orienteringen er omvendt af det ønskede, og skift så fortegnet på fluxen til slut. Her har vi valgt at fikse det med en ombytning af vektorernes rækkefølge.

Nu til lågets normalvektor:

n_laag = r_laag.diff(u).cross(r_laag.diff(v))
simplify(n_laag)

\[\begin{split}\displaystyle \left[\begin{matrix}0\\0\\u\end{matrix}\right]\end{split}\]

Det peger opad, hvilket stemmer med en ønsket udadrettet orientering. Intet skal ændres her.

Fluxen

Lad os sammen med den paraboloideformede væg plotte kun de vektorer i vektorfeltet, der starter på vægfladen. Vi kan gøre dette med slice-argumentet på plot_vector:

flade = dtuplot.plot3d_parametric_surface(*r_vaeg, (u, 0, 1), (v, -pi, pi), n=8, show=False)

p_felt2 = dtuplot.plot_vector(V(x, y, z),(x, -1, 1),(y, -1, 1),(z, 0, 1),
    slice=flade[0],
    use_cm=False,
    quiver_kw={"alpha": 0.5, "length": 0.1, "color": "red"},
    show=False,
)
p_vaeg.camera = {"elev": 60, "azim": -50}

(p_vaeg + p_felt2).show()

Bemærk, hvordan vektorfeltets vektorer tegnes fra alle hjørner. Det bliver en smule overvældende med så mange, når p_vaeg vises samtidig, så vi har reduceret n en smule for at vise færre pile. På plottet ser det ud til, at flowet er opad gennem containeren. Med udadrettede normalvektorer vil vi derfor forvente en samlet set negativ flux gennem paraboloidevæggen og en samlet set positiv flux gennem låget.

Summen af de to fluxer giver den totale flux gennem containeren. Vi beregner dem hver for sig - først paraboloidevæggen:

integrand_vaeg = n_vaeg.dot(V(*r_vaeg))
Flux_vaeg = integrate(integrand_vaeg, (u, 0, 1), (v, -pi, pi))
Flux_vaeg

Fluxen ud gennem paraboloidevæggen er \(-\pi\). Ganske som forventet er den negativ.

Så låget:

integrand_laag = n_laag.dot(V(*r_laag))
Flux_laag = integrate(integrand_laag, (u, 0, 1), (v, -pi, pi))
Flux_laag

Fluxen gennem låget er \(2\pi\); positiv som forventet. Samlet flux:

Flux = Flux_laag + Flux_vaeg

Den samlede flux ud gennem den lukkede overflade er dermed \(\mathrm{Flux}(\pmb{V},F_{væg}\cup F_{låg})=-\pi+2\pi=\pi\). Den samlede flux er positiv, så der er et samlet udadrettet flow indefra containeren.

Flux gennem stykke af kugleskal

Vi får givet vektorfeltet:

x, y, z = symbols("x y z")
V = Matrix([-y + x, x, z])
V

\[\begin{split}\displaystyle \left[\begin{matrix}x - y\\x\\z\end{matrix}\right]\end{split}\]

og et stykke af en kugleflade:

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

\[\begin{split}\displaystyle \left[\begin{matrix}2 \sin{\left(u \right)} \cos{\left(v \right)}\\2 \sin{\left(u \right)} \sin{\left(v \right)}\\2 \cos{\left(u \right)}\end{matrix}\right]\end{split}\]

hvor \(u\in[\pi/6,\pi/2],v\in[0,\pi]\).

p_flade = dtuplot.plot3d_spherical(
    radius,
    (u, pi / 6, pi / 2),(v, 0, pi),
    aspect="equal",
    camera={"elev": 25, "azim": 55},
    show=False,
)
dummy_flade = dtuplot.plot3d_spherical(
    radius, (u, pi / 6, pi / 2), (v, 0, pi), show=False, n=8
)
p_felt = dtuplot.plot_vector(
    V,
    (x, -1, 1),(y, -1, 1),(z, 0, 1),
    slice=dummy_flade[0],
    use_cm=False,
    quiver_kw={"alpha": 0.5, "length": 0.2, "color": "red"},
    show=False,
)

(p_flade + p_felt).show()

Normalvektoren:

n = r.diff(u).cross(r.diff(v))
n

\[\begin{split}\displaystyle \left[\begin{matrix}4 \sin^{2}{\left(u \right)} \cos{\left(v \right)}\\4 \sin^{2}{\left(u \right)} \sin{\left(v \right)}\\4 \sin{\left(u \right)} \sin^{2}{\left(v \right)} \cos{\left(u \right)} + 4 \sin{\left(u \right)} \cos{\left(u \right)} \cos^{2}{\left(v \right)}\end{matrix}\right]\end{split}\]

Da vi ikke har nogen ønsket orientering, bliver justering af normalvektoren ikke relevant. Vi skal dog stadig vide, hvilken retning, den peger, så vi kan fortolke fortegnet på den resulterende flux korrekt. Normalvektoren peger her fra centrum og ud a kuglestykket.

Fluxen beregnes til:

integrand = n.dot(V.subs({x: r[0], y: r[1], z: r[2]}))
integrate(integrand, (u, pi / 6, pi / 2), (v, 0, pi))

En positiv værdi, så det samlede flow er dermed langs med normalvektorretningen, altså udad/væk fra kuglens centrum.

Flux gennem massiv halvkugle

Betragt i \(\Bbb R^3\) en massiv halvkugle med radius \(a\) parametriseret som følger, for \(u\in[0,a]\), \(v\in[0,\pi/2]\) og \(w\in[-\pi,\pi]\), samt et vektorfelt \(\boldsymbol V\):

x,y,z,u,v,w=symbols('x y z u v w')
a = symbols("a", real=True, positive=True)
V = Matrix([8 * x, 8, 4 * z**3])
r = u * Matrix([sin(v) * cos(w), sin(v) * sin(w), cos(v)])
V, r

\[\begin{split}\displaystyle \left( \left[\begin{matrix}8 x\\8\\4 z^{3}\end{matrix}\right], \ \left[\begin{matrix}u \sin{\left(v \right)} \cos{\left(w \right)}\\u \sin{\left(v \right)} \sin{\left(w \right)}\\u \cos{\left(v \right)}\end{matrix}\right]\right)\end{split}\]

Vi skal finde fluxen gennem dette massive legemes randflade, altså gennem dets lukkede overflade. Overfladen består af to stykker, som vi håndterer individuelt. Første fladestykke opnås ved, at \(u\) fikseres til radiusværdien \(a\):

r1 = r.subs(u, a)
n1 = r1.diff(v).cross(r1.diff(w))
integrand1 = n1.dot(V.subs({x: r1[0], y: r1[1], z: r1[2]}))
flux1 = integrate(integrand1, (v, 0, pi / 2), (w, -pi, pi))
r1,n1,flux1

\[\begin{split}\displaystyle \left( \left[\begin{matrix}a \sin{\left(v \right)} \cos{\left(w \right)}\\a \sin{\left(v \right)} \sin{\left(w \right)}\\a \cos{\left(v \right)}\end{matrix}\right], \ \left[\begin{matrix}a^{2} \sin^{2}{\left(v \right)} \cos{\left(w \right)}\\a^{2} \sin^{2}{\left(v \right)} \sin{\left(w \right)}\\a^{2} \sin{\left(v \right)} \sin^{2}{\left(w \right)} \cos{\left(v \right)} + a^{2} \sin{\left(v \right)} \cos{\left(v \right)} \cos^{2}{\left(w \right)}\end{matrix}\right], \ \frac{8 \pi a^{5}}{5} + \frac{16 \pi a^{3}}{3}\right)\end{split}\]

Det andet fladestykke svarer til en fiksering af \(v\) til værdien \(\pi/2\), hvilket resulterer i en skive med radius \(a\):

r2 = r.subs(v, pi / 2)
n2 = r2.diff(u).cross(r2.diff(w))
integrand2 = n2.dot(V.subs({x: r2[0], y: r2[1], z: r2[2]}))
flux2 = integrate(integrand2, (u, 0, a), (w, -pi, pi))
r2,n2,flux2

\[\begin{split}\displaystyle \left( \left[\begin{matrix}u \cos{\left(w \right)}\\u \sin{\left(w \right)}\\0\end{matrix}\right], \ \left[\begin{matrix}0\\0\\u \sin^{2}{\left(w \right)} + u \cos^{2}{\left(w \right)}\end{matrix}\right], \ 0\right)\end{split}\]

Er vi nu klar til at lægge de to fluxer sammen for at få den samlede flux gennem denne massive halvkugle? Pas på! Vi skal først sikre os, at deres normalvektorer stemmer overens - ellers vil de to fluxers fortegn ikke passe sammen.

Et hurtigt tjek afgør, at normalvektorerne for første fladestykke peger udad, mens de for det andet fladestykke peger indad. Altså stemmer fluxernes fortegn ikke. Det er dog blot et hurtigt fiks at ændre fortegnet på den ene af dem, fx fluxen gennem det andet fladestykke:

flux_gennem_halvkugle = flux1 + (-flux2)
flux_gennem_halvkugle

Dette resultat skal nu forstås med udadrettet positiv orientering. Bemærk, at det ikke er vigtigt, om normalvektorernes orientering er indad eller udad for en sådan lukket flade, så længe de stemmer overens!

Gauss’ divergenssætning (ikke pensum)

Resten af denne demo indeholder materiale, der kan være relevant for visse studerende i projektperioden i Matematik 1b, men som ikke er en del af kursets pensum.

Gauss’ sætning er en brugbar læresætning i forbindelse med fluxudregninger gennem lukkede flader. Den lader os finde fluxen gennem en lukket flade ved at “skubbe beregningen ind på det indre”, hvilket i visse tilfælde er lettere. Hvis vi mere specifikt betragter en lukket flade som overfladen/randen af et massivt legeme \(\Omega\), så siger Gauss’ sætning, at fluxen ud igennem denne randflade \(\partial \Omega\) er lig med rumintegralet af divergensen af vektorfeltet gennem \(\Omega\):

\[\mathrm{Flux}(\boldsymbol V,\partial \Omega)=\int_{\partial \Omega}\boldsymbol V\cdot \mathrm d\boldsymbol S=\int_{\Omega}\mathrm{div}(\boldsymbol V) \mathrm dS.\]

Flux gennem massiv halvkugle med Gauss’ sætning (ikke pensum)

Betragt igen vektorfeltet \(\boldsymbol V\) og den massive halvkugle parametriseret af \(\boldsymbol r\), som vi fandt fluxen igennem i forrige afsnit:

V, r

\[\begin{split}\displaystyle \left( \left[\begin{matrix}8 x\\8\\4 z^{3}\end{matrix}\right], \ \left[\begin{matrix}u \sin{\left(v \right)} \cos{\left(w \right)}\\u \sin{\left(v \right)} \sin{\left(w \right)}\\u \cos{\left(v \right)}\end{matrix}\right]\right)\end{split}\]

Her er \(u\in[0,a],v\in[0,\pi/2],w\in[-\pi,\pi]\). Vi fandt fluxen til at være:

flux_gennem_halvkugle

Da dette er en lukket flade, gælder Gauss’ sætning. Lad os se, om vi kan få samme resultat ved brug af Gauss’ sætning. Vi vil derfor nu udføre et rumintegral af vektorfeltets divergens gennem halvkuglens indre.

Til rumintegralet skal vi bruge Jacobimatricens determinant:

M = Matrix.hstack(r.diff(u), r.diff(v), r.diff(w))
jacobifunk = M.det()
jacobifunk

Absolutværdien skal findes, men vi ser, at determinanten altid er positiv, da det eneste led, der ikke er kvadreret, er \(\sin(v)\), som er positiv på \(v\)s interval. Nu vektorfeltets divergens:

divV = dtutools.div(V, [x, y, z])
divV_r = divV.subs({x: r[0], y: r[1], z: r[2]})
divV,divV_r

hvor vi i samme omgang lavede restriktionen af \(\mathrm{div}\boldsymbol V\) til parametriseringen. Nu kan rumintegralet udregnes, her som et tripelintegral:

integrate(divV_r * jacobifunk, (u, 0, a), (v, 0, pi / 2), (w, -pi, pi))

Vi får samme resultat, som forventet!

Bemærk, at vi her ikke har koncentreret os om orientering. Det skyldes, at Gauss’ sætning automatisk antager en udadrettet positiv orientering. At resultatet stemmer med det, vi fandt tidligere, skyldes, at vi netop endte med at vælge en udadrettet orientering i den tidligere flux-beregning.

Fluxen gennem en rampe (ikke pensum)

Vi får givet vektorfeltet

V = Matrix([-x + cos(z), -y * x, 3 * z + exp(y)])
V

\[\begin{split}\displaystyle \left[\begin{matrix}- x + \cos{\left(z \right)}\\- x y\\3 z + e^{y}\end{matrix}\right]\end{split}\]

og vi betragter det rumlige legeme \(\Omega\), som vi parametriserede i demo 8 på følgende måde, for \(u\in[0,3],v\in[0,2],w\in[0,1]\):

u, v, w = symbols("u v w")
r = Matrix([u, v, w * v**2])
r

\[\begin{split}\displaystyle \left[\begin{matrix}u\\v\\v^{2} w\end{matrix}\right]\end{split}\]

a = symbols("a")
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_felt = 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_felt).show()

Vi vil nu bestemme fluxen af vektorfeltet \(\pmb{V}\) ud gennem rampen \(\Omega\)s randoverflade. Rampens overflade består af mange mindre stykker, hvilket vil gøre fluxberegningen omstændig og tidskrævende ved brug af den sædvanlige fluxformel. Men Gauss’ sætning kan bruges på lukkede flader! Og Gauss’ sætning er langt enklere at bruge her. Altså, vi vil i stedet finde rumintegralet af vektorfeltets divergens gennem rampens indre.

Jacobimatricen og dens determinant:

M = Matrix.hstack(r.diff(u), r.diff(v), r.diff(w))
M, det(M)

\[\begin{split}\displaystyle \left( \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 2 v w & v^{2}\end{matrix}\right], \ v^{2}\right)\end{split}\]

Vi bemærker, at \(v>0\), hvilket betyder, at determinanten altid er positiv, så Jacobifunktionen er:

jacobifunk = M.det()
jacobifunk

Divergensen:

divV = dtutools.div(V, (x, y, z))
divV

Restriktion af divergensen til det rumlige parameterområde:

divV_r = divV.subs({x: r[0], y: r[1], z: r[2]})
divV_r

Rumintegralet:

integrate(divV_r * jacobifunk, (u, 0, 3), (v, 0, 2), (w, 0, 1))

Gauss sætning sikrer, at dette resultat er lig med fluxen:

(6)\[\begin{equation} \int_{\delta\Omega}\pmb{V}\cdot \mathrm{d}\pmb{S} = \int_0^1\int_0^2\int_0^3(2-u)\cdot v^2\,\mathrm{d}u\,\mathrm{d}v\,\mathrm{d}w = 4. \end{equation}\]

Dette eksempel tydeliggør, at Gauss’ sætning kan spare os en masse arbejde og tid: Vi undgik fem fluxberegninger gennem overfladens fem fladestykker, der hver ville have krævet fem parametriseringer med fem Jacobifunktioner, som så ville skulle lægges sammen for den samlede flux. I stedet dannede vi blot én parametrisering og fandt dens ene Jacobifunktion. Vi skulle dog udføre et tripelintegral i stedet for kun dobbeltintegraler, og vi skulle også finde divergensen. Om Gauss’ sætning kan spare dig tid, eller om du hellere vil holde dig til den sædvanlige metode til flux-udregninger, afhænger af situationen.