• No results found

2.9. Noggrannhet och precision vid numeriska ber¨ akningar

N/A
N/A
Protected

Academic year: 2021

Share "2.9. Noggrannhet och precision vid numeriska ber¨ akningar"

Copied!
25
0
0

Loading.... (view fulltext now)

Full text

(1)

forts. p˚a f¨oreg˚aende f¨orel¨asning:

I allm¨anhet l¨onar det sig att s¨atta ut decimalpunkterna f¨or reella tal. Om man inte enkelt kan skriva talen under varandra s˚a kan det l¨att uppst˚a fel vid inmatningen, vilket visas av f¨oljande exempel. Antag, att vi anv¨ander READ-satsen

READ (*, 100) X,Y,Z,T,U,V 100 FORMAT (6F5.2)

f¨or att l¨asa in talen

1.8 6.5 14.9 5.1 1.2 0.4 17.4 25. 4. -6.1 4.8 -7.8

Vid l¨asningen av den andra raden uppst˚ar d˚a ett fel (kan du se var?). Om inmatningen sker fr˚an en terminal,

¨ar det d¨arf¨or oftast bekv¨amare att anv¨anda ett fritt format (som anges med en asterisk (*)).

N¨ar man l¨aser in teckenstr¨angar, r¨acker det med att anv¨anda A-formatet som s˚adant, utan att ange antalet tecken, eftersom teckenstr¨angens l¨angd anges i CHARACTER-deklarationen. Inga extra tecken kommer s˚aledes att l¨asas in. Ett mindre antal tecken kan man d¨aremot nog l¨asa in, s˚asom framg˚ar av nedanst˚aende exempel:

(2)

CHARACTER (10) :: A ...

READ (10, 100) A 100 FORMAT (A1) ...

H¨ar l¨ases endast f¨orsta tecknet in i A, de ¨ovriga 9 tecknen fylls ut med blanka tecken.

Om det finns f¨alt som man inte vill l¨asa in, kan man anv¨anda X-specifikationen (eller TR-specifikationen), som framg˚ar av exemplet nedan:

READ (5,100) X, Y

100 FORMAT (F5.2, 7X, F8.4)

Man kan ocks˚a hoppa ¨over hela poster vid inmatningen, genom att anv¨anda formatspecifikationen /.

Observera i detta fall, att t.ex. instruktionen READ (5, 100) X, Y

100 FORMAT (F5.2/F10.4)

kommer l¨asa X fr˚an en post och Y fr˚an f¨oljande.

(3)

