Ohjelmointi ja numeeriset menetelmät, luento 6, 22.2.

Omat muuttujatyypit

Muuttujatyypin äärittely:

        type kaupunki
          real :: pituus, leveys
          character (len=20) :: nimi
        end type

        type (kaupunki) :: turku, k1, k2
        type (kaupunki), dimension(1000) :: luettelo

        real :: x

        turku%pituus = 22.3
        turku%leveys = 60.45
        turku%nimi = 'Turku'

        x = luettelo(i)%pituus

        k1 = k2
        k2 = luettelo(327)

        k1 = kaupunki(25.0, 60.0, 'Helsinki')

Moduuli

        module m
        ! tyyppien ja muuttujien maarittelyt                
        ...
        contains
        subroutine sub 
        ! aliohjelma

        end
        ...
        end module

Moduulin muuttujamäärittelyt ja aliohjelmat saadaan käyttöön use-lauseella.

Moduuli käännetään erikseen ja linkitetään sitä käyttävään ohjelmaan.

Esimerkki: vakioiden määrittely

     module realnum
       implicit none
       integer, parameter :: single=selected_real_kind(5)
       integer, parameter :: double=selected_real_kind(10)
     end module

     program realtest
        use realnum
        implicit none
        real (kind=single) :: x
        real (kind=double) :: y
        ...
    end program

Use-lause

Moduulin tunnuksiin voidaan viitata myös eri nimillä:

     program realtest
        use realnum, myreal=>double
        implicit none
        real (kind=myreal) :: x
        ...
    end program

     program realtest
        use realnum, only: double
         ...
     end program

     program realtest
        use realnum, only: myreal=>double
         ...
     end program

Moduulia käyttävä ohjelma voisi olla esimerkiksi:

     program koe
        use ratpack
        implicit none
        type (rational) :: p, q, r
        p=rat(1,2)
        q=rat(2,3)
        r=p+q
        q=q+1
        write(*,*) nomin(r),denom(r)
        write(*,*) nomin(q),denom(q)
        p=5
        p= .inv. p
        write(*,*) nomin(p),denom(p)
     end program

        7 6
        5 3
        1 5

        module ratpack
        type rational
           integer :: nominator, denominator
        end type
        contains

        function mul (p, q)
        ! lasketaan rationaalilukujen p ja q tulo
          type (rational) mul
          type (rational), intent(in):: p, q
          integer a,b,c,d
          type (rational) r
          a=p%nominator; b=p%denominator
          c=q%nominator; d=q%denominator
          r%nominator = a*c
          r%denominator = b*d
          mul = simplify (r)
       end function

       function div (p, q)
       ! lasketaan rationaalilukujen p ja q osamaara
       ..
       end function

       function add (p, q)
       ! lasketaan rationaalilukujen p ja q summa
       ..
       end function

       function sub (p, q)
       ! lasketaan rationaalilukujen p ja q erotus
       ..
       end function

       function simplify (q)
       ! supistetaan rationaaliluku q yksinkertaisimpaan 
       ! muotoon Eukleideen algoritmilla
       ! funktion arvo on supistettu rationaaliluku
       type (rational) simplify
         type (rational), intent(in):: q
         integer  a, b
         a=q%nominator
         b=q%denominator
         do
           if (a < b) then
              c=a; a=b; b=c; 
           end if
           do while (a >= b) 
             a = a-b 
           end do
           if (a == 0) exit
         end do
         ! b on nyt osoittajan ja nimittajan
         ! suurin yhteinen tekija
         simplify%nominator = q%nominator/b
         simplify%denominator = q%denominator/b
      end function

      end module

     program koe
        use ratpack
        implicit none
        type (rational) :: p, q, r

        p%nominator=1;  p%denominator=2
        q%nominator=3;  q%denominator=4
        r=add(p,q)
        write(*,*) r%nominator, r%denominator
     end program

  f95 -o koe koe.f90 ratpack.o 

Moduulissa voi olla myös yksityisiä (private) määrittelyjä, jotka eivät näy käyttäjälle.

Myös julkisen tietotyypin komponentit voivat olla yksityisiä.

        real, private, dimension(100) :: table
        
        type rational
           integer, private :: nominator, denominator
        end type

