06 - numerikus integrálás

Láttuk előadáson, hogy hogyan lehet olyan függvényt írni, mely egy másik függvényt vár paraméterben. Most azt nézzük meg, hogy hogyan tudunk ilyennel dolgozni.

függvénygörbe alatti terület

Alapvetően a határozott integrál nem más, mint területszámítás: a lenti képen pl. azt látjuk, hog ha az x^2 - 2x + 1 függvény határozott integrálját akarjuk kiszámolni a [0,3] intervallumon, akkor arra vagyunk kíváncsiak, hogy mekkora a függvény görbéje alatti, a képen narancssárgával sraffozott terület.

(Azzal, hogy ami a 0 alatt van, az mínuszban számít, ilyen most a képen nincs.)

Ilyesmire szükség lehet pl. fizikai mennyiségek számolásakor, pl. ha az x tengely az idő (fizikában sokszor az szokott lenni), a függvény maga meg valami holminak a pillanatnyi sebességét írja le (feltéve, hogy az irány állandó), akkor két időpont közti határozott integrálja ennek a függvénynek megadja a megtett utat.

numerikusan

Kalkulusból tanultunk komplikált módszereket és trükköket arra, hogy ha adott a függvényünk valami zárt képlettel (a képen pl egy polinom van, ami jónak, könnyen integrálhatónak számít), akkor szimbolikusan is ki tudjuk számolni az integrált, megkapjuk a primitív függvényt, stb.

A valóság ennyire nem szokott a kezünkre játszani: inkább az a jellemző, hogy az integrálandó f függvényre nem ismerünk zárt képletet (vagy ha igen, akkor azt nem tudjuk kiintegrálni szimbolikusan), csak mérési pontjaink vannak néhány értékkel, vagy behelyettesíteni tudunk a függvénybe, kapunk számokat, de ennél több adatunk nincs.

Ha ,,elég sűrűn'' mintavételezzük a függvényt, akkor a határozott integrál értékét egy [a,b] szakaszon (azaz ott a függvénygörbe alatti területet) tudjuk közelíteni pl. az úgynevezett trapéz módszerrel:

• választunk (vagy kapunk egy n egész számot)
• ennyi egyforma kis részre osztjuk az [a,b] intervallumot
• minden részen egy-egy trapézzal közelítjük a területet:
  • ha a kis intervallum az [x,y], akkor
  • f(x) lesz a trapéz egyik oldalának a hossza,
  • f(y) a másik,
  • az intervallum hossza meg a magassága
• és ezeket a trapéz területeket összeadjuk.

Pl. ha a fenti kép függvényében n=3 egyforma részre vágjuk az intervallumot, akkor ilyesmit kapunk:

A képen van három trapéz, mind ,,oldalra van fordítva'', a jobb oldalinak pl azt egyik oldala 4 hosszú, a másik 1, a magassága pedig 1 (ez a [2,3] intervallumon levő trapéz). A másik kettő meg egy-egy háromszög, olyan trapéz, aminek a rövid oldala 0 hosszú.

a feladat

Implementáljunk egy trapez függvényt, ami paraméterben megkapja az integrálandó függvényt, az intervallum két végpontját, a-t és b-t, meg egy n intet, ami azt mondja meg, hogy az intervallumot hány részre osszuk fel, és visszaadja a határozott integrál trapéz módszeres közelítését!

a megoldás

Az első kérdés, hogy mi legyen a függvény fejléce?

def trapez( f: Double=>Double, a: Double, b: Double, n: Int ): Double = ???

Az intervallum két végpontjának ugyan nem mondta a feladat, hogy milyen függvényeket akarunk integrálni, de ha valami mérési pontokról van szó, akkor a double jó választásnak hangzik, az osztások száma viszont egész kell legyen; az érdekes rész itt az, hogy egy double-t paraméterben váró és double-ra kiértékelődő függvényt így adunk át: `f: Double => Double`, ez az utóbbi az ilyen szignatúrájú függvény típus.

Válasz mutatása

A függvényen belül implementáljunk egy olyan függvényt, ami egy [x,y] részintervallumon számolja ki egy trapéz területét és ezt adja vissza!

def trapez( f: Double => Double, a: Double, b: Double, n: Int ): Double = {
  def trapez( x: Double, y: Double ) =  //trapéz terület: két oldal összege, szorozva  
    ( f(x) + f(y) ) * (y-x) / 2.0       //magassággal, felezve

  ??? //ez csak azért, hogy forduljon: ez a ??? kifejezés elmegy Double-nek is
}

Válasz mutatása

Itt azt érdemes megfigyelni, hogy ha kapunk egy f függvényt paraméterben, abba a megfelelő típusú értékeket simán behelyettesíthetjük és a kifejezés értéke persze olyan típusú lesz, mint f kimenetének a típusa. Így pl. most hogy f: Double => Double, és a: Double, azaz megfelel f inputjának, írhatjuk, hogy f(a).

Érdemes lehet először is kiszámolni egy kicsi intervallum hosszát. Tegyük ezt meg, de hibakezeléssel: ha 0 vagy kisebb az input n-ünk, akkor számoljunk 1-gyel! Tegyük el a kis intervallumok hosszát egy (mondjuk) h nevű értékbe.

