Ohjelmointi ja numeeriset menetelmät, luento 7, 1.3.

Numeerinen integrointi

Integrointi on yhteenlaskua, joka on tasoittava operaatio: lähtötietojen satunnaiset virheet kumoavat toisensa.

Derivointiin liittyy lähes yhtäsuurten suureiden vähennyslasku, joka on virheitä vahvistava operaatio.

Jos käsitellään havaintodataa, se voi olla syytä tasoittaa tai esittää aineistoon sovitettavalla funktiolla. Ei useinkaan tarpeen integroitaessa, välttämätöntä derivoitaessa.

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.

Puolisuunnikasmenetelmä

\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.

Simpsonin menetelmä

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.

program simpson
  implicit none
  integer i
  do i=4,10
    write(*,*) i, simpsonint(0.0, 1.0, i)
  end do
contains
  
real function f(x) 
  ! integroitava funktio
  real, intent(in) :: x
  f=x**4
end function

real function simpsonint(a, b, n)
  ! integroidaan funktio f(x) valin [a, b] yli Simpsonin
  ! menetelmalla kayttamalla n jakovalia (n oltava parillinen)
  real, intent(in) :: a, b
  integer, intent(in) :: n

  real :: h, &    ! askelpituus
          s2, s4  ! osasummat
  integer i

  if (2*(n/2) /= n) then
     write(*,*) 'jakovalien maaran oltava parillinen:',n
     simpsonint=0.0
     return
  end if

  h = (b-a)/n
  s2=0.0
  s4=0.0
  do i=1,n-1,2
     s4=s4+f(a+i*h)
  end do
  do i=2,n-2,2
     s2=s2+f(a+i*h)
  end do
  simpsonint = (h/3)*(f(a)+f(b)+2*s2+4*s4)
end function

end program

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⁴.

Rombergin menetelmä

Olkoon I on integraalin tarkka arvo. Voidaan osoittaa, että puolisuunnikasmenetelmän arvo askelpituuden funktiona R₀(h) on

R₀ (h ) = I + C₂h² + C₄h⁴ + ...,

missä kertoimet C_i eivät riipu askelesta h.

R₀ (h / 2) = I + C₂ h² / 4 + C₄ h⁴ / 16 + ... ,

Lasketaan näistä lineaarikombinaatio

R₁ (h ) = (1 / 3) ( 4 R₀ (h / 2) - R₀ (h ) )
= I + C'₄h⁴ + ... ,

Alkuperäinen menetelmä on kertalukua h², mutta nistä muodostettu kombinaatio kertalukua h⁴.

Jos aineistona askelpituudella h taulukoituja arvoja, lasketaan integraalit askelilla 2h ja h.

Menetelmää voidaan jatkaa edelleen korkeampiin kertalukuihin.

Newtonin-Cotesin menetelmät

I = \sum w_i f( x_i),

missä (x_i), i = 0, ... , n on jokin sopivasti valittu pistejoukko.

Jos pisteet x_i valitaan tasavälisesti, saadaan Newtonin-Cotesin menetelminä tunnetut integrointimenetelmät.

Kaikki edellä esiintyneet menetelmät kuuluvat tähän ryhmään.

Edellä esitetyt menetelmät ovat yksinkertaisia ja helposti ohjelmoitavia. Moniin tarkoituksiin ne ovat täysin käyttökelpoisia. Suuri tarkkuus vaatii lyhyttä jakoväliä, mikä tekee menetelmistä melko hitaita.

Gaussin kvadratuuri

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

Useampiulotteiset integraalit

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 Gauss tai jokin Newton-Cotes.

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

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.

Havaintodatan integrointi

Jos aineisto ei tasavälistä, korjataan menetelmiä niin, että askelpituus muuttuu; esimerkiksi puolisuunnikassääntö:

\int f(x)dx = (1/2) ( (x₁-x₀)f(x₀) + (x₂-x₁)f(x₁) + ... + (x_n-x_n-1)f(x_n))

Toinen mahdollisuus: sovitetaan aineistoon jokin funktio, jolloin voidaan käyttää mitä tahansa menetelmää tai laskea integraali analyyttisesti.

Jos aineisto saadaan esimerkiksi simuloinnista tai on muuten raskas lakea, mutta voidaan laskea mielivaltaiselle pisteelle, kannattaa käyttää Gaussin menetelmää.

Numeerinen derivointi

f(x_s) \approx P_n(x_s) = f₀ + s \Delta f₀ + ( s\choose 2 ) \Delta²f₀ + ... + (s \choose n ) \Deltaⁿ f₀.

Tämän avulla saadaan

f'(x_s) \approx P'_n(x_s) = (d/ds) P_n(x_s) (ds/dx) = (d/ds) P_n(x_s) (1/h)
= (1/h) ( \Delta f₀ + ((2s -1)/2) \Delta²f₀ + ... )

Kun s=0, saadaan

f'(x₀) = (1/h) ( \Delta f₀ - (1/2) \Delta²f₀ + (1/3) \Delta³f₀ - ... +- (1/n) \Deltaⁿ f₀ ) $$

Tämän virhe on luokkaa

1/(n+1) hⁿ f⁽ⁿ⁺¹⁾(\xi).

Symmetrisempiä versioita:

f'(x₀) = (1/2h) (f₁ - f_-1).

f'(x₀) = (1/12h) (f_-2 - 8f_-1 + 8f₁ - f₂).

Jos kyseessä mittausdata, jossa kohinaa, se on ensin tasoitettava esimerkiksi sovittamalla siihen pienimmän neliösumman käyrä.

Satunnaisluvuista

Press, Teukolsky, Vetterling, Flannery: Numerical Recipes, Cambridge University Press.

"We guarantee that each number is random individually, but we don't guarantee that more than one of them is random."

Satunnaislukujen (oikeammin pseudosatunnaislukujen) jono on deterministinen. Samoilla alkuarvoilla saadaan aina sama lukujono. Alkuarvo voidaan yleensä valita.

Lineaarinen kongruenssimenetelmä

missä a, b ja m ovat luonnollisia lukuja.

Luvut toistuvat viimeistään m numeron kuluttua. Jos vakiot on valittu huonosti, jakso voi olla paljon lyhempi.

Peräkkäisten lukujen korrelaatio: Jos valitaan n lukua kerrallaan esittämään pisteen paikkaa n-ulotteisessa avaruudessa, pisteet eivät täytä koko avaruutta, vaan sijoittuvat n-1 -ulotteisille hypertasoille, joita on korkeintaan m^1/k kappaletta, usein paljon vähemmän.

Satunnaisluvun vähiten merkitsevät bitit ovat vähiten satunnaisia.

Satunnaislukua ei pidä paloitella osiin; osat eivät ole yhtä satunnaisia.

Fortran 90:ssä satunnaislukuja voi tuottaa aliohjelmalla

        call random_number (v)

Satunnaislukugeneraattori voidaan alustaa kutsulla

        call random_seed (n)

     subroutine zz(x,n)
     integer, intent(in) :: n
     real, dimension (n) :: x
       ...
     real, dimension(n,n) :: matrix
        ...

     real, allocatable, dimension (:,:) :: matrix
       ...
     read(*,*) n
     allocate(matrix(n,n))
       ...

     if (.not. allocated(matrix)) &
        allocate(matrix(1:10, -10:10))

     deallocate(matrix)