Numerično integriranje

Uvod¶

V okviru tega poglavja bomo za dano funkcijo $f(x)$ izračunali določen integral:

kjer sta $a$ in $b$ meji integriranja, $f(x)$ pa so vrednosti funkcije, ki jih pridobimo iz tabele vrednosti ali s pomočjo analitične funkcije.

Pri numeričnem integriranju integral ocenimo z $I$ in velja

kjer je $E$ napaka ocene integrala.

Numerični integral bomo računali na podlagi diskretne vsote:

kjer so $A_i$ uteži, $x_i$ pa vozlišča na intervalu $[a, b]$ in je $m+1$ število vozlišč.

Ogledali si bomo dva različna pristopa k numerični integraciji:

1. Newton-Cotesov pristop, ki temelji na ekvidistantnih vozliščih (konstanten korak integracije) in

2. Gaussov integracijski pristop, kjer so vozlišča postavljena tako, da se doseže natančnost za polinome.

Motivacijski primer¶

Pri numeričnem integriranju si bomo pomagali s konkretnim primerom:

Pripravimo si vozlišča. Osnovni korak naj bo $h=0.25$, v tem primeru imamo štiri podintervale in pet vozlišč, pri koraku $2h$ so tri vozliščne točke in pri koraku $4h$ samo dve (skrajni):

Pripravimo še funkcijske vrednosti:

Pripravimo prikaz podatkov:

Prikažimo podatke:

Analitično izračunajmo točen rezultat:

Newton-Cotesov pristop¶

Na sliki je prikazan splošen primer, kjer je razdalja med vozlišči $x_i$ enaka $h$ (gre za ekvidistantno delitev).

V okviru tega poglavja si bomo najprej pogledali trapezno ter sestavljeno trapezno pravilo, pozneje pa bosta sledili še Simpsonova ter Rombergova metoda.

Trapezno pravilo¶

Trapezno pravilo vrednosti podintegralske funkcije $f(x)$ na (pod)intervalu $[x_0, x_1]$ interpolira z linearno funkcijo. Za dve vozliščni točki to pomeni, da površino pod grafom funkcije približno izračunamo kot:

In so uteži:

Numerična implementacija¶

Numerična implementacija je:

Numerični zgled¶

V konkretnem primeru to pomeni, da prvo in zadnjo funkcijsko vrednost utežimo s $h/2$. V našem primeru je $h=1$:

Pripravimo sliko:

Prikažemo:

Napaka trapeznega pravila¶

Razlika med analitično vrednostjo integrala in numeričnim približkom $I$ je napaka metode:

Če je funkcija $f(x)$ vsaj dvakrat odvedljiva, se lahko (glejte npr. vir: Burden, Faires, Burden: Numerical Analysis 10th Ed) izpelje ocena napake, ki velja samo za trapezni približek prek celega integracijskega intervala:

kjer je $h=b-a$ in $\xi$ neznana vrednost na intervalu $[a,b]$.

Sestavljeno trapezno pravilo¶

Če razdelimo interval $[a, b]$ na $n$ podintervalov in na vsakem uporabimo trapezno pravilo integriranja, govorimo o sestavljenem trapeznem pravilu (angl. composite trapezoidal rule).

V tem primeru za vsak podinterval $i$ uporabimo trapezno pravilo in torej za meje podinterval $x_i$ in $x_{i+1}$ uporabimo uteži $A_i=A_{i+1}=h/2$. Ker so notranja vozlišča podvojena, sledi:

Pri tem smo predpostavili podintervale enake širine:

Sledi torej:

Numerična implementacija¶

Numerična implementacija je:

Numerični zgled¶

Zgoraj smo že pripravili podatke za dva podintervala (tri vozlišča):

Izračunajmo oceno integrala s sestavljenim trapeznim pravilom:

Pripravimo sliko:

Prikažemo:

Napaka sestavljenega trapeznega pravila¶

Napaka sestavljenega trapeznega pravila izhaja iz napake trapeznega pravila; pri tem tako napako naredimo $n$-krat.

Ker velja $h\cdot n=b-a$, izpeljemo napako sestavljenega trapeznega pravila kot:

$\eta$ je vrednost na intervalu $[a,b]$. Napaka je drugega reda $\mathcal{O}(h^2)$.

Boljši približek integrala¶

