12.1 Terület, összeg és numerikus integrál

12.1.1 Cél és konkrét helyzet: egy teljes mennyiség területként

A Northstar rendeléseihez tartozó szimulációban egy kiszállító cég azt vizsgálja, hogy egy új útvonalon mennyi időt tölt a furgon, és mekkora utat tesz meg. A sebesség nem állandó: a furgon először lassan halad, majd a forgalom ritkulásával gyorsul. Legyen a sebesség időfüggvénye szemléltetésként

$$v(t)=0{,}5+0{,}02\,t,$$

ahol $t$ a percben mért idő, $v(t)$ pedig km/percben mért sebesség. A teljes megtett utat a sebesség alatti terület adja a 0 és 60 perc közötti intervallumon. Ebben az alfejezetben a kérdés: hogyan lehet ezt a területet számítógéppel, rácsfinomítással közelíteni, és honnan tudjuk, hogy a közelítés milyen hibájú?

Az induló gyakorlás egy még egyszerűbb, kézzel ellenőrizhető matematikai modellen futtatjuk le:

• bemenet: $f(x)=x^2$ és az $[a,b]=[0,1]$ intervallum;
• megfigyelési egység: egy kiértékelési pont $x_i\in[0,1]$;
• cél: az $\int_0^1 x^2\,dx$ terület közelítése három rácsfinomsággal;
• ellenőrzési feltétel: a pontos érték $1/3$, a trapezoid-közelítés N=4 rácsra $\frac{11}{32}=0{,}34375$, és a hiba $N$ növelésével másodfokban csökken.

A megoldás előtt becsüld meg a területet egyetlen téglalappal: a függvény mindenütt legfeljebb 1 magas, ezért $\int_0^1 x^2 dx$ 0 és 1 között van. Egy jobb geometriai érvelés szerint az $x^2$ görbe a $(0,0),(1,0),(1,1)$ háromszög alatt fut, ezért a terület kisebb, mint $\frac12$; a pontos $\frac13$ ehhez közeli, de a numerikus rácsnak meg kell mutatnia, mennyire.

12.1.2 Riemann-szumma, trapezoid és a terület mértékegysége

A Riemann-szumma egy véges összeg, amely a területet téglalapokkal közelíti. Osszuk az $[a,b]$ intervallumot $n$ darab, nem feltétlenül egyenlő részintervallumra:

$$a=x_0<x_1<\cdots<x_{n-1}<x_n=b.$$

A $i$-edik részintervallum szélessége

$$\Delta x_i=x_i-x_{i-1},$$

és ebből a részintervallumból választunk egy $ξ_i\in[x_{i-1},x_i]$ pontot, ahol a téglalap magasságát leolvassuk. A Riemann-szumma

$$S_n=\sum_{i=1}^{n}f(ξ_i)\Delta x_i.$$

A jelölések:

• $x_i$: a rács csomópontja; a kódban az x tömb egy eleme, float64 szám, mértékegysége az $x$ tengely egysége (itt például perc vagy dimenziótlan);
• $f(ξ_i)$: a függvény értéke a mintapontban; ez magasság, mértékegysége a függvény egysége (itt például km/perc);
• $\Delta x_i$: részintervallum szélessége, mértékegysége az $x$ tengely egysége;
• egy tag $f(ξ_i)\Delta x_i$: kis téglalap területe, mértékegysége a két egység szorzata.

Ha minden intervallum azonos $h=(b-a)/n$ szélességű, és a mintapont a bal végpont, akkor a bal oldali téglalapösszeget kapjuk; ha a jobb végpont, akkor a jobb oldalit. Az x = np.linspace(a,b,n+1) tömb alakja $(n+1,)$, mert $n$ intervallumhoz $n+1$ csomópont kell.

A trapezoid-szabály téglalapok helyett egyenes szakaszokkal köti össze a pontokat, tehát minden részintervallumon egy kis trapéz területét számítja:

$$T_n=\sum_{i=1}^{n}\frac{f(x_{i-1})+f(x_i)}{2}\Delta x_i.$$

Egyenletes rácson ez éppen a bal és jobb oldali összeg átlaga:

$$T_n=\frac{L_n+R_n}{2}.$$