Vill man verkligen hoppa ¨over en post mellan X och Y, m˚aste man anv¨anda formatet 100 FORMAT (F5.2//F10.4). Man kan ocks˚a hoppa ¨over en post vid inl¨asningen med ett tomt format:

READ (5, 100) 100 FORMAT ()

F¨or att l¨asa och skriva p˚a andra enheter ¨an en sk¨armterminal, anv¨ands logiska enhetsnummer. I Fortran reserveras vanligen 5 f¨or standard-in och 6 f¨or standard-ut. F¨or att l¨asa och skriva p˚a skivfiler (t.ex.) kan man anv¨anda andra enhetsnummer.

En skivfil ¨oppnas i allm¨anhet f¨or inmatning eller utskrift med en OPEN-sats. N¨ar man ¨ar f¨ardig, st¨anger man den med en CLOSE-sats (observera, att huvudprogrammets END-sats ocks˚a automatiskt st¨anger alla ¨oppna filer). En OPEN-sats kan t.ex. se ut s˚a h¨ar: OPEN (10, FILE=’fil.dat’, STATUS=’NEW’). Denna sats kommer att ¨oppna en skivfil i enhet 10 med namnet fil.dat, som inte existerar f¨orut (om den existerar, anv¨ander man specifikationen STATUS=’OLD’). Om man anv¨ander parametern STATUS=’SCRATCH’ kom- mer filen att strykas, d˚a programmet avslutas. OPEN-satsen inneh˚aller som vi ser en serie specifikationer,

˚atskiljda av kommatecken (observera apostroferna, som omger namnet som st˚ar efter likhetstecknet). En annan viktig specifikation i OPEN-satsen ¨ar ACCESS, som antingen kan vara ACCESS=’SEQUENTIAL’ eller ACCESS=’DIRECT’. Denna specifikation anger filtypen. I det f¨orra fallet ¨ar filen sekventiell, dvs den kan l¨asas bara i en riktning (utg˚angsantagandet ifall man inte anger filtypen). I det senare fallet kan den l¨asas i vilken riktning som helst (direktaccess).

(4)

Varje fil t¨ankes uppdelad p˚a poster (eng. ’record’), som f¨or textfiler kan tolkas som rader. Om filtypen till˚ater direkt access, m˚aste ocks˚a storleken av varje post anges med specifikationen RECL=num, d¨ar num anger antalet bytes (t.ex. 80).

CLOSE-satsen som st¨anger en fil, har ocks˚a flere specifikationer, men vanligen anv¨ands den i formen CLOSE (num) (eller CLOSE (UNIT = num)), d¨ar num anger enhetsnummern. En sats som ¨ar nyttig f¨or sekventiella filer, ¨ar REWIND, som anv¨ands f¨or att b¨orja behandla filen igen fr˚an b¨orjan. Som ett exempel kan vi studera en del av ett program, som f¨orst skriver data p˚a en fil, och sedan l¨aser samma fil p˚a nytt:

OPEN (10, FILE=’test.dat’, STATUS=’NEW’) WRITE (10, 10) x,y, (z(i), i=1,n)

REWIND 10

READ (10,10) a,b, (c(k), k=1,n) 10 FORMAT (10F8.4)

Observera, att data i en fil alltid l¨ases p˚a samma s¨att som de skrivits! Ist¨allet f¨or att skriva formaterade data, som i fallet ovan, kan man ocks˚a skriva oformateradedata. Is˚afall b¨or OPEN-satsen inneh˚alla specifikationen FORM=’UNFORMATTED’, och formatnummern utel¨amnas d˚a vid l¨asning och skrivning. I detta fall kommer data att skrivas i bin¨ar form, vilket leder till t¨atare packning (men g¨or det sv˚arare att anv¨anda datafilen p˚a en annan dator).

(5)

Som ett exempel skall vi studera ett program, som l¨aser en fil, som skrivits oformaterat, och skriver ut formaterade data p˚a en annan fil.

PROGRAM bin2text

! Program som konverterar en bin¨ar fil till text IMPLICIT NONE

REAL :: x(100), y(100), a(100,100) INTEGER :: i,j,k,n

OPEN (10, FILE=’bin.dat’, FORM=’UNFORMATTED’, STATUS=’OLD’) READ (10) n

READ (10) (x(i), i=1,n),(y(i), i=1,n) READ (10) ((a(i,j),j=1,n), i=1,n) CLOSE (10)

OPEN (20, FILE=’asc.dat’, STATUS=’NEW’) WRITE (20, *) ’Elementen i vektorn x:’

WRITE (20, 100) DO i=1,n,8

WRITE (20, 100) (x(j), j=i,min(i+7,n)) END DO

WRITE (20, ’(/A)’) ’ Elementen i vektorn y:’

WRITE (20, 100) DO i=1,n,8

WRITE (20, 100) (y(j), j=i,min(i+7,n)) END DO

(6)

WRITE (20, ’(/A)’) ’ Elementen i matrisen A:’

DO i=1,n

DO j=1,n,8

WRITE (20, 100) (a(i,k), k=j, min(j+7,n)) END DO

END DO STOP

100 FORMAT (T1,8F9.4) END PROGRAM bin2text

Som vi ser, kommer programmet att skriva (h¨ogst) ˚atta tal per rad (observera anv¨andningen av funktionen min!). Vid utskriften anv¨ands implicerade DO-satser, som vi tidigare har n¨amnt. Formatbeteckningen

’(/A)’ leder till att texten f¨oreg˚as av en tom rad (anges av snedstrecket). Den bin¨ara filen bin.dat i ovanst˚aende exempel skulle man t.ex. kunna alstra med programmet

PROGRAM bintest

! Program som genererar en bin¨ar testfil IMPLICIT NONE

REAL :: x(100), y(100), a(100,100) INTEGER :: n, i, j

! L¨as in en l¨amplig dimension (h¨ogst 100):

(7)

DO

PRINT *, ’Matrisdimensionen (h¨ogst 100):’

READ *,n

IF (n > 0 .and. n <= 100) EXIT END DO

! Konstruera matriserna:

DO i=1,n

x(i) = 10.*i y(i) = i/10.

DO j=1,n

a(i,j) = i*j END DO

END DO

! ¨Oppna filen och skriv ut matriserna

OPEN (10, FILE=’bin.dat’, FORM=’UNFORMATTED’, STATUS=’NEW’) WRITE (10) n

WRITE (10) (x(i), i=1,n), (y(i), i=1,n) WRITE (10) ((a(i,j), j=1,n), i=1,n)

STOP

END PROGRAM bintest

(8)

Matrisdimensionen testas h¨ar med en ¨oppen DO-konstruktion f¨or att kontrollera, att den h˚aller sig innanf¨or stipulerade gr¨anser. Med n = 5 ser utskriften fr˚an bin2text ut p˚a f¨oljande s¨att:

Elementen i vektorn x:

10.0000 20.0000 30.0000 40.0000 50.0000 Elementen i vektorn y:

0.1000 0.2000 0.3000 0.4000 0.5000 Elementen i matrisen A:

1.0000 2.0000 3.0000 4.0000 5.0000 2.0000 4.0000 6.0000 8.0000 10.0000 3.0000 6.0000 9.0000 12.0000 15.0000 4.0000 8.0000 12.0000 16.0000 20.0000 5.0000 10.0000 15.0000 20.0000 25.0000

Vi har tidigare beskrivit hur man ¨oppnar filer f¨or in- och utmatning i Fortran. Vi skall nu diskutera ytterligare n˚agra parametrar, som ing˚ar i OPEN-satsen.

(9)

En av dem ¨ar IOSTAT, som anger en heltalsvariabel, vars v¨arde ¨ar lika med 0, om filen ¨oppnats utan problem, men annars olika 0. Ett exempel p˚a anv¨andningen ges av f¨oljande programfragment, som fr˚agar efter namnet p˚a en existerande fil, som ¨oppnas. Ifall n˚agot problem uppst˚ar, ges ett felmeddelande.

PRINT *, "Ange filnamnet"

READ ’(A)’, filnamn

OPEN (UNIT=1, FILE=filnamn,STATUS="OLD",IOSTAT=ios) IF (ios /= 0) THEN

PRINT *, "Kan inte ¨oppna ",filnamn ...

END IF

Parametern IOSTAT kan ocks˚a anv¨andas i in- och utmatningssatser. V¨ardet 0 anger d˚a att satsen utf¨orts utan problem, medan ett positivt v¨arde anger att ett fel uppst˚att. Vid in- och utmatning anger ett negativt v¨arde (oftast) att filen tagit slut. Filens slut anges av ett s¨arskilt sluttecken i filen, som skrivs med satsen ENDFILE n, d¨ar n anger enhetens nummer. Detta kan t.ex. utnyttjas om man vill l¨agga till data i slutet av en sekventiell fil, s˚asom visas av f¨oljande exempel:

(10)

...

! Hoppa till slutet DO

READ (UNIT=1, IOSTAT=ios) z

IF (ios<0) EXIT ! negativ status betyder slutet END DO

! G˚a ett steg tillbaka BACKSPACE 1

! Skriv ny information WRITE (1) x

...

Man kan ocks˚a l¨agga till data i slutet av en fil genom att ¨oppna den med parametern POSITION =

’APPEND’. BACKSPACE-satsen beh¨ovs alltid, om filen inneh˚aller ett sluttecken (beror p˚a systemet). Ett dylikt sluttecken kan ocks˚a explicit skrivas in i filen med ENDFILE-satsen. Ist¨allet f¨or IOSTAT kan man i I/O-satser ocks˚a anv¨anda parametern END = num, som ¨overflyttar kontrollen till en instruktion med sats- nummern num, ifall filen tar slut. Anv¨andningen av IOSTAT rekommenderas dock, eftersom den m¨ojligg¨or b¨attre kontroll ¨over vad som verkligen h¨ander.

Ytterligare exempel p˚a filhantering ges i f¨oljande program, som s¨oker upp ett namn i en fil som inneh˚aller t.ex. namn p˚a anv¨andare och anv¨andarkoder, och skriver ut anv¨andarkoden, ifall namnet kan hittas (f¨or att avsluta s¨okningen, anv¨ands * *).

(11)

PROGRAM search

! Program som s¨oker upp anv¨andarnamn i en sekventiell fil IMPLICIT NONE

CHARACTER(LEN=20) :: fil_namn, f_namn, t_namn, f_n, t_n, a_n INTEGER :: ios

LOGICAL :: funnet

! Ange filnamn och ¨oppna filen:

DO

PRINT *, "Ange filnamn:"

READ *, fil_namn

OPEN (UNIT=1, FILE=fil_namn, STATUS="OLD", IOSTAT=ios) IF (ios == 0) EXIT

PRINT *, "Fel p˚a filen!"

END DO

! Ange f¨ornamn och tillnamn DO

funnet = .FALSE.

PRINT *, "F¨ornamn och tillnamn:"

READ *, f_namn, t_namn

IF (f_namn == "*" .AND. t_namn == "*") EXIT REWIND (1)

! S¨ok upp den r¨atta posten DO

READ (UNIT=1, FMT=*, IOSTAT=ios) f_n, t_n, a_n

(12)

IF (ios<0) THEN

! Filen tog slut EXIT

ELSE IF (ios>0) THEN

! Fel i filen

PRINT ’("Fel nr ",I3," vid l¨asning av fil ",A)’,ios,fil_namn STOP

END IF

IF (f_namn == f_n .AND. t_namn == t_n) THEN PRINT *, "Anv¨andarkoden = ",a_n

funnet = .TRUE.

END IF END DO

IF (.NOT. funnet) PRINT *, "Namnet ej funnet!"

END DO STOP

END PROGRAM search

Vid in- och utmatning i Fortran brukar satserna READ, WRITE och PRINT vanligen inleda behandlingen av en ny post. Det finns fall d˚a man skulle vilja l¨asa endast en del av en post, och resten senare. Detta, som inte var m¨ojligt med FORTRAN 77, kan g¨oras med Fortran 90.

(13)

Metoden kallas icke-framskridande (”non-advancing”) in- och utmatning, och kan anv¨andas endast vid formaterad behandling av sekventiella filer. I detta fall beh¨over man inte behandla en hel post s˚asom fallet ¨ar vid normal filhantering. Dessutom f˚ar man reda p˚a antalet tecken som posten inneh˚aller, samt underr¨attas om filen tagit slut under l¨asningen.

Icke-framskridande in- och utmatning anges med parametern ADVANCE i en READ eller WRITE sats, som har formen ADVANCE = "YES" eller ADVANCE = "NO". I det f¨orstn¨amnda fallet anv¨ands den normala, framskridande metoden, medan man i det senare fallet anv¨ander den nya, icke-framskridande metoden.

Om ADVANCE ="NO" ing˚ar i en READ-sats, finns det tre m¨ojligheter:

1) Posten l¨astes inte i sin helhet, och det uppstod inga fel under l¨asningen. I detta fall kommer n¨asta READ-sats att forts¨atta med behandlingen av samma post.

