Demo 9: Plan-, rum-, kurve- og fladeintegraler med transformationssætningen

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()

Når skalarfunktioner skal integreres over ikke-akseparallelle områder som cirkelskiver, kugler og kurver, kan vi ikke direkte udføre sædvanlig Riemann-integration med hensyn til variablerne \(x_1,\cdots,x_n\). Lærebogen giver os dog en teknik kendt som transformationssætningen eller variabelsubstitution, der lader os skubbe integrationen over på parametervariablerne i stedet - hvis vi altså er i stand til at lave en regulær parameterfremstilling af området.

Denne demo behandler først integration over plane områder i \((x,y)\)-planet og gennem rumlige legemer i \((x,y,z)\)-rummet, altså hhv. planintegralet og rumintegralet. Dernæst forsætter demoen med integration langs kurver og over flader, altså hhv. kurveintegralet og fladeintegralet, i \(n\)-dimensionelle rum. Men før vi kan komme i gang med disse integrationer, tager vi et kig på en vigtig komponent i dem alle: restriktionen af funktionen til parameterområdet.

Restriktioner til parametriseret geometri

For at kunne integrere over vores parametriserede område ved brug af transformationssætningen skal vi først have helt styr på såkaldt restriktioner. Betragt følgende skalarfunktion \(f\) og vektorfunktion \(\boldsymbol r\):

\[\begin{equation*} f: \mathbb{R}^3 \to \mathbb{R}, \quad f(x,y,z) = x+y+z \end{equation*}\]

\[\begin{equation*} \boldsymbol{r}(u,v,w) = \begin{bmatrix} u -v \\ w^2 \\ v(u+w) \end{bmatrix}, \quad u,v,w \in \mathbb{R}. \end{equation*}\]

Vi vil gerne bestemme den sammensatte funktion \(f(\boldsymbol{r}(u,v,w))\). Hvis \(\boldsymbol r\) repræsenterer et parametriseret geometrisk område, kalder vi dette en restriktion af \(f\) til \(\boldsymbol r\). Det, vi vil have til at ske, er, at første element i \(\boldsymbol r\) erstatter \(x\), dets andet element erstatter \(y\) og så videre.

Lad os definere \(\boldsymbol r\) som en simpel vektor:

u,v,w = symbols('u v w')
r = Matrix([u -v, w**2, v*(u+w)])
r

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

\(f\) kan defineres som enten et Sympy-udtryk (dvs. lagret som en simpel variabel), en sædvanlig Python-funktion, eller som en lambda-funktion:

x,y,z = symbols('x,y,z', real =True)

# Sympy-udtryk
f_sym = x+y+z 

# Sædvanlig Python-funktion
def f_funktion(x,y,z):
    return x+y+z

# Lambda-funktion
f_lam = lambda x,y,z : x+y+z

En sammensætning som \(f(\boldsymbol{r}(u,v,w))\) bestemmes ved brug af forskellig syntaks alt efter dit valg blandt ovenstående:

fr_sym = f_sym.subs(zip((x,y,z), r)).simplify()
fr_funktion = f_funktion(*r).simplify()
fr_lam = f_lam(*r).simplify()

fr_sym, fr_funktion, fr_lam

Sympy-udtrykket gør brug af zip, der erstatter variablerne med vektorelementerne én for én - som en “lynlås” (eng: “zipper”) - således, at første vektorelement r[0] erstatter x, r[1] erstatter y, og r[2] erstatter z.

Dette er ikke nødvendigt ved Python-funktioner eller lambda-funktioner. Her ser vi i stedet asterisken * igen, placeret foran vektorsymbolet. Denne fortæller Python, at koordinaterne i r skal gives til funktionen som individuelle argumenter.

Planintegralet

Lad os nu gennemgå et eksempel på integration over et geometrisk område. Vi får givet en funktion \(f:\Bbb R^2\to\Bbb R\) ved:

\[\begin{equation*} f(x,y) = 2xy. \end{equation*}\]

f = lambda x,y: 2 * x * y
f(x,y)

Lad os som geometrisk område genbruge den plane cirkelskive \(C\), som vi fik parametriseret i forrige demo:

u,v = symbols('u v')
rC = Matrix([2,1]) + v * Matrix([cos(u),sin(u)])
rC

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