A terület tehát egy mértékegységgel rendelkező mennyiség: ha $f$ sebesség km/percben és $x$ idő percben, akkor a terület km. Ha később $f$ valószínűségi sűrűség lesz, akkor a mértékegysége $1/x$ egység, ezért a szorzat dimenziótlan valószínűség — de erről majd a Sűrűség és eloszlásfüggvény alfejezetben.

A határ fogalma itt ugyanaz, mint amit a deriváltnál láttunk: rögzített szabályú közelítési folyamatot vizsgálunk. Az integrált úgy definiáljuk, hogy a legnagyobb részintervallum-szélesség nullához tart:

$$\int_a^b f(x)\,dx=\lim_{\max_i \Delta x_i\to 0}\sum_{i=1}^{n}f(ξ_i)\Delta x_i,$$

feltéve, hogy ez a határérték létezik, és nem függ a felosztás módjától vagy a $ξ_i$ mintapont választásától. A folytonos, szelíd függvényekre ez teljesül; a későbbi szakaszos vagy szinguláris függvényeknél elővigyázatosság kell.

12.1.3 Teljes minimális munkapélda: az $x^2$ területe

A kiinduló feladat pontos értéke a hatványfüggvény elemi integrálási szabályából:

$$\int_0^1 x^2\,dx=\left[\frac{x^3}{3}\right]_0^1=\frac13.$$

A Riemann-szumma finomodását az alábbi NumPy-kóddal mutathatjuk meg. A kódban a trapezoid-összeget számítjuk, és összehasonlítjuk a zárt képletből adódó $T_n=\frac{2n^2+1}{6n^2}$ értékkel. Az $x^2$ esetén ez a képlet pontosan megadja a trapezoid-közelítést egyenletes rácson:

import numpy as np

def f(x):
    return x * x

a, b = 0.0, 1.0
for n in (4, 16, 64):
    x = np.linspace(a, b, n + 1)          # shape (n+1,), float64
    h = (b - a) / n
    trap = 0.5 * np.sum((f(x[:-1]) + f(x[1:])) * h)
    expected = (2 * n * n + 1) / (6 * n * n)
    print(f"n={n:2d}  trap={trap:.12f}  expected={expected:.12f}  diff={trap - expected:.2e}")

A kód kimenete ebben az írásmódban nem lefuttatott eredményjegyzékből származik; az alábbi táblázat a fenti képletekből közvetlenül számított, elméletileg várt érték. A $n=4$ sort rögtön a táblázat után kézzel is ellenőrizzük:

Az $x^2$ görbe szemléltetése:

A kézi ellenőrzés N=4-re. A rács: $x=[0,\ 0{,}25,\ 0{,}5,\ 0{,}75,\ 1]$, $h=0{,}25$. A bal oldali téglalapok magassága: $f(0)=0,\ f(0{,}25)=0{,}0625,\ f(0{,}5)=0{,}25,\ f(0{,}75)=0{,}5625$; a szélesség mindig $0{,}25$. A bal összeg:

$$L_4=h(0+0{,}0625+0{,}25+0{,}5625)=0{,}25\cdot 0{,}875=0{,}21875.$$

A jobb oldali összeg a jobb végpontokat veszi: $f(0{,}25)+f(0{,}5)+f(0{,}75)+f(1)=0{,}0625+0{,}25+0{,}5625+1=1{,}875$, szorozva $h$-val: $R_4=0{,}46875$. A trapezoid az átlaguk:

$$T_4=\frac{0{,}21875+0{,}46875}{2}=0{,}34375=\frac{11}{32}.$$

A pontos értékhez képesti hiba tehát

$$\left|\frac{11}{32}-\frac13\right|=\frac{1}{96}\approx 0{,}0104167.$$

Ez a $T_n=\frac{2n^2+1}{6n^2}$ képletből is látszik: $T_4=\frac{33}{96}=\frac{11}{32}$, illetve $T_n=\frac13+\frac{1}{6n^2}$, tehát a hiba elméletileg $\frac{1}{6n^2}$. A fenti táblázat diff oszlopa azért 0, mert a kódban a kitevőhatvány-rácson a választott rácson a numerikus összeg és az analitikus képlet lebegőpontos értéke az adott végrehajtásban egybeeshet; általános feladathoz inkább np.isclose vagy 1e-12 tolerancia javasolt.