V nadaljevanju si bomo pogledali t. i. Richardsonovo ekstrapolacijo, pri kateri na podlagi ocene integrala s korakom $h$ in $2h$ izračunamo boljši približek.

V kolikor integral $I$ numerično izračunamo pri dveh različnih korakih $h$ in $2\,h$, velja:

kjer sta $I_h$ in $I_{2h}$ približka integrala s korakom $h$ in $2h$ ter $E_h$ in $E_{2h}$ oceni napake pri koraku $h$ in $2h$. Izpeljemo $I_{2h}-I_{h} = E_{h}-E_{2h}$

Naprej zapišemo:

Ob predpostavki, da je $f''\left (\eta \right )$ pri koraku $h$ in $2h$ enak ($\eta$ je pri koraku $h$ in $2h$ dejansko različen), zapišemo:

Sledi:

Sedaj lahko ocenimo napako pri koraku $h$:

Na podlagi ocene napake lahko izračunamo boljši numerični približek $I_{h}^*$:

ali

Numerični zgled¶

Predhodno smo s trapeznim pravilom že izračunali integral pri koraku $h=1$ in pri koraku $h=0,5$, rezultata sta bila:

S pomočjo zgornje formule izračunamo boljši približek:

Sestavljeno trapezno pravilo je implementirano tudi v paketu numpy, s funkcijo numpy.trapezoid (v starejših različicah numpy.trapz):

trapezoid(y, x=None, dx=1.0, axis=-1)

• y predstavlja tabelo funkcijskih vrednosti,

• x je opcijski parameter in definira vozlišča; če parameter ni definiran, se privzame ekvidistančna vozlišča na razdalji dx,

• dx definira (konstanten) korak integracije, ima privzeto vrednost 1,

• axis definira os po kateri se integrira (v primeru, da je y večdimenzijsko numerično polje).

Funkcija vrne izračunani integral po sestavljenem trapeznem pravilu. Več informacij lahko najdete v dokumentaciji.

numpy implementacija sestavljenega trapeznega pravila¶

Poglejmo si primer:

Simpsonova in druge metode¶

Zgoraj smo si pogledali trapezno pravilo, ki temelji na linearni interpolacijski funkciji na posameznem podintervalu. Z interpolacijo višjega reda lahko izpeljemo še druge integracijske metode.

Izračunati moramo:

Tabeliramo podintegralsko funkcijo $f(x)$ in tabelo interpoliramo z Lagrangevim interpolacijskim polinomom $P_n(x)$ stopnje $n$:

kjer so $y_i=f(x_i)$ funkcijske vrednosti v vozliščih $x_i$ in je Lagrangev polinom $l_i$ definiran kot:

Za numerični izračun integrala $\int_{a}^b f(x)\,dx$ (meje so: $a=x_0$, $b=x_n$) namesto funkcije $f(x)$ vstavimo v integral interpolacijski polinom $P_n(x)$:

Ker je integriranje linearna operacija, lahko zamenjamo integriranje in vsoto:

Lagrangev polinom integriramo in dobimo uteži $A_i$:

Izpeljava trapeznega pravila z uporabo Lagrangevih polinomov¶

Poglejmo si kako z Lagrangevim interpolacijskim polinomom prve stopnje strojno izpeljemo uteži $A_i$ v primeru trapeznega pravila.

Najprej v simbolni obliki definirajmo spremenljivko x, vozlišči x0 in x1 ter korak h:

Pripravimo Python funkcijo, ki v simbolni obliki vrne seznam $n$ koeficientov Lagrangevih polinomov $[l_0(x), l_1(x),\dots, l_{n-1}(x)]$ stopnje $n-1$:

Najprej poglejmo Lagrangeva polinoma za linearno interpolacijo ($n=2$):

Sedaj Lagrangev polinom $l_0(x)$ integriramo čez celotni interval:

Izraz poenostavimo in dobimo:

Ker je širina podintervala konstantna je $x_1=h_0+h$, izvedemo zamenjavo:

Zgornje korake za Lagrangev polinom $l_0(x)$ lahko posplošimo za seznam Lagrangevih polinomov:

Izpeljali smo uteži, ki smo jih uporabili pri trapezni metodi:

Trapezno pravilo je:

Ocena napake je (vir: Burden, Faires, Burden: Numerical Analysis 10th Ed):

Izračun uteži za Simpsonovo pravilo¶

Potem ko smo zgoraj pokazali strnjen izračun za trapezno pravilo, lahko podobno izvedemo za kvadratno interpolacijo čez tri točke ($n=3$).