hvor \(u\in [0, 2\pi] ,v\in[0, 1]\). Den akseparallelle firkant \([0, 2\pi] \times [0, 1]\) i \((u,v)\)-planet afbildes til denne cirkelskive i \((x,y)\)-planet:

dtuplot.plot_implicit(((x < 2*pi) & (x > 0) & (y > 0) & (y < 1)),
                      (x, -0.5, 2*pi + 0.5), (y, -0.5,1.5),
                      title="(u,v)-planet",
                      xlabel="u", ylabel="v",
                      axis_center='auto')

/builds/pgcs/pg1/venv/lib/python3.11/site-packages/spb/series.py:2255: UserWarning: The provided expression contains Boolean functions. In order to plot the expression, the algorithm automatically switched to an adaptive sampling.

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b05f7c3d0>

dtuplot.plot_parametric_region(
    *rC, (u, 0, 2*pi), (v, 0, 1),
    aspect="equal",
    title="(x,y)-planet"
)

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b03f2c610>

Vores opgave er nu at bestemme planintegralet af \(f\) over \(C\). I matematisk notation skal vi udregne:

\[\int_C f(x,y)\; \mathrm{d}\pmb{x}.\]

Cirkelskiven er ikke akseparallel, så vi får brug for transformationssætningen, som lærebogen giver os følgende formel for, når der er tale om et planintegral:

\[\begin{equation*} \int_C f(x,y)\; \mathrm{d}\pmb{x} = \int_{0}^{1} \int_{0}^{2\pi} f(\boldsymbol r_C(u,v))\cdot|\det (\boldsymbol J_{\boldsymbol r_C}) |\;\mathrm{d}u\,\mathrm{d}v. \end{equation*}\]

Jacobimatricen \(\boldsymbol J_{\boldsymbol r_C}\) er defineret ved:

\[\begin{equation*} \boldsymbol{J}_{\boldsymbol r_C} = \left[\begin{matrix} \vdots & \vdots \\ \frac{\partial \boldsymbol r(u,v)}{\partial u} & \frac{\partial \boldsymbol r(u,v)}{\partial v} \\ \vdots & \vdots \end{matrix}\right] = \left[\begin{matrix}\nabla(2 + v\cos\left(u \right))^T\\ \nabla(1 + v\cos\left(u \right))^T\end{matrix}\right]. \end{equation*}\]

JacobiM = Matrix.hstack(rC.diff(u), rC.diff(v))
JacobiM

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

Absolutværdien af Jacobideterminanten \(|\det(\boldsymbol J_{\boldsymbol r_C})|\) er

Jacobifunk = abs(det(JacobiM)).simplify()
Jacobifunk

Denne værdi kan opfattes som et mål for, hvor meget \((u,v)\)-parameterområdet deformeres lokalt, som det gennemgår afbildningen \(\boldsymbol r_C\). Vores planintegral er hermed:

\[\begin{equation*} \int_C f(x,y)\; \mathrm{d}\pmb{x} = \int_{0}^{1} \int_{0}^{2\pi} f(2+v\cos(u),1+v\sin(u))\cdot|v|\; \mathrm{d}u \,\mathrm{d}v, \end{equation*}\]

og vi kan færdiggøre berengingen i én kodelinje:

integrate(f(*rC) * Jacobifunk,(u,0,2*pi),(v,0,1))

Rumintegralet

I forrige demo parametriserede vi den massive enhedskugle \(\Omega\) centreret i origo, for \(u\in[0,2\pi], v\in[0,\pi], w\in[0,1]\):

u,v,w = symbols('u v w')
r_omega = Matrix([[w*cos(u)*cos(v)],[w*sin(v)*cos(u)],[w*sin(u)]])
r_omega

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

Vi vil i det følgende beregne rumintegraler gennem denne kugle. Da dette ikke er et akseparallelt område i \((x,y,z)\)-rummet, har vi igen brug for, at transformationssætningen skubber integrationen over på det akseparallelle parameterområde i \((u,v,w)\)-rummet:

\[\begin{equation*} \int_{\Omega} f\; \mathrm{d}\pmb{x}= \int_{0}^{1} \int_{0}^{\pi}\int_{0}^{2\pi} f(\boldsymbol{r}_{\Omega})\,\,|\mathrm{det}(\boldsymbol J_{\boldsymbol r_{\Omega}})| \;\mathrm du\,\mathrm dv\,\mathrm dw. \end{equation*}\]