A gyakorlati munkához a scipy.integrate modul adaptív, nem egyenletes rácsot használó megoldókat is kínál; ezek a fajta egyszerű esetben sokkal kisebb hibával dolgoznak, mint a rögzített trapezoid. A pontos hívásformát és a lokális verzió számát a laborban a SciPy integrációs dokumentációjával érdemes ellenőrizni.

12.1.4 Eredmény értelmezése: hiba és a bal-jobb keretezés

Négy tanulságot érdemes kiemelni:

1. A bal és jobb oldali összeg keretezi az integrált egy monoton függvénynél. Mivel $x^2$ növekvő $[0,1]$-en, a bal téglalapok alábecsülik, a jobb téglalapok felülbecsülik a területet: $0{,}21875<0{,}3333<0{,}46875$. A trapezoid ezek átlaga, ezért a hiba részben kiegyenlítődik.

2. A trapezoid hibája $O(h^2)$. Az $x^2$-példán a hiba pontosan $\frac{1}{6n^2}=O(h^2)$, mert $h=1/n$. Amikor $n$-t négyszerezzük, a hiba durván tizenhatodolódik. Ez a sima függvények esetén szokásos viselkedés; a szakaszos függvényekre csak $O(h)$ marad, ahogy az egyik feladat mutatja.

3. A görbe magassága és a görbe alatti terület két külön mennyiség. Az $f(1)=1$ magasság nem azt jelenti, hogy a terület 1; a terület $1/3$, mert a magasság az intervallum nagy részén kisebb 1-nél. Az egységek is különböznek: a magasság mértékegysége az $f$ egysége, a területé az $f$ egysége szorozva az $x$ egységével.

4. A határátmenet finomabb és finomabb rácson ellenőrizhető. A $T_{64}=1/3+1/(6\cdot64^2)\approx0{,}3333740$ a valódi $1/3\approx0{,}3333333$ értéktől körülbelül $0{,}0000407$-tel tér el; négy tizedesjegyre még nem azonos a szokásos kerekítésnél; a határozott integrál ezt a folyamatot jelenti, akkor is, ha a gyakorlatban véges rácson állunk meg.

A kompozit trapézszabály másodrendű hibabecslése elegendő simasági feltétel mellett, például korlátos második derivált esetén adható. Szakadásoknál vagy nem sima integrandusnál a konvergenciarend eltérhet; a rácsfinomítás önmagában nem bizonyít helyességet.

12.1.5 Feladatok

12.1.5.1 1. feladat: mértékegység-interpretáció és lineáris sebesség

A kiszállító furgon sebessége $v(t)=0{,}5+0{,}02\,t$ km/perc, ahol $t\in[0,60]$ perc. Indokold saját szóval, miért adja $\int_0^{60}v(t)\,dt$ a megtett utat, és mi a kapott terület mértékegysége! Ezután számítsd ki a trapezoid-közelítést egyetlen részintervallummal, $n=1$-gyel, vagyis a két végpontot összekötő egyenes alatti területként. Hasonlítsd össze a pontos $66$ km értékkel! Miért kapható itt $n=1$-gyel is pontos eredmény?

12.1.5.2 2. feladat: implementálás másik függvényen

Módosítsd a példakódot úgy, hogy az $f(x)=\sin(x)$ függvényt integrálja a $[0,\pi]$ intervallumon, $n=4$, $n=8$ és $n=16$ mellett! A pontos érték $2$. Számítsd ki a trapezoid-közelítéseket, ellenőrizd, hogy a hibák közelítőleg négydelődnek-e! Hogyan változik a hibaelőjel?

12.1.5.3 3. feladat: hibadiagnózis egy elfelejtett szorzónál

Egy tanuló az $f(x)=e^{-x}$ függvényt integrálja a $[0,1]$ intervallumon. Az alábbi numerikus adatokkal rendelkezünk:

x = np.linspace(0, 1, 17)
h = (1 - 0) / 16
wrong = 0.5 * np.sum(np.exp(-x[:-1]) + np.exp(-x[1:]))  # hibás

A wrong kimenet körülbelül 10,1, a pontos terület pedig $1-e^{-1}\approx 0{,}6321$. Mi a hiba? Add meg a javított egy sort, és becsüld meg, miért ad a javított verzió 0,6321-hez közeli értéket!