2) Om posten blivit l¨ast f¨orbi sluttecknet (vilket kan testas med IOSTAT), kommer positionen i filen att vara omedelbart efter ifr˚agavarande post. Detta uppfattas dock inte som ett feltillst˚and.

3) Om n˚agot fel intr¨affade under l¨asningen (t.ex. filen tog slut) kommer positionen i filen att vara omedel- bart efter senast l¨asta post.

F¨or att f˚a reda p˚a hur m˚anga tecken som blivit l¨asta, anv¨ands parametern SIZE. S˚alunda kommer t.ex.

den icke-framskridande READ-satsen

(14)

READ (UNIT=1, FMT=100, ADVANCE="NO", SIZE =nc, &

IOSTAT=ios) x,y,z

att lagra antalet l¨asta tecken i variabeln nc, och heltalsvariabeln ios kommer att ange hur l¨asningen har lyckats.

Vid icke-framskridande utmatning finns det bara tv˚a m¨ojligheter, antingen uppstod inga fel, och d˚a kommer n¨asta WRITE-sats att forts¨atta med samma post, eller ocks˚a uppstod det n˚agot fel, och positionen i filen kommer att vara efter senast skrivna post.

De flesta filer, som vi hittills haft att g¨ora med, har varit sekventiella, dvs de kan l¨asas och skrivas bara i en riktning. Men man kan ocks˚a behandla posterna i en fil i en godtycklig ordning, om man ¨oppnar filen med parametern ACCESS="DIRECT". I detta fall m˚aste alla poster i filen ha samma l¨angd, vilket anges genom specifikationen RECL =n, d¨ar n ¨ar ett heltal, som anger den konstanta postl¨angden. Poster- na i en direktaccess-fil uppfattas vara numrerade. S˚alunda kan man t.ex. skriva post 100 oformaterat med satsen WRITE (UNIT=1, REC=100) x,y,z. En formaterad utskrift av formen WRITE (UNIT=1, REC=100, FMT=’(10I5)’) (m(i), i=1,100) anger att tio poster kommer att skrivas, utg˚aende fr˚an post 100. Observera ocks˚a d˚a man ¨oppnar en fil f¨or direkt access, att parametern FORM =’UNFORMATTED’