Tyypin sisäistä määrittelyä voidaan muuttaa ilman, että se vaikuttaa käyttäjän ohjelmiin. Nimiä voidaan muuttaa tai koko esitystapa vaihtaa.

Käyttäjä ei voi vahingossa sotkea tietorakenteita.

Rationaalilukumoduuliin täytyy lisätä proseduuri, jolla muodostetaan rationaalilukuja:

        function rat (n, d)
        ! muodostetaan rationaaliluku, jonka 
        ! osoittaja=n nimittaja=d
          type (rational) rat
          integer, intent(in):: n, d
          rat%nominator=n
          rat%denominator=d
       end function

        function nomin (p)
        ! palautetaan rationaaliluvun osoittaja
          integer nomin
          type (rational), intent(in):: p
          nomin=p%nominator
        end function

        function denom (p)
        ! palautetaan rationaaliluvun nimittaja
          integer denom
          type (rational), intent(in):: p
          denom=p%denominator
       end function

     program koe
        use ratpack
        implicit none
        type (rational) :: p, q, r

        p=rat(1, 2)
        p=rat(3, 4)
        r=add(p,q)
        write(*,*) nomin(p), denom(p)
    end program

        type (rational) :: p, q, r
        ...
        r=p+q

Pienimmän neliösumman sovitus

R² = \sum | y_i - f(x_i) | ²

Sovitettava funktio f sisältää muuttujan lisäksi joukon vakioita (parametreja) a_k, k=1, ..., K.

Residuaali R on sovitettavan funktion parametrien funktio: R=R(a₁, ... , a_K). Etsitään sellaiset parametrien arvot, joilla R tulee mahdollisimman pieneksi.

Oletetaan, että f on derivoituva parametrien suhteen. Residuaalin R minimissä on

\d R/ \d a₁=0, ... \d R/ \d a_K=0.

Tästä saadaan yhtälöryhmä sovitettavassa funktiossa esiintyville parametreille.

Jos funktio f voidaan esittää joidenkin kantafunktioiden f_i lineaarikombinaationa

f(x)= \sum a_i f_i(x),

saadaan lineaarinen yhtälöryhmä riippumatta siitä, mitä muotoa kantafunktiot f_i ovat.

Esimerkiksi

f(x) = a₀ + a₁ x + a₂ x²,
f(x) = a₁ sin x + a₂ cos x,
f(x) = a₀ + a₁ e^x + a₂ e^3x.

Johdetaan ratkaisu, kun sovitettava funktio on muotoa

f(x) = a + b x,

eli kantafunktiot ovat 1 ja x.

Minimoitava suure on

R=\sum (y_i - a - b x_i)².

Osittaisderivaatat ovat

\dR / \da = -2\sum (y_i - a - bx_i),
\dR / \db = -2\sum (y_i - a - bx_i) x_i.

Saadaan normaaliyhtälöt:

\sum y_i - a \sum 1 - b \sum x_i = 0,
\sum x_i y_i - a \sum x_i - b \sum x_i² = 0.

eli

a N + b S_x = S_y,
a S_x + b S_xx = S_xy

missä

S_x = \sum x_i,
S_y = \sum y_i,
S_xx = \sum x_i²,
S_xy = \sum x_i y_i.

Datapisteiden koordinaatit ovat yleensä mittaustuloksia, joilla kullakin on oma virheensä.

Kirjoitetaan minimoitava suure muotoon

