( 1 ) /(*) = 0 .

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXXV (1992)

M. Sz y s z k o w i c z

Ottawa

Wyznaczanie zera funkcji

( Praca wpłynęła do Redakcji 2Ą.03.1992)

1. Wstęp Dane jest równanie nieliniowe postaci

( 1 ^{) /(*) =} 0 ^.

Liczba r, która spełnia powyższe równanie jest nazywana pierwiastkiem rów- nania bądź zerem funkcji / . Zadanie polega na znalezieniu wartości r. Rów- nanie może nie mieć rozwiązań rzeczywistych, jak na przykład równanie

sin(ar) — 3 = 0

lub może mieć nieskończenie wiele rozwiązań, jak równanie sin(ai) — 1 = 0.

Znane są wzory na pierwiastki równania (1) w sytuacji, gdy funkcja / jest wielomianem drugiego, trzeciego lub czwartego stopnia. Około roku 1830 Galois wykazał, że nie da się uzyskać podobnych wzorów dla wielomianów stopnia piątego oraz wyższych. W ogólnym przypadku nie istnieją formuły, które określają rozwiązanie równania nieliniowego. Trzeba posłużyć się me- todami iteracyjnymi i komputerem, celem otrzymania zera ciągłej funkcji rzeczywistej.

Przez „rozwiązanie” równania (1) na maszynie cyfrowej rozumie się otrzymanie takiej liczby r*, dla której zachodzą warunki:

|/(r*)| « 0 lub |r — r*| « 0.

Symbol « zwykle zależy od dokładności użytej arytmetyki. Liczba r* pra- wie spełnia równanie bądź liczba r* leży blisko prawdziwego rozwiązania r.

Istnieje wiele metod i algorytmów numerycznych (patrz [1], [4], [5]), które wyznaczają wartość r*.

34 M. Szyszkowicz

Celem tej pracy jest przedstawienie efektywnej metody wyznaczania jed- nego zera funkcji. Podany jest algorytm oraz jego realizacja w postaci proce- dury o nazwie BISEC. Zakłada się, że znany jest przedział [a, 6], na którym funkcja / jest ciągła i zmienia znak, tj. f{a)f(b) < 0. Zero bądź zera leżą w przedziale [a, 6]. Fakt zmiany znaku przez funkcję, daje możliwość genero- wania ciągu przybliżeń, który dąży do rozwiązania. Przedstawiona metoda jest zawsze zbieżna i korzysta jedynie z wartości funkcji.

Najprostszą metodą numeryczną wyznaczania rozwiązania równania (1), jest metoda połowienia przedziału. Metoda ta znana jest również jako me-

toda bisekcji lub metoda Bolzano. Przedział [a, 6] dzieli się na pół i wybiera się tę połowę, na której funkcja zmienia znak. Tak wybrana połowa jest nowym przedziałem [a, 6]. Proces ten się kontynuuje, aż do spełnienia przy- jętego kryterium zatrzymania.

Drugą, dobrze znaną metodą wykorzystującą zmianę znaku funkcji na przedziale [a, 6], jest metoda reguła falsi. Metoda ta opiera się na interpolacji linowej, zbudowanej na węzłach (a, /(a )) oraz (6, /(&)). Punkt przecięcia, tak wyznaczonej prostej z osią x-ów, jest jednym z końców nowego przedziału [a, 6]. Stale wybiera się ten przedział, na którym f ( a) f ( b ) < 0. Istnieje kilka modyfikacji metody reguła falsi [1], jak również kompilacji tej metody z metodą bisekcji. Zwykle zmodyfikowane wersje reguła falsi lub hybrydowe połączenia działają lepiej niż metoda oryginalna [3].

W tej pracy wykorzystuje się metodę połowienia przedziału, reguła falsi oraz metodę siecznych. Metoda siecznych jest realizowana jedynie w przy- padku, gdy wyznaczone nią przybliżenie nie wyprowadza poza przedział [a, 6].

2. Algorytm W każdym etapie algorytmu znane są trzy punkty a, b oraz c, dla których zachodzi

f( a) f ( b) <

0

^, /(a )/(c ) >

0

^{, /(c ) /}

0

^.

Pierwszy warunek orzeka o zmianie znaku na przedziale [a, 6]. Drugi mówi o tym, że a oraz c leżą po tej samej stronie pierwiastka r. Podobnie postępuje się, gdy /(6 ) /( c ) > 0. Za pierwszym razem, wartość c jest wyznaczona z połowienia przedziału [a, 6]. Wykonuje się liniową interpolację z węzłami (c ,/( c ) ) oraz (6, /(&)) i otrzymuje punkt p z metody reguła falsi. Dokonuje się liniowej ekstrapolacji z węzłami (a ,/(a )) oraz (c ,/(c )) i wyznacza się punkt q z metody siecznych. Nowe przybliżenie d dla pierwiastka r, oblicza się następująco

d = ( p+q) /2, gdy q G (c, 6).