Izračun uteži je analogen zgornjemu:

Simpsonovo pravilo (to pravilo se imenuje tudi Simpsonovo 1/3 pravilo) je:

Ocena napake je (vir: Burden, Faires, Burden: Numerical Analysis 10th Ed):

Pri tem je treba izpostaviti, da je napaka lokalno 5. reda $\mathcal{O}(h^5)$ in definirana z neznano vrednostjo četrtega odvoda $f^{(4)}(\xi)$, posledično je to pravilo točno za polinome stopnje 3 ali manj.

Primer uporabe:

Pripravimo sliko. Ker Simpsonovo pravilo temelji na kvadratni interpolaciji, moramo najprej pripraviti interpolacijski polinom (pomagamo si z numpy.polynomial.Polynomial):

Prikažemo:

scipy.integrate.newton_cotes¶

Koeficiente integracijskega pristopa Newton-Cotes pridobimo tudi s pomočjo scipy.integrate.newton_cotes() dokumentacija:

newton_cotes(rn, equal=0)

kjer sta parametra:

• rn, ki definira število podintervalov (mogoč je tudi nekonstanten korak, glejte dokumentacijo),

• equal, ki definira ali se zahteva konstantno široke podintervale.

Funkcija vrne terko, pri čemer prvi element predstavlja numerično polje uteži in drugi člen oceno napake.

Poglejmo si primer:

Izračun uteži za Simpsonovo 3/8 pravilo¶

Nadaljujemo lahko s kubično interpolacijo čez štiri točke ($n=4$):

Simpsonovo 3/8 pravilo je:

Ocena napake je (vir: Burden, Faires, Burden: Numerical Analysis 10th Ed):

Poglejmo si primer uporabe. Uporabimo pripravljeno tabelo vrednosti funkcije v štirih točkah:

Pripravimo še prikaz:

Prikažemo:

Sestavljeno Simpsonovo pravilo¶

Interval $[a, b]$ razdelimo na sodo število $n$ podintervalov enake širine $h=(b - a)/n$, s čimer so definirana vozlišča $x_i=a+i\,h$ za $i=0,1,\dots,n$. V tem primeru zapišemo sestavljeno Simpsonovo pravilo:

kjer je $\eta$ neznana vrednost na intervalu $[a, b]$. Napaka je četrtega reda $\mathcal{O}(h^4)$.

Numerična implementacija:

Pripravimo sliko:

Prikažemo:

Boljša ocena integrala¶

Napaka sestavljene Simpsonove metode je definirana z:

kjer je $\eta$ neznana vrednost na intervalu $[a, b]$.

Izboljšano oceno integrala določimo na podoben način kakor pri sestavljeni trapezni metodi; integral $I$ ocenjujemo pri dveh različnih korakih $h$ in $2\,h$, velja natančno:

kjer je $I_h$ približek integrala s korakom $h$ in $E_h$ ocena napake pri koraku $h$; analogno velja za $I_{2h}$ in $E_{2h}$.

Če predpostavimo, da je $f^{(4)}\left (\eta \right )$ v obeh primerih enak, lahko določimo razliko $I_{2h}-I_{h} = E_{h}-E_{2h}$.

Naprej zapišemo:

Ob predpostavki, da je $f^{(4)}\left (\eta \right )$ pri koraku $h$ in $2h$ enak ($\eta$ je pri koraku $h$ in $2h$ dejansko različen), zapišemo:

Sledi:

Sedaj lahko ocenimo napako pri koraku $h$:

Na podlagi ocene napake lahko izračunamo boljši približek $I_{h}^*$:

ali

Numerični zgled¶

Predhodno smo s Simpsonovim pravilom že izračunali integral pri koraku $h=0,5$ in pri koraku $h=0,25$, rezultata sta bila:

S pomočjo zgornje formule izračunamo boljši približek:

Pridobimo boljši numerični približek, izgubimo pa oceno napake!

Simpsonova metoda v scipy.integrate.simpson¶

V scipy je implementirano sestavljeno Simpsonovo pravilo v scipy.integrate.simpson() (dokumentacija):

simpson(y, x=None, *, dx=1.0, axis=-1)

kjer so parametri:

• y tabela funkcijskih vrednosti, ki jih integriramo,

• x vozlišča, gre za opcijski parameter, če je x=None, se predpostavi ekvidistantne podintervale širine dx,

