---

Turbo Pascal. Programowanie Autor: Tomasz M. Sadowski Format: B5, stron: 136 Copyright © 1996 by Wydawnictwo Helion wersja spakowana

---

Rozwiązujemy dowolne równanie

Problem rozwiązania dowolnego równania nie należy bynajmniej do zadań typu kamienia filozoficznego czy uniwersalnego rozpuszczalnika. Chociaż dla pewnych klas równań istnieją dawno opracowane wzory i metody, spora grupa równań nieliniowych nie daje się rozwiązać "na papierze" w sposób dokładny. Nie oznacza to jednak, że nie można ich w ogóle rozwiązać: kwestii przybliżonego rozwiązywania równań nieliniowych poświęcono z niezłym skutkiem spory dział matematyki znany jako metody numeryczne. W niniejszym rozdziale spróbujemy zademonstrować jedną z najprostszych metod rozwiązywania równań nieliniowych - tak zwaną metodę bisekcji. Jest ona na tyle prosta, że praktycznie nie wymaga znajomości skomplikowanej matematyki (a jedynie logicznego myślenia), a jednocześnie całkiem skuteczna.

Rysunek 9. Zasada działania metody bisekcji

Załóżmy, że chcemy znaleźć miejsce zerowe funkcji f(x) przedstawionej powyżej. Szukamy więc takiego punktu c na osi x, dla którego f(c) = 0. Wybierając na osi x dwa punkty a i b takie, by wartości f(a) i f(b) miały przeciwne znaki, możemy z pewnością stwierdzić, że funkcja ma pomiędzy a i b co najmniej jedno miejsce zerowe (gdyż musi gdzieś przeciąć oś x). Podzielmy przedział pomiędzy a i b na dwie równe części. Wartość funkcji w nowo uzyskanym punkcie c₁ ma nadal znak przeciwny do f(a), a zatem miejsce zerowe leży gdzieś pomiędzy a i c₁. Podzielmy nowo uzyskany przedział na połowy: tym razem okazuje się, że wartość funkcji w punkcie c₂ ma znak przeciwny do f(c₁), a więc miejsce zerowe leży pomiędzy c₁ a c₂. Wykonując kolejne podziały dojdziemy w końcu do punktu... no, właśnie. Zauważ, że z każdym kolejnym podziałem zbliżamy się do miejsca zerowego coraz wolniej: za każdym razem przedział poszukiwania zawężany jest dwukrotnie, tak więc (teoretycznie) dokładne położenie miejsca zerowego (odpowiadające przedziałowi o zerowej szerokości, czyli pojedynczemu punktowi) osiągniemy po nieskończonej liczbie podziałów. Chyba nie masz zamiaru czekać aż tak długo!

Rozwiązaniem tego problemu jest właściwe sformułowanie tzw. kryterium zakończenia obliczeń lub kryterium stopu. W naszym przypadku kolejne cykle podziałów (zwane fachowo iteracjami) będziemy prowadzili tak długo, aż bezwzględna wartość funkcji w punkcie c_n zmaleje poniżej zadanego progu.

Spróbujmy zapisać w bardziej formalny sposób pojedynczy podział:

	początek
		wyznacz punkt c w połowie odległości pomiędzy a i b
		jeżeli f(a) ma znak różny od f(c), to przenieś punkt b do punktu c
		w przeciwnym przypadku przenieś punkt a do punktu c
	koniec

Ale jak zabrać się za wykonywanie kolejnych podziałów? Oczywiście nie będziesz musiał wielokrotnie wpisywać odpowiednich instrukcji: zamiast tego wykorzystasz do ich cyklicznego wykonywania instrukcję pętli.

Spośród trzech dostępnych w Pascalu struktur pętli do realizacji naszego zadania odpowiednie są dwie - while i repeat (trzecią pętlą, for, zajmiemy się nieco później). Należą one do grupy pętli sterowanych warunkiem, co oznacza, że wykonują się tak długo, jak długo spełnione będzie odpowiednie kryterium (while - powtarzaj, dopóki...) lub do momentu, kiedy zostanie ono spełnione (repeat - powtarzaj, aż...). Nie trzeba chyba dodawać, że owym kryterium będzie właśnie nasze kryterium stopu:

zakończ iteracje, jeśli wartość bezwzględna f(c) jest mniejsza od zadanego progu.

Której instrukcji użyć? W zasadzie w naszym przypadku jest to obojętne. Działanie pętli while i repeat jest bardzo zbliżone, zaś różnica sprowadza się do faktu, że w pierwszej z nich warunek przerwania pętli sprawdzany jest na początku:

while warunek
do instrukcja

zaś w drugiej - na końcu pętli:

repeat instrukcja
until warunek

W konsekwencji instrukcje tworzące zawartość pętli repeat muszą wykonać się co najmniej raz, zaś w przypadku pętli while mogą nie wykonać się ani razu (jeśli warunek będzie od razu spełniony). Ponadto musisz pamiętać, że sformułowanie warunku jest dla obu pętli dokładnie przeciwne (pętla while wykonuje się tak długo, jak długo warunek jest spełniony, repeat - tak długo, jak długo jest niespełniony). Zapamiętaj również, że umieszczając kilka instrukcji w pętli while musisz połączyć je w instrukcję złożoną słowami begin i end, zaś w przypadku pętli repeat konieczność ta nie występuje (ogranicznikami instrukcji złożonej są tu słowa repeat i until).

Spróbujmy zapisać nasz algorytm z wykorzystaniem pętli repeat. Będzie on wyglądał mniej więcej tak:

	początek
		wprowadź wartości krańców przedziału i dokładności
		powtarzaj
			wyznacz punkt c w połowie odległości pomiędzy a i b
			jeżeli f(a) ma znak różny od f(c), to przenieś punkt b do punktu c
			w przeciwnym przypadku przenieś punkt a do punktu c
		aż do momentu, gdy wartość bezwzględna f(c) jest mniejsza od zadanego progu
		wypisz znalezione miejsce zerowe
	koniec

Załóżmy, że chcemy rozwiązać równanie

1 - e ^{sin x * cos x}

posiadające (jak nietrudno obliczyć na kartce) pierwiastek w punkcie x = 0. Tłumacząc powyższy schemat na angielski (i Pascal) otrzymamy następujący program:

	program Bisekcja;
	{ Program rozwiązuje równania nieliniowe metodą bisekcji }
	var
		a, b, c : real;			{ granice przedziału i punkt podziału }
		eps : real;			{ dokładność }
	begin
		writeln('Program znajduje miejsce zerowe funkcji')
		writeln('w przedziale [a;b]');
		write('Podaj wartosc a: ');			{ wprowadź granice przedziału }
		readln(a);
		write('Podaj wartosc b: ');
		readln(b);
		write('Podaj dokladnosc: ');
		readln(eps);
		repeat
			c := (a + b)/2;		{ podziel przedział na pół }
			if (1 - exp(sin(a)*cos(a)))*(1 - exp(sin(c)*cos(c))) < 0
			then
				b := c	{ funkcja ma przeciwne znaki w a i c }
			else
				a := c;	{ funkcja ma przeciwne znaki w b i c }
		writeln(c);
		until abs(1 - exp(sin(c)*cos(c))) < eps;
		{ badamy wartość bezwzględną! }
		writeln('Miejsce zerowe: c = ',c:12:8);
		readln;
	end.

Nasz program dość dokładnie odpowiada schematowi podanemu powyżej i w zasadzie nie wymaga dodatkowych komentarzy (mimochodem poznałeś kilka nowych funkcji arytmetycznych - exp, sin, cos i abs). Warto jednak zwrócić uwagę na kryterium stopu:

until abs(1 - exp(sin(c)*cos(c))) < eps;

Jak wspomnieliśmy, zapis ten pozwala uniknąć nieskończenie długiego dochodzenia kolejnych podziałów do punktu stanowiącego właściwe rozwiązanie. Ale czy tylko? Spróbuj wykonać poniższy program:

	program Dodawanie;
	{ ilustracja kumulacji błędów zaokrągleń }
	{ zmiennoprzecinkowych }
	var x : real;
	begin
		x := 0.0
		repeat
			x := x + 0.1;	{ po 100 dodawaniach powinniśmy... }
			writeln(x:16:12);
		until x = 10.0;		{ ... otrzymać dokładnie 10 }
	end.

