Ohjelmointi ja numeeriset menetelmät, luento 3

Ohjelmointi ja numeeriset menetelmät, luento 3, 1.2.

Funktiot, jatkoa ...

Rekursiiviset funktiot

        recursive function fibonacci(n) result(f)
        integer :: f
        integer, intent(in) :: n
        if (n<=2) then
           f=1
        else
           f=fibonacci(n-2)+fibonacci(n-1)
        endif
        end

Melko vähän käyttöä numeerisissa tehtävissä.

Sopii esimerkiksi puiden käsittelyyn:

        program tree
        integer, parameter :: n=5
        integer left(n), right(n), data(n)
        left(:)  = (/  2,  4,  5,  0,  0 /)
        right(:) = (/  3,  0,  0,  0,  0 /)
        data(:)  = (/ 10, 20, 30, 40, 50 /)
        writetree (1)

        contains
        recursive subroutine writetree(n)
          integer n
          if (left(n)  > 0) writetree(left(n))
          write (6,*) data(n)
          if (right(n)  > 0) writetree(right(n))
        end
        end program

Yhtälön ratkaiseminen

Aikaisemmin esiintyi suora iterointi, jossa yhtälö kirjoitetaan muotoon x=f(x).


program solve
! ratkaistaan yhtalo x=f(x) suoralla iteroinnilla
  implicit none
  real :: x, limit=1.0e-7
  x = iter(0.5)
  write (*,*) x, f(x)

  contains

  real function f(x)
  real, intent(in) :: x
    f = (1.0+x) ** 0.2
  end function

  real function iter (x0)
  real, intent(in) :: x0
  real :: x1, x2

  x1 = x0
  x2 = f(x1)
  ! iteroidaan, kunnes tulos ei muutu
  do while (abs(x2-x1) > limit)
     x1 = x2
     x2 = f(x1)
  write (*,*) x1,x2
  end do
  iter = x1
  end function
end program solve


>f95 -o solve solve.f90
Compiling file solve.f90.

>./solve
 1.08447182 1.15824175
 1.15824175 1.16632617
 1.16632617 1.16719866
 1.16719866 1.16729259
 1.16729259 1.16730273
 1.16730273 1.16730380
 1.16730380 1.16730392
 1.16730392 1.16730404
 1.16730404 1.16730404
 1.16730404 1.16730404

Varsinainen ratkaisija iter ei riipu funktiosta f.
Samalla ratkaisijalla voidaan ratkaista erilaisia yhtälöitä.
Tieto funktiosta välitettävä parametrina ratkaisjalle.

Parametrina välitettävän funktion oltava eri tiedostossa.

Tiedosto func.f90:


  real function f(x)
  real, intent(in) :: x
    f = (1.0+x) ** 0.2
  end function

Itse ohjelma tiedostossa solve.f90:


program solve
! ratkaistaan yhtalo x=f(x) suoralla iteroinnilla
  implicit none
  real :: x, limit=1.0e-7
  real, external :: f
  x = iter(f, 0.5)
  write (*,*) x, f(x)

  contains

  real function iter (f, x0)
  real, intent(in) :: x0
  real :: x1, x2
  real, external :: f

  x1 = x0
  x2 = f(x1)
  ! iteroidaan, kunnes tulos ei muutu
  do while (abs(x2-x1) > limit)
     x1 = x2
     x2 = f(x1)
  write (*,*) x1,x2
  end do
  iter = x1
  end function
end program solve

Käännös ja suoritus:


>f95 -c func.f90

>f95 -o solve solve.f90 func.o

>./solve
 1.08447182 1.15824175
 1.15824175 1.16632617
 1.16632617 1.16719866
 1.16719866 1.16729259
 1.16729259 1.16730273
 1.16730273 1.16730380
 1.16730380 1.16730392
 1.16730392 1.16730404
 1.16730404 1.16730404
 1.16730404 1.16730404

Ulkoisen funktion parametreja ei pystytä tarkistamaan.

Ratkaisu interface-lohko:


  real function iter (f, x0)
  real, intent(in) :: x0
  real :: x1, x2
  interface
     real function f(x)
       real, intent(in) :: x
     end function
  end interface
  ...
  end function