Kuglens rumfang

Jf. lærebogen kan vi finde “størrelsen på geometrien”, hvilket i dette tilfælde med en tre-dimensionel kugle vil være dens rumfang (eller volumen) ved at lade funktionen \(f:\mathbb{R}^3 \to \mathbb{R}\) i transformationssætningen antage den konstante værdi \(f(\boldsymbol{x}) = 1\). Dette forenkler rumintegralformlen til:

\[\begin{equation*} \mathrm{Rumfang}=\int_{\Omega} 1\; \mathrm{d}\pmb{x}= \int_{0}^{1} \int_{0}^{\pi}\int_{0}^{2\pi} |\mathrm{det}(\boldsymbol J_{\boldsymbol r_{\Omega}})| \;\mathrm du\,\mathrm dv\,\mathrm dw. \end{equation*}\]

Jacobimatricen af vores \(\Omega\)-parameterfremstilling er

\[\begin{equation*} \boldsymbol{J}_{\boldsymbol r_{\Omega}} = \left[\begin{matrix} \vdots & \vdots & \vdots \\ \frac{\partial \boldsymbol r_{\Omega}(u,v,w)}{\partial u} & \frac{\partial \boldsymbol r_{\Omega}(u,v,w)}{\partial v} & \frac{\partial \boldsymbol r_{\Omega}(u,v,w)}{\partial w} \\ \vdots & \vdots & \vdots \end{matrix}\right] = \left[\begin{matrix}\nabla(w\cos{\left(u \right)} \cos{\left(v \right)})^T\\ \nabla(w \cos{\left(u \right)}\sin{\left(v \right)})^T\\ \nabla(w\sin{\left(u \right)})^T\end{matrix}\right], \end{equation*}\]

og den findes nemt med dtumathtoolss jacobian-kommando:

J_omega = r_omega.jacobian([u,v,w])
J_omega

\[\begin{split}\displaystyle \left[\begin{matrix}- w \sin{\left(u \right)} \cos{\left(v \right)} & - w \sin{\left(v \right)} \cos{\left(u \right)} & \cos{\left(u \right)} \cos{\left(v \right)}\\- w \sin{\left(u \right)} \sin{\left(v \right)} & w \cos{\left(u \right)} \cos{\left(v \right)} & \sin{\left(v \right)} \cos{\left(u \right)}\\w \cos{\left(u \right)} & 0 & \sin{\left(u \right)}\end{matrix}\right]\end{split}\]

Absoluteværdien af Jacobideterminanten er:

Jacobifunk = abs(det(J_omega)).simplify()
Jacobifunk

og vi kan nu beregne rumintegralet for at finde rumfanget:

V = integrate(Jacobifunk, (u, 0, 2*pi), (v, 0, pi), (w, 0, 1))
V

Sammenlign dette resultat med rumfangsformlen for en enhedskugle, som kendes fra grundlæggende geometri.

Massen af kuglen

Vi vil nu fortolke funktionen \(f:\mathbb{R}^3 \to \mathbb{R}\) som en massetæthedsfunktion (også blot kaldet en densitetsfunktion) for vores kuglelegeme \(\Omega\). En sådan massetæthed vil have en passende fysisk enhed så som \(\mathrm{kg/m^3}\), der ikke er vigtige for os i matematisk sammenhæng. Med denne \(f\) er vi i stand til at beregne \(\Omega\)s samlede masse. Dette er blot endnu et rumintegral af \(f\) over geometrien, som igen udføres vha. transformationssætningen.

Massetæthedsfunktionen vi vælge at gøre brug af i dette eksempel er:

\[\begin{equation*} f : \mathbb{R}^3 \to \mathbb{R}, \quad f(\boldsymbol{x}) = 2||\boldsymbol x||. \end{equation*}\]

f = lambda x,y,z: 2*sqrt(x**2 + y**2 + z**2)

Jacobimatricen for \(\boldsymbol r_\Omega\) fra før, og dermed også absolutværdien af dens Jacobideterminant, kan genbruges her, da den knytter sig til parametriseringen og ikke til funktionen. Denne gang er \(f\) ikke lig \(1\), så vi skal huske at inkludere restriktionen af \(f\) til \(\boldsymbol r_\Omega\) i integranden:

