write (*,*) x
write (6,*) x
write (*,100) x
100 format('x=',F8.3)
write(*,'("x=",F8.3)') x
write(*,"('x=',F8.3)") x
write(*,'(''x='',F8.3)') x
character (len=80) :: form
character (len=2) :: fi
...
fi='I2'; if (i>99) fi='I3'
form='('//fi//',F5.2)'
write(6,form) i,x(i)
Muotoinosa tulkitaan vasta suoritushtkellä.
Syötön ja tulostuksen ohjaus Unixissa:
read(*,...), read(5,...) lukee tiedostosta input:
prog < inputwrite(*,...), write(6,..) tulostaa tiedostoon output:
prog > output prog < input > outputJos output olemassa, edellinen versio tuhoutuu. Useamman ajon tulostus peräkkäin samaan tiedostoon:
prog > output prog >> outputSyöttötietojen lukeminen komentotiedostosta (here document)
prog < < EOF ... ... EOF
Syöttö- ja tulostuslauseessa voi olla kaksi osoitetta virhetilanteita varten:
read(5,1, ERR=100, END=900) x
1 format(f5.2)
...
100 write(6,"('pieleen meni')")
...
900 write(6,"('tiedosto loppui')")
F90:ssä myös:
integer status
do
read(5,form, iostat=status) x
if (status > 0) then
write(6,"('pieleen meni')")
exit
else if (status < 0) then
write(6,"('tiedosto loppui')")
exit
end if
...
end do
integer status
do
read(5,form, iostat=status) x
if (status /= 0) then
call error(status)
exit
end if
...
end do
...
subroutine error(code)
integer code
write(6,"('error',I4)") code
end subroutine
Laitenumeroihin 5 ja 6 liittyviä tiedostoja ei tarvitse avata tai sulkea.
Muihin laitenumeroihin voi liittyä jokin oletusnimi (fort0001.dat tms.)
Tiedosto avataan open-lauseella:
open (1, 'fyle')
! luodaan uusi tiedosto
! tiedostoa ei saa olla ennestaan olemassa
open (1, file='fyle', status='new')
! avataan olemassaoleva tiedosto
open (1, file='fyle', status='old', ERR=900)
! jos tiedosto olemassa ja siihen kirjoitetaan,
! vanha tiedosto tuhotaan ja luodaan uusi
open (1, file='fyle', status='unknown')
Tiedosto suljetaan close-lauseella
close(1)
Binääritulostus:
open(1, file='bin.dat', form='unformatted')
write (1) x
Normaalisti tiedostoa luetaan ja kirjoitetaan
tietue kerrallaan alusta alkaen.
Suorasaantitiedoston määrittely
open(1, file='lista', access='direct', recl=128)
read(1, rec=i) x
Suorasaantitiedosto oletusarvoisesti binääritiedosto.
Suorasaantitiedostolle on aina ilmoitettava tietueen pituus recl tavuina. Luku- ja kirjoituslauseessa ilmoitettava tietueen numero (rec)
Olkoon g tunnettu funktio ja f sen approksimaatio.
Funktioiden välistä "etäisyyttä" eli normia || f-g || pyritään minimoimaan. Normi voidaan määritellä monella tavalla (jos pisteitä äärellinen määrä, integraali korvataan summalla):
L1-normi
|| f-g ||1 = \int | f-g | dx
Sallii suuret poikkeamat, jos ne rajoittuvat lyhyelle välille.
L2-normi vastaa euklidisen avaruuden etäisyyttä
|| f-g ||2 = \sqrt( \int | f-g |2 dx).
Yleisemmin Lp-normi
|| f-g ||p = ( \int | f-g |p dx )1/p.
L\infty-normi eli maksiminormi
|| f-g || \infty = max | f-g |.
Estää suuret poikkeamat, mutta sovitus ei välttämättä missään kovin hyvä.
y = f(x0) + f'(x0) (x-x0),
missä f'(x0) on funktion f derivaatan df/dx arvo pisteessä x0. Jos nyt x on lähellä x0:aa, ovat itse funktion ja sen tangentin kuvaajat lähellä toisiaan pisteessä x. Funktion arvoa tässä pisteessä voidaan niin ollen arvioida tangentin avulla:
f(x) = f(x0) + f'(x0) (x - x0).
Tämä arvio on sitä huonompi, mitä enemmän derivaatta f' muuttuu välillä [x0, x]. Tätä muutosta puolestaan kuvaa toinen derivaatta f'' jne. Mitä tarkempi tulos halutaan sitä korkeampia derivaattoja on otettava mukaan. Voidaan osoittaa, että funktion arvo pisteessä x on
f(x) =
f(x0) +
f'(x0) (x-x0) +
(1/2) f''(x0) (x-x0)2 +
... +
(1/n!) f(n) (x0)
(x-x0)n + ...
missä f(n)(x0) on funktion n:nnen derivaatan arvo pisteessä x0 ja n! on n-kertoma = 1×2×3×...×n. Tätä sarjaa sanotaan funktion f Taylorin sarjaksi pisteessä x0.
Seuraavassa eräitä useimmin käytettyjä sarjoja. Kaikissa näissä on x0 = 0; tällaista sarjaa kutsutaan joskus myös MacLaurinin sarjaksi.
| 1/(1 + x) | = 1 - x + x2 - x3 + ... | suppenee, kun |x| < 1 |
| 1/(1 - x) | = 1 + x + x2 + x3 + ... | |
| (1 + x) | = 1 + x/2 - x2/8 + x3/16 - ... | |
| (1 - x)1/2 | = 1 - x/2 - x2/8 - x3/16 - ... | |
| 1/(1 + x)1/2 | = 1 - x/2 + 3x2/8 - x3/16 - ... | |
| 1/(1 - x)1/2 | = 1 + x/2 + 3x2/8 + 5x3/16 + ... | |
| ex | = 1 + x + x2/2! + x3/3! + ... + xn/n! + ... | kaikilla x |
| ln(1 + x) | = x- x2/2 + x3/3 - x4/4 +... | -1 < x <= 1 |
| sin x | = x - x3/3! + x5/5! - ... | kaikilla x |
| cos x | = 1 - x2/2! + x4/4! - ... | kaikilla x |
| tan x | = x + x3/3 + 2x5/15 + ... | |x| < pi/2 |
Monissa käytännön tapauksissa nämä sarjat suppenevat varsin nopeasti, jolloin riittää ottaa mukaan vain muutama ensimmäinen termi. Suurimpana hyötynä tästä on mahdollisuus korvata monimutkaisissa lausekkeissa esiintyvät hankalat funktiot kuten neliöjuuret yksinkertaisilla polynomeilla, jolloin lausekkeita pystytään usein sieventämään huomattavasti. Tavallisimmin käytettyjä ovatkin esimerkiksi seuraavat lineaariapproksimaatiot:
(1 + x)1/2 = 1 + x/2,
1 / (1 + x)1/2 = 1 - x,
jne.
(a0 + a1x + ... + anxn ) / (1 + b1x + ... + bmxm )
Kertoimien etsiminen johtaa yleensä optimointitehtävään.
Usein käytetty yksinkertainen menetelmä on Pade-approksimaatio.
Esimerkiksi eksponenttifunktion Taylorin sarja pisteessä x=0 on
f(x) = ex = 1 + x + (1/2)x2 + (1/6)x3 + ...
Yritetään muodostaa rationaaliapproksimaatio, joka on muotoa
R(x) = (a + bx + cx2) / (1 + dx).
Vaadimme, että tämä antaa pisteessä x=0 saman arvon kuin Taylorin sarja, ja että myös molempien derivaatat ovat samoja.
Tarkastellaan erotusta
f(x)-R(x)
= 1 + x + (1/2)x2 + (1/6)x3
- (a + bx + cx2) / (1 + dx)
= [(1 + x + (1/2)x2
+ (1/6)x3)(1 + dx)
- (a + bx + cx2) ] / [1 + dx].
Vaadimme, että origossa tämä erotus ja sen derivaatat ovat nollia. Osoittajan on oltava identtisesti nolla.
Kirjoitetaan osoittaja auki ja ryhmitellään termit muuttujan x potenssien mukaan:
[1-a] + [1 + d - b]x + [(1/2) + d - c]x2 + [(1/6) + (1/2)d]x3 +(1/6)x4.
Jotta tämä olisi identtisesti nolla, on x:n jokaisen potenssin kertoimen erikseen oltava nolla. Viimeistä termiä ei saada häviämään, mutta kyseessä on kolmannen asteen approksimaatiota, joten termi voidaan heittää pois. Saadaan yhälöryhmä
1 - a = 0,
1 + d - b = 0,
(1/2) + d - c = 0,
(1/6) + (1/2)d = 0,
josta
a = 1,
b = 2/3,
c = 1/6,
d = -1/3.
Approksimaatio on
f(x) = [1 + (2/3)x + (1/6)x2] / [1 - (1/3)x].
Esimerkiksi pisteessä x=1 neljän termin Taylorin sarja antaisi arvoksi 8/3 \approx 2.667 (suht. virhe 0.019). Rationaaliapproksimaatiosta saadaan f(1)=11/4=2.750 (suht. virhe 0.012).
Sama vika kuin Taylorin sarjalla: se antaa yhdessä pisteessä funktion arvon täsmälleen oikein, mutta ei pyri minimoimaan virhettä muissa pisteissä.
i x f(x) \Delta f \Delta^2 f \Delta^3 f
-1 x_{-1} y_{i-1}
\Delta f_{-1}
0 x_0 y_0 \Delta^2 f_{-1}
\Delta f_0 \Delta^3 f_{-1}
1 x_1 y_1 \Delta^2 f_0
\Delta f_1
2 x_2 y_2
Merkitään h = x1 - x0.
Askeloperaattori Ef(x) = f(x + h).
Eteenpäinen differenssi \Delta f(x) = f(x + h) - f(x).
Taaksepäinen differenssi \nabla f(x) = f(x) - f(x - h).
Kaikki ovat lineaarisia operaattoreita.
Operaattoreille voidaan johtaa mm. seuraavat ominaisuudet:
E2 f(x) = E(Ef(x)) = E(f(x + h)) = f(x+2h).
En f(x) = f(x + nh).
\Delta f(x) = f(x + h) - f(x) = Ef(x) - f(x) = (E-1) f(x),
josta
\Delta = E-1.
Lineaarinen interpolointi:
P1 (x0 + sh)
= y0 +
(x - x0)
[y1 - y0] /
[x1 - x0]
= y0 +
[(x - x0) / h]
(y1 - y0)
= y0 +
s [y1 - y0]
= f0 + s\Delta f0.
Newtonin-Gregoryn interpolointipolynomi:
P(x0+sh)
= Es f(x0)
= (1+\Delta)s f(x0)
= [ 1 + s\Delta + {s \choose 2 }\Delta2
+ ... ] f(x0)
= f0 + s\Delta f0
+ {s \choose 2} \Delta2 f0 + ...
Esimerkiksi toisen asteen interpolointipolynomi:
P2(x0 + sh) = f0 + s\Delta f0 + (s(s-1)/2) \Delta2 f0.
Kun s = 0, 1, 2, tämä kulkee pisteiden (x0, y0), (x0 + h, y1) ja (x0+2h, y2) kautta:
P2(x0)
= P2(x0+0h)
= y0,
P2(x1)
= P2(x0+1h)
= y0 + \Delta f0
= y1,
P2(x2)
= P2(x0+2h)
= y0 + 2\Delta f0
+ \Delta2 f0
= y2
P2 voidaan käsittää myös s:n funktioksi, kun s = (x-x0)/h on mielivaltainen reaaliluku.
Samalla tavoin voidaan johtaa myös korkeamman asteisia interpolointipolynomeja.
Eteenpäisten differenssien sijasta voidaan käyttää myös taaksepäisiä tai molempia yhdessä.
x y
10 2.0
30 3.0
50 3.8
75 4.8
100 5.2
Aineistossa 5 pistettä, joten se voidaan esittää
neljännen asteen polynomin avulla:
y = P4(n) = a0 + a1x + a2x2 + a3x3 + a4x4.
Sijoittamalla annetut arvot saadaan yhtälöryhmä
a0 + 10 a1
+ 100 a2 + 1000 a3
+ 10000 a4 = 2.0,
a0 + 30 a1
+ 900 a2 + 27000 a3
+ 810000 a4 = 3.0,
a0 + 50 a1
+ 2500 a2 + 125000 a3
+ 6 250000 a4 = 3.8,
a0 + 75 a1
+ 5625 a2 + 421875 a3
+ 31 640625 a4 = 4.8,
a0 + 100 a1
+ 10000 a2 + 1000000 a3
+ 100000000 a4 = 5.2.
Tämä on lineaarinen yhtälöryhmä, jonka ratkaisu on
a0 =1.234,
a1 =0.09115,
a2 =-0.001672,
a3 =0.00002347,
a4 =-0.0000001189.
Työlästä!
Jos esitettäviä pisteitä on n kappaletta, kardinaalifunktioita tarvitaan n kappaletta, Li, i = 1, ... ,n. Valitaan ne siten, että Li(xi) = 1 ja Li(xj)=0, kun i ja j ovat eri suuria.
Edellisessä esimerkissä n=4, joten kardinaalifunktiot ovat neljännen asteen polynomeja. Algebran peruslauseen mukaan ne voidaan kirjoittaa neljän tekijän tulona:
L1(x) =
A1 (x-x2)
(x-x3) (x-x4)
(x-x5),
...
L5(x) =
A5 (x-x1)
(x-x2) (x-x3)
(x-x4).
Vakiot Ai saadaan ehdosta Li(xi)=1; esim.
A (x1-x2) (x1-x3) (x1-x4) (x1-x5)=1.
Ensimmäinen kardinaalifunktio on siten
L1(x) = [(x - x2) / (x1 - x2)] [(x - x3) / (x1 - x3)] [(x - x4) / (x1 - x4)] [(x - x5) / (x1 - x5)]
Esimerkkitapauksen kardinaalifunktiot
L1(x)
= (1/4680000)(x-30)(x-50)(x-75)(x-100),
L2(x)
= (1/-800000)(x-10)(x-50)(x-75)(x-100),
L3(x)
= (1/700000)(x-10)(x-30)(x-75)(x-100),
L4(x)
= (1/-1421875)(x-10)(x-30)(x-50)(x-100),
L5(x)
= (1/7875000)(x-10)(x-30)(x-50)(x-75).
Interpolointipolynomi voidaan muodostaa kardinaalifunktioiden lineaarikombinaationa:
P(x) = y1 L1(x) + y2 L2(x) + ... + ynn Ln (x).
Esimerkissä
P(x) = 2 L1(x) + 3 L2(x) + 3.8 L3(x) + 4.8 L4(x) + 5.2 L(x).
Polynomeja voi käyttää interpolointiin, mutta ei ekstrapolointiin.
x y
0.0 0
1.2 6
2.0 11
3.5 9
4.1 17
5.0 24
Aineistoa ei voi kunnolla kuvata
yhdellä funktiolla.
Käytetään paloittaista sovitusta: aineisto jaetaan sopiviin väleihin, ja kuhunkin väliin sovitetaan eri funktio.
Spline-funktiot ovat kolmannen asteen polynomeja, joiden kertoimet valitaan siten, että osavälien rajoilla toiset derivaatat ovat jatkuvia.
Kertoimille saadaan lineaarinen yhtälöryhmä.
Esimerkin aineistossa on 6 pistettä 5 väliä. Kuvataan käyrää viidellä kolmannen asteen polynomilla
Si(x)=ai+bi ((x-xi) / hi)
+ci ((x-xi) / hi )2
+di ((x-xi) / hi )3,
i=1, ... ,5,
missä
hi = x_i+1 - xi.
Ensimmäinen ja toinen derivaatta:
S'i(x) = 1/hi
(bi + 2ci ((x-xi) / hi)
+ 3di ((x-xi / hi)2 ),
S''i(x) = 1/(hi)2)
(2ci + 6di ((x-xi) / hi).
Ensimmäistä väliä esittävän polynomin täytyy kulkea molempien päätepisteiden kautta:
S1(0) = a1 = 0,
S1(1.2) =
a1 + b1 + c1 + d1 = 6.
Ensimmäisen välin päätepisteessä polynomin S1 ensimmäinen ja toinen derivaatta ovat samoja kuin polynomin S2 vastaavat derivaatat toisen välin alussa:
S'1(1.2) = S'2(1.2),
S''1(1.2) = S''2(1.2),
eli
(1 / 1.2) ( b1 + 2c1 + 3d1 ) = (1 / 0.8) b2,
(1 / 1.22) (2c1 + 6d1) = (1 / 0.82) 2c2.
Esimerkissä välejä 5 kappaletta, ja kutakin väliä kohti saadaan neljä yhtälöä, joten periaatteessa saadaan 20 yhtälöä. Derivaattoja koskeva yhtälö ei päde viimeiselle välille: Yhtälöitä on 18, mutta määrättäviä vakioita 20.
Kaksi vakiota voidaan valita halutulla tavalla.
Luonnollinen splini: toiset derivaatat nollia ensimmäisessä ja viimeisessä pisteessä.
S''1(0) = 2c1 = 0,
S''5(5) = 2c5 + 6d5 = 0.
Saadaan lineaarinen yhtälöryhmä
a1 = 0,
a1 + b1 + c1 + d1 = 6,
b1 + 2c1 + 3d1 - 1.500b2 = 0,
2c1 + 6d1 - 4.5002c2= 0,
2c1 = 0,
a2 = 6,
a2 + b2 + c2 + d2 = 11,
b2 + 2c2 + 3d2 - 0.533b3 = 0,
2c2 + 6d2 - 0.5692c3= 0,
.
.
.
a5 = 17,
a5 + b5 + c5 + d5 = 24,
2c5 + 6d5 = 0,
Kertoimet ai saadaan suoraan ja niiden arvot voidaan sijoittaa yhtälöihin, joissa ne esiintyvät.
Tämän ratkaisu on
i ai bi ci di 1 0.0 4.54 0.00 1.46 2 6.0 5.94 1.94 -2.89 3 11.0 2.19 -23.68 19.50 4 9.0 5.32 5.57 -2.89 5 17.0 11.67 -7.00 2.33
\pspic{10cm}{pict/splin.ps}
Edellä saatu yhtälöryhmä ei ole kätevin mahdollinen.
y = ai + bi (x-xi) + ci (x-xi)2 + di(x-xi)3,
yi = ai,
yi+1 = ai + bihi + cihi2 + dihi3,
y' = bi + 2ci(x-xi) + 3di(x-xi)2,
y'i= bi,
y'i+1 = bi + 2cihi +3di hi2,
y'' = 2ci +6di(x-xi),
y''i = 2ci,
y''i+1 = 2ci + 6dihi
Otetaan uusiksi muuttujiksi toisen derivaatan arvot Di=y''i. Tuntemattomille vakioille saadaan lausekkeet
ai = yi,
ci = Di /2,
di = (Di+1 - Di)/6hi,
bi = (yi+1-yi) / hi - (2hiDi + hiDi+1) / 6.
Välin i alussa y'i=bi. Edellisen välin avulla laskettuna derivaatta on
y'i = bi-1 + 2ci-1(xi-xi-1) + 3di-1(xi-xi-1)2
= bi-1 + 2ci-1 hi-1 + 3di-1hi-12.
Asetetaan nämä lausekkeet yhtä suuriksi ja lausutaan vakiot derivaattojen Di ja y-arvojen avulla:
yi'= (yi+1-yi / hi) - (2hiDi + hiDi+1 / 6)
= 3 ( {Di-Di+1 / 6hi-1 )hi-12
+2 (( Di-1) / 2 ) hi-1
+ ( yi-yi-1) / hi-1
- ( 2hi-1Di-1 + hi-1Di) / 6,
joka sievenee muotoon
hi-1Di-1 + (2hi-1+2hi)Di + hiDi+1
= 6 [ ( yi+1-yi / hi)
- (yi - yi-1 / hi-1) ],
i=2, ... ,n-1.
Tässä on n-2 yhtälöä ja n tuntematonta Di. Lisäksi voidaan valita esimerkiksi D1 = Dn=0, jolloin kerroinmatriisi on
2(h1+h2) h2 0 0 ... 0 0 h2 2(h2+h3) h3 0 ... 0 0 0 h3 2(h3+h_4) ... 0 0 . . . 0 0 0 ... hn-2 2(hn-2+hn-1)
Tämä on tridiagonaalinen yhtälöryhmä, joka on helppo ratkaista.
subroutine cubicspline(n, x, y, a, b, c, d)
integer, intent(in) :: n
real, intent(in), dimension(maxpoint) :: x, y
real, intent(out), dimension(maxpoint) :: a, b, c, d
integer i
real, dimension(maxpoint,4) :: u
real, dimension(maxpoint) :: s, h
real t
do i=1,n-1 ! h = askel x-suunnassa
h(i) = x(i+1)-x(i)
end do
do i=1,n-2 ! lasketaan kerroinmatriisi
u(i, 1) = h(i)
u(i, 2) = 2*(h(i)+h(i+1))
u(i, 3) = h(i+1) ! oikean puolen vektori
u(i, 4) = 6.0* ((y(i+2)-y(i+1))/h(i+1) &
- (y(i+1)-y(i))/h(i))
end do
u(1,1) = 0.0
u(n-2,3) = 0.0
do i=2,n-2 ! eliminointi
t = u(i, 1)/u(i-1, 2)
u(i, 2) = u(i,2) - t * u(i-1, 3)
u(i, 4) = u(i,4) - t * u(i-1, 4)
end do
u(n-2,4) = u(n-2,4) / u(n-2,2) ! takaisinsijoitus
do i=n-3,1,-1
u(i,4) = (u(i, 4) - u(i, 3) * u(i+1, 4))/u(i,2)
end do
s(1) = 0.0 ! toiset derivaatat
do i=2,n-1
s(i)=u(i-1,4)
end do
s(n) = 0.0
do i=1,n-1 ! splinien kertoimet
a(i) = y(i)
b(i) = (y(i+1)-y(i))/h(i) - (2.0*h(i)*s(i) + h(i)*s(i+1))/6.0
c(i) = s(i) / 2.0
d(i) = (s(i+1)-s(i))/(6*h(i))
end do
end subroutine
Splinejä voidaan jäykistää, jolloin ne mutkittelevat vähemmän, mutta eivät enää kuvaa aineistoa aivan täsmällisesti.
Splinit sopivat interpolointiin, mutta eivät ekstrapolointiin.
Jos aineistossa on jyrkkiä vaihteluita tai pitkiä tyhjiä välejä, käyrään voi tulla ylimääräisiä mutkia ja heilahteluja.
Splinien ongelmia voi korjata lisäämällä uusia havaintopisteitä. Tärkeätä on, että pisteitä on riittävän tiheässä jyrkkien mutkien ympäristössä.
Jos aineiston yhtä pistettä muutetaan, koko ratkaisu on laskettava uudestaan. Splinien kerroinmatriisi on nauhamatriisi, ja yhden pisteen muutoksesta aiheutuva häiriö vaimenee yleensä nopeasti kaueammas siirryttäessä.
Esimerkki: Rungen funktio
y = 1 / (1 + 25x2).
Jos aineisto on yleistä käyrää esittävä pistejoukko, x-koordinaattia ei voi käyttää riippumattomana muuttujana.
Jos molemmat koordinaatit riippuvat jostakin parametrista t, se voidaan valita riippumattomaksi muuttujaksi. Saadaan kaksi funktiota x=x(t) ja y=y(t), jolloin vastaava piste tasossa on (x(t), y(t)).
Muussa tapauksessa valitaan jokin parametrointi. Esimerkiksi
xi(t)
= ai + bit
+ cit2
+ dit3,
yi(t)
= a'i + b'it
+ c'it2
+ d'it3.
missä 0 <= t <= 1.
z(t) = (1-t)3z1 + 3 t(1-t)2z2 + 3 t2(1-t)z3 + t3z4,
missä 0 <= t <= 1 ja pisteet zi ovat käyrän neljä kontrollipistettä.
Pierre Bezier alkoi käyttää näitä 1960-luvulla tietokoneavusteisessa suunnittelussa.
Bezier-käyrä määritellään kahden päätepisteen ja kahden kontrollipisteen avulla. Käyrä kulkee aina päätepisteiden kautta. Kontrollipisteillä ilmoitetaan tangenttien suunnat päätepisteissä.
Käyrä lähtee päätepistestä kontrollipisteen suuntaan sitä suorempana mitä kauempana kontrollipiste on.
Sopii hyvin juuri interaktiiviseen käyttöön.