Muutetaan yhtälöä ja alkuarvoa:


  x = iter(f, 1.5)

  real function f(x)
  real, intent(in) :: x
    f = x ** 2
  end function

 2.25000000 5.06250000
 5.06250000 25.6289063
 25.6289063 656.840820
 656.840820 431439.875
 431439.875 1.86140361E+11
 1.86140361E+11 3.46482343E+22
A floating overflow exception was detected.

Ratkaisu hajaantuu!

Ratkaisu voi myös jäädä heilahtelemaan eri arvojen välille.

Muutetaan yhtälö muotoon x=f^-1(x) ja yritetään uudestaan:


  x = iter(f, 1.5)

  real function f(x)
  real, intent(in) :: x
    f = sqrt(x)
  end function

 1.22474492 1.10668194
 1.10668194 1.05198956
 1.05198956 1.02566540
 1.02566540 1.01275146
 1.01275146 1.00635552
 1.00635552 1.00317276
 1.00317276 1.00158513
 1.00158513 1.00079226
 1.00079226 1.00039601
 1.00039601 1.00019801
 1.00019801 1.00009894
 1.00009894 1.00004947
 1.00004947 1.00002468
 1.00002468 1.00001228
 1.00001228 1.00000608
 1.00000608 1.00000298
 1.00000298 1.00000143
 1.00000143 1.00000072
 1.00000072 1.00000036
 1.00000036 1.00000012
 1.00000012 1.00000000
 1.00000000 1.00000000

Välinpuolitusmenetelmä

f(x) = x⁵-x-1 = 0.

f(1)=-1 < 0 ja f(2)=29 > 0, joten yhtälöllä on varmasti ainakin yksi ratkaisu välillä 1 Puolitetaan tätä väliä, ja katsotaan kullakin askeleella, kummassa puoliskossa ratkaisun täytyy olla.

f(1.5)=5.09 > 0, joten polynomi vaihtaa merkkiä välillä 1 < x < 1.5.

f(1.25)=0.802 > 0, joten merkin täytyy vaihtua välillä 1 < x < 1.25.

Koska f(1.125)=-0.323 < 0, ratkaisun täytyy olla välillä 1.125 < x <1.25.

Välin puolittumisen vuoksi kullakin iteraatiokierroksella saadaan yksi bitti lisäinformaatiota. Yksi desimaaliluku vastaa suunnilleen kolmea bittiä, joten ratkaisun tarkentamiseen yhdellä merkitsevällä numerolla tarvitaan noin kolme iteraatiota.

Menetelmä suppenee varmasti, mutta hitaasti.

   program halve
   implicit none
   real x
   x = solve (1.0, 2.0)
   write(*,*) x, f(x)

   contains

   real function f(x)
   real, intent(in) :: x
      f = x**5 - x - 1
   end function f

   real function solve (xmin, xmax)
     real, intent(in) :: xmin, xmax
     real :: x1, y1, &       ! valin alkupaa
             x2, y2, &       ! valin loppupaa
             x0, y0          ! keskikohta

     x1=xmin; y1=f(x1) 
     x2=xmax; y2=f(x2)

     x0=(x1+x2)/2 ;  y0=f(x0)

     do
       if (y1 * y0 < 0) then ! ratkaisu valilla (x1,x0)
          x2=x0
          y2=f(x2)
       else                  ! ratkaisu valilla (x0,x2)
          x1=x0
          y1=f(x1)
       end if
       x0=(x1+x2)/2
       y0=f(x0)
       if (abs(y0) < 0.0001) then ! tarkkuus saavutettu
          solve=x0
          exit
       end if
     end do
   end function solve
   end program halve

Regula falsi

Yhtälön f(x) = x⁵-x-1 = 0 tapauksessa f(1) on itseisarvoltaan paljon pienempi kuin f(2). Tuntuu luonnolliselta, että ratkaisun on oltava lähempänä ykköstä kuin kakkosta.

Jos väli on lyhyt, funktiota voidaan approksimoida suoralla. Jos välin päätepisteet ovat a ja b, suoran yhtälö on

y = f(a) + [f(b)-f(a)]/[b-a] (x-a).

Tämä leikkaa x-akselin pisteessä

