73 
UPORABNI PROGRAMI 
Prilagodljivo numerično integriranje 
% Informatica UP 18 
% Adaptive Numerical Integration 
% oktober 1984 
% pripravil Vladimir Batagelj 
% sistem DEC-10 
% 
% 
% 
% 
POSTOPKI ZA PRIUGODLJIVO NUMERIČNO INTEGRIRANJE 
Vladimir Batagelj 
VTOZD Matematika in mehanika 
Univerza E. Kardelja v LJubljani 
Namen tega sestavka Je opozoriti na postopke za 
prilagodljivo (adaptivno) numerično integriranje in na 
vire I 2, ^, 5, 6j , kjer so podrobneje opisani. 
Postopki za prilagodljivo numerično integriranje so 
področje, na katerem se uspešno prepletata numerična 
analiza in računalništvo, natančneje, načrtovanje in 
analiza algoritmov. Tako je na primer v numeričnl 
knjižnici NAG podprogram D01AGA/F, ki temelji na 
Oliverjevi prilagodljivi metodi, priporočen kot najboljši 
splošni integracijski podprogram. 
Če želimo pri' običajnih postopkih za numerično 
integriranje C 1 3 (»večati natančnost Izračuna, 
povečujemo število n točk ali vozlov, ki jih upoštevamo 
pri Izračunu. S tem sicer zmanjšamo napako metode, žal pa 
se nam pri tem zaradi večanja števila operacij večata 
zaokrožltvena napaka in čas (cena) računanja. Zato se 
obnaša celotna napaka tako, kot prikazuje slika: 
napaka 
Pri prilagodljivih postopkih pa drobimo interval samo tam, 
kjer Je treba. Zato običajno za. dosego željene natačnostl 
potrebujejo manj korakov - pri (skoraj) isti napaki metode 
bo manjša zaokrožltvena napaka ps ' tudi postopek je 
hitrejši in s tem bolj dconcolčen. 
Osnovna zamisel prilagodljivih postopkov numeričnega 
integriranja je znana v računalništvu kot pristop "deli 
in vladaj". Recimo, da želimo izračunati vrednost 
P(a,b,eps) Integrala: 
I(a,b) = j f(x) dx 
pri čemer dopuščamo napako eps . Naj bo še S(a,b) ena 
izmed kvadraturnih formul, ki jih poznamo iz numerične 
analize [17 , In dajejo oceno za vrednost I(a,b) . 
Določimo oceno S. = S(a,b) za neznani I(a,b) . Nato 
razpolovimo Interval Cajb"] a točko m = (a+b)/2 in 
izračunamo še oceno S- = S(a,m) * S(m,b) . Če se oceni 
razlikujeta za manj kot eps , vzamemo Sp za 
P(a,b,eps) . Sicer isti postopek rekurzivno ponovimo na 
obeh podlntervalih £a,m] in [m,b] , pri čemer pa na 
vsakem podintervalu dopuščamo le pol manjšo napako 
eps/2 . Ko enkrat poznamo P(a,m,eps/2) in 
P(m,b,eps/2) , lahko izračunamo tudi 
P(a,b,eps) = P(a,m,eps/2) + P(m,b,eps/2) 
V praksi se Izkaže, da Je zmanjšanje dovoljene napake na 
eps/2 pri rekurzivnl ponovitvi postopka preveč 
pesimistično. Zato v postopek vpeljemo nov parameter q , 
1 ^q <2 in dovoljujemo napako eps/q . Izkaže se £33, 
da je čisto v redu že q = 1.4 . 
V tem sestavku je opisani postopek razdelan za primer, ko 
oceno vrednosti Integrala na intervalu [Ja,b] računamo 
po Simpsonovem obrazcu: 
S(a,b) = ^ ( f(aO + i).f(^) + f(b) ) 
Treba pa Je takoj pripomniti, da podprogrami iz numeričnlh 
knjižnic uporabljajo veliko bolj "zvite" ocene in tudi sam 
postopek Je še nadalje izpopolnjen. 
Za primerjavo Je dodan še navadni postopek za numerično 
integrirčinje po Slmpsonu. 
V obeh programih se nahajajo spremenljivke, za merjenje 
posameznih značilnih količin: 
i - zaporedna številka 
countf - število izračunanih funkcijskih vrednosti 
depth - globina rekurzije 
Te spremenljivke in stavki, ki Jih vsebujejo, ne vplivajo 
na samo računanje, zato Jih lethko tudi izločimo. 
Za to, da bodo prednosti prilagodljivega postopka očitne. 
Je bila primerjava narejena na integralu (ploščina 
I četrtine enotskega. kroga): 
1 
if 
2 
yr dx 
ki Je za navadni postopek neugoden. Njegova štirikratna 
vrednost Je število pi . 
Priloženi tabeli govorita sami zase. Omenimo le pomen 
, posameznih stolpcev: 
NAVADNO: zaporedna številka. Izračunana vrednost 
integrala, dejanska napaka, število izračunanih 
funkcijskih vrednosti; 
PRILAGODLJIVO: zaporedna številka, izračunana vrednost 
integrala, dovoljena napaka, de.lanska napaka, število 

7« 
izračunanih Ainkoijskih vrednosti, največja globina 
rekurzije. 
Izračuni so bili opravljeni na računalniku DEC-10 , 
Univerze v LJubljani. 
PROGRAM slmpsonCoutput); 
CONST pi : 3.111592653589793 ; 
ernnln z 0.00000005 { 
VAR error int : real; 1 n, oountf : integer; 
FUNCTION f( x: real ); real ( 
BEGIN 
oountf ;s oountf • 1 
f := sqrt(ab3(1-sqr(x 
.END (• r «) j 
FUNCTION integral 
( FUNCTION f:real 
VAR h, p 
BEGIN 
X : real ; o. 
h := (b - a)/(2«n) j | 
X := a 
POR i 
X :: 
p ••: 
0 : 
END; 
; 0 := 1 ; 
:= 1 TO 2«n-
X + h ; 
: p + C«f(x) j 
: 6 - 0 
integral ;= hip/3 
END (I integral «] ; 
BEGIN 
)) 
a,b:real; n:integer ): real; 
i : integer ; 
J := f(a) + f(b) 
DO BEGIN 
» 
viritelnC ':5,'NAVADNO INTEGRIRANJE PO SIMPSONU •); 
writeln; 
n := 1 
REPEAT 
i := 0 J 
oountf := 0 ; 1 :3 
int 
1 + 1 ; 
:= 1»integral(f,0,1,n) ; 
error := pl - intj 
writeln( 1:7, int, 
n :: 2m 
UNTIL abs( error ) <= 
END. 
NAVADNO INTEGRIRANJE 
1 
2 
3 
. 1* 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
2.976068E+00 
3.083595E+00 
3.121ia9E+00 
3.131398E+00 
3.139052E+00 
3. iioegsE+oo 
3. l'tl276EtOO 
3.11118lEt00 
3.111553E+O0 
3. l'(1579E+00 
3.111588E+00 
3.111591E+00 
3.111592EtOO 
3.111592EtOO 
3.111592E+00 
3.l'(1592E+00 
3.111593E+00 
3.1l<1591EtOO 
3.111593E+00 
3.111601E+00 
error, countf:8 ); 
errmln 
PO SIMPSONU 
1.655219E-01 
5.799752E-02 
2.0')0317E-02 
7.19'(966E-03 
2.510'l39E-03 
8.975565E-01 
3.I7186IE-OI 
1.121163E-01 
3.978610E-05 
1.385808E-05 
1.917383E-06 
1.758337E-06 
5.662111E-07 
2.9fl0232E-07 
1.172325E-07 
5.066395E-07 
-8.9'l0697E-08 
1.817912E-06 
-8.9't0697E-08 
-8.II31057E-06 
3 
5 
9 
17 
33 
65 
129 
257 
513 
1025 
. 2019 
1097 
8193 
16385 
32769 
65537 
131073 
262115 
521289 
1018577 
PROGRAM sinpsonCoutput); 
CONST pi = 3.111592653589793 i 
errrnin s O.OOOOOOOV | 
VAR error, int : real; i, oountf, depth : integer; 
FUNCTION f( x: real ): real ; 
BEGIN 
countf := oountf -» 1 ; 
f := sqrt(abs(1-3qr(x))) 
END (1 f «) ; 
FUNCTION integral 
( FUNCTION f:real; a,b,error: 
VAR m, C, fa, fB, fb : real ; 
FUNCTION Irec 
real ): real; 
d : Integer ; 
( a, m, b, fa, tm. Vb, estimate, error:real 
CONST q = 1.5 ; 
VAR C, na, mb, ftna, fmb, esta. 
BEGIN 
estb, nevest 
d := d + 1 ; IF d > depth THEN depth := d ; 
ma := (a -f (n)/2 ; fpa := f(ma) ; 
Bib := (m + b)/2 ; JlJib :r f(nib) ; 
' C := (b - a)/12 ; 
esta :s C • ( f a + lifma + foi ) ; 
estb :: C • ( fm + lafmb + fb ) ; 
nevest := esta * estb ; 
IF absCestimate-neviest) > error THEN 
Irec := lreo(a,iiia,ro,fa,ftiia 
lrec(m,nib,b, fm, fmb 
ELSE 
Irec := nevest ; 
d:= d - 1 
END (« Irec •) ; 
BEGIN 
d := 0 ; m := (a + b)/2 ; 0 : 
fa := f(a) ; fin := fdn) ; fb 
,flii, esta, error 
,fb,estb,error 
= (b - a)/6 ; 
:= f(b) ; 
): real; 
real 
'q> + 
'q) 
i 
integral := irec(a,m,b,fa,fte,fb,c»(fa+1afin+fb),error) 
END (K Integral *) ; 
BEGIN 
vriteC ':15); 
vrltelnCPRILAGODLJIVO INTEGRIRANJE PO SIMPSONU')J 
vriteln ; error ;= 1 ; i := 0 
WHILE error > errmln DO BEGIN 
oountf := 0 ; depth := 0 ; 
J 
1 := 1 • 1 ; 
int := 1»inteKral(f,0,1,error) ; 
vrltelnCi,int,error,pl-lnt,oountf;6,depth:6); 
error := error / 2 
END 
END. 
PRILAGODLJIVO INTEGRIRANJE PO SIMPSONU 
1 3.083595E+00 1.000000E+00 
2 3.O83595E+OO 5.OOOOOOE-OI 
3 3.083595E+00 2.500000E-01 
1 3.083595E+00 1.250000E-01 
5 3.O83595E+OO 6.250000E-02 
6 3.O83595E+O0 3.125000E-02 
7 3.121189E+00 1.5fa2500E-02 
8 3.131383E+00 7.81250OE-O3 
9 3.139032E+O0 3.9O625OE-O3 
10 3.111253E+00 1.953125E-03 
11 3. ini58EtOO 9.765625E-01 
12 3.II153OE+OO 1.882813E-01 
13 3.111556E*00 2.111106E-01 
11 3.111565E+00 1.220703E-01 
15 3.111583E+00 6.IO3516E-O5 
16 3.111588E+00 3.051758E-05 
17 3.111590E+00 1.52587 9E-05 
18 3.111691E+00 7.629395E-06 
19 3.111592E+00 3.811697E-06 
20 3.111592E+0O 1.907319E-O6 
21 3.111592E+00 9.536713E-07 
22 3.111593E+00 1.768372E-07 
23 3.111593E+00 2.3811d6E-07 
21 3.111593E+00 1.192093E-07 
25 3.111593E+00 5.960161E-08 
26 3.111593E+00 2.980232E-08 
27 3.111593E+00 1.1901I6B-O8 
5.799752E-02 
5.799752E-02 
5.799752E-02 
5.799752E-02 
5.799752E-02 
5.799752E-02 
2.010350E-02 
7.2O9718E-O3 
2.560167E-03 
3.397763E-01 
1.318853E-01 
6.216567E-05 
3.686517E-05 
2.783537E-05 
9.891371E-06 
I.IIO7IIE-O6 
2.533197E-06 
l.ei771lE-06 
9.23872OE-07 
3.576279E-07 
2.682209E-07 
1.190116E-07 
8.9IO697E-O8 
2.980232E-08 
O.OOOOOOEtOO 
O.OOOOOOE>00 
O.OOOOOOE^O 
5 
5 
5 
5 
5 
5 
9 
13 
17 
25 
29 
33 
37 
11 
19 
57 
65 
73 
85 
105 
117 
133 
153 
185 
213 
253 
285 
2 
3 
1 
6 
7 
8 
9 
10 
11 
12 
13 
11 
15 
17 
18 
19 
20 
21 
22 
23 
21 

75 
LITERATURA: 
[ll Bohte Z.: Numerične metode, (v I. Vldav: Višja 
matematika III). LJubljana, 1973 
[2I De Boor C: CADRE: an algorlthm for numerloal 
guadrature. (v Mathematloal software, J.R.Rlce, ed.). 
Aoademlo Preaa, New york, Il17-1t9 
[3].Dennls J.E. Jr., More J.J.: Computer solutlon of 
matheaatical problema. ( v Conway R., Griea D.: An 
Intraductlon to Prograinnlng). Klnthrop, Cambrldge, Mass., 
1975, 358-366 
[ij] Kahaner D.K.: Comparlslon of numerloal guadrature 
formulas. (v Mathematloal software, J.R.Bloe, ed.). 
Aoademlo Press, New Kork, 229-259 
[5] NAG fortran llbrary manual, mark 6. NAG Ltd., 1977 
[6] Oliver J.: A doubly-adaptlve Clenshaw-Curtls quadrature 
method. The Computer Journal, 15(1972), 141-117 
Vsa literatura Je dosegljiva v Matematični knjižnici. 
Jadranska 19, LJubljana. 
Procedura magicodd v listi 1 oblikuje magični 
kvadrat podobno kot procedura maglceven, le da 
Je stopnja n matrike x liha In ni manjša od 3. 
Funkcija maglcterra v listi 1 izračuna element 
s(i,J) magičnega kvadrata s(l ,, n, 1 .. n) za 
liho vrednost h, ki nI.manjša od 3. 
Lista 1 vsebuje še procedur za izpis magičnega 
kvadrata kvadr_lzpi8 In proceduro za preizkuša­
nje magičnosti kvadrata magic_teat, ki najprej 
določi magično vsoto in nato preizkusi na to 
vsoto vse vrstice, vse stolpce in vse diagona­
le. -
Z glavnim preizkusnim programom liste 1 se upo­
rabijo zadevne procedure in funkcija (maglce- ' 
ven, magicodd in magioterm) z izpisom magičnih 
matrik In njihovim testiranjem. 
3. Izvajanje programa 
MAGIČNI KVADRATI SODE IN LIHE STOPNJE 
S = = 3I£SSS£ 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
% 
% 
% 
% 
% 
Informatica UP 19 
Magic Squares of Even and Odd Degree 
september 1984 
priredil Anton P, Železnikar 
sistem CP-M, Delta Partner 
prevajalnik Janus-Ada, verzija 1.5.0 
% 
% 
% 
% 
% 
1. Področje uporabe 
Problematika magičnih kvadratov sodi v podro­
čje matematične in programirne rekreacije. Tu 
al Je mogoče izmišljati različne algoritme za 
določevanje magičnih kvadratov sode in lihe 
stopnje. 
2. Opis programa 
Procedura maglceven v listi 1 določa: "za kva­
dratno matriko n krat n' elemente z vrednostmi 
od> 1 do n krat n in Jih razmeŠča v lekalko-
grafskem zaporedju v polje x(l .. n; 1 .'. n), 
ko Je n sodo celo število, ki nI manjše od 4 in 
oblikuje na ta način tklm. magični kvadrat. 
V magičnem kvadratu so vsote elementov poljubne 
vrstice, poljubnega stolpca in poljubne 
diagonale medseboj enake. 
Lista 2 prikazuje Izvajanje programa. Magični 
kvadrata 1 in 2 sta posledici uporabe procedure 
maglceven za n = 4 in za n = 6. Magična kvadra­
ta 3 in 4 nastaneta z uporabo procedure 
magicodd za n = 3 in za n = 5.; Magični kvadrat 
5 se pojavi z uporabo funkcije magioterm za n = 
5. Magični kvadrat 6 je rezultat uporabe proce­
dure maglceven za n = 16 in magični kvadrat 7 
Je rezultat uporabe procedure magicodd za n = 
17. . 
MaslCni kvadrati sodo In 
in lihe stopnje 
UITH util» 
PACKAGE BODY maaicl 18 
n« CONSTANT (»17» " 
SUBTVPE d IS inteser RANOE i .. n> 
TYFE Btolpec IS ARRAV <d) OF Inteaer« . 
TYPE POlJo IS ARRAV <d ) OF stolpec» 
-T Naslednja procedura Je predmet na<le poior-
noatl (maslAni kvadrat sodo stopnJe) 
F-ROCEIiURE masicčven < n« IN inteiiert 
m OUT pol Je > 18 
Ar b, n2> nn> i« integerl^ 
p» qf r• Boolean) - / 
FROCEIiURE alpha ( p, q. a« IN inteser« 
h« IN Boplean > 18 
hh« Boolean« , 
EEGIN -
hh «" h» , 
FOR i IN p .. q LOOP 
hh 1- NOT hh» 
.IF hh THEN 
x<iXa) •» nn - «t»n ••• 1 + n -it 
ELSE' 
x< i >(a > •- a*n - n ••• il 
END 1F> 
ENI) LOOPI 
END alpha* 
Lista 1, Procedure in programi za generiranje 
magičnih kvadratov (nadaljevanje na.drugi str.)