Demo 7: Riemann-integraler

Demo af Karl Johan Funch Måstrup Kristensen, Magnus Troen og Jakob Lemvig. Revideret marts 2026 af shsp.

from sympy import *
from dtumathtools import *

init_printing()

I integrationsarbejde er der ingen vej uden om Sympys integrate()-kommando. De første to afsnit i denne demo vil vise brugen af denne kommando til at finde stamfunktioner og bestemte integraler.

I mere avancerede tilfælde af integration kan integrate() dog have svært ved at udregne et integral. Du vil da skulle anvende din viden fra kurset til at gøre udtrykket mere “spiseligt” for Sympy; af og til skal Sympy have det skåret ud i mindre bidder for ikke at kvæles og give op, andre gange skal en helt anden og alternativ metode benyttes, måske endda en numerisk tilgang. De sidste to afsnit vil se på tilgange til at håndtere sådanne integraler, der er for avancerede til at blive løst med integrate() direkte.

Integration i én dimension

Det ubestemte integral (stamfunktioner)

Som lærebogen beskriver, er stamfunktioner til en funktion \(f\) enhver funktion \(F\), der opfylder \(\mathrm dF(x)/\mathrm dx=f(x)\), skrevet med integrationsnotation som:

\[\begin{equation*} \int f(x)\; \mathrm dx=F(x). \end{equation*}\]

Glem ikke, at hvis man lægger en konstant til en hvilken som helst stamfunktion, så får man en anden stamfunktion. Lægges en vilkårlig konstant til en hvilken som helst stamfunktion, så fås et udtryk for alle stamfunktioner (det, vi typisk kalder det ubestemte integral), \(F(x)+c\,,\,c\in\Bbb R\).

Lad os som et eksempel integrere funktionen \(f:\mathbb{R} \to \mathbb{R}\) givet ved

\[\begin{equation*} f(x) = \mathrm e^x \sin(x). \end{equation*}\]

x = symbols("x", real=True)
f = exp(x)*sin(x)
f

Sympy giver altid den mindste stamfunktion (den uden nogen konstant lagt til) med kommandoen integrate(integrand, variabel). Bemærk, at variabel-argumentet kan udelades, når funktionen kun har én variabel.

F = integrate(f, x)
F

plot(f, F, (x,-5,5), legend=True)

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7f3acc756f50>

Som nævnt tilføjer Sympy ikke konstanter efter en integration. Du kan altid manuelt tillægge en vilkårlig konstant \(C\in\Bbb R\), hvis du skal bruge det fulde ubestemte integral, der repræsenterer alle stamfunktioner:

C = symbols('C')
F = integrate(f, x) + C
F

Et bestemt integral

Et bestemt integral over et interval \([a,b]\subset \Bbb R\) er jf. integralregningens fundamentalsætning (se lærebogen) forskellen i en stamfunktions værdier, når den evalueres i de to punkter \(a,b\):

\[\begin{equation*} \int_a^b f(x)\; \mathrm dx=F(b)-F(a). \end{equation*}\]

Vi kalder typisk denne værdi for “arealet under grafen” i intervallet. Lad os sige, at vi ønsker at udregne:

\[\begin{equation*} \int_{-1}^1 \mathrm e^x \sin(x)\; \mathrm dx. \end{equation*}\]

Med Sympy tilføjes blot interval-endepunkterne til kommandoen: integrate(integrand, (variabel, a, b)):

integrate(f, (x, -1, 1))  ,  integrate(f, (x, -1, 1)).evalf()

Alternativt kan vi manuelt udføre subtraktionen, da vi allerede har bestemt en stamfunktion:

F.subs(x, 1) - F.subs(x, -1)

Integration i flere dimensioner

Har man med funktioner af adskillige variabler at gøre, ser man ofte gentagne integrationer (når der er tale om akseparallelle områder fx). En dobbeltintegration,

\[\begin{equation*} \int\int f(x_1, x_2)\; \mathrm dx_1 \; \mathrm dx_2, \end{equation*}\]

kan udføres i Sympy ved integrate(integrand, x1, x2). Søges det bestemte integral over givne intervaller for variablerne,

\[\begin{equation*} \int_{a_1}^{b_1}\int_{a_2}^{b_2} f(x_1, x_2)\; \mathrm dx_1 \; \mathrm dx_2, \end{equation*}\]

skal der simpelthen bare tilføjes flere variabel-tupler i kommandoen med syntaksen integrate(integrand, (x1, a1, b1), (x2, a2, b2)).

Lad os betragte funktionen \(f:\Bbb R^2\to\Bbb R\) givet ved \(f(x_1,x_2)=(x_1x_2)^3 - x_1^2 + \sin(3 x_2) + x_1x_2^2 + 4\):

x1, x2 = symbols('x1 x2')
f = (x1*x2)**3 - x1**2 + sin(3 * x2) + x1*x2**2 + 4
f