x = a - f(a) [b-a] / [f(b)-f(a)].

Iterointia voidaan jatkaa valitsemalla uudeksi väliksi [a,x] tai [x,b] riippuen siitä, kummalla välillä funktio vaihtaa merkkiään.

Sekanttimenetelmä

Edellisen menetelmän vikana on hitaanpuoleinen suppeneminen, kun on päästy lähelle juurta.

Jos kaksi viimeksi löydettyä approksimaatiota ovat x_n ja x_n-1, niiden määrittämän suoran yhtälö on

y = f(x_n)+ [f(x_n-1)-f(x_n)] / [x_n-1-x_n] (x - x_n).

Tästä ratkaistu uusi nollakohta on

x=x_n+1 = x_n - f(x_n) [x_n-1-x_n] / [f(x_n-1)-f(x_n)].

Tämä menetelmä suppenee usein aikaisempia nopeammin.

f(x) = 1 /x⁴-1 = 0.

Koska f(1/2) > 0 ja f(2) < 0, ratkaisun täytyy olla välillä 1/2 < x < 2. Oletetaan, että tämän välin päätepisteet ovat kaksi ensimmäistä approksimaatiota.

x₀ = 0.5000, f(x₀)=15,
x₁ = 2.0000, f(x₁)=-0.9375,
x₂ = 1.9118, f(x₂)=-0.9251,
x₃ =-8.482.

Peräkkäisten approksimaatioiden jono ei suppene ja sekanttimenetelmä ei siis löydä ratkaisua lainkaan

Välinpuolitus- ja regula falsi -menetelmä takaavat, että iteroidut arvot pysyvät koko ajan lyhenevällä välillä. Sekanttimenetelmässä mitään tällaista väliä ei ole, ja iteroitujen arvojen jono saattaa jopa ruveta hajaantumaan.

Newtonin menetelmä

Edellä esiintyi kulmakerroin

[f(x_n-1)-f(x_n)] / [x_n-1-x_n].

Kun väli lyhenee, tämä kulmakerroin lähenee funktion derivaattaa. Arvioidaan tätä lauseketta derivaatan avulla, jolloin iteraatioaskel tulee muotoon

x_n+1 = x_n - f(x_n) / f'(x_n).

Tämä iteraatio, jossa käytetään myös funktion derivaatan arvoja, on Newtonin tai Newtonin ja Raphsonin menetelmä.

Derivaatan on oltava laskettavissa. Käytännössä tämä tarkoittaa, että funktiolla on analyyttinen lauseke.

Sekanttimenetelmän ongelma koskee myös Newtonin menetelmää. Jos funktion derivaatta on pieni, menetelmä ei enää suppene. Seuraava lause antaa ehdon menetelmän suppenemiselle:

Olkoon f derivoituva välillä [a,b], f(a)f(b) < 0, f'(x) ei ole 0 ja f''(x) on samanmerkkinen koko välillä. Newtonin menetelmä suppenee silloin mielivaltaiselle alkuarvolle x₀ välillä [a,b], jos

| f(a) / f'(a) | < b-a, ja
| f(b) / f'(b) | < b-a.

Yhtälön ratkaisu Newtonin menetelmällä


program newton
   implicit none
   real :: x
   integer :: i
   real, external :: f, df   

   do i=1,10
     x = solve (f, df, 1.5, maxiter=i)
     write(*,*) i, x, f(x)
   end do

contains

   real function solve (f, df, x0, maxiter, limit)
   ! rartkaistaan yhtalo f(x)=0 Newtonin iteraatiolla
   ! df(x)=f'(x) funktion derivaatta
   ! x0 = iteraation alkuarvo
   ! maxiter = itraatioiden maksimimaara, oletus 100
   ! limit = ratkaisun tarkkuus, oletus 1E-7
     real, intent(in) :: x0  
     integer, intent(in), optional :: maxiter
     real, intent(in), optional :: limit

     interface
       real function f(x)
         real, intent(in) :: x
       end function f
       real function df(x)
         real, intent(in) :: x
       end function df
     end interface

     real :: x, y,   &       ! uusi arvo
               xx            ! edellinen arvo
     integer :: n            ! kierroslaskuri

     real :: xlimit=1E-7, ylimit=1E-7
     integer :: nlimit=100

     
     if (present(maxiter)) nlimit=maxiter
     if (present(limit)) then 
        xlimit=limit ; ylimit=limit
     end if

     n=1
     xx=x0

     do
       y=f(xx)
       x = xx - y / df(xx) 
       if (abs(y) < ylimit) exit
       if (abs(x-xx) < xlimit) exit
       n=n+1
       if (n > nlimit) exit
       xx = x 
     end do
     solve = x 
   end function solve