Jak łatwo sprawdzić, nasza pętla wcale nie ma zamiaru zakończyć działania po osiągnięciu wartości 10. Czyżby więc dziesięciokrotne dodanie 0.1 do 0.0 nie dawało 1.0? Istotnie: liczby rzeczywiste reprezentowane są w komputerze ze skończoną dokładnością i operacje na większości z nich obarczone są pewnym błędem, który kumuluje się w przypadku wykonywania większej ilości obliczeń. Dlatego też porównania zmiennoprzecinkowe wykorzystujące operator = (sprawdzający dokładną równość) zwykle okazują się nieskuteczne. Rozwiązaniem jest użycie operatorów >= i <=, czyli nierówności nieostrych. Poprawne kryterium zakończenia pętli miałoby u nas postać

until x >= 10.0;

Wielu początkujących programistów zapomina o niedoskonałości komputera, próbując przełożyć "idealny" zapis matematyczny wprost na język programowania. Efekty mogą być niezbyt przyjemne, zwłaszcza w przypadku wykorzystania takich zapisów w pętlach (dodatkowo występuje wówczas kumulacja błędów) i instrukcjach warunkowych.

Programując operacje na liczbach rzeczywistych należy zawsze pamiętać o ich niedokładnej reprezentacji.

Omówione wyżej pętle while i repeat należą do grupy pętli sterowanych warunkiem, co oznacza, że liczba przebiegów pętli nie jest z góry określona. Typowym zastosowaniem pętli while i repeat jest wykonywanie obliczeń (lub innych działań) aż do spełnienia określonego kryterium, wprowadzanie danych z klawiatury lub pliku i wszelkie inne operacje, w których nie jesteśmy w stanie wcześniej ustalić, jak wiele cykli trzeba będzie wykonać. Czasem jednak okazuje się, że musimy zrealizować ściśle określoną liczbę operacji, np. wypisać na ekranie wszystkie litery od "A" do "Z" wraz z ich kodami ASCII. Oczywiście i to można zrobić za pomocą pętli while lub repeat, jednak znacznie poręczniejsza okazuje się tu pętla for, sterowana nie warunkiem, lecz licznikiem. A oto i sam program:

	program ASCII;
	{ wyprowadza znaki ASCII od 'A' do 'Z' i ich kody }
	var c : char;			{ wypisywany znak }
	begin
		for c := 'A' to 'Z' do
			write(c:4, ord(c):4)
	end.

Jak widać, składnia pętli for jest następująca:

	for licznik := początek to koniec do
		instrukcja

co tłumaczy się "na nasze" jako "dla licznika przebiegającego wartości od początek do koniec wykonuj instrukcję". Pętlę "w dół" (wartość końcowa licznika jest mniejsza od początkowej) zapisuje się nieco inaczej:

	for licznik := początek downto koniec do
		instrukcja

Warto o tym pamiętać, gdyż użycie do tego celu poprzedniej postaci spowoduje, że pętla nie wykona się ani razu. Jeśli wreszcie wartość początkowa licznika jest równa końcowej, pętla wykona się dokładnie jeden raz (co ma niewielki sens) niezależnie od tego, której z powyższych form użyjemy.

Programując pętlę for musisz dokładnie wiedzieć, ile razy będzie się ona wykonywała, czyli znać wartość początkową i końcową licznika. Nie jest oczywiście powiedziane, że wartości te musisz podać w postaci stałych: równie dobrze możesz użyć do tego celu zmiennych lub wyrażeń, byle tylko dawały one w wyniku wartości typu porządkowego, np:

for i := 1 to MaxIndex div 2 do { ... }

Sam licznik również musi być zmienną typu porządkowego (zwykle integer lub char) i w trakcie wykonywania pętli zmienia się z krokiem 1 od wartości początek do koniec włącznie. Zaprogramowanie w Pascalu pętli z krokiem różnym od 1 (w szczególności niecałkowitym) jest nieco bardziej skomplikowane; zwykle lepiej w tym celu użyć pętli while lub repeat, chociaż można to zrobić tak:

	for i := 1 to 20 do
		begin
			writeln(i);
			Inc(i)	{ dodaj 1 do licznika }
		end;