Jeśli punkt q leży poza przedziałem (c, 6), to wartość d oblicza się jako środek przedziału [c,6]. W następnym kroku przedział [a, 6] jest zastąpiony przez [c,6]. Nowy punkt c jest określony przez wartość d. Wyznacza się nowe p

Wyznaczanie zera funkcji 35 oraz q i proces się powtarza. W sytuacji, gdy wyprowadzone przybliżenie d wypadło z przedziału (a, 6) bądź /(c ) = 0, to wykonanie algorytmu się przerywa. Za rozwiązanie bierze się ostatnią obliczoną wartość c.

Algorytm można krótko scharakteryzować, jako metodę polegającą na generowaniu średnich arytmetycznych z przybliżeń zera metodą siecznych i reguła falsi. Autor tej pracy po opracowaniu swojego algorytmu stwier- dził, że podobne postępowanie zostało zaproponowane w [2]. Obie metody asymptotycznie zachowują się tak samo i mają ten sam rząd zbieżności 1.839.

3. Procedura BISEC Algorytm został zrealizowany w języku Fortran w postaci procedury (real function) BISEC. Procedura BISEC była testo- wana dla wielu równań nieliniowych. Tabela 1 podaje otrzymane wyniki dla niektórych zadań testowych z [3].

Tabela 1. Liczba obliczeń wartości funkcji / celem wyznaczenia pierwia- stka równania f ( x) = 0 przez procedurę BISEC.

Funkcja f ( x) = (x — l)(z — 2)(x — 3)(z — A){x — 5)(x — 6) Przedział Liczba obliczeń /

[0,1.5] 10

[1.1.2.5] 9

[2.1.3.5] 9

[3.1.4.5] 8

[4.1.5.5] 8

[5.1.6.5] 9

Funkcja f(x) = (nx — l ) / ( ( n — l)x) na przedziale [.01,1].

(wartość n, liczba obliczeń / ) : (2, 7), (5,11), (15,12), (20,13).

Funkcja f(x) = x n na przedziale [—1,10].

(wartość n, liczba obliczeń / ) : (3,116), (5,100), (7,75), (9,49), (19,19), (25,13).

Liczba c (rozwiązanie) zawarta jest w zredukowanym przedziale (a, 6), gdzie f{a)f{b) < 0, c 6 (a, b) oraz /(c ) « 0. Punkty a, b i c są sąsiednimi licz- bami zmiennopozycyjnymi bądź wartości funkcji na przedziale (a, 6) są tak małe, że algorytm powtarza przybliżenia lub wyprowadza poza (a, b). Skoń- czona długość mantysy gwarantuje zatrzymanie. Procedura BISEC również zakańcza działanie, jeśli f ( c) = 0. W opracowaniu procedury zwrócono baczną uwagę na poprawne zaprogramowanie poszczególnych etapów algo- rytmu. Metoda bisekcji jest wykonana instrukcją

D = A + ( B - A ) * 0.5, a nie, jak często się zdarza widzieć

D = (A + B)/

2.

36 M. Szyszkowicz

Matematycznie instrukcje są równoważne, jednak w arytmetyce zmiennopo- zycyjnej jest między nimi różnica. Warunek f(a)f(b) < 0 testuje się badając wyrażenie

SFB * F(A) <

0

^,

gdzie SFB ma wartość 1 lub —1 zgodnie ze znakiem /(&). Wartość SFB wylicza się raz (jedno wywołanie funkcji SIGN). Nie jest poprawne badanie

F(A) * F(B) <

0

^.

Ten test zawodzi dość szybko, i tak na przykład dla f ( A) = 10-30 oraz f ( B) = —

10-30

daje wartość —

10 ^~60,

co może być na niektórych kompute- rach zerem. Test na użycie metody siecznych, przeprowadza się bez wyzna- czania punktu q oraz bez obliczania modułów wartości funkcji / . Porównuje się dwie wielkości, w procedurze Z oraz ZC. Pozwala to także uniknąć dzie- lenia przez zero w trakcie realizacji wzoru metody siecznych. Jeśli w trzech kolejnych krokach (ilość ustalona doświadczalnie), aproksymacje‘występują po tej samej stronie pierwiastka, to w następnym, czwartym kroku stosuje się metodę bisekcji. Działanie algorytmu przenosi się do innej części przedziału [a, 6], gdzie algorytm podejmuje nową akcję.

Autor dziękuje recenzentowi za konstruktywne uwagi.

Literatura.

1. A. B jó r c k , G. D a h lq u is t, Metody numeryczne, Warszawa, PW N 1983.

2. D. Le, Three new rapidly convergent algorithms for finding a zero of a function, SIAM J. Sci. Stat. Comput. 2, (1985), 193-208.

3. D. N e r in c k x , A . H a e g e m a n s, A comparision of non-linear equation solvers, J.

Comput. Appl. Math. 2 (1976), 145-148.

4. A . R a ls t o n , W stęp do analizy numerycznej, Warszawa, PW N 1971.

5. J.F. T r a u b , Iterative methods of the solution of equations, Englewood Cliffs, N.J., 1964.