end program newton

Funktio ja sen derivaatta määritellään ulkoisina funktioina:

   real function f(x)
   real, intent(in) :: x
      f = cos(x) - x
   end function f

   real function df(x)
   real, intent(in) :: x
      df = -sin(x) - 1
   end function df

Käännös ja suoritus:

  >f95 -c func.f90 
  >f95 -o newton newton.f90 func.o
  >newton

 1 0.784472406 -7.67112970E-02
 2 0.739518702 -7.25686550E-04
 3 0.739085197 -1.19209290E-07
 4 0.739085138 0.00000000E+00
 5 0.739085138 0.00000000E+00
 6 0.739085138 0.00000000E+00
 7 0.739085138 0.00000000E+00
 8 0.739085138 0.00000000E+00
 9 0.739085138 0.00000000E+00
 10 0.739085138 0.00000000E+00

Yhtälöryhmät

Yhtälöitä ja tuntemattomia voi olla useampia. Menetelmiin tämä ei aiheuta oleellisia muutoksia.


program solve2
! ratkaistaan yhtalopari x=f(x,y), y=g(x,y) suoralla iteroinnilla
  implicit none
  real :: x, y, limit=1.0e-7
  call iter(1.5,2.5,x,y)
  write (*,*) x,y, f(x,y)

  contains

  real function f(x,y)
  real, intent(in) :: x,y
    f = (1+x)/y
  end function

  real function g(x,y)
  real, intent(in) :: x,y
    g = sqrt(3*x+1)
  end function

  subroutine iter (x0, y0, xs, ys)
  real, intent(in) :: x0, y0
  real, intent(out) :: xs, ys
  real :: x1, x2, y1, y2

  x1 = x0       ; y1 = y0
  x2 = f(x1,y1) ; y2 = g(x1,y1)

  ! iteroidaan, kunnes tulos ei muutu
  do while (abs(x2-x1)>limit .or. abs(y2-y1)> limit)
     x1 = x2    ;  y1 = y2 ;
     x2 = f(x1,y1) ; y2 = g(x1,y1)
     write(*,*) x1,y1,'  ', x2,y2
  end do
  xs = x1 ; ys = y1
  end subroutine
end program solve2

Ratkaisun tarkkuus

Tuntuu luonnolliselta vaatia, että jos x~ on yhtälön f(x)=0 likimääräinen ratkaisu, täytyy olla

| f(x~) | < epsilon,

missä epsilon on jokin ennalta asetettu tarkkuusvaatimus.

Jos f'(x~) on lähellä nollaa, |f(x~)| voi olla pieni laajalla alueella todellisen ratkaisun ympäristössä.

Ratkaisun tulisi olla sellainen, että myös |x_0-x~| on pieni, missä x_0 on funktion todellinen nollakohta. Koska x_0 on tuntematon, tätä ei tietenkään pystytä laskemaan. Jos funktion derivaatta tunnetaan, tätä virhettä voidaan arvioida.

Nollakohdan x_0 ympäristössä funktio on likimain

f(x) \approx f'(x_0)(x-x_0),

josta

x-x_0 \approx f(x)/f'(x_0).

Kun x=x~, on

x~ - x_0 \approx f(x~)/f'(x_0).

Jos funktion derivaatta ei muutu kovin rajusti nollakohdan ympäristössä, on likimain

x~ - x_0 \approx f(x~)/f'(x~).

Ratkaisun likiarvo on siten lähellä todellista nollakohtaa, jos

| f(x~)/f'(x~)| < epsilon.

Tämän lisäksi vaaditaan tietenkin, että myös itse funktion arvo f(x~) on pieni.