Powyższa konstrukcja wypisze na ekranie wszystkie liczby nieparzyste zawarte między 0 a 20. Sztuczka, której tu użyliśmy, polega na zmianie licznika w trakcie wykonywania pętli (normalnie zajmuje się tym sama pętla, a my nie powinniśmy się do tego wtrącać) i - chociaż legalna i poprawna - jest niezbyt bezpieczna. Spróbuj np. zmienić wartość końcową licznika na 21 (program "przeskoczy" zakończenie pętli) lub instrukcję Inc(i) na Dec(i) (zmniejszenie licznika o 1 zamiast zwiększenia - pętla "zapętli się" na wartości 1).

Programując pętlę for należy stanowczo unikać jawnego modyfikowania licznika wewnątrz pętli. Aby zaprogramować pętlę z krokiem różnym od 1, lepiej użyć konstrukcji while lub repeat.

Wzmiankowane wyżej instrukcje Inc i Dec pozwalają odpowiednio na zwiększenie (inkrementację) lub zmniejszenie (dekrementację) wartości argumentu typu porządkowego o zadaną wartość (również typu porządkowego) i odpowiadają instrukcjom dodania lub odjęcia wartości, np.

Dec(k, 30); { to samo, co k := k - 30; }

Na tym zakończymy rozważania na temat pętli w Pascalu (i kilku drobnych dodatków), choć do zastosowań pętli będziemy wracać jeszcze wielokrotnie. W kolejnym rozdziale zajmiemy się pewną niedogodnością programu Bisekcja, polegającą na kłopotliwym zadawaniu postaci rozwiązywanego równania. Zauważ, że nasza funkcja nieliniowa

f(x) = 1 - e ^{sin x * cos x}

jest w programie wykorzystywana dwa razy: raz podczas właściwej iteracji, drugi raz w kryterium stopu. Jest to niewygodne i nieeleganckie, gdyż chcąc znaleźć miejsce zerowe innej funkcji musimy wpisywać ją dwa razy. Jednocześnie sam algorytm pozwala nam znaleźć miejsce zerowe dowolnej funkcji f(x). Czy zatem nie dałoby się zamknąć naszej funkcji w jakieś pudełko, do którego wkładalibyśmy tylko wartość argumentu i odbierali wynik? Oczywiście, że tak.

W ramach ćwiczeń praktycznych proponuję Ci następujące zadania:

• zmodyfikuj program Bisekcja tak, by wykorzystywał pętlę while;
• napisz program wypisujący wartości funkcji f(x) = 1 - e ^{sin x * cos x}
w przedziale od x = 0 do x = 360 stopni (2) z krokiem 10 stopni. Pamiętaj, że argumenty funkcji trygonometrycznych podawane są w radianach. Spróbuj wykorzystać wszystkie poznane pętle.

Zapamiętaj

• Do cyklicznego (iteracyjnego) wykonywania instrukcji służą w Pascalu pętle.
• Instrukcje pętli mogą być sterowane warunkiem (while-do, repeat-until) lub licznikiem (for-to/downto).
• Pętle while i repeat używane są wówczas, gdy nie znamy z góry liczby przebiegów pętli, możemy natomiast określić warunek jej zakończenia.
• Programując kryterium zakończenia pętli while i repeat (i nie tylko) należy pamiętać o niedokładności reprezentacji liczb rzeczywistych i kumulacji błędów obliczeń.
• Pętla for wykorzystywana jest w sytuacjach, gdy możemy dokładnie określić liczbę przebiegów.
• Licznik pętli for jest zawsze typu porządkowego i zmienia się z krokiem 1.
• Modyfikacja licznika wewnątrz pętli for jest dopuszczalna, ale może prowadzić do trudnych do wykrycia błędów i należy jej unikać.
• Do inkrementacji lub dekrementacji zmiennych typu porządkowego można wykorzystać instrukcje Inc i Dec.

Poprzedni | Spis treści | Następny | Wersja spakowana |

---