p = dtuplot.plot3d(f, (x1, -1.2, 1.2), (x2, -1.2, 1.2), zlim=(-2, 6), use_cm=True, colorbar=False,
                   camera = {'azim': 25, 'elev': 20}, wireframe=True)

Det ubestemte integral er, for \(C\in\Bbb R\):

integrate(f, x1, x2) + C

Et bestemt integral ville i dette tilfælde med to variable svare til “rumfanget under grafen”. Det bestemte integral i området \((x_1,x_2)\in[-1.2,1.2]^2\) som et eksempel er:

integrate(f, (x1, -1.2, 1.2), (x2, -1.2, 1.2))

Tips og tricks til at få Sympy til at integrere

Når Sympy løber ind i problemer med et besværligt integral, følger her nogle tips, der ofte kan hjælpe. Problemerne kan resultere i, at Sympy ikke giver et brugbart resultat eller bruger alt for lang tid på en udregning, der muligvis aldrig ender.

Definér variabler med antagelser/afgrænsninger

Fortæl Sympy, hvad du ved om dine variabler, når du definerer dem som symboler. Er det fx en positiv, reel variabel, så fortæl det til Sympy. De mest brugbare afgrænsninger er følgende:

• real = True

• positive = True

• negative = True

• nonnegative = True

• nonzero = True

Nogle eksempler:

x1 = symbols('x1', real = True, nonnegative = True) # x1 er reel og ikke negativ
x2 = symbols('x2', nonzero = True) # x2 er forskellig fra nul

Forenkl før der integreres

Af og til kan det hjælpe at kalde simplify på integranden, før der integreres. Det kan nogle gange reducere tiden, Sympy kræver for at gennemføre integrationen.

Workflow’et vil være noget à la dette:

x1, x2 = symbols('x1 x2',real=True)
f = 1 / (x1 + x2)

integrand = simplify(f)
F = integrate(integrand, x1, x2)
F

Udnyt integralernes linearitet

For at kunne integrere udtrykket

\[\begin{equation*} \int_{c}^{d}\int_{a}^{b} \alpha f_1(x,y) + \beta f_2(x,y) \:\mathrm dx\,\mathrm dy, \end{equation*}\]

vil det ofte være behjælpeligt først at omskrive til:

\[\begin{equation*} \alpha \int_{c}^d\int_{a}^{b} f_1(x,y) \,\mathrm dx\,\mathrm dy + \beta \int_{c}^{d}\int_{a}^{b} f_2(x,y)\,\mathrm dx\,\mathrm dy. \end{equation*}\]

Integrationen udføres nu i to dele, hvorefter man manuelt skal samle det igen. Som et eksempel kan vi forsøge at udregne:

\[\begin{equation*} \int_{2}^{5}\int_{3}^{10} \frac{27}3 \cos(x) - 4\, \mathrm{atan}(y -x) \:\mathrm dx\,\mathrm dy. \end{equation*}\]

x,y = symbols('x y')
f = S(27)/3 * cos(x) - 4 * atan(y -x)

# Integralernes grænser
a = 3
b = 10
c = 2
d = 5

Forsøger du at lade Sympy beregne hele integralet i én omgang som i følgende kodelinje, vil du nok hurtigt løbe tør for tålmodighed, og du vil skulle afbryde ved at stoppe Python-kernen i VS Code. Eller måske nægter Sympy simpelthen at give et resultat.

#integrate(f, (x, a, b), (y, c, d)) # Kør ikke denne linje med mindre du gerne vil opleve problemet

Så lad os i stedet opdele integrationen ved at gøre brug af lineariteten, hvilket vil resultere i to mindre integraler:

\[\begin{equation*} \frac{27}3\int_{2}^{5}\int_{3}^{10} \cos(x) \:\mathrm dx\,\mathrm dy - 4 \int_{2}^{5}\int_{3}^{10} \mathrm{atan}(y -x) \:\mathrm dx\,\mathrm dy. \end{equation*}\]

integrand1 = simplify(cos(x)) # Her er simplify faktisk ikke nødvendig
integrand2 = simplify(atan(y - x))

F = S(27)/3 * integrate(integrand1, (x, a,b), (y, c,d)) - 4 * integrate(integrand2, (x, a,b), (y, c,d))
F

Lad os evaluere dette som decimaltal for et brugbart resultat:

F.evalf()

Find en stamfunktion først og indsæt grænser bagefter

Nogle gange vil Sympy ikke evaluere et bestemt integral:

f = tan(ln(x**2))**3 / x

integrate(simplify(f), (x,3,10))

Men Sympy kan godt finde en stamfunktion:

F = integrate(f, x)
F

Derfor kan vi omgå problemet ved manuelt at udregne det bestemte integral vha. integralregningens fundamentalsætning:

F_definite = F.subs(x, 10) - F.subs(x, 3)
F_definite

F_definite.evalf()

Brug .doit()-kommandoen