12.1.5.4 4. feladat: szakaszos függvény — érvényességi határ

Legyen

$$f(x)=\begin{cases} 0, & x<0{,}3,\\ 1, & x\ge 0{,}3, \end{cases}$$

és integráljuk $[0,1]$-en, a pontos terület $0{,}7$. Számítsd ki a trapezoid-közelítést $n=10$ és $n=20$ esetén! Mit tapasztalsz a hibával kapcsolatban, ha az $n$-t megduplázod? Miért nem érvényes itt az $O(h^2)$ osztályozás?

12.1.6 Megoldások és indoklás

12.1.6.1 1. feladat megoldása

A sebesség az út idő szerinti változási rátája: ha rövid $\Delta t$ időt nézünk, a megtett út közelítőleg $v(t)\Delta t$ km. Az összegzés a kis lépéseket teszi teljessé, ezért $\int v\,dt$ a teljes út. A mértékegység: (km/perc) · perc = km.

Az $n=1$ trapezoid a két végpont közötti szakasz alatti terület:

$$T_1=\frac{v(0)+v(60)}{2}\cdot 60=\frac{0{,}5+1{,}7}{2}\cdot 60=1{,}1\cdot 60=66\,\text{km}.$$

Pontos eredményt kapunk, mert a trapezoid két pont között egyenest feltételez, és a sebesség itt lineáris: az egyenes alatti terület valóban a két végpontmagasság átlaga szorozva a szélességgel. Ez egy fontos érvényességi korlát: a trapezoid lineáris függvényen a rácstól függetlenül egzakt, görbe függvényen viszont rácsfüggő hibát hagy.

12.1.6.2 2. feladat megoldása

A kód módosítása:

def g(x):
    return np.sin(x)

a, b, exact = 0.0, np.pi, 2.0
for n in (4, 8, 16):
    x = np.linspace(a, b, n + 1)
    h = (b - a) / n
    trap = 0.5 * np.sum((g(x[:-1]) + g(x[1:])) * h)
    print(f"n={n:2d}  trap={trap:.6f}  error={trap - exact:.6f}")

Az elméletileg várt értékek:

• n=4-nél a trapéz magasságok a $0,\ \frac{\pi}{4},\ \frac{\pi}{2},\ \frac{3\pi}{4},\ \pi$ pontokon: $0$, $0{,}7071068$, $1$, $0{,}7071068$, $0$. Az összeg belső tagja $0{,}5\cdot 0+0{,}7071068+1+0{,}7071068+0{,}5\cdot 0=2{,}4142136$, szorozva $h=\pi/4$-gyel: $T_4\approx 1{,}8961$, hiba $\approx -0{,}1039$.
• n=8-nál $T_8\approx 1{,}9742$, hiba $\approx -0{,}0258$.
• n=16-nál $T_{16}\approx 1{,}9936$, hiba $\approx -0{,}0064$.

Az egymás utáni hibák aránya közelítőleg $\frac{0{,}1039}{0{,}0258}\approx 4{,}0$ és $\frac{0{,}0258}{0{,}0064}\approx 4{,}0$, tehát az $O(h^2)$ viselkedés igazolódik. A hiba negatív: a szinusz $[0,\pi]$-n konvex módon görbül lefelé, a trapéz a húr alatti területet adja, amely kisebb a görbe alatti területnél.

12.1.6.3 3. feladat megoldása

A hibás kód elfelejtette a $h$ szorzót. A wrong összeg a trapézmagasságok összege, vagyis egy hosszúság jellegű mennyiség, nem terület; mert $\sum$ tagjai $\frac{f(x_{i-1})+f(x_i)}{2}$ mértékegységűek, hiányzik a · h. Emiatt az eredmény körülbelül a 16-szorosa a területnek, hiszen $h=1/16$. A javítás:

correct = 0.5 * np.sum((np.exp(-x[:-1]) + np.exp(-x[1:])) * h)

A várt eredmény a $0{,}6321205588\ldots$ értékhez közel lesz, a maradék hiba $O(h^2)$ nagyságú. A hibaforrás azonos azzal a hibával, amikor hisztogrammagasságokat összegzünk oszlopszélesség nélkül: az összeg nem terület.