¨ar underf¨orst˚add.

(15)

Posterna i en fil med direkt access kan skrivas i vilken ordning som helst. En post som har blivit skriven, kan skrivas ¨over p˚a nytt, men inte strykas. Om postens l¨angd ¨ar kortare ¨an specifikationen RECL, kommer den att fyllas ut till r¨att l¨angd.

Ett exempel p˚a ett program som anv¨ander direkt access f¨or att l¨asa en fil, skriven med direkt access, visas nedan.

PROGRAM direct

! Programmet testar direkt access IMPLICIT NONE

CHARACTER :: post*80, fil_namn*30 INTEGER :: ios, post_nr

! Ange filnamn och ¨oppna filen:

DO

PRINT *, "Ange filnamn:"

READ *, fil_namn

OPEN (UNIT=1, FILE=fil_namn, ACCESS="DIRECT", &

FORM="FORMATTED", STATUS="OLD", RECL=50, IOSTAT=ios) IF (ios == 0) EXIT

PRINT *, "Fel p˚a filen!"

END DO

! Ange postnummer DO

PRINT *, ’Ange postnummer (0 avslutar):’

(16)

READ *, post_nr

IF (post_nr <= 0) EXIT

! S¨ok upp posten, och skriv ut den

READ (UNIT=1, FMT = ’(A)’, REC = post_nr, IOSTAT = ios) post IF (ios == 0) THEN

