************* tu 3 krótkie listingi i 7 ilustracji ********* CAÎECZKI Utarî sië poglâd, ûe do obliczeï, komputerowej symulacji zjawisk fizycznych czy innych zastosowaï zwiâzanych z matematykâ lub fizykâ jedynym sîusznym komputerem jest pecet. No cóû, dyskusje nad tym, jaki komputer jest lepszy, nie prowadzâ z reguîy do niczego. Faktem jest, ûe pecet w tych zastosowaniach jest popularniejszy. Ale maluch teû jest najpopularniejszym samochodem w Polsce, choê na pewno nie najlepszym. Bolesîaw Szczerba Ûeby jasno wyraziê swoje zdanie, powiem tylko, ûe jestem studentem cholernie pecetycznej fizyki jâdrowej i jakoô ûyjë bez peceta w domu. Snobizm? Nie, wygoda. Do skomplikowanych obliczeï, wymagajâcych duûej mocy obliczeniowej, mam na uczelni RISC-a czy od biedy 486 DX2 66 MHz. Do napisania i skompilowania programu w C czy zaîadowania programu do niewielkiej obróbki danych wystarczy mi w zupeînoôci amigowy PC-Task 2.03, a raczej uruchomione pod nim odpowiednie programy. Do wiëkszej obróbki danych i 486 w domu na nic by mi sië nie przydaî, bo po pierwsze, przez duûâ obróbkë rozumiem pliki z danymi rzëdu 20 MB, a takich plików mogâ byê dziesiâtki i setki, wiëc nie przeniósîbym ich do domu; po drugie, nawet gdybym je przeniósî, to po co mam traciê na nie miejsce na twardym dysku; po trzecie, za pracë po godzinach nikt mi nie pîaci. Zostawmy jednak të dyskusjë na kiedy indziej i przejdúmy do rzeczy. Pokaûë, jak w prosty sposób uûywaê naszej Amy do obliczeï przydatnych w matematyce i fizyce. Zacznë od rzeczy najprostszej, choê bardzo czësto stosowanej w fizyce -- od liczenia caîek oznaczonych. Wbrew pozorom to nie jest taka banalna sprawa, gdyû wiele caîek nie da sië policzyê analitycznie. Takie wypadki sâ najczëstsze i liczy sië je numerycznie. Poniewaû nie wszyscy Czytelnicy zetknëli sië z caîkami w liceum, a jest to niesîychanie fajna sprawa i absolutnie dla ludzi, przypomnë pokrótce, o co chodzi w tej zabawie. Osoby znajâce rachunek caîkowy mogâ pominâê poniûszy akapit. Trochë teorii Bodajûe w trzeciej klasie liceum uczy sië o pochodnych. I ten poziom wiedzy zakîadam u Czytelnika. W razie czego moûna sië douczyê z podrëcznika do matematyki dla tejûe klasy. Caîki przychodzâ chyba w czwartej klasie, ale, o ile wiem, tylko w klasach o profilu matematyczno-fizycznym. W ogromnym uproszczeniu: caîka to taka odwrotna pochodnej. Funkcja, z której dana nam funkcja jest pochodnâ, nosi nazwë funkcji pierwotnej wzglëdem danej. Podawanie funkcji pierwotnych do danych to wîaônie caîkowanie nie oznaczone. I tak np. pochodnâ funkcji sinus jest cosinus. Oznacza to, ûe funkcjâ pierwotnâ do cosinusa jest wîaônie sinus, bo tenûe musimy zróûniczkowaê, aby otrzymaê cosinus. Mówimy, ûe caîkâ nie oznaczonâ funkcji cosinus jest funkcja sinus, ale z dokîadnoôciâ do staîej. Z dokîadnoôciâ do staîej, bo caîkâ z cosinusa moûe byê zarówno sin(x)+3, sin(x)-3.14, jak i sin(x) (=sin(x)+0). Dlaczego tak? Policz pochodne podanych trzech funkcji: wszystkie one sâ równe cos(x). Zapisujemy to tak, jak pokazano na rysunku 1. Zapis dx oznacza, ûe caîkowanie przebiega po zmiennej x. Caîkowanie oznaczone polega na liczeniu caîki nie oznaczonej w pewnych granicach. Popatrzmy na rysunek 2. To caîka nie oznaczona. Na rysunku 3. ta sama caîka, tylko w granicach od 2 do 3. Co robimy? Odgadujemy postaê funkcji pierwotnej do danej funkcji (tu chodzi o pierwotnâ do x^2, czyli x^3/3), obliczamy wartoôê funkcji pierwotnej w górnej granicy caîkowania (liczba nad znakiem caîki -- tu 3), obliczamy jej wartoôê w dolnej granicy (liczba pod znakiem caîki -- tu 2) i odejmujemy pierwsze od drugiego. Warto zaznaczyê, ûe caîka nie oznaczona jest funkcjâ, podczas gdy caîka oznaczona liczbâ. Teraz dopiero, dla tych, którzy przeûyli powyûsze proste wywody, zaczyna sië zabawa. Skâdinâd wiadomo (tego juû nie bëdë dowodziî -- jest to w kaûdej ksiâûce do analizy matematycznej), ûe caîka oznaczona z jakiejô funkcji to dokîadnie fragment pola pod krzywâ, wyznaczonâ przez të funkcjë. A to juû nie byle co. Przedstawia to rysunek 4. Czy wyobraûacie sobie, ûe w ten sposób moûna liczyê pola pod dowolnymi (no, caîkowalnymi, tzn. majâcymi funkcjë pierwotnâ) krzywymi? Moûemy policzyê, ile np. wynosi pole jednego "garbu" funkcji sinus czy cosinus -- wystarczy caîkë, podanâ na rysunku 1., policzyê w granicach -pi/2 do pi/2, i gotowe. Moûna teû wyprowadziê wzór na pole koîa (czy ktokolwiek z Was sië zastanawiaî, skâd sië wziâî wzór na pole koîa?). Co ciekawsze, granice caîek mogâ byê nieskoïczone, tzn. moûemy liczyê caîki np. od zera do plus nieskoïczonoôci. I okazuje sië, ûe pola pod niektórymi krzywymi, mimo ûe przedstawiajâ geometrycznie nieograniczone figury, majâ SKOÏCZONE pola. Trochë praktyki Podany sposób liczenia caîek oznaczonych to tzw. sposób analityczny. Wymyôlamy funkcjë pierwotnâ, liczymy jej wartoôci w granicach i mamy pole. Jednak co zrobiê, gdy dana funkcja NIE MA funkcji pierwotnej? A takich wypadków jest sporo -- przykîadem moûe byê funkcja exp(-x^2). Albo co zrobiê, gdy mamy BARDZO skomplikowanâ funkcjë i nawet nie wiemy, czy funkcja pierwotna istnieje (nie mówiâc o jej policzeniu)? Jak na zîoôê akurat z takimi caîkami mamy najczëôciej do czynienia w fizyce. I tu z pomocâ przychodzâ komputery. Pokaûë, w jaki sposób liczyê caîki oznaczone (i tym samym odpowiednie pola) numerycznie. Tu juû nie ma najmniejszego znaczenia, czy funkcja pierwotna istnieje, czy nie. Nie musimy jej W OGÓLE szukaê. Przy liczeniu numerycznym caîek oznaczonych skorzystamy z omówionego przed chwilâ zwiâzku jej wartoôci z polem pod krzywâ (rysunek 4.). Sprawa jest banalna: znamy dokîadnâ postaê funkcji, znamy a i b, cóû nam wiëcej trzeba? Metoda prostokâtów Spójrzmy na rysunek 5. Caîe pole dzielimy na N prostokâtów o równych podstawach. Pole prostokâta to chyba kaûdy umie policzyê (jak nie umie, to raczej nie doszedî do tego fragmentu...) Zakîadamy sobie jakieô N, znamy a oraz b (przy zaîoûeniu, ûe b>a -- patrz rysunek). Dîugoôê podstawy kaûdego z prostokâtów wynosi (b-a)/N, to chyba oczywiste. Natomiast wysokoôê danego prostokâta to po prostu wartoôê funkcji f(x) w dowolnym punkcie naleûâcym do jego podstawy -- dla porzâdku umówmy sië, ûe bierzemy skrajny prawy punkt. Wspóîrzëdna x pierwszego takiego punktu to a+(b-a)/N, drugiego a+2(b-a)/N, trzeciego a+3(b-a)/N itd. Przedostatni z nich ma wspóîrzëdnâ x równâ a+(N-1)(b-a)/N, a ostatni a+N(b-a)/N, czyli b. Bierze sië to z dodawania do pierwszego punktu przedziaîu, czyli punktu (a,0), kolejnych wielokrotnoôci dîugoôci podstawy prostokâta. Jedyne, co trzeba zrobiê, to policzyê wartoôci funkcji w podanych powyûej punktach. Teraz trzeba policzyê pola wszystkich prostokâtów (to juû îatwizna) i dodaê do siebie. Otrzymana liczba to wartoôê caîki oznaczonej w granicach od a do b. Jasne jak drut i proste jak sîoïce. Metoda trapezów Podana powyûej metoda prostokâtów ma jednak pewnâ wadë. Popatrzcie na rysunek. Ile wolnego miejsca zostaje miëdzy krzywâ a prostokâtami? Oszustwo to widaê tym bardziej, im bardziej stroma jest krzywa. Jest to nieekonomiczne. Po prostu mamy sakramencko wielki niedomiar lub nadmiar, który moûe nam zepsuê caîe obliczenia przy niedostatecznie gëstym podziale, gdyû moûe je obarczyê duûym bîëdem. W praktyce stosuje sië raczej innâ metodë -- metodë trapezów. Jak nietrudno sië domyôliê, polega ona na tym, ûe pole pod krzywâ jest przybliûane trapezami (patrz rysunek 6.). Wysokoôê trapezu przejmuje rolë podstawy prostokâta i równieû wynosi (b-a)/N, gdzie N jest liczbâ trapezów, na jakie dzielimy pole. Róûnica polega jedynie na tym, ûe do policzenia pola trapezu musimy znaê trzy liczby: obie podstawy i wysokoôê. Wysokoôê juû znamy, a podstawy? Teû! Pierwsza to wartoôê funkcji w skrajnym lewym punkcie podstawy trapezu, a druga to wartoôê funkcji w jej skrajnym prawym punkcie. Wspóîrzëdne punktów pozostajâ bez zmian, z tym ûe startujemy tym razem od punktu (a,0), a nie jak poprzednio od (a+(b-a)/N,0). Po zsumowaniu pól wszystkich trapezów mamy o wiele dokîadniejszâ wartoôê caîki. Daje to o wiele lepsze rezultaty i nie zaleûy juû tak silnie od stromoôci funkcji, jak metoda prostokâtów. Poza tym jest to metoda szybsza -- aby obliczyê wartoôê caîki z tâ samâ dokîadnoôciâ metodâ prostokâtów, musielibyômy zaîoûyê gëstszy podziaî (tzn. dzieliê pole na wiëcej czëôci), w zwiâzku z tym obliczenia trwaîyby dîuûej. A fizycy potrafiâ doceniê nawet niewielki przyrost prëdkoôci algorytmu, zwîaszcza gdy obliczenia idâ w dni, tygodnie, miesiâce... W praniu Pora pokazaê, co z tym wszystkim moûna zrobiê -- przedstawiam trzy listingi programów do samodzielnego wklepania. Jak zwykle zakîadam ôladowâ znajomoôê C oraz chëê programowania, poîâczonâ z podstawowymi pojëciami w tym zakresie. A wiëc przyjmujë, ûe Czytelnik wie, co to pëtle czy instrukcje warunkowe. Nie bëdë teû wyjaôniaî, po co przy instrukcji 'if' sâ nawiasy ani co robiâ dwa plusy koîo siebie. Podane listingi powinny sië daê skompilowaê dowolnym kompilatorem C, ja uûywaîem SAS/C 6.5. Dla porzâdku podajë teû zawartoôê pliku SCOptions (jednakowâ dla kaûdego listingu, choê przypominam, ûe kaûdy z nich to ODDZIELNY projekt): MATH=STANDARD ERRORREXX OPTIMIZE LINK OPTIMIZERTIME Jeôli masz lepszy procesor, moûesz przyspieszyê program, dopisujâc np. linië: CPU=68020 lub przy innym procesorze odpowiednio innâ. Pierwsze dwa programy pokazujâ, jak policzyê liczbë PI. Na rysunku 7. widzimy wykres funkcji 1/(1+x^2). Caîka oznaczona z tej funkcji w granicach od zera do jeden to dokîadnie PI/4 -- jest to policzone analitycznie. Funkcjâ pierwotnâ do danej jest tu funkcja odwrotna do tangensa -- arcus tangens. Jej wartoôê w zerze to zero, jej wartoôê w jedynce to PI/4, bo tg(PI/4) to jeden. Korzystamy z faktu, ûe pole pod krzywâ jest równe odpowiedniej caîce oznaczonej; wiemy, ûe ta caîka ma wartoôê PI/4, wiëc co nam szkodzi policzyê jâ niezaleûnie w sposób numeryczny? Wystarczy policzone numerycznie pole pomnoûyê przez 4 i mamy wartoôê PI! Na listingu 1. moûna zobaczyê, jak policzyê to pole metodâ prostokâtów, a na listingu 2. metodâ trapezów. Moûna samemu stwierdziê, która metoda jest lepsza, porównujâc czas obliczeï i ich dokîadnoôê. Programiki po skompilowaniu trzeba uruchamiaê z CLI lub Shella, podajâc jako argument gëstoôê podziaîu, czyli N. Na tyle czëôci program podzieli pole pod krzywâ przedstawionâ na rysunku 7. Na poczâtku nie radzë przekraczaê miliona, bo to moûe potrwaê. I jeszcze jedno: uwaûni Czytelnicy nie bëdâ czekaê ze stoperem w rëku na koniec obliczeï... Jak policzyê np. caîkë oznaczonâ z funkcji f(x)=x^2 w DOWOLNYM przedziale , pokazaîem na listingu 3. Przy uruchomieniu programu z CLI lub Shella trzeba pamiëtaê o podaniu trzech argumentów: podziaîu oraz granic. Specjalnie wybraîem takâ prostâ funkcjë, ûeby kaûdy mógî "sprawdziê", czy program nie oszukuje -- caîkë z podanej funkcji moûna policzyê analitycznie (patrz rysunek 3.). Funkcja f(x) jest zadeklarowana poza funkcjâ main(); nic nie stoi na przeszkodzie, aby w to miejsce wpisaê innâ funkcjë -- choêby wspomnianâ przedtem exp(-x^2). Z jej pomocâ równieû moûna policzyê liczbë PI (pytanie do studentów: jak?). W jakiej postaci podawaê funkcje, aby byîy zrozumiane przez kompilator, pokazano w pliku math.h. Moûna sië pobawiê, ale to nie jest rozwiâzanie uniwersalne. Fajnie byîoby napisaê program, który liczy caîkë z dowolnej podanej przez uûytkownika funkcji w zadanym obszarze. Teoretycznie jest to bardzo proste. Diabeî jak zwykle siedzi w szczegóîach. Przede wszystkim trzeba zaûâdaê od uûytkownika czterech argumentów: funkcji, obu granic oraz gëstoôci podziaîu. W tym wypadku mamy funkcjë jako tekst pod argv[1]. Trzeba tylko napisaê procedurkë dekodujâcâ lub, jak kto woli, "wyîawiajâcâ" funkcje ze zmiennej tekstowej. Jako lekturë uzupeîniajâcâ polecam plik "string.h". Moûna to traktowaê jako zadanie domowe. Ps. Chciaîbym w tym miejscu, korzystajâc z okazji, pozdrowiê mojego kolegë ze studiów -- Arka Wâsiïskiego, który nadal jeszcze uûywa Amigi.