12.1.6.4 4. feladat megoldása

$n=10$ esetén $h=0{,}1$, a rács: $0, 0{,}1,\ldots,1$. A függvényértékek: $f(0)=f(0{,}1)=f(0{,}2)=0$, és $f(0{,}3)=f(0{,}4)=\cdots=f(1)=1$. A trapézösszeg:

$$T_{10}=0{,}1\left(0{,}5\cdot0+0+0+1+1+1+1+1+1+1+0{,}5\cdot1\right)=0{,}1\cdot7{,}5=0{,}75.$$

Tehát a hiba $0{,}75-0{,}7=0{,}05$.

$n=20$ esetén $h=0{,}05$, a rács $0, 0{,}05,\ldots,1$. A belső csomópontok közül az $x<0{,}3$ pontok: $0{,}05,\ldots,0{,}25$, ez 5 darab nulla; a maradék belső pontok $0{,}3$–tól $0{,}95$-ig 14 darab egyes. A trapézösszeg:

$$T_{20}=0{,}05\left(0+5\cdot0+14\cdot1+0{,}5\cdot1\right)=0{,}05\cdot14{,}5=0{,}725.$$

A hiba $0{,}025$. A rács finomításakor a hiba feleződött, nem negyedelődött: a szakadás körüli két részintervallumon a trapéz mindig a valódi ugrás felét adja, és ez a járulék $h$-val arányos. Ezért a trapezoid $O(h^2)$ konvergenciája csak kellően sima függvényekre támaszkodhat; ugrás vagy csúcs esetén a hiba $O(h)$-ra romlik.

12.1.7 Gyakori tévedés, visszautalás és önellenőrzés

Gyakori tévedés: a görbe magassága helyett a görbe alatti területet tévesen a végponti magassággal azonosítják. Ellenpélda: $f(x)=x^2$-re $f(1)=1$, tehát valaki azt mondhatná, hogy a terület $1$; a fenti számítás mutatja, hogy a terület $1/3$. A magasság pillanatnyi érték, a terület az összes magasság súlyozott összege az intervallumon.

Visszautalás Meredekség és numerikus derivált alfejezetre: a deriváltnál a differenciahányados a változási ráta volt, és a mértékegység a kimenet/bemenet hányadosa. Az integrál ennek a fordított iránya: a rátából visszanyert teljes mennyiség, amelynek mértékegysége a kimenet·bemenet szorzat. A két fogalom nem ugyanaz a szám, hanem egymás inverz műveletei.

Visszautalás CSV-től Parquet-ig: séma és partíció alfejezetre: a Parquet-fájlok ellenőrzésénél a sorösszeg fejezte ki a teljes minor unit értéket, és a duplikálódás miatt a join után a teljes összeg megnőtt. Az integrál ugyanennek az összegzési gondolatnak a folytonos megfelelője: véges sorösszeg helyett végtelenül sok, nullához tartó szélességű járulékot adunk össze.

Önellenőrzés:

1. Mi a különbség $\Delta x_i$ és $f(x_i)$ között?
2. Miért keretezi a bal és jobb oldali összeg az integrált egy növekvő függvényen?
3. Milyen mértékegysége van $\int v(t)\,dt$-nek, ha $v$ km/h és $t$ óra?
4. Mikor romlik a trapezoid $O(h^2)$ hibája $O(h)$-rá?

A válaszok a fenti szövegben ellenőrizhetők; ha bármelyik nem világos, érdemes újra elvégezni a kézi N=4 számítást.

12.1.8 Rövid összefoglaló

A Riemann-szumma téglalapokkal, a trapezoid-szabály trapézokkal közelíti a görbe alatti területet; mindkettő véges összeg, amelynek határértéke a finomabb rácson a határozott integrál. A terület mértékegysége az integrandus és a tengely mértékegységének szorzata, ezért például a sebesség-idő integrál út hosszúságot ad. Az $x^2$ területét három rácsfinomságon kiszámolva láttuk: a bal és jobb oldali összeg keretezi a pontos $1/3$ értéket, a trapezoid hibája pedig $O(h^2)$ sima függvényre, míg szakadásos függvényre csak $O(h)$ marad. A következő alfejezet ezt a területfogalmat a sűrűségre és eloszlásfüggvényre alkalmazza.