• dx širina ekvidistantnih podintervalov oz korak integriranja,

• axis os integriranja (pomembno v primeru večdimenzijskega numeričnega polja).

V starejših različicah SciPy se je funkcija imenovala simps in je imela še parameter even; oboje je bilo odstranjeno v SciPy 1.14.

Poglejmo si primer, najprej uvozimo funkcijo simpson:

Rombergova metoda¶

Rombergova metoda temelji na Richardsonovi ekstrapolaciji. Predpostavimo, da integral $\int_a^b f(x)\textrm{d}x$ integriramo na intervalu $[a,b]$, ki ga razdelimo na $2^{n-1}$ podintervalov ($n=1,2,4,8, \dots$).

Rezultat trapeznega pravila pri $n=1$ označimo z $R_{\underbrace{1}_{n},\underbrace{1}_{j}}$, pri čemer $j$ označuje natančnost pridobljenega rezultata $\mathcal{O}(h^{2j})$.

Če uporabimo trapezno integracijsko pravilo pri $n=1,2,4 \dots$ podintervalih, izračunamo:

• $R_{1,1}$, korak $h_1=h$, red natančnosti $\mathcal{O}(h_1^2)$,

• $R_{2,1}$, korak $h_2=h/2$, red natančnosti $\mathcal{O}(h_2^2)$,

• $R_{3,1}$, korak $h_3=h/4$, red natančnosti $\mathcal{O}(h_3^2)$,

• $R_{4,1}$, korak $h_4=h/8$, red natančnosti $\mathcal{O}(h_4^2)$,

• $\dots$

• $R_{n,1}$, korak $h_4=h/(2^{n-1})$, red natančnosti $\mathcal{O}(h_n^2)$.

Na podlagi Richardsonove ekstrapolacije izračunamo boljši približek

• $R_{2,2} = R_{2,1} + \frac{1}{3}\left(R_{2,1}-R_{1,1}\right)$, korak $h_2=h/2$, red natančnosti $\mathcal{O}(h_2^4)$,

• $R_{3,2} = R_{3,1} + \frac{1}{3}\left(R_{3,1}-R_{2,1}\right)$, korak $h_3=h/4$, red natančnosti $\mathcal{O}(h_3^4)$,

• $\dots$

• $R_{n,2} = R_{n,1} + \frac{1}{3}\left(R_{n,1}-R_{n-1,1}\right)$, korak $h_n=h/(2^{n-1})$, red natančnosti $\mathcal{O}(h_n^4)$.

Nato nadaljujemo z Richardsonovo ekstrapolacijo:

• $R_{3,3} = R_{3,2} + \frac{1}{15}\left(R_{3,2}-R_{2,2}\right)$, korak $h_3=h/4$, red natančnosti $\mathcal{O}(h_3^6)$,

• $\dots$

• $R_{n,3} = R_{n,2} + \frac{1}{15}\left(R_{n,2}-R_{n-1,2}\right)$, korak $h_n=h/(2^{n-1})$, red natančnosti $\mathcal{O}(h_n^6)$.

Richardsonovo extrapolacijo lahko posplošimo:

• $R_{n,j} = R_{n,j-1} + \frac{1}{4^{j-1}-1}\left(R_{n,j-1}-R_{n-1,j-1}\right)$, korak $h_n=h/(2^{n-1})$, red natančnosti $\mathcal{O}(h_n^{2j})$

Pri tem je boljši približek $R_{2,2}$ pri koraku $h/2$ enak rezultatu, ki ga dobimo po Simpsonovi metodi pri koraku $h/2$. Podobno je boljši približek $R_{3,2}$ pri koraku $h/4$ enak numeričnemu integralu Simpsonove metode pri koraku $h/4$. In je dalje $R_{3,3}$ enak popravljenemu približku Simpsonove metode pri koraku $h/4$.

Pripravimo numerične podatke:

Poglejmo si primer od zgoraj. Najprej s sestavljeno trapezno metodo izračunamo integral pri različnih korakih (drugi red napake):

Nato izračunamo boljše približke (dobimo četrti red napake):

Rezultati predstavljajo rezultat Simpsonove metode pri koraku $h=0,5$, $h=0,25$ in $h=0.125$:

Ponovno izračunamo boljše približke (dobimo šesti red napake):

Rezultat predstavlja boljši rezultat Simpsonove pri koraku $h=0,25$ in $h=0.125$:

Ponovno izračunamo boljše približke (dobimo osmi red napake):

Rombergova metoda nam torej ponuja visoko natančnost rezultata za relativno majhno numerično ceno. Napako pa ocenimo:

Rombergova metoda v scipy.integrate.romb¶

V scipy je implementirana Rombergova metoda v scipy.integrate.romb() (dokumentacija):

romb(y, dx=1.0, axis=-1, show=False)

kjer so parametri:

• y tabela funkcijskih vrednosti, ki jih integriramo,

• dx širina ekvidistantnih podintervalov oz korak integriranja,

• axis os integriranja (pomembno v primeru večdimenzijskega numeričnega polja),

• show če je True prikaže elemente Richardsonove ekstrapolacije.

Poglejmo si primer od zgoraj:

Gaussov integracijski pristop¶

Newton-Cotesov pristop temelji na polinomu $n$-te stopnje in napaka je $n+1$ stopnje. To pomeni, da integracija daje točen rezultat, če je integrirana funkcija polinom stopnje $n$ ali manj; pri sodem $n$ zaradi simetrije celo stopnje $n+1$ (Simpsonovo pravilo, $n=2$, je npr. točno tudi za kubične polinome).

Gaussov integracijski pristop je v principu drugačen: cilj je integral funkcije $f(x)$ nadomestiti z uteženo vsoto vrednosti funkcije pri diskretnih vrednostih $f(x_i)$:

Pri tem je neznana utež $A_i$ in tudi lega vozlišča $x_i$. Za stopnje polinoma $n$ bomo potrebovali tudi več točk $(x_i, f(x_i))$.

V nadaljevanju bomo spoznali, da lahko zelo učinkovito izračunamo numerično točen integral. Prednost Gaussove integracije je tudi, da lahko izračuna integral funkcij s singularnostmi (npr: $\int_0^1\sin(x)/\sqrt{(x)}\,dx$).

Gaussova kvadratura z enim vozliščem¶

Predpostavimo, da želimo integrirati polinom stopnje $n=1$ (linearna funkcija):

Izračunajmo:

Po drugi strani pa želimo integral izračunati glede na ustrezno uteženo $A_0$ vrednost funkcije $f(x_0)$ v neznanem vozlišču $x_0$ (samo eno vozlišče!):

Z enačenjem zgornjih izrazov izpeljemo:

$a_0$ in $a_1$ sta koeficienta linearne funkcije, ki lahko zavzame poljubne vrednosti, zato velja:

Gre za sistem linearnih enačb z rešitvijo:

Če je integrirana funkcija linearna, bomo samo na podlagi vrednosti v eni točki izračunali pravo vrednost!

Da je Gaussov integracijski pristop neodvisen od mej integriranja $x_L$, $x_D$, pa uvedemo standardne meje.

Standardne meje: $x_L=-1$ in $x_D=1$¶

Zaradi splošnosti meje $x\in[x_L, x_D]$ transformiramo v $\xi\in[-1, +1]$ s pomočjo:

in

Velja:

kjer je:

V primeru standardnih mej, je pri eni Gaussovi točki utež $A_0=2$ in $x_0=0$ vrednost, pri kateri moramo izračunati funkcijo $f(x_0)$.

Strojno izpeljevanje uteži in Gaussove točke¶

Definirajmo simbole in nastavimo enačbo:

Pripravimo dve enačbi (za prvo predpostavimo $a_0=0, a_1=1$, za drugo predpostavimo $a_0=1, a_1=0$) ter rešimo sistem za A_0 in x_0:

Za dodatno razlago priporočam video posnetek.

Gaussova integracijska metoda z več vozlišči¶

Izpeljati želimo Gaussovo integracijsko metodo, ki bo upoštevala na intervalu $[a,b]$ $v$ vozlišč in bo točno izračunala integral polinomov do stopnje $n=2v-1$. Veljati mora:

Pri izpeljavi se bomo omejili na standardne meje $x_L=-1$, $x_D=1$,

kjer je polinom stopnje $n=2v-1$ definiran kot:

Z dvema Gaussovima točkama/vozliščema točno izračunamo integral polinoma do tretjega reda, s tremi Gaussovimi vozlišči pa točno izračunamo integral polinoma do petega reda!

Strojno izpeljevanje¶

Pripravimo si najprej simbolni zapis polinoma in ustreznih spremenljivk:

Sedaj pa poiščimo uteži $A_i$ in vozlišča $x_i$ za primer dveh Gaussovih vozlišč; polinom je torej tretje stopnje.

Polinom:

Podobno kakor zgoraj za eno vozlišče, tukaj definirajmo enačbe:

Rešimo jih za neznane $x_i$ in $A_i$:

Določili smo seznam dveh (enakih) rešitev.

Najprej sta definirani vozlišči: $x_0=-\sqrt{3}/3$ in $x_1=\sqrt{3}/3$, katerima pripadata uteži $A_0=A_1=1$.

Koda zgoraj je izpeljana v splošnem - število vozlišč v lahko povečate ter izračunate vozlišča ter pripadajoče uteži.

Spodaj je podana tabela vozlišč in uteži za eno, dve in tri vozlišča (za meje $a=-1$, $b=1$):

Število točk 1:

Število točk 2:

Število točk 3:

Za več vozlišč in tudi oceno napake, glejte Mathworld Legendre-Gauss Quadrature.

Za primer treh Gaussovih točk numerični integral izračunamo (standardne meje):

Numerična implementacija¶

Numerična implementacija (vključno s transformacijo mej) za eno, dve ali tri vozlišča:

Poglejmo si zgled. Najprej definirajmo funkcijo, ki jo želimo integrirati:

Sedaj pa funkcijo (v konkretnem primeru f) in ne vrednosti (npr. f(0.)) posredujemo funkciji Gaussova. Najprej za eno vozlišče, nato dve in tri:

scipy.integrate.fixed_quad¶

Znotraj scipy je Gaussova integracijska metoda implementirana v okviru funkcije scipy.integrate.fixed_quad() (dokumentacija):

fixed_quad(func, a, b, args=(), n=5)

kjer so parametri:

• func je ime funkcije, ki jo kličemo,

• a je spodnja meja,

• b je zgornja meja,

• args je terka morebitnih dodatnih argumentov funkcije func,

• n je število vozlišč Gaussove integracije, privzeto n=5.

Funkcija vrne terko z rezultatom integriranja val in vrednost None: (val, None)

Poglejmo primer od zgoraj:

Rezultat je enak predhodnemu:

scipy.integrate¶

scipy.integrate je močno orodje za numerično integriranje (glejte dokumentacijo). V nadaljevanju si bomo pogledali nekatere funkcije.

Integracijske funkcije, ki zahtevajo definicijsko funkcijo:¶

• quad(func, a, b[, args, full_output, ...]) izračuna določeni integral func(x) v mejah [a, b],

• dblquad(func, a, b, gfun, hfun[, args, ...]) izračuna določeni integral func(x,y),

• tplquad(func, a, b, gfun, hfun, qfun, rfun) izračuna določeni integral func(x,y,z),

• nquad(func, ranges[, args, opts, full_output]) izračuna določeni integral $n$ dimenzijske funkcije func( ...),

scipy.integrate.quad¶

Poglejmo si zelo uporabno funkcijo za integriranje, scipy.integrate.quad() (dokumentacija):

quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, 
     points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)

Izbrani parametri so:

• func Python funkcija, ki jo integriramo,

• a spodnja meja integriranja (lahko se uporabi np.inf za mejo v neskončnosti),

• b zgornja meja integriranja (lahko se uporabi np.inf za mejo v neskončnosti),

• full_output za prikaz vseh rezultatov, privzeto 0,

• epsabs dovoljena absolutna napaka,

• epsrel dovoljena relativna napaka.

Poglejmo primer od zgoraj:

scipy.integrate.dblquad¶

Gre za podobno funkcijo kot quad, vendar za integriranje po dveh spremenljivkah (dokumentacija). Funkcija scipy.integrate.dblquad() izračuna dvojni integral func(y, x) v mejah od x = [a, b] in y = [gfun(x), hfun(x)].

dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

Izbrani parametri so:

• func je Python funkcija, ki jo integriramo,

• a je spodnja meja integriranja x (lahko se uporabi np.inf za mejo v neskončnosti),

• b je zgornja meja integriranja x (lahko se uporabi np.inf za mejo v neskončnosti),

• gfun je Python funkcija, ki definira spodnjo mejo y v odvisnosti od x,

• hfun je Python funkcija, ki definira zgornjo mejo y v odvisnosti od x.

Poglejmo primer izračuna površine polkroga s polmerom 1:

Definirajmo ustrezne Python funkcije in izračunajmo rezultat:

Integracijske funkcije, ki zahtevajo tabelo vrednosti:¶

• trapezoid(y[, x, dx, axis]) sestavljeno trapezno pravilo,

• cumulative_trapezoid(y[, x, dx, axis, initial]) kumulativni integral podintegralske funkcije (vrne rezultat v vsakem vozlišču),

• simpson(y[, x, dx, axis]) Simpsonova metoda,

• romb(y[, dx, axis, show]) Rombergova metoda.

Tukaj si bomo na primeru ogledali funkcijo cumulative_trapezoid. Oglejmo si primer, ko na maso $m=1$ (enote izpustimo) deluje sila $F=\sin(t)$. Iz 2. Newtonovega zakona sledi: $F=m\,\ddot x$. Pospešek integriramo, da izračunamo hitrost; nato integriramo še enkrat za pomik. Definirajmo najprej funkcijo za pospešek:

Zanima nas dogajanje v času 3 sekund:

Izračunajmo tabelo pospeškov ter nato integrirajmo za hitrost (pri tem je pomembno, da definiramo začetno vrednost):

Hitrost sedaj še enkrat integrirajmo, da izračunamo pot:

Prikažimo rezultat:

Vprašanja za vaje¶

Vprašanje 1: Z uporabo orodij paketa sympy simbolno določite vrednost statičnega momenta prereza v obliki četrtine kroga na sliki:

Pripravite tudi numerično funkcijo odvisnosti integranda $f(\varphi) = \Big(\frac{r^3}{3}\, \sin (\varphi) \Big)$ pri podani vrednosti polmera $R$.

Podatki:

$R = 7.5$ mm

Vprašanje 2: Uporabite numerično funkcijo, pripravljeno pri prejšnji nalogi, in izračunajte ter grafično prikažite vrednosti integranda iz prve naloge $\Big(f(\varphi) = \frac{r^3}{3}\sin (\varphi) \Big)$ pri 100 diskretnih točkah na intervalu $\varphi \in [0, \pi/2]$.

Izračunajte vrednost določenega integrala $S$ še numerično, s pomočjo trapezne in Simpsonove 1/3 metode iz paketov numpy oziroma scipy.

Trapezno pravilo¶

Vprašanje 3: Z uporabo lastne implementacije trapeznega pravila (osnovnega, ne sestavljenega) izračunajte določeni integral

Dobljeno vrednost primerjajte z rezultatom funkcije numpy.trapezoid (ali scipy.integrate.trapezoid), kjer opazovan interval razdelite na 10 ekvidistantnih odsekov.

Izboljšan približek z Richardsonovo ekstrapolacijo:¶

Vprašanje 4 : Uporabite Richardsonovo ekstrapolacijo in določite natančnejšo vrednost integrala iz prejšnje naloge tako, da opazovan interval $x \in [0, 2]$ najprej razdelite na 2, nato še na 4 ekvidistantne segmente in uporabite sestavljeno trapezno metodo iz paketa numpy.

Simpsonovo 1/3 pravilo¶

Vprašanje 5: Določen integral iz prejšnje naloge izračunajte z lastno kodo Simpsonove 1/3 metode (osnovne, ne sestavljene).

Dobljeno vrednost primerjajte z rezultatom funkcije scipy.integrate.simpson, kjer opazovan interval razdelite na 7 ekvidistantnih točk.

Vprašanje 6: Podane so izmerjene vrednosti pospeška avtomobila v smeri vožnje v času. Z uporabo kumulativne trapezne metode (scipy.integrate.cumulative_trapezoid) določite vrednosti hitrosti in opravljene poti avtomobila v opazovanem časovnem intervalu.

Avto pri času $t = 0$ miruje (namig: argument initial).

Gaussova integracijska metoda (Gaussova kvadratura)¶

Vprašanje 7: Z uporabo sestavljene trapezne metode iz scipy in metode Gaussove integracije (scipy.integrate.quad) izračunajte vrednosti integralov:

Pri uporabi trapezne metode razdelite opazovan interval na 10 točk.

Vprašanje 8: Integral:

izračunajte še z uporabo funkcije scipy.integrate.fixed_quad, ki ji lahko podamo tudi parameter števila vozlišč pri integraciji, n. Najprej uporabite integracijo z enim vozliščem, nato z dvema in nazadnje še s tremi vozlišči.