PRINT *, post ELSE

PRINT *, "L¨asfel!"

END IF END DO

PRINT *, "Filen behandlad"

STOP

END PROGRAM direct

Om filen inte kunde ¨oppnas, kommer heltalsvariabeln ios att f˚a ett v¨arde olika noll. Detta sker b˚ade ifall filnamnet ¨ar or¨att skrivet, eller filen har fel format (t.ex. variabel postl¨angd, etc.). Om det angivna postnumret inte ¨ar positivt, avbryts programmet, medan postnumret i annat fall anges p˚a nytt. Om man anger ett st¨orre postnummer ¨an vad som finns i filen, kommer heltalsvariabeln ios att f˚a ett v¨arde olika noll, vilket anges som ”l¨asfel”.

(17)

Filer av det slag som vi tillsvidare diskuterat, brukar ocks˚a kallas yttre filer f¨or att skilja dem fr˚an s.k. inre filer, som egentligen inte ¨ar n˚agra filer, men beter sig som s˚adana. En inre fil ¨ar egentligen bara ett s¨att att konvertera mellan olika format. En dylik fil ¨ar vanligen en variabel av typ CHARACTER. B¨ast f¨orklarar man anv¨andningen med ett exempel.

CHARACTER (LEN=50) :: rad

WRITE (UNIT=rad, FMT=’(2F10.5,3I10)’) x,y,i,j,k