Hvis Sympy giver et output med en ufærdig integration, kan vi i nogle tilfælde tvinge Sympy til at evaluere integratet med integrate(integrand, (variabel,a,b)).doit(). Dette risikerer dog at give fejl, men det er et forsøg værd.

Numeriske metoder

I nogle tilfælde bruger Sympy alligevel meget lang tid eller kan overhovedet ikke løse en integration, selv efter du har hjulpet den alt, hvad du kan. Da har du ikke andet valg end at løse integralet numerisk. Sådanne numeriske metoder kræver dertilegnede numeriske biblioteker så som Numpy og Scipy. Har du dem ikke allerede installeret fx fra andre kurser, så gør det nu via en terminal (bemærk, at dtumathtools automatisk installerer numpy-pakken):

# conda install scipy  # Hvis du bruger pip3, så skriv pip i stedet for conda

Scipy-biblioteket indeholder kommandoerne quad, som er en numerisk integrator til funktioner af én variabel, og nquad, som er til funktioner af flere variable. Lad os hente disse ind og anvende dem i de følgende delafsnit:

from scipy.integrate import quad, nquad
import numpy as np

Numerisk integration i én dimension

Lad os tage et kig på integralet

\[\begin{equation*} \int_{\pi}^{10} \frac{\tan^3(\ln(x^2))}{x} \,\mathrm dx. \end{equation*}\]

f = tan(ln(x**2))**3 / x

integrate(f, (x,pi,10))

Sympy vil ikke evaluere integralet, så vi vil forsøge at finde en numerisk løsning.

Først konverterer vi \(f\) til en Python-funktion, specifikt en lambda-funktion. Dette gøres med kommandoen lambdify, der groft sagt er en genvej til at definere \(f\) med def f(x)::

f_num = lambdify(x, f, 'numpy')
f_num(3)

En sådan numerisk defineret funktion kan ikke vises/printes på skærmen med sin forskrift (fx vil det give en fejl at taste f_num, selv hvis du først definerer x som symbol). Men vi kan godt evaluere den i udvalgte punkter som gjort ovenfor med 3.

Nu kan quad bruges til at bestemme integralet numerisk:

F_num, fejl = quad(f_num, np.pi, 10)
F_num, fejl

Bemærk, at det er vigtigt, at vi indsætter Numpys version af \(\pi\), hvorfor vi skrev np.pi.

Funktionen quad giver to outputs: Først den numeriske approksimation af integralet og dernæst en estimeret fejl (afvigelse) fra den eksakte løsning.

Denne funktion er faktisk en af dem, som Sympy kan finde en analytisk løsning til, hvis vi blot hjælper lidt til. Lad os sammenligne den numeriske løsning med en analytisk løsning. Vi lægger ud med at forsøge at tvinge Sympy til at evaluere integralet med .doit() - det returnerer dog bare en fejl:

# integrate(f, (x, pi, 10)).doit() # Fjern udkommenteringen for at se fejlen

Lad os derfor i stedet se, om ikke Sympy kan give os en stamfunktion, som vi kan bruge til manuelt at udregne det bestemte integral:

F_stam = integrate(f, x)
F_stam

F_analytisk = trigsimp(F_stam.subs(x,10) - F_stam.subs(x,pi))
F_analytisk

Det virkede! (trigsimp-kommandoen er en trigonometric simplifier, der tvinger trigonometriske led til at gå ud med hinanden og reducere mest muligt vha. trigonometriske relationer.)

Lad os sammenligne de to resultater:

F_analytisk.evalf() - F_num

Den numeriske løsning giver faktisk en ganske god approksimation, men de to resultater er ikke helt ens. Bemærk, at den nøjagtige fejl, som quad indfører, er markant mindre end den estimerede fejl på \(3.4 \cdot 10^{−11}\), så det er ikke så slemt igen.

Numerisk integration i to dimensioner

Som et sidste eksempel vil vi arbejde på integralet:

\[\begin{equation*} \int_{0}^{\pi}\int_{-10}^{20} \sin(x)\cdot y^2\: \mathrm dy \,\mathrm dx. \end{equation*}\]

Her har Sympy ingen problemer med at finde resultatet, så lad os sammenligne det følgende analytiske resultat med en numeriske metode:

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

f = sin(x)*y**2

F_analytisk = integrate(f, (y,-10,20), (x,0,pi))
F_analytisk

For en numerisk evaluering af integralet definerer vi igen \(f\) som en Nympy-funktion vha. lambdify:

f_num = lambdify((x,y), f, 'numpy')   # lambdify kan tage flere variabler
f_num(3,4)

Integralet kan nu bestemmes med nquad. For en funktion med \(n\) variable, kaldes nquad med syntaksen nquad(function, [[interval til variabel_1], [interval til variabel_2], ..., [interval til variabel_n]]). I vores tilfælde med to variabler skriver vi:

F_num, fejl = nquad(f_num, [[0, np.pi], [-10,20]])
F_num, fejl

En sammenligning af de to metoder:

F_num - F_analytisk