R² = \sum ((y_i - f(x_i) \sigma_i)²,

missä \sigma_i on havainnon y_i virhe.

Jos kaikki virheet \sigma_i ovat samoja, tämä poikkeaa aikaisemmasta versiosta vain vakiokertoimella, jolloin normaaliyhtälöiden ratkaisut eivät muutu.

Kun y:n virheet otetaan huomioon, pienimmän neliösumman suoran normaaliyhtälöt ovat:

a S + b S_x = b S_y,
a S_x + b S_xx = S_xy,

missä

S = \sum (1 / \sigma_i²),
S_x = \sum (x_i / \sigma_i²),
S_y= \sum (y_i / \sigma_i²),
S_xx = \sum (x_i² / \sigma_i²),
S_xy = \sum (x_i y_i / \sigma_i²).

Yleisessä tapauksessa normaaliyhtälöt ovat

\sum (y_i - f(x_i) / \sigma_i²) (\d f(x_i) / \d a_k=0, k=1,...,K.

Suoran tapauksessa normaaliyhtälöiden ratkaisu on

a = (S_xx S_y - S_x S_xy / D),
b = (S S_xy - S_x S_y / D),

missä

D = S S_xy - (S_x)².

Suoran parametrien hajonnat ovat

\sigma_a = \sqrt( S_xx / D),
\sigma_b = \sqrt( S / D).

Matriisiformalismi

y(x) = a₁\phi₁(\x) + ... + a_K\phi_K(x).

Selitettävän muuttujan y arvoista muodostetaan pystyvektori

y = ( y₁, y₂, ... y_n) ^T

Muuttujan x arvojen x_i avulla lasketaan matriisi

   A =  \phi1(x1) \phi2(x1)  ... \phiK(x1)
        \phi1(x2) \phi2(x2)  ... \phiK(x2)
        ...         ...               ...
        \phi1(xn) \phi2(xn)  ... \phiK(xn)

a = (a₁, a₂, ... a_K)^T.

Jotta sovitettava funktio kuvaisi dataa, täytyy olla

A a = y.

Tämä yhtälö voidaan ratkaista täsmällisesti vain, jos A on neliömatriisi, eli havaintoja on yhtä monta kuin kantafunktioita.

Yleisessä tapauksessa voidaan etsiä vektori a, joka antaa minimiarvon normille || A a - y ||.

Pienimmän neliösumman ratkaisu on

Ca = d,

missä

C = A^T A

ja

d= A^T y

Tämä K:n yhtälön ryhmä on normaaliyhtälöiden matriisimuoto.

Jos mittausten virheet otetaan huomioon

C = A^T \Sigma^-1 A

d = A^T \Sigma^-1 y,

missä \Sigma^-1 on kovarianssimatriisin käänteismatriisi.

Jos mittaukset ovat riippumattomia, kovarianssimatriisi on

  \sigma12     0               ...    0 
  0              \sigma22      ...    0 
  ...            ...             ...    0
  0              0               ...    \sigman2

  1/\sigma12     0               ...    0 
  0                1/\sigma22    ...    0 
  ...              ...             ...    0
  0                0               ...    1/\sigman2

C_ij = \sum_l (\phi_i(x_l) \phi_j(x_l) / \sigma_l²).

d_i = \sum_l (\phi_i (x_l) y_l / \sigma_l ²).

Käänteismatriisi C^-1 on kovarianssimatriisi, jonka lävistäjältä löytyvät kertoimien varianssit:

\sigma_a_i = \sqrt(C^-1_ii).

Jos C^-1 on lävistäjämatriisi, parametrit ovat riippumattomia.

 1  1  0.5
 3  1  1.0
 4  3  0.5
 6  4  0.5

  A  =    1  1 
          1  3 
          1  4 
          1  6 

  b  =    1 
          1 
          3 
          4

 
 \Sigma =       0.25  0    0     0 
                0     1    0     0 
                0     0    0.25  0 
                0     0    0     0.25

 \Sigma-1 =  4     0    0     0 
                0     1    0     0 
                0     0    4     0 
                0     0    0     4

C = A^T \Sigma^-1 A =

     13  
     47 
     47 
    221

d = A^T \Sigma^-1 y =

     33 
    151

    | 13   47 |  | a |   |  33 |
    |         |  |   | = |     |
    | 47  221 |  | b |   | 151 |

Kerroinmatriisin käänteismatriisi on

C^-1 =

    0.333,  -0.071 
   -0.071    0.020

\sigma_a = \sqrt(0.333) = 0.577,
\sigma_b = \sqrt(0.020) = 0.140.

On etsittävä matriisi X, jolle

AX = I,

missä I on yksikkömatriisi. Tämä voidaan kirjoittaa n² lineaarisen yhtälön ryhmäksi.

Jos esimerkiksi n=2, on

 | a11  a12 |  | x11  x12 |     | 1  0 |
 |                |  |                |  =  |      |
 | a21  a22 |  | x21  x22 |     | 0  1 |

  a11 x11 + a12 x21 = 1,
  a11 x12 + a12 x22 = 0,
  a21 x11 + a22 x21 = 0,
  a21 x12 + a22 x22 = 1.

Polynomit

f(x) = a₀ + a₁x + a₂x + ... + a_nxⁿ.

Polynomi voi käyttäytyä huonosti havaintopisteiden ulkopuolella.

Kun astelukua lisätään, polynomi saadaan kulkemaan kaikkien haluttujen pisteiden kautta, mutta niiden välillä se voi heilahdella voimakkaasti. Usein korkein järkevä asteluku on noin 4--5.

Jos aineistossa on pitkiä aukkoja, sovitettuun käyrään saattaa ilmaantua asiaankuulumattomia mutkia.

Regularisointimenetelmät "jäykistävät" sovitettavaa funktiota ja estävät sen liialliset heilahtelut. Esimerkiksi voidaan rajoittaa korkeimmanasteisten termien kertoimia, mikä merkitsee korkeampien derivaattojen pitämistä pieninä.

Regularisointi voidaan toteuttaa esimerkiksi korvaamalla minimoitava residuaali lausekkeella

\sum (y_i - P(x_i))² + \lambda\sum P''(x_i)²,

missä \lambda on jokin vakio. Tämän lausekkeen minimointi johtaa ratkaisuun, jossa polynomien toiset derivaatat ovat pieniä, joten käyrässä ei esiinny jyrkkiä mutkia.

Esimerkki: Kuvataan Rungen funktiota 9 pisteellä ja sovitetaan niihin eri asteisia polynomeja:

Fourier'n sarjat

f(x) = A₀ + \sum A_k cos (2\pi kx/P) + \sum B_k sin (2\pi kx/P),

missä P on jakson pituus.

Jos aineisto ei ole tasavälistä, kertoimien laskeminen on hankalaa. Helpoin tapa on ratkaista ongelma pienimmän neliösumman sovituksena. Kantafunktiona ovat sin kx ja cos kx.

Mahdollisia ongelmia:

1) Lineaarisessa sovituksessa jakso on tunnettava etukäteen. Mikäli sekin halutaan sovittaa samanaikaisesti, joudutaan ratkaisemaan epälineaarinen tehtävä.

2) Jos aineisto ei ole tasavälinen, kantafunktiot eivät ole ortogonaalisia. Myös sarjan alkupään kertoimet muuttuvat, jos sarjaan lisätään uusia termejä. Jos sarjalle käytetään aineistosta riippuvaa katkaisukriteeriä (esimerkiksi lisätään termejä, kunnes R ei enää oleellisesti pienene), eri aineistoja ei enää voi luotettavasti vertailla keskenään. Mikäli eri aineistoja halutaan verrata, on jokainen aineisto esitettävä yhtä monella termillä.

3) Aineiston tulisi kattaa ainakin yksi kokonainen jakso, eikä siinä saisi olla pitkiä katkoja. Muuten kertoimien virheet voivat tulla hyvin suuriksi.

4) Jos mitattavassa suureessa esiintyy Nyquistin taajuutta korkeampia taajuuksia, ne voivat aiheuttaa mitattuihin arvoihin matalia taajuuksia, joita todellisuudessa ei ole olemassakaan.

5) Jos aineistossa esiintyy jyrkkiä hyppäyksiä, sarja suppenee niiden lähellä varsin hitaasti.

Epälineaariset sovitukset

f(x) = a e^-bx,

Sovitetaankin funktion f logaritmi

ln f(x) = a' - bx,

missä a'=ln a. Tuloksena on lineaarinen tehtävä.

Mikäli sovitettavan funktion derivaatat parametrien suhteen pystytään laskemaan, voidaan parametreille johtaa yhtälöryhmä, kuten lineaarisessakin tapauksessa. Yhtälöryhmä on kuitenkin epälineaarinen.

Jos derivaattoja ei pystytä laskemaan analyyttisesti, tai jos tuloksena on kovin mutkikas yhtälöryhmä, on helpompaa käyttää jotakin minimointiohjelmaa.