READ (UNIT=rad, FMT=’(2(F8.3,2X),3(I8,2X))’) x,y,i,j,k

WRITE-satsen kommer h¨ar att alstra en teckenstr¨ang i variabeln rad som best˚ar av v¨ardena av de reella variablerna x och y med fem decimaler och tio tecken vardera, samt heltalsvariablerna i, j och k ocks˚a med tio tecken var. Dessa variabler matas in i samma teckenstr¨ang och f¨orkortas till ˚atta tecken var.

(18)

2.9. Noggrannhet och precision vid numeriska ber¨ akningar

Nogrannhet i datorframst¨allningen av ett flyttal best¨ams av antalet bitar i mantissan, och vilka tal som kan framst¨allas p˚a detta s¨att best¨ams av vilka v¨arden exponenten kan anta. Detta beror givetvis p˚a datorns ordl¨angd, vilket st¨aller till med problem d˚a man ¨overf¨or program mellan olika datorer. F¨or att avhj¨alpa detta, har man inf¨ort parametriserade variabler i Fortran 90 vilket g¨or det m¨ojligt att kontrollera precisionen och exponentens r¨ackvidd. Parametriseringen sker med en s¨arskild parameter KIND.

Vi skall belysa detta med n˚agra exempel:

REAL :: x, y, z

REAL, DIMENSION(100) :: a REAL(KIND=4) :: b, c

REAL(KIND=2) :: v(5)

Variablerna x, y, z och vektorn a ¨ar alla ”vanliga” variabler av typ REAL, medan b och c har KIND typ 4 och vektorn v har KIND typ 2. F¨or att ta reda p˚a KIND-typen, kan man anv¨anda en inbyggd funktion KIND. S˚alunda kommer t.ex. satsen i = KIND(b), om b deklarerats som ovan, att ge v¨ardet 4 ˚at variabeln i. Det ¨ar inte n¨odv¨andigt att uts¨atta parameternamnet KIND i en deklaration, det r¨acker att skriva t.ex.

REAL(4), men det blir tydligare om man anv¨ander den fullst¨andigare formen REAL(KIND=4).

(19)

I allm¨anhet ¨ar KIND-typerna datorberoende, men det finns en m¨ojlighet i Fortran 90 att konstruera KIND- parametrar, som ¨ar helt datoroberoende. Detta sker med den inbyggda funktionen SELECTED REAL KIND, som kan ha tv˚a heltaliga argument, P och R (det finns ¨aven SELECTED INT KIND). Det f¨orst- n¨amnda argumentet P anger minimiantalet decimaler som beh¨ovs, och det andra argumentet R anger den minsta r¨ackvidden f¨or exponenten som beh¨ovs. Om man t.ex. deklarerar en variabel med satsen REAL(KIND=SELECTED REAL KIND (P=7,R=30)) :: a, s˚a betyder det att man f¨or att framst¨alla a beh¨over minst 7 decimaler och att exponentens r¨ackvidd b¨or vara minst 30. Om samma noggrannhet an- v¨ands i flera deklarationer, l¨onar det sig att definiera en s¨arskild konstant f¨or KIND-parametern, t.ex. p˚a f¨oljande s¨att:

INTEGER, PARAMETER :: real_7_30 = SELECTED_REAL_KIND(P=7,R=30) ...