def trapez( f: Double => Double, a: Double, b: Double, n: Int ): Double = {

  def trapez( x: Double, y: Double ) =  
    ( f(x) + f(y) ) * (y-x) / 2.0       

  val h = (b-a) / ( if ( n <= 0 ) 1 else n )  //az if is csak egy kifejezés, ér osztani az értékével

  ??? //ez csak azért, hogy forduljon: ez a ??? kifejezés elmegy Double-nek is
}

Válasz mutatása

Így is lehet, de ismerkedjünk meg a max függvénnyel, ami a Math objektumon belül van deklarálva, ezért Math.max-ként tudjuk meghívni (ahogy korábban a tesztekben Gyakorlat.fib(..)-ként hívtuk a Fibonacci függvényt): ő kap két (mondjuk) Intet és visszaadja a maximumukat.

Az a belső if ha belegondolunk, nem más, mint Math.max( n, 1 ): ha n legalább 1, akkor annyi marad, ha pedig nem, akkor 1 lesz. Tehát kicserélhetjük:

def trapez( f: Double => Double, a: Double, b: Double, n: Int ): Double = {

  def trapez( x: Double, y: Double ) =  
    ( f(x) + f(y) ) * (y-x) / 2.0       

  val h = (b-a) / Math.max(n,1)  //h a kis trapézok magassága lesz

  ??? //ez csak azért, hogy forduljon: ez a ??? kifejezés elmegy Double-nek is
}

Válasz mutatása

Most írjuk meg az ,,adjuk össze a kis trapézok területét'' függvényt. Ha nem megy elsőre, gondoljuk át: hogy írnánk meg imperatív nyelven ezt?

double sum = 0.0;
for (int i = 0; i < n; i++ ) {
  sum += trapez( a + i*h, a+(i+1)*h );
}
return sum;

Válasz mutatása

Ebből kiindulva, a for ciklus tailrecbe való átírásával milyen tailrec kódot kapunk?

def trapez( f: Double => Double, a: Double, b: Double, n: Int ): Double = {

  def trapez( x: Double, y: Double ) =
    ( f(x) + f(y) ) * (y-x) / 2.0

  val h = (b-a) / Math.max( n, 1 )

  @tailrec
  def forLoop( i: Int, sum: Double ): Double = {
    if ( i >= n ) sum
    else forLoop( i + 1, sum + trapez( a+h*i, a+h*(i+1) ) )
  }

  forLoop(0, 0.0)

}

A for ciklusból csak az `i` ciklusváltozót és a `sum` változót tartjuk meg, mert az `n` és `f` még mindig ugyanazok, mint amit argumentumként kaptunk; kezdünk a ciklus kilépési feltételével, ekkor visszaadjuk az összeget, amit kiszámoltunk, egyébként meg aktualizálunk; `i` értéke nő eggyel, `sum` értéke meg az aktuális trapézunk területével.

Válasz mutatása

hogy hívjuk meg?

Tesztelni is kéne a kódunkat, hogyan integráljuk ki pl. a rajzolt példánkat a [0,3] intervallumon, azt 3 részre osztva?

println( trapez( x => x*x - 2*x + 1, 0.0, 3.0, 3 ) )  //prints 3.5

Válasz mutatása

Amit innen látunk: ha ismert, hogy mi a függvény paraméterünk szignatúrája (most igen: a trapez függvény paraméterként egy Double => Double függvényt vár), akkor egyszerűen input => érték formában (ha a számítás összetettebb, akkor persze az érték részt kapcsosban) meg tudunk adni egy függvényt.

(egyébként ha n=100-zal hívjuk meg, akkor a visszaadott érték 3.00045, egész közel lesz a tényleges értékhez, ami 3 - ez a függvény ,,konvex'', ami azt jelenti, hogy ha szakaszokat huzigálunk két pontja közé, az nem fog lemenni a függvény görbéje alá, ezért a trapéz módszer mindig fölélövi valamennyivel egy konvex függvény integrálját. De ha elég sok a pontunk, ahol kiértékelünk, akkor elég közel tudunk menni a tényleges integrálhoz.)

Futtassuk le a teszteseteket is, hogy lássuk, valóban stimmel az előre elvárt függvények és osztás számok esetén a területszámítás.

mini upgrade

Kicsit kevesebb lehet a kerekítési hiba és pontosabb az eredmény, ha nem azt csináljuk, hogy a kis trapézok területeit osztogatjuk kettővel, és ezt adjuk össze, hanem csak a végén osztjuk kettővel az összeget. Hogy néz ki úgy a kód, ha ezt megtesszük?

def trapez( f: Double => Double, a: Double, b: Double, n: Int ): Double = {

  def trapez( x: Double, y: Double ) =
    ( f(x) + f(y) ) * (y-x) // itt nem osztunk kettővel

  val h = (b-a) / Math.max( n, 1 )

  @tailrec
  def forLoop( i: Int, sum: Double ): Double = {
    if ( i >= n ) sum
    else forLoop( i + 1, sum + trapez( a+h*i, a+h*(i+1) ) )
  }

  forLoop( 0, 0.0 ) / 2.0  //így adjuk vissza az összeg felét

}

Válasz mutatása

---

Utolsó frissítés: 2021-02-07 23:06:47