integrand = (f(*r_omega) * Jacobifunk).simplify()
integrand

Det nye rumintegral, som denne gang repræsenterer kuglens samlede masse, er:

M = integrate(integrand, (u, 0, 2*pi), (v, 0, pi), (w, 0, 1))
M

Kurveintegralet

Et kurveintegral, \(\int_K\;f\, \mathrm{d}\pmb{s}\), af en skalarfunktion \(f\) langs en kurve \(K\) kan udføres vha. transformationssætningen i lærebogen, der lader os skubbe integrationen over på parameterintervallet, hvis vi blot har en ordentlig parameterfremstilling af kurven klar. Da bliver formlen for kurveintegralet følgende:

\[\begin{equation*} \int_K\;f\, \mathrm{d}\pmb{s} = \int_{a}^{b} \,f(\boldsymbol r(u))\,\Vert \boldsymbol r'(u) \Vert \mathrm{d}u. \end{equation*}\]

Kurvelængde

Vi får givet en parameterkurve for \(u \in [0, 2\pi]\):

u = symbols('u', real=True)
r = Matrix([sin(u), sin(u)*cos(u)])
r

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

p_kurve = dtuplot.plot_parametric(*r, (u,0,2*pi), use_cm=False, label="r(u)",axis_center="auto")

Vi vil nu gerne finde længden af denne usædvanlige kurve, som vi ikke har nogen elementær geometrisk formel for. Jf. lærebogen kan længden bestemmes som et kurveintegral, hvis vi sætter \(f(x)=1\).

\[\begin{equation*} L=\int_K\; \mathrm{d}\pmb{s} = \int_{0}^{2\pi} \Vert \boldsymbol r'(u) \Vert \mathrm{d}u. \end{equation*}\]

Tangentvektoren \(\boldsymbol r'(u)\) er:

dr = r.diff(u)
dr

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

og dens norm findes med denne passende kommando fra dtutools:

Jacobifunk = dtutools.l2_norm(dr)

Kurvens længde er dermed:

L=integrate(Jacobifunk, (u,0,2*pi)).n()

hvor .n() tvinger resultatet til at blive vist som decimaltal. Altså er kurven ca. 6 enheder lang.

Tangentlinjer

Da vi nu har fundet tangentvektoren, kan vi også let konstruere tangentlinjen fra ethvert punkt på kurven, hvis vi skulle få brug for fx til et flot plot. En sådan tangentlinje konstrueres ganske enkelt som enhver anden ret linje ved et punkt derpå plus en retningsvektor gange en parameter:

\[\text{tangentlinje}=\boldsymbol r(u_0)+t\cdot \boldsymbol r'(u_0)\,,\quad t\in\Bbb R.\]

Tangentlinjen med berøringspunktet \(x_0=r(\pi/3)\) som et eksempel parametriseres dermed som følger, for \(t\in \Bbb R\):

t = symbols("t")
r_tan = r.subs(u,pi/3) + t*dr.subs(u,pi/3)
r_tan

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{t}{2} + \frac{\sqrt{3}}{2}\\- \frac{t}{2} + \frac{\sqrt{3}}{4}\end{matrix}\right]\end{split}\]

p_punkt = dtuplot.scatter(r.subs(u,pi/3), show=False)
p_tan = dtuplot.plot_parametric(*r_tan, (t,-1,1), use_cm=False, label="r'(pi/3)", show=False)

(p_kurve + p_punkt + p_tan).show()

Kurveintegral langs rumkurve

Vi får nu givet en funktion:

x,y,z = symbols("x y z")
f = lambda x,y,z: sqrt(x**2 + y**2 + z**2)
f(x,y,z)

og en parameterkurve \(K\) for \(u\in[0,5]\) - da vi nu er i et 3D rum, kaldes en sådan kurve også for en rumkurve:

u = symbols('u', real=True)
r = Matrix([u*cos(u), u*sin(u), u])
r

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

p_rumkurve = dtuplot.plot3d_parametric_line(*r, (u,0,2*pi), use_cm=False, label="r(u)",aspect="equal", legend=True)

Lad os beregne kurveintegralet:

\[\int_K f(x,y,z)\, \mathrm{d}\pmb{s}= \int_{0}^{5} \,f(\boldsymbol r(u))\,\Vert \boldsymbol r'(u) \Vert \mathrm{d}u.\]

Tangentvektoren og dens norm (altså vores Jacobifunktion) er:

dr = r.diff(u)
Jacobifunkn = dr.norm().simplify()
Jacobifunk

Bemærk: Hvis ikke u var defineret som en reel variabel med symbols('u', real=True), så ville man ikke kunne bruge dr.norm(), og man ville i stedet skulle bruge dtutools.l2_norm(dr).

Restriktionen:

restriktion = f(*r).simplify()
restriktion

Alle komponenter er nu klar, og vi beregner kurveintegralet langs kurven:

integrate( restriktion * Jacobifunk ,(u,0,5)).evalf()

Skulle vi få brug for kurvens længde, så udelades blot funktionen fra formlen:

integrate(Jacobifunk,(u,0,5)).evalf()

Antagelser

Bemærk, at vores restriktion i ovenstående kurveintegraludregning indeholder et absolutværdi-symbol, \(|u|\):

restriktion

Sympy vil oftest være ligeglad og vil stadig kunne udføre integrationen uden problemer, som vi så ovenfor. Men hvis integrationen er mere avanceret, og Sympy kæmper med den i længere tid, så kan vi give Sympy en hjælpende hånd ved manuelt at håndtere denne absolutværdi først.

I dette tilfælde har vi \(u\in [0, 5]\), hvorfor \(u\) altid vil være positiv. Absolutværdi-symbolet kan derfor fjernes uden problemer. Det ved Sympy ikke fra definition u = symbols('u', real=True), så vi redefinerer til

u = symbols('u', real=True, nonnegative=True)

restriktion

Hvis Sympy ikke automatisk reducerer absolutværdi-linjerne væk, så brug refine():

refine(restriktion, Q.nonnegative(u))

Q.nonnegative(u) gennemtvinger antagelsen om ikke-negativitet. For andre antagelser erstattes .nonnegative()-delen med en af mulighederne i denne tabel.

Fladeintegralet

Betragt funktionen \(f: \mathbb{R}^3 \to \mathbb{R}\) givet ved

\[\begin{equation*} f(x,y,z) = 8 z. \end{equation*}\]

x,y,z = symbols('x y z', real=True)
def f(x,y,z):
    return 8*z

f(x,y,z)

Betragt også en flade \(S\) givet ved følgende parameterfremstilling, hvor \(u \in [0,\frac{\pi}{2}],v \in [0,1]\):

u,v = symbols('u v', real=True, nonnegative=True)
r = Matrix([u*cos(u),u*sin(u),u*v])
r

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

Denne gang huskede vi nonnegative=True, da vi endnu engang har en situation, hvor hverken \(u\) eller \(v\) kan være negative. Denne geometriske form kaldes cylinderflade, da den strækker sig opad som et stykke af en cylindrisk “væg”:

dtuplot.plot3d_parametric_surface(*r,(u,0,pi/2),(v,0,1), aspect='equal')

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b51ccf550>

Vores opgave er at bestemme fladeintegralet af \(f\) over fladen \(S\). Til det formål anvender vi os igen af transformationssætningen, som jf. lærebogen i tilfældet af et fladeintegral giver os følgende formel:

\[\begin{equation*} \int_S f\; \mathrm{d}\pmb{A} = \int_{0}^{1} \int_{0}^{\pi/2} f(\boldsymbol r_S(u,v))\cdot||\frac{\partial \boldsymbol r_C}{\partial u}\times \frac{\partial \boldsymbol r_C}{\partial v}||\;\mathrm{d}u\,\mathrm{d}v. \end{equation*}\]

Jacobifunktionen er denne gang normen af krydsproduktet mellem de to partielt afledte af parametriseringen:

krydsprodukt = r.diff(u).cross(r.diff(v))
Jacobifunk = sqrt((krydsprodukt.T * krydsprodukt)[0]).simplify()
Jacobifunk

Integranden:

integrand = f(*r) * Jacobifunk
integrand

Fladeintegralet:

integrate(integrand,(v,0,1),(u,0,pi/2)).evalf()

Havde vi ønsket at finde arealet af fladen, så skulle vi igen blot udelade funktionen:

A=integrate(Jacobifunk,(v,0,1),(u,0,pi/2)).evalf()
A