REAL(KIND=real_7_30) :: a ...

I de flesta datorer ¨ar det m¨ojligt att lagra flyttal i enkel eller dubbel precision, och de har h˚ardvara f¨or att utf¨ora aritmetiska operationer med dem. Om datorn r¨aknar med sex signifikanta siffor och har en exponentr¨ackvidd av 40, kommer den i ovanst˚aende fall att r¨akna med dubbel precisions noggrannhet, men om ordl¨angden ˚a andra sidan ¨ar s˚a l˚ang, att enkel precision r¨acker till f¨or att framst¨alla tal med 15 signifikanta siffror, och exponentens r¨ackvidd ¨ar 300, s˚a r¨acker enkel precision till i detta fall.

(20)

F¨or att b¨attre inse betydelsen av precision vid numeriska ber¨akningar skall vi studera ett program, som ber¨aknar uttrycket

s

 1 2

2

3 . . . 2n − 1 2n

2

·  2 1

3

2 . . . 2n 2n − 1



f¨or olika v¨arden av n, och med olika stor precision (som man f˚ar om man ¨andrar argumentet P i SELECTED REAL KIND:

PROGRAM prec_test IMPLICIT NONE

INTEGER, PARAMETER :: real6 = SELECTED_REAL_KIND(P=6,R=32) REAL(KIND=real6) :: x, p, q, s

INTEGER :: k, n DO

PRINT *, ’Ge v¨ardet av n=’

READ *,n

IF (n <= 0) EXIT x = 1.

p = 1.

DO k=1,2*n-1 x = REAL(k)

p = (x/(x+1.))*p END DO

x = 1.

(21)

q = 1.

DO k=1,2*n-1 x = REAL(k)

q = ((x+1.)/x)*q END DO

s = SQRT(p*p)*q PRINT *, s

END DO STOP

END PROGRAM prec_test

Det matematiskt korrekta resultatet ¨ar 1, men avrundningsfel kommer att leda till fel, som ¨ar s¨arskilt markanta vid l˚ag precision (t.ex. 6 decimaler), vilket framg˚ar av nedanst˚aende exempel:

OUTO...ADB_KURS>r prec_test2 Ge v¨ardet av n=

1000

1.000008

Ge v¨ardet av n=

100000 1.000280

Ge v¨ardet av n=

1000000 1.004621

(22)

Om man ¨okar precisionen (t.ex. till 14 decimaler), minskar avvikelsen betydligt, vilket man ser genom att testk¨ora en modifikation av programmet ovan:

OUTO...ADB_KURS>r prec_test3 Ge v¨ardet av n=

1000

0.999999999999998 Ge v¨ardet av n=

100000

0.999999999999947 Ge v¨ardet av n=

1000000

1.00000000000001

Vi skall ocks˚a studera hur KIND–typerna kan anv¨andas f¨or ange precision i numeriska metoder. En av de enklaste metoderna att numeriskt ber¨akna en rot till en ekvation brukar kallas funktionsiteration. Om man p˚a n˚agot s¨att kan skriva ekvationen f (x) = 0 i formen x = g(x), s˚a kan man iterera fram en rot genom utg˚a fr˚an en approximation till roten, substituera den i h¨ogra membrum samt r¨akna ut en ny approximation till roten. Detta betyder att man ber¨aknar roten genom successiva approximationer: xn+1 = g(xn), n = 0, 1, 2, . . .. Ber¨akningen avbryts, d˚a skillnaden mellan tv˚a p˚a varandra f¨oljande approximationer ¨ar mindre ¨an en given toleransgr¨ans: xn+1 − xn < .

(23)

Metoden har till¨ampats i f¨oljande program, som ber¨aknar den reella roten till ekvationen x3−2x−5 = 0.

Denna ekvation, som studerades redan av Newton 1669, kan skrivas i formen x = (2x + 5)1/3, vilket g¨or det m¨ojligt att till¨ampa funktionsiteration.

MODULE precision

! Modul som anger variablernas precison IMPLICIT NONE

INTEGER, PARAMETER :: r=SELECTED_REAL_KIND(P=14,R=30) END MODULE precision

MODULE functions IMPLICIT NONE CONTAINS

FUNCTION fn(x) USE precision IMPLICIT NONE

REAL(KIND=r) :: fn

REAL(KIND=r), INTENT(IN) :: x fn = (2._r*x + 5._r)**(1._r/3._r) RETURN

END FUNCTION fn END MODULE functions PROGRAM fun_it

! Program f¨or iterativ rotbest¨amning: x = fn(x) USE precision

(24)

USE functions IMPLICIT NONE

REAL(KIND=r) :: x, x1, dx INTEGER :: i, maxit

PRINT *, "Ge en approximation till roten"

READ *, x1

PRINT *, "Ge maximiantalet iterationer"

READ *, maxit

PRINT *, "Ge ¨onskad precision"

READ *, dx

! Iterera:

DO i=1,maxit x = fn(x1)

IF (ABS(x-x1) < dx) EXIT x1 = x

PRINT *, " x = ",x END DO

IF (i > maxit) PRINT *, "Proceduren har inte konvergerat"

PRINT *, "Antal iterationer: ",i PRINT *, " x = ",x

STOP

END PROGRAM fun_it

(25)

I detta program, liksom i det f¨oreg˚aende programmet, anv¨andes 14 decimalers precision, som definierats i en s¨arskild modul, kallad ”precision”. Variabler med denna precision har deklarerats som REAL(KIND=r).

P˚a ett liknande s¨att har ocks˚a konstanter angetts med samma precision, t.ex. 2. r.

Konvergensen ¨ar r¨att snabb, med utg˚angsv¨ardet x = 2 f˚as ett resultat med tio korrekta decimaler efter 12 iterationer, och med utg˚angsv¨ardet x = 100 efter 15 iterationer:

Ge en approximation till roten 2 Ge maximiantalet iterationer 30 Ge ¨onskad precision 1.E-10

x = 2.08008382305190 x = 2.09235067779758 x = 2.09421699601252 x = 2.09450065219465 x = 2.09454375753281 x = 2.09455030780827 x = 2.09455130318276 x = 2.09455145443897 x = 2.09455147742373 x = 2.09455148091647 x = 2.09455148144722

Antal iterationer: 12 x = 2.09455148152787

References

Related documents

Resonemang, inf¨ orda beteck- ningar och utr¨ akningar f˚ ar inte vara s˚ a knapph¨ andigt presenterade att de blir sv˚ ara att f¨ olja.. ¨ Aven endast delvis l¨ osta problem kan

Tillräckliga data saknas ännu för TAF men i en liten studie på hivnegativa kvinnor kunde den aktiva metaboliten av tenofovir (tenofovir diphos- phate) inte detekteras

Rutinen som anv¨ands f¨ or att definiera operatorn, kan ha antingen ett eller tv˚ a argument, men eftersom funktionen normalt definieras i samma modul som inneh˚

F¨ or att g¨ ora det litet noggrannare kan du ber¨ akna normen av skillnaden mellan resultatvektorerna med MATLAB–funktionen norm, och j¨ amf¨ ora normerna f¨ or olika v¨ arden

F¨ ors¨ ok f¨ orb¨ attra programmet genom pivotering, dvs upps¨ ok elementet med det st¨ orsta absoluta v¨ ardet i A(k:n,k) (observera, att k-1 m˚ aste adderas f¨ or att f˚ a

Pr¨ ova olika summeringsordning, och f¨ ors¨ ok ocks˚ a summera positiva och negativa termer skilt f¨

Skriv en funktion som ber¨ aknar utbrytningshastigheten utg˚ aende fr˚ an planetens massa och radie4. Skriv ut resultatet i tabellform, med fem v¨ arden p˚ a

Skriv ett Fortran-program, som l¨ aser filens namn, kontrollerar dess existens, l¨ aser in spektret i minnet, byter om intensiteternas ordningsf¨ oljd, och skriver ut det