Numeerinen integrointi

Funktio f on integroituva välillä [x₁=a, x_n=b], jos Riemannin summalla

R = \sum_i=1^n-1 f (z_i) h_i,

missä

h_i = x_i+1 - x_i, x_i <= z_i <= x_i+1,

on raja-arvo, kun jaon silmä h = max_i { h_i } -> 0.

Muuttamalla jakovälejä ja pisteiden z_i valintaa saadaan erilaisia integrointimenetelmiä.

Jaetaan integroimisväli yhtä leveisiin viipaleisiin. Olkoon kunkin leveys h.

1) Lasketaan integroitava kunkin välin alkupisteessä:

\int_x0^x_n f (x) dx = h( f (x₀) + f (x₀ + h) + ... + f (x₀ + (n - 1) h) ).

Esimerkki:

I = \int₀¹ x² d x.

Oikea arvo on 1/3. Lasketaan nyt integraalin likiarvo jakamalla integroimisväli neljään osaan.

I₄ = (1 / 4) (0 + (1 / 16) + (1 / 4) + (9 / 16)) = (14 / 64) \approx 0.219.

Koska integroitava funktio on koko välillä kasvava ja funktio lasketaan välin alkupäässä, saadaan liian pieni arvo.

2) Lasketaan arvo välin keskikohdalla

\int_x0^x_n f (x) d x = h (f (x₀ + h/2) + f (x₀ + 3 h / 2) + ... ).

I₄ = (1 / 4) ( (1 / 64) + (9 / 64) + (25 / 64) + (49 / 64) ) = (21 / 64) \approx 0.328.

\int f (x ),d x
= h [ (f₀ + f₁)/2 + (f₁ + f₂)/2 + ... (f_n-1 + f_n)/2 ]  
= (h / 2) (f₀ + 2f₁ + 2f₂ + ... + 2f_n-1 + f_n)
= (h / 2) [f (x₀ ) + 2 f (x₀ + h ) + 2 f (x₀ + 2h ) + ... 2 f (x₀ + (n-1)h ) + f (x₀ + nh) ]

Esimerkkitapaus

I₄
= (1 / 8) [ 0 + 2 (1 / 16) + 2 (4 / 16) + 2 (9 / 16) + 1 ]
= (1 / 8) (44 / 16) \approx 0.344.

Lasketaan ensin integraali yhden osavälin yli. Newtonin-Gregoryn interpolointikaavan mukaan on

f (x) \approx f₀ + s \Delta f₀ + (s (s-1) / 2) \Delta² f₀.

\int_x0 ^x₂ f (x ) d x
\approx \int_x0 ^x₂ (f₀ + s\Delta f₀ + (s(s-1)/ 2) \Delta² f₀) dx
= h \int₀2 (f₀ + s\Delta f₀ + (s(s-1)/ 2) \Delta² f₀) d s
= h (2 f₀ + 2 \Delta f₀ + (1 / 3) \Delta²f₀)
= (h / 3) (6 f₀ + 6 \Delta f₀ + \Delta² f₀ )

Sijoitetaan tähän

\Delta f₀ = f₁ - f₀
\Delta² f₀ = \Delta f₁-\Delta f₀ = (f₂ - f₁) - (f₁ - f₀)
= f₀ - 2 f₁ + f_2.

Integraali kahden jakovälin yli on

\int_x0 ^x₂ f (x ) d x \approx (h / 3) (f₀ + 4 f₁ + f₂).

Integraali koko välin yli on

\int_x0 ^x_n f (x ) d x
\approx (h / 3) (f (x₀) + 4 f (x₀+h ) + 2 f (x₀+2h ) + 4 f (x₀+3h ) + 2 f (x₀+4h ) + ... 4f (x₀+(n-1)h) + f (x₀+nh).

Jakovälejä oltava parillinen määrä.

Esimerkkitapaus:

I₄ = (1 / 12) [0 + 4 (1 / 16) + 2 (1 / 4) + 4 (9 / 16) + 1]
= 16 / (12 x 4) = 1 / 3.

Tässä integroitavaa approksimoitiin toisen asteen polynomilla. Koska integroitava itse on toista astetta, tulos on tarkka.

Virhearvio

\Delta f₀ = f₁ - f₀ \approx h (df / d x)
\Delta² f₀ = \Delta f₁ - \Delta f₀ \approx h² (d² f / dx²)
...

Lineaarinen interpolaatio: virhe luokkaa h² f" (x ).

Kun tämä integroidaan yhden osavälin yli saadaan lokaali virhe. Välin leveys on h, joten integraali on luokkaa h³ f" (z ), missä z on jokin välin piste.

Välejä on 1/h kappaletta, joten globaali virhe on luokkaa h² f" (z ).

Simpsonin menetelmässä lokaali virhe on kertalukua h³ ja globaali virhe h⁴.

Oletetaan, että funktio on korkeintaan n-asteinen polynomi. Yritetään löytää sellaiset pisteet x_i ja kertoimet w_i,että

\sum w_i f ( x_i)

antaa polynomien integraaleille täsmälleen oikean arvon.

Esimerkki: käytetään vain kahta pistettä. Valitaan vielä integroimisväli symmetriseksi, [-1, 1].

Jotta kaava toimisi oikein kaikille astetta n astetta oleville polynomeille, sen on annettava oikeat arvot myös funktioiden 1, x, x², ... , xⁿ integraaleille:

\int_-1¹ 1 dx = 2 = w₁ + w₂,
\int_-1¹ x dx = 0 = w₁x₁ + w₂x₂,
\int_-1¹ x² d x = (2 / 3) = w₁x₁² + w₂x₂²,
\int_-1¹ x³ d x = 0 = w₁x₁³ + w₂x₂³.

Tässä on neljä yhtälöä ja neljä tuntematonta, w₁, w₂, x₁ ja x₂. Toinen ja neljäs yhtälö toteutuvat, jos valitaan x₂ = -x₁ ja w₁ = w₂. Silloin ensimmäisen yhtälön perusteella on w₁ = w₂ = 1, ja kolmannesta yhtälöstä saadaan x₁ = -x₂=1/\sqrt 3.

Yleiset kolmannen asteen polynomit saadaan edellä esiintyneiden funktioiden lineaarikombinaatioina. Siten mielivaltaiselle korkeintaan kolmatta astetta olevalle polynomille p₃ pätee

\int_-1¹ p₃(x ) d x = p₃(1 / \sqrt 3) + p₃(- 1 / \sqrt 3).

Kokeillaan

\int_-1¹ (x³ + 2x² + 1) d x =
\subst_-1¹ x⁴ / 4 + 2x³ / 3 + x =
= (1 / 4) + (2 / 3) + 1 - ( (1 / 4) - (2 / 3) - 1)
= 10 / 3.

Sama Gaussin kahden pisteen kvadratuurilla:

\int_-1¹ (x³ + 2x² + 1) d x =
= 1 / (3 \sqrt 3) + (2 / 3) + 1 - ( 1 / (3 \sqrt 3) - (2 / 3) - 1)
= 10 / 3.

Korkeampaa astetta oleville polynomeille saadaan vastaavat yhtälöt, joiden ratkaiseminen on työlästä. Pisteiden x_i ratkaiseminen suoraan yhtälöryhmistä ei ole kovin käytännöllistä. Siksi polynomit esitetäänkin ensin jonkin ortogonaalin kantafunktiojoukon avulla.

Jos pisteitä on n kappaletta, niiden paikat ovat Legendren polynomin P_n (x) nollakohtia.

P₀ (x ) = 1,
P₁ (x ) = x,
P₂ (x ) = (1 / 2) ( 3x² - 1),
P₃ (x ) = (1 / 2) ( 5x³ - 3x),
...
(2n+1) x P_n (x ) = (n+1) P_n+1(x ) + n P_n-1 (x ).

Esim. kolmen pisteen menetelmän pisteet saadaan yhtälöstä 5x³-3x=0, josta x₁=-\sqrt(3/5)=-0.7746, x₂=0, x₁=\sqrt(3/5)=0.7746$.

Näitä ei kannata ratkaista joka kerta uudestaan, vaan käytetään valmiiksi taulukoituja arvoja.

Pistettä x_i vastaava paino on

w_i = [2 / (1-x_i²) (P' (x_i )) ² ].

Polynomien derivaatat saadaan kaavoista

P'₀ (x ) = 0,
P'₁ (x ) = 1,
P'₂ (x ) = 3 x,
P'₃ (x ) = (1 / 2) (15 x² - 3),
...
P'_n+1 (x ) = P'_n-1 (x ) + (2n + 1) P_n (x ).

Mielivaltainen integroimisväli voidaan muuntaa väliksi [-1, 1] sijoituksella

y = (b - a) t / 2 + (b + a) / 2,

d y = (b - a)/2 d t,

jolloin integraali on

\int_a^b f (y ) d y = (b - a)/2 \sum w_i f( y_i),

missä

y_i = ((b - a)/2) x_i + (b + a)/2.

Esimerkki:

\int₀^\pi/2 sin x d x

Lasketaan tämä Gaussin kolmen pisteen menetelmällä. Integroimisvälin muunnos on

y_i = (\pi / 4) x_i + \pi / 4.

Integraalin laskemiseen tarvitaan seuraavat suureet:

   i w_i            x_i            y_i            w_i sin y_i
  -1 0.55555 55556 -0.77459 66692  0.17703 13620  0.09783 78406
   0 0.88888 88889  0.00000 00000  0.78539 81634  0.62853 93611
   1 0.55555 55556  0.77459 66692  1.39376 49648  0.54687 26838
 \sum                                             1.27324 98855

 n    x_i                 w_i
 2 
    0.57735 02691 89626   1.00000 00000 00000

 3 
    0.00000 00000 00000   0.88888 88888 88889
    0.77459 66692 41484   0.55555 55555 55556

 4 		       
    0.33998 10435 84856   0.65214 51548 62546
    0.86113 63115 94053   0.34785 48451 37453

 5 
    0.00000 00000 00000   0.56888 88888 88889
    0.53846 93101 05684   0.47862 86704 99366    
    0.90617 98459 38664   0.23692 68850 56189

 6 
    0.23861 91860 83197   0.46791 39345 72691
    0.66120 93864 66265   0.36076 15730 48138    
    0.93246 95142 03152   0.17132 44923 79171

 7 
    0.00000 00000 00000   0.41795 91836 73469
    0.40584 51513 77398   0.38183 00505 05118   
    0.74153 11855 99395   0.27970 53914 89277  
    0.94910 79123 42759   0.12948 49661 68868

 8 
    0.18343 46424 95650   0.36268 37833 78362
    0.52553 24099 16329   0.31370 66458 77887
    0.79666 64774 13627   0.22238 10344 53375
    0.96028 98564 97536   0.10122 85362 90376

 9 		       
    0.00000 00000 00000   0.33023 93550 01260
    0.32425 34234 03809   0.31234 70770 40003    
    0.61337 14327 00590   0.26061 06964 02935    
    0.83603 11073 26636   0.18064 81606 94857    
    0.96816 02395 07626   0.08127 43883 61575

Esimerkiksi kaksiulotteinen integraali:

\int_a^b \int_c^d f( x, y) d x dy
= \sum w_i (\int_c^d f (x_i, y) d y)
= \sum w_i \sum w'_j f (x_i, y_j )

Painoista ja pisteistä riippuen menetelmä voi tässä olla mikä tahansa edellisistä.

Kun dimensio suurempi kuin noin 5, kannattaa käyttää Monte Carlo -menetelmää.

I = \int₀¹ f (x ) dx

suorakaidemenetelmällä laskettu arvo on

\sum_i=0^n-1 f (x_i ) / n.

Tämä on funktion keskiarvo välillä [0,1]. Yleisessä tapauksessa integraali on funktion keskiarvo kerrottuna välin pituudella.

Periaatteessa sama tulos saadaan, jos integraali lasketaan satunnaisissa pisteissä t_i, jotka jakautuvat tasaisesti integrointivälille:

I' = \sum_i=0^n-1 f (t_i ) / n.

Tämä on satunnaismuuttuja, jonka odotusarvo on integraalin oikea arvo. Kun N kasvaa rajatta I' lähestyy integraalin todellista arvoa.

Jos on laskettava integraali

\int_A f dV

jonkin mutkikkaan alueen A yli, lasketaan

\int_B g d V,

missä B on jokin yksinkertainen alue, joka sisältää A :n ja g (x ) = f (x ), jos x on alueen A sisällä, ja 0 muuten. Integraalin likiarvo on

I = \sum g (x_i ),

missä satunnaisluvut x_i ovat jakautuneet tasaisesti alueeseen B.

Esimerkiksi n-ulotteisen pallon tilavuuden laskeminen. Tuotetaan kuution sisälle tasaisesti jakautuneita pisteitä

r_i = (X₁, X₂, ... , X_n ).

Nyt g(r_i ) = 1, jos | r_i | on pienempi kuin pallon säde ja 0 muuten. Summa \sum g (r_i ) lähenee pallon ja ympäröivän kuution tilavuuksien suhdetta.

Tuloksen virhe on verrannollinen lausekkeeseen 1 / \sqrt N. Laskentaa voidaan jatkaa tarpeen mukaan. Jos käytetään tavallista kiinteän laskentahilan menetelmää, koko homma on uusittava, jos tarkkuutta halutaan parantaa.

Tarkkuus paranee hitaasti, joten menetelmä ei sovellu kaikkiin tehtäviin.

Virhe ei riipu avaruuden dimensiosta, joten menetelmää kannattaa käyttää, kun on laskettava moniulotteisia integraaleja.