plik


ÿþWprowadzenie do Scilaba wersja 0.9999 ± Bruno Pinçon Institut Elie Cartan Nancy E.S.I.A.L. Université Henri Poincaré Email : Bruno.Pincon@iecn.u-nancy.fr PrzekBad z jzyka francuskiego : Piotr FulmaDski Katarzyna Szulc WydziaB Matematyki i Informatyki Institut Elie Cartan Nancy 1 Uniwerstytetu Aódzkiego Université Henri Poincaré Email : fulmanp@math.uni.lodz.pl Email : Katarzyna.Szulc@iecn.u-nancy.fr Skrypt ten byB pocztkowo opracowany przez studentów E.S.I.A.L. (Ecole Supérieure d Informatique et Application de Lorraine). Opisuje on niewielk cz[ mo|liwo[ci Scilaba, w szczególno[ci te, które pozwalaj zastosowa notacje analizy numerycznej i maBych symulacji stochastycznych takich jak: " operacje na macierzach i wektorach o wspóBrzdnych rzeczywistych; " programowanie w Scilabie; " prosta grafika; " niektóre funkcje dla dwóch wymienionych powy|ej (generowanie liczb losowych, rozwizywanie równaD, ...) Scilab umo|liwia wykonanie wielu innych operacji, w szczególno[ci w dziedzinie automatyki, obróbki sygnaBów dzwiko- wych, symulacji systemów dynamicznych (przy pomocy scicos)... Jako|e zamierzam systematycznie uzupeBnia ten dokument, jestem w peBni otwarty na wszelkie uwagi, sugestie i krytyk pozwalajce mi na jego uleprzenie (równie| w przypadku bBdów ortograficznych), piszcie do mnie1. MaBa historia tego dokumentu: " wersja 0.999: modyfikacje rozdziaBu o grafice oraz uzupeBnienie tre[ci rozdziaBu o programowaniu; wersja relatywna do Scilab-2.7 ; " wersja 0.9999 (ten dokument): dostosowanie rozdziaBu o grafice do Dowej grafigi obiektowejZcilaba; wersja relatywna do Scilab-4.0. W wyniku dopisania kilku paragrafów, dokument straciB na spójno[ci. Istnieje jednak obecnie inny podrcznik, ktory jest dostpny na stronie internetowej Scilaba (czytaj dalej). 1 T sam uwage chca przekaza czytelnikom autorzy tBumaczenia polskiego (przyp. tBum.) Podzikowania " dla Doc Scilab, który czsto pomagaB mi na forum u|ytkowiników; " dla Bertrand Guiheneuf, który dostarczyB mi magiczn "[cieszk"do skompilowania Scilaba 2.3.1 pod linuxem (kompilacja pod linuxem nowszych wersji nie powoduje problemów); " dla moich kolegów i przyjacióB, Stéphane Mottelet2 , Antoine Grall, Christine Bernier-Katzentsev oraz Didier Schmitt ; " dla Patrice Moreaux za uwa|ne przeczytanie i poprawienie bBdów; " dla Helmut Jarausch, który przetBumaczyB ten dokument na jzyk niemiecki i który zwróciB moj uwag na kilka niejasno- [ci; " i dla wszystkich czytelników za ich wsparcie, uwagi i korekty. 2 dziki za ten .pdf Stéphane! 2 Spis tre[ci 1 Wiadomosci wstepne 2 1.1 Co to jest Scilab? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Jak korzysta z tego dokumentu? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Podstawy pracy z Scilab-ie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Gdzie znalez informacje na temat Scilab-a? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.5 Jaki jest statut programu Scilab? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Operacje na macierzach i wektorach 4 2.1 Wprowadzanie macierzy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Typowe wektory i macierze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Wyra|enia w Scilab-ie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3.1 Kilka podstawowych przykBadów wyra|eD macierzowych . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3.2 DziaBania na elementach macierzy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.3 Rozwizywanie ukBadów równaD liniowych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.4 Indeksowanie, wydobywanie podmacierzy, konkatenacji macierzy i wektorów . . . . . . . . . . . . . . 11 2.4 Informacje na temat [rodowiska pracy(*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 WywoBanie pomocy z linii poleceD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6 Generator prostych wykresów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.7 Pisanie i wykonywanie skryptów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.8 Dodatkowe informacje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.8.1 Skracanie instrukcji przy zapisie macierzowym . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.8.2 PozostaBe uwagi dotyczce rozwizywania ukBadów równaD liniowych (*) . . . . . . . . . . . . . . . . . 16 2.8.3 Kilka prostych macierzy (*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8.4 Funkcjesizeilength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 wiczenia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 Programowanie w Scilabie 24 3.1 Ptle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.1 Ptlafor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.2 Ptlawhile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Instrukcja warunkowa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.1 Instrukcjaif-then-else . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.2 Instrukcjaselect-case(*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 Inne rodzaje zmiennych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.1 AaDcuchy zanków . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.2 Listy (*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.3 Kilka przykBadów z wektorami i macierzami logicznymi . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Funkcje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4.1 Przekazywanie parametrów (*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4.2 Wstrzymywanie funkcji . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4 Grafika 36 4.1 Okna graficzne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Wprowadzenie doplot2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3 plot2dz argumentami opcjonalnymi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.4 Inne wersjeplot2d:plot2d2,plot2d3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5 Rysowanie wikszej ilo[ci krzywych zBo|onych z ró|nej ilo[ci punktów . . . . . . . . . . . . . . . . . . . . . . 43 4.6 Zabawa z kontekstem graficznym . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.7 Tworzenie histogramów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.8 Zapisywanie grafiki w ró|nych formatach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.9 Prosta animacja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.10 Powierzchnie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.10.1 Wprowadzenie doplot3d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.10.2 Kolory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.10.3 plot3dz des facettes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.10.4 Rysowanie powierzchni opisanej przy pomocy x = x(u, v), y = y(u, v), z = z(u, v) . . . . . . . . . . . 53 4.10.5 plot3dz interpolacj kolorów . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.11 Krzywe w przestrzeni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 1 4.12 Ró|no[ci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 Zastosowania i uzupeBnienia 58 5.1 Równania ró|niczkowe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1.1 Podstawowe u|ycieode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1.2 Van der Pol jeszcze raz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.1.3 Troche wicej oode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Generowniw liczb losowych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2.1 Funkcjarand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2.2 Funkcjagrand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.3 Dystrubuanty i ich odwrotno[ci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 Proste symulacje stochastyczne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4.1 Wprowadzenie i notacja . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4.2 PrzedziaBy ufno[ci . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4.3 Wykres dystrubuanty empirycznej . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.4.4 Test Ç2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.4.5 Test KoBmogorowa-Smirnova . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.4.6 wiczenia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 6 Zbiór drobiazgów 72 6.1 Definiowanie wektora i macierzy wspóBczynnik po wspóBczynniku . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.2 Na temat warto[ci zwracanych przez funkcj . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.3 Funkcja zostaBa zmidyfikowana ale... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.4 Problem zrand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.5 Wektory wierszowe, wektory kolumnowe... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.6 Operator porównania . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.7 Liczby zespolone a liczby rzeczywiste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.8 Proste instrukcje a funkcje Scilaba . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.9 Obliczanie wyra|eD logicznych . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 A Odpowiedzi do wiczeD z rozdziaBu 2 75 B Rozwizania wiczeD z rozdziaBu 3 76 C Rozwizania wiczeD z rozdziaBu 4 79 2 1 Wiadomosci wstepne 1.1 Co to jest Scilab? Dla osób znajcych ju| MATLAB-a odpowiedz jest prosta  Scilab jest jego darmowym (wicej szczeguBów zwizanych z tym tematem w dalszej cz[ci) odpowiednikiem, powstaBym3 w I.N.R.I.A. (Institut National de Recherche en Informatique et Auto- matique). SkBadnia, z wyjtkiem nielicznych poleceD, jest taka sama (istotne ró|nice wystepuj w przypadku grafiki). Osobom nie znajcym MATLAB-a chc powiedzie, |e Scilab jest przystepnym [rodowiskiem do wykonywania obliczeD numerycznych dziki zaimplementowanym niezbednym metodom: " rozwiazywanie ukladów liniowych, " wyznaczanie warto[ci wBasnych, wektorów wBasnych, " dekompozycja dla warto[ci osobliwych i pseBdo-osobliwych, " szybka transformacja Fouriera, " rozwiazywanie równaD ró|niczkowych, " algorytmy optymalizacji, " rozwiazywanie równan nieliniowych, " generowanie liczb losowych, " dla wielu niezaawansowanych zastosowaD algebry liniowej w automatyce. Ponadto Scilab wyposa|ony jest w funkcje sBu|ce do tworzenia grafiki, zarówno niskopoziomowej (wielokty, odczytywanie wspóBrzdnych poBo|enia kursora, itp.) jak i wysokopoziomowej (krzywe, powierzchnie itp.). Wprowadzony jzyk programo- wania, dziki operowaniu notacj macierzow, jest prostym, ale pot|nym i efektywnym narzdziem. Wyszukiwanie bBdów w programach jest szybkie dziki Batwemu operowaniu zmiennymi. W przypadku, gdy obliczenia bd zbyt czasochBonne (jzyk jest interpretowany. . . ) mo|liwe jest Batwe poBczenie programu Scilab-a z podprogramami napisanymi w C czy FORTRAN-ie. 1.2 Jak korzysta z tego dokumentu? W rozdziale drugim wyja[niono podstawy pracy z Scilab-em jako narzdziem do obliczeD macierzowych. Wystarczy prze[ledzi proponowane przykBady. Sekcje oznaczone gwiazdk podczas pierwszego czytania mo|na pomin. Osoby zainteresowane zagadnieniami dotyczcymi grafiki, moga przeanalizowac pierwsze przykBady z rozdziaBu czwartego. RozdziaB trzeci wyja[nia podstawy programowania w Scilab-ie. RozdziaB szósty pokazuje jak unikn najczstrzych bBdów popeBnianych podczas pracy w Scilab-ie lub nieporozumieD wynikajcych z jego specyfikacji (prosz pode[lij mi tak|e i swoje!). Istniej nieznaczne ró|nice pomidzy [rodowiskami graficznymi przeznaczonymi dla Unix i Windows polegajce na odmiennym sposobie rozmieszczenia przycisków i menu. W tym dokumencie opieram sie na wersji Unix aczkolwiek u|ytkownicy wersji przeznaczonej dla systemu Windows nie powinni natrafi na problemy ze znalezieniem odpowiednich opcji / przyciskow. 1.3 Podstawy pracy z Scilab-ie W najprostszym przypadku, Scilab mo|e by wykorzystywany jako  kalkulator zdolny wykonywa obliczenia na wektorach i macierzach liczb rzeczywistych i/lub zespolonych (ale tak|e na zwykBych skalarach) oraz do wizualizacji krzywych i po- wierzchni. W najprostrzych zadaniech raczej nie ma potrzeby pisania programów. Dosy szybko zaczniemy jednak korzysta ze skryptów (zbiorów instrukcji, poleceD Scilab-a) a nastpnie funkcji. Oczywi[cie niezbdne w takich sytuacjach staje sie u|ycie edytora tekstu, na przykBad emacs (Unix, Windows), wordpad (Windows), vi lub vim (Unix a tak|e Windows)... 1.4 Gdzie znalez informacje na temat Scilab-a? W dlalszej cz[ci dokumentu zakBda si, |e czytelnik dysponuje wersj 2.7 programu. W celu uzyskania dodatkowych informacji odsyBam do  Scilab home page : http://scilabsoft.inria.fr na której znalez mo|na ró|norodn dokumentacj, oraz efekty pracy innych u|ytkowników.  Scilab group wydaB (w latach 1999  2001) okoBo dwudziestu artykuBów w czasopi[mie  Linux magazie . Polecam je szczególnie, gdy| poruszaj wikszo[ zagadnieD zwizanych ze Scilab-em (wikszo[ci z nich w tym opracowaniu nawet si nie porusza). Uzyska je mo|na pod adresem 3 Scilab wyko|ystuje du| ilo[ funkcji pochodzcych z ró|nych miejsc, dostpnych czsto przez Netlib 3 http://www.saphir-control.fr/articles/ Na temat Scilab-a prowadzone jest równie| forum dyskusyjne, w ramach którego istnieje mo|liwo[ zadawania pytaD, doko- nywania uwag, udzielania odpowiedzi na wcze[niej postawione pytania, itd: comp.sys.math.scilab Wszystkie zamieszczone tam wiadomo[ci s archiwizowane i dostpne ze strony domowej Scilab-a po wybraniu Scilab newsgroup archive. Wybierajc jeden z dwóch odsyBaczy Books and Articles on Scilab lub Scilab Related Links umieszczonych na stronie gBównej uzyskujemy dostp do wielu ró|nch dokumentów. W szczególno[ci: " wprowadzenie B. Ycart (Démarrer en Scilab et statistiques en Scilab); "  Scilab Bag Of Tricks autorstwa Lydia E. van Dijk i Christoph L. Spiel, który jest raczej przeznaczony dla osób dobrze znajcych ju| Scilab-a (rozwój tej ksi|ki zostaB niestety brutalnie przerwany kilka lat temu); " Travaux Pratiques sur Scilab classes par themes umo|liwia dostp do projektów realizowanych przez studentów ENPC; " wprowadzenie do informatyki z u|ciem Scilab-a (verb+http://kiwi.emse.fr/SCILAB/+). Oczywi[cie w zale|no[ci od potrzeb znalez mo|na du|o innych opracowaD traktujcych problem w nieco odmienny ni| za- mieszczony tutaj sposób. 1.5 Jaki jest statut programu Scilab? Osoby dobrze znajce oprogramowanie (na licencji GPL) z pewno[ci interesuje statut Scilab-a jako programu darmowego. Czstym pytaniem na forum jest : Scilab: is it really free? Yes it is. Scilab is not distributed under GPL or other standard free software copyrights (because of historical reasons), but Scilab is an Open Source Software and is free for academic and industrial use, without any restrictions. There are of course the usual restrictions concerning its redistribution; the only specific requirement is that we ask Scilab users to send us a notice (email is enough). For more details see Notice.ps or Notice.tex in the Scilab package. Answers to two frequently asked questions: Yes, Scilab can be included a commercial package (provided proper copyright notice is included). Yes, Scilab can be placed on commercial CD s (such as various Linux distributions). 2 Operacje na macierzach i wektorach Ta cz[ daje podstawy do poznania zastosowaD Scilab-a jako narzdzia do operacji macierzowych. Aby rozpocz prac z Scilab-em wystarczy wpisa scilab w terminalu4. Je[li wszystko przebiega prawidBowo, na ekranie uka|e si okno Scilab-a z gBównym menu zawierajcym w szczególno[ci przysciskHelp,Demosa w oknie wprowadzania poleceD uka|e si ========== scilab-2.7 Copyright (C) 1989-2003 INRIA/ENPC ========== Startup execution: loading initial environment --> gdzie-->jest znakiem zachty. 4 Albo klikn odpowiedniaÛ ikon. 4 2.1 Wprowadzanie macierzy Podstawowym typem danych w Scilab-ie jest macierz liczb rzeczywistych lub zespolonych. Najprostszym sposobem definiowa- nia macierzy (wektora, skalara bdcych w istocie szczególnymi przypadkami macierzy) w [rodowisku Scilab jest wprowadzenie z klawiatury listy elementów macierzy, stosujc nastpujc konwencj: " elementy tego samego wiersza oddzielone s spacj lub przecinkiem; " lista elementów musi by ujta w nawias kwadratowy []; " w przypadku wprowadzania skalarów nawias kwadratowy mo|na pomin; " ka|dy wiersz macierzy, z wyjtkiem ostatniego, musi by zakoDczony [rednikiem. Dla przykBadu komenda: -->A=[1 1 1;2 4 8;3 9 27] da na wyj[ciu A = ! 1. 1. 1. ! ! 2. 4. 8. ! ! 3. 9. 27. ! W przypadku, gdy instrukcja zostanie zakoDczona [rednikiem, wynik nie pojawi si na ekranie. Wpiszmy na przykBad -->b=[2 10 44 190]; Aby zobaczy wspóBrzdne wprowadzonego wektora, wystarczy napisa -->b a odpowiedzi bdzie b = ! 2. 10. 44. 190. ! Bardzo dBuga instrukcja mo|e by napisana w kilku linijkach, przy czym przechodzc do linii nastpnej, lini poprzedni nale|y zakoDczy trzema kropkami, jak w poni|szym przykBadzie: -->T = [ 1 0 0 0 0 0 ;... --> 1 2 0 0 0 0 ;... --> 1 2 3 0 0 0 ;... --> 1 2 3 0 0 0 ;... --> 1 2 3 4 0 0 ;... --> 1 2 3 4 5 0 ;... --> 1 2 3 4 5 6 ] co daje T = ! 1. 0. 0. 0. 0. 0. ! ! 1. 2. 0. 0. 0. 0. ! ! 1. 2. 3. 0. 0. 0. ! ! 1. 2. 3. 4. 0. 0. ! ! 1. 2. 3. 4. 5. 0. ! ! 1. 2. 3. 4. 5. 6. ! W przypadku wprowadzania liczb zespolonych stosuje si nastpujc skBadnie: -->c=1 + 2*%i c = 1. + 2.i -->Y = [ 1 + %i , -2 + 3*%i ; -1 , %i] Y = ! 1. + i - 2. + 3.i ! ! - 1. i ! 5 2.2 Typowe wektory i macierze Istniej funkcj do konstrukcji typowych macierzy i wektorów. Macierz jednostkowa. Aby otrzyma macierz jednostkow o wymiarach 4 na 4 stosujemy instrukcj: -->I=eye(4,4) I = ! 1. 0. 0. 0. ! ! 0. 1. 0. 0. ! ! 0. 0. 1. 0. ! ! 0. 0. 0. 1. ! Argumentami funkcjieye(n,m)jest liczba wierszy (n) oraz kolumn (m) (Uwaga : je|eli n < m (odpowiednio n > m wówczas otrzymamy macierz odwzorowania wzajemnie jednoznacznego, surjekcja (odpowiednio injekcja) przestrzeni kanonicznej Km na Kn). Macierz diagonalna, diagonalna macierzy. Aby otrzyma macierz diagonaln, w której elementy na gBównej przektnej pochodz z wcze[niej zdefiniowanego wektorabwpisujemy -->B=diag(b) B = ! 2. 0. 0. 0. ! ! 0. 10. 0. 0. ! ! 0. 0. 44. 0. ! ! 0. 0. 0. 190. ! Uwaga: Piszcbnadal mamy dostp do wcze[niej zdefiniowanego wektora, co pokazuje, |e Scilab rozró|nia wielko[ liter. Zastosowanie funkcjidiagna dowolnej macierzy pozwala uzyska gBówn przektn macierzy jako wektor kolumnowy -->b=diag(B) b = ! 2. ! ! 10. ! ! 44. ! ! 190. ! (Funkcja ta przyjmuje tak|e drugi, opcjonalny, argument  porównaj wiczenia.) Macierze zerowa i jedynkowa. Funkcje zeros i ones pozwalaj odpowiednio stworzy macierze zerowe i macierze skBadajce si z jedynek. Podobnie jak dla funkcjieyeich argumentami s liczba wierszy i liczba kolumn. -->C = ones(3,4) C = ! 1. 1. 1. 1. ! ! 1. 1. 1. 1. ! ! 1. 1. 1. 1. ! Mo|na tak|e u|y jako argument nazw macierzy ju| zdefiniowanej w [rodowisku. W efekcie otrzymujemy macierz o takich wymiarach macierzy bdcej argumentem -->O = zeros(C) O = ! 0. 0. 0. 0. ! ! 0. 0. 0. 0. ! ! 0. 0. 0. 0. ! Macierz trójktna. Funkcjetriuitrilpozwalaj otrzyma macierz trójktn górn i doln: -->U = triu(C) U = ! 1. 1. 1. 1. ! ! 0. 1. 1. 1. ! ! 0. 0. 1. 1. ! Macierze o elementach losowych. Funkcjarandpozwala utworzy macierz o elementach peudolosowych (pochodzcych z przedziaBu [0,1); mo|liwe jest tak|e u|ycie rozkBadu normalnego jak i podanie zarodka dla generatora liczb pseudolosowych): 6 -->M = rand(2, 5) M = ! 0.2113249 0.0002211 0.6653811 0.8497452 0.8782165 ! ! 0.7560439 0.3303271 0.6283918 0.6857310 0.0683740 ! n-elementowy wektor o staBej ró|nicy midzy elementami. Aby wprowadzi wektor (wierszowy)x o n wspóBrzdnych xn-x1 równomiernie rozmieszczonych pomidzy x1 i xn (innymi sBowy tak aby xi+1 -xi = , n wzBów, zatem n-1), u|ywamy n-1 funkcjilinspace. -->x = linspace(0,1,11) x = ! 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ! Instrukcj analogiczn pozwalajc utworzy wektor o zadanej warto[ci pierwszej wspóBrzdnej, ustalonej ró|nicy pomidzy wspóBrzdnymi i ostatniej wspóBrzdnej nie wikszej ni| zadana warto[ jest -->y = 0:0.3:1 y = ! 0. 0.3 0.6 0.9 ! SkBadnia jest nastpujca y = wartPocz:przyrost:wartGranicz. Tak dBugo jak pracujemy z liczbami caBkowitymi nie ma problemu z ustaleniem warto[ci granicznej odpowiadajcej ostatniej skBadowej wektora. -->i = 0:2:12 i = ! 0. 2. 4. 6. 8. 10. 12. ! Dla liczb rzeczywistych jest to zdecydowanie trudniejsze do okre[lenia: (i) przyrost mo|e posiad nieskoDczone rozwinicie w reprezentacji binarnej lub jego skoDczone rozwinicie mo|e wybiega poza zakres reprezentacji maszynowej liczb rzeczywistych powodujc ich zaokrglenia; (ii) bBdy zaokrgleD numerycznych kumuluj si w miar obliczania kolejnych skBadowych wektora. -->xx = 0:0.05:0.60 ans = ! 0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 ! Uwaga: w zale|no[ci od komputera na jakim uruchomiony zostanie Scilab powy|szy przykBad mo|e da ró|nie wyniki, to znaczy 0.6 mo|e pojawi si jako ostatnia skBadowa. Czsto przyrost równy jest 1, mo|na go w takiej sytuacji pomin: -->ind = 1:5 ind = ! 1. 2. 3. 4. 5. ! W przypadku gdy przyrost jest dodatni (ujemny) orazwartPocz>wartGrancz(wartPocz<wartGrancz) otrzymujemy wektor bez wspóBrzdnych nazywany w Scilab-ie macierz pust: -->i=3:-1:4 i = [] -->i=1:0 i = [] 2.3 Wyra|enia w Scilab-ie Scilab jest jzykiem posiadajcym bardzo prost skBadni, w której instrukcje przypisania maj posta zmienna = wyrazenie lub pro[ciej wyrazenie 7 gdzie w ostatnim przypadku warto[wyrazeniajest przypisana do domy[lnej zmiennej o nazwieans. Wyra|enia w Scilab- ie mog by tak proste (je[li chodzi o zapis) jak wyra|enia skalarne w innych jzykach programowania, ale mog skBada si z macierzy i wektorów co czsto sprawia trudno[ pocztkujacym u|ytkownikom tego jzyka. Dla wyra|eD  skalarnych mamy standardowe operatory+, -, /i ^i najcz[ciej stosowane funkcje przedstawione w tabeli 1. Zauwa|my, |e funkcje *, wymienione w tabeli poni|ej funkcji “ s dostpne od wersji 2.4 oraz, |e Scilab proponuje tak|e inne funkcje specjalne, wsród których mo|na znalez funkcje Bessela, eliptyczne, itd. abs warto[ bezwzgldna, moduB exp eksponent log logarytm naturalny log10 logarytm o podstawie 10 cos cosinus (argument w radianach) sin sinus (argument w radianach) sin(x) sinc x tan tangente (argument w radianach) cotg cotangente (argument w radianach) acos arccos asin arcsin atan arctg cosh cosinus hiperboliczny sinh sinus hiperboliczny tanh tangens hiperboliczny acosh argch asinh argsh atanh argth sqrt pierwiastek kwadratowy floor E(x) = ( #x #) = n Ô! n d" x < n + 1, x " N ceil #x # = n Ô! n - 1 < x d" n, x " N int int(x) = #x # je[li x > 0 oraz = #x # dla x d" 0 2 x 2 " erf funkcja bBdu erf(x) = e-t dt À 0 2 +" 2 " erfc dopeBnienie funkcji bBdu okre[lone przez ercf(x) = 1 - erf(x) = e-t dt À x +" gamma “(x) = tx-1e-tdt 0 lngamma ln(“(x)) d dlgamma ln(“(x)) dx Tablica 1: Wybrane funkcj u|ywane przez Scilab-a. tab:fonc Definiujc zmienn (bdc skalarem, wektorem, macierz) jej warto[ (warto[ci) nie musi by wyra|ona przez konkretn liczb, ale tak|e przez wyra|enie którego warto[ zostanie jej przypisana. -->M = [sin(%pi/3) sqrt(2) 5^(3/2); exp(-1) cosh(3.7) (1-sqrt(-3))/2] M = ! 0.8660254 1.4142136 11.18034 ! ! 0.3678794 20.236014 0.5 - 0.8660254i ! Uwaga: Powy|szy przykBad ilustruje potencjalne niebezpieczeDstwo podczas obliczania pierwiastka kwadratowego z liczy ujem- nej. Scilab rozwa|a czy ma do czynienia z liczbami zespolonymi i zwraca jeden z pierwiastów jako rezultat. 2.3.1 Kilka podstawowych przykBadów wyra|eD macierzowych Dostpne s wszystkie proste dziaBania wykonywane na macierzach: suma dwóch macierzy, iloczyn macierzy, iloczyn skalarny i macierzowy itd. Oto kilka przykBadów (w których wykorzystujemy wcze[niej zdefiniowane macierze). Uwaga: Tekst wystpu- jacy w danym wierszu po znaku//oznacza dla Scilab-a komentarz. Nie jest interpretowany a jedynie dostarcza pewnych uwag i wyja[nieD. -->D = A + ones(A) // napisz A aby zobaczyc jej zawartosc D = ! 2. 2. 2. ! ! 3. 5. 9. ! 8 ! 4. 10. 28. ! -->A + M // nie mozna wykonac dzialania dodawania (niezgodnosc wymiarow), // jaka jest odpowiedz Scilab-a? !--error 8 inconsistent addition -->E = A*C // C jest macierza jedynkowa o wymiarach (3,4) E = ! 3. 3. 3. 3. ! ! 14. 14. 14. 14. ! ! 39. 39. 39. 39. ! --> C*A // nie mozna wykonac mnozenia (niezgodnosc wymiarow), // jaka jest odpowiedz Scilab-a? !--error 10 inconsistent multiplication --> At = A // transpozycje otrzymuje sie stawiajac za nazwa macierzy znak apostrofu At = ! 1. 2. 3. ! ! 1. 4. 9. ! ! 1. 8. 27. ! --> Ac = A + %i*eye(3,3) // tworzymy macierz o elementach zespolonych Ac = ! 1. + i 1. 1. ! ! 2. 4. + i 8. ! ! 3. 9. 27. + i ! --> Ac_adj = Ac // w ten sposob otrzymujemy macierz transponowana o elementach // zespolonych zastapionych ich sprzezeniami Ac_adj = ! 1. - i 2. 3. ! ! 1. 4. - i 9. ! ! 1. 8. 27. - i ! -->x = linspace(0,1,5) // konstrukcja wektora kolumnowego x = ! 0. ! ! 0.25 ! ! 0.5 ! ! 0.75 ! ! 1. ! -->y = (1:5) // inny wektor kolumnowy y = ! 1. ! ! 2. ! ! 3. ! ! 4. ! ! 5. ! -->p = y *x // iloczyn skalarny wektorow x i y p = 10. -->Pext = y*x // otrzymujemy macierz 5x5 rzedu 1, dlaczego? Pext = 9 ! 0. 0.25 0.5 0.75 1. ! ! 0. 0.5 1. 1.5 2. ! ! 0. 0.75 1.5 2.25 3. ! ! 0. 1. 2. 3. 4. ! ! 0. 1.25 2.5 3.75 5. ! --> Pext / 0.25 // macierz mozna podzielic przez skalar ans = ! 0. 1. 2. 3. 4. ! ! 0. 2. 4. 6. 8. ! ! 0. 3. 6. 9. 12. ! ! 0. 4. 8. 12. 16. ! ! 0. 5. 10. 15. 20. ! --> A^2 // podniesienie do potegi drugiej macierzy ans = ! 6. 14. 36. ! ! 34. 90. 250. ! ! 102. 282. 804. ! --> [0 1 0] ans // mozna uzyc zmiennej ans, ktora zawiera wynik ostatniego * // dzialania gdy nie zostal on przypisany do zadnej zmiennej ans = ! 34. 90. 250. ! --> Pext*x - y + rand(5,2)*rand(2,5)*ones(x) + triu(Pext)*tril(Pext)*y; --> // wpisz ans aby zobaczyc wynik Inn, bardzo interestujc cech charakterystyczn, jest mo|liwo[ podania jako argumentu dla funkcji (z tabeli 1) macierzy zamiast kolejnych jej elementów. Innymi sBowy wpisanie instrukcjif(A)oznacza obliczenie warto[ci funkcjifna kolejnych elementach macierzyA; otrzymamy w ten sposób macierz [f(aij)]. PrzykBady: -->sqrt(A) ans = ! 1. 1. 1. ! ! 1.4142136 2. 2.8284271 ! ! 1.7320508 3. 5.1961524 ! -->exp(A) ans = ! 2.7182818 2.7182818 2.7182818 ! ! 7.3890561 54.59815 2980.958 ! ! 20.085537 8103.0839 5.320D+11 ! Uwaga: dla funkcji, które maj sens dla macierzy (co innego dla funkcji które stosuje sie dla ka|dego elementu macierzy. . . ) na przykBad funkcja eksponent, nazwa funkcji jest poprzedzona litermw ten sposób aby otrzyma eksponent macierzy A wystarczy wprowadzi komendexpm 2.3.2 DziaBania na elementach macierzy Aby pomno|y lub podzieli dwie macierz, A i B, o tych samych wymiarach, w taki spsób aby wynikiem byBa macierz, równie| o tych samych wymiarach, w której ka|dy element jest iloczynem (ilorazem) odpowiednich elementów macierzy A i B nale|y u|yc operatorów .* lub ./. A.*B jest macierz o elementach [aijbij] natomiast A./B jest macierz o lementach [aij/bij]. Podobnie mo|na podnie[ do potgi ka|dy z elementów macierzy wpisujc operator .^: A.^p pozwoli otrzyma macierz o wyrazach [ap ]. Rozwa|my przykBad: ij -->A./A ans = ! 1. 1. 1. ! ! 1. 1. 1. ! ! 1. 1. 1. ! 10 Uwagi: " W przypadku gdyAnie jest macierz kwadratow dziaBanieA^nbdzie wykonywane na kolejnych elementach macierzy A. Zaleca si jednak stosowanie zapisuA.^nponiewa| jest on bardziej czytelny. ij " Je[lisjest skalarem orazAjest macierz wówczass.^Adaje macierz o wyrazach sa . 2.3.3 Rozwizywanie ukBadów równaD liniowych Aby rozwiza ukBad równaD liniowych gdzie macierz wspóBczynników jest kwadratowa, Scilab stosuje rozkBad LU z cz[ciow zamian wierszy prowadzc do rozwiazania dwóch trójktnych ukBadów równaD. Jest to jednak operacja niewidoczna dla u|ytkownika dziki wyko|ystaniu operatora\: -->b=(1:3) // tworzymy wektor b b = ! 1. ! ! 2. ! ! 3. ! -->x=A\b // rozwiazujemy Ax=b x = ! 1. ! ! 0. ! ! 0. ! -->A*x - b // sprawdzamy poprawnosc wyniku ans = ! 0. ! ! 0. ! ! 0. ! Aby zapamita ten sposób postpowania, nale|y mie na uwadze ukBad pocztkowy Ax = b a nastpnie pomno|y ukBad lewostronnie przez A-1 (co oznacza podzielenie przez macierz A). Sposób ten daje dokBadny wynik, ale w ogólno[ci wystpuj bBdy zaokrglenia spowodowane arytmetyk liczb zmiennoprzecinkowych. -->R = rand(100,100); // stawiamy srednik na koncu aby uniknac zalewu ekranu liczbami -->y = rand(100,1); // srednik, jak wyzej -->x=R\y; // rozwiazanie ukladu Rx=y -->norm(R*x-y) // funkcja norm pozwala obliczyc norme wektorow (macierzy), // (obliczyc mozemy dwie normy - euklidesowa i hermitea) ans = 1.134D-13 Uwaga: Nie otrzymacie identycznego z moim je|eli funkcja rand nie zostanie u|yta tak jak w powy|szym przykBadzie. . . W momencie gdy rozwizanie ukBadu liniowego jest wtpliwe, Scilab wy[wietla informacje ostrzegajce i pozwalajce podj od- powiednie w takiej sytuacji dziaBania. 2.3.4 Indeksowanie, wydobywanie podmacierzy, konkatenacji macierzy i wektorów Aby odniesc si do konkretnego elemetnu macierzy wystarczy przy nazwie poda w nawiasie jego indeksy. Na przykBad: -->A33=A(3,3) A33 = 27. -->x_30 = x(30,1) x_30 = - 1.2935412 11 -->x(1,30) !--error 21 invalid index -->x(30) ans = - 1.2935412 Uwaga: Je|eli macierz jest wektorem kolumnowym wystarczy jednynie wpisa numer linii, w której znajduje si szukany element; analogicznie postpujemy w przypadku wektora wierszowego. Zalet jzyka Scilab jest mo|liwo[ Batwego wydobywanie podmacierzy z macierzy wyj[ciowej. -->A(:,2) // aby uzyskac 2 kolumne,... ans = ! 1. ! ! 4. ! ! 9. ! -->A(3,:) // ... 3 wiersz ans = ! 3. 9. 27. ! -->A(1:2,1:2) // podmacierz glowna rzedu 2 ans = ! 1. 1. ! ! 2. 4. ! Omówmy teraz ogóln skBadni. Niech macierz A ma wymiary (n.m), niech dalej v1 = (i1, . . . , ip) oraz v2 = (j1, . . . , jq) oznaczaj wektory (wierszowe lub kolumnowe), w których warto[ci s takie, |e 1 d" ik d" n et 1 d" jk d" m, wówczas A(v1,v2) oznacza macierz o wymiarach (p, q) utowrzon z wyrazów macierzy A odpowiadajcych wierszom i1, i2, . . . , ip oraz kolumnom j1, j2, . . . , jq . -->A([1 3],[2 3]) ans = ! 1. 1. ! ! 9. 27. ! -->A([3 1],[2 1]) ans = ! 9. 3. ! ! 1. 1. ! W praktyce dokonujemy prostrzych ekstrakcji, wydobywa si elementy umieszczone w przylegajcych blokach na przykBad w kolumnach lub wierszach. W takim przypadku u|yjemy wyra|enia i_poczatkowe:przyrost:i_koncowe w celu wygenerowania wektora wskazników. Natomiast aby wygenerowa peBny obszar odpowiadajcy wymiarowi u|yjemy operatora :(jak wida to w pierwszym przykBadzie). Zatem aby otrzyma podmacierz zBo|on z 1 i 3 wiersza zastosujemy -->A(1:2:3,:) // lub inaczej A([1 3],:) ans = ! 1. 1. 1. ! ! 3. 9. 27. ! Przejdzmy teraz do operacji konkatencaji macierzy, która umo|liwia poBczenie (ustawiajc obok siebie) wiele macierzy w celu otrzymania jednej zwanej macierz blokow. Dla przykBadu rozwa|my nastpujc macierz podzielon na bloki: ëø öø 1 2 3 4 ìø ÷ø 1 4 9 16 A11 A12 ìø ÷ø A = = . íø øø 1 8 27 64 A21 A22 1 16 81 256 Nale|y zatem zdefiniowa podmacierze A11, A12, A21, A22: 12 -->A11=1; -->A12=[2 3 4]; -->A21=[1;1;1]; -->A22=[4 9 16;8 27 64;16 81 256]; ostatecznie otrzymujemy macierz A powstaB z poBczenia 4 bloków: -->A=[A11 A12; A21 A22] A = ! 1. 2. 3. 4. ! ! 1. 4. 9. 16. ! ! 1. 8. 27. 64. ! ! 1. 16. 81. 256. ! Z punktu widzenia syntaktyki, macierze blokowe traktowane s jak zwykBe skalary (nale|y przy tym oczywi[cie pamita o zgodno[ci wymiarów odpowidnich macierzy blokowych). Istnieje odrbna skBadnia sBu|ca do jednoczesnego usunicia z macierzy wierszy lub kolumn: niech v = (k1, k2, . . . , kp) bdzie wektorem skBadajcym si z numerów wierszy lub kolumn macierzy M. PolecenieM(v,:)=[] spowoduje usunicie wierszy o numerach k1, k2, . . . , kp z macierzy M, natomiastM(:,v)=[]usunie kolumny k1, k2, . . . , kp. Ponadto je[li u jest wektorem (wierszowym lub kolumnowym),u(v)=[]usunie odpowiednie skBadowe. Kolejno[ elementów macierzy. Macierze w Scilab-ie s skBadowane kolumna za kolumn i ta kolejno[ elementów wyko- rzystywana jest w wielu funkcjach (porównaj dla przykBadu poleceniematrix, która umo|liwia zmian wymiarów macierzy). W szczególno[ci dla operacji wstawiania i wydobywania mo|liwe jest u|ycie domy[lnego porzdku przy wykorzystaniu jedny- nie pojedyDczego wektora indeksów (w miejsce dwóch, jako wskaznik kolumn lub wierszy). Oto kilka przykBadów opartych na wcze[niej zdefiniowanej macierzy A: -->A(5) ans = 2. -->A(5:9) ans = ! 2. ! ! 4. ! ! 8. ! ! 16. ! ! 3. ! -->A(5:9) = -1 // spowoduje wstawienie elementu -1 A = ! 1. - 1. - 1. 4. ! ! 1. - 1. 9. 16. ! ! 1. - 1. 27. 64. ! ! 1. - 1. 81. 256. ! 2.4 Informacje na temat [rodowiska pracy(*) Wystarczy wpisawhoa otrzymamy w ten sposób nastpujce informacje -->who your variables are... Anew A A22 A21 A12 A11 x_30 A33 x y R b Pext p Ac_adj Ac At E D cosh ind xx i linspace M U O zeros C B I Y c T startup ierr scicos_pal home PWD TMPDIR percentlib fraclablib soundlib xdesslib utillib tdcslib siglib s2flib roblib optlib metalib 13 elemlib commlib polylib autolib armalib alglib mtlblib SCI %F %T %z %s %nan %inf old newstacksize $ %t %f %eps %io %i %e %pi using 14875 elements out of 1000000. and 75 variables out of 1023 " Zmienne, które zostaBy wprowadzoneAnew, A,A22, A21, ...,bw porzdku odwrotnym do ich wczytywania; " Nazwy bibliotek Scilab-a (posiadajcych rozszerzenie lib) i funkcji:cosh. To znaczy funkcje (te opisane w jzyku Scilab- a) oraz biblioteki traktowane s przez Scilab-a jak zmienne. Uwaga: Procedury w Scilab-ie zaprogramowane w Fortran 77 i C nazywane s  prymitywami Scilab-a i nie s uwa|ane za zmienne w Scilab-ie. W dalszej cz[ci tego dokumentu u|ywa si czasem przesadnie sformuBowania  prymityw aby zaznaczy funkcje Scilab-a (programów w jzyku Scilab), które s zawarte w standardowo dostpnym [rodowisku. " StaBe predefiniowane, takie jak À, e i i epsilon maszynowy eps oraz dwie inne staBe, klasyczne w arytmetyce zmienno- przecinkowej: nan - ang. not a number, inf  ". Zmienne, których nazwa poprzedzona jest znakiem % nie mog zosta usunite. " Bardzo istotna zmiennanewstacksizeokre[lajca rozmiar stosu (to znaczy ilo[ dostpnej pamici). " Na koniec wy[wietlona zostaje informacja o ilo[ci u|ytych sBów 8 bitowych w odniesieniu do ich caBkowitej (dostpnej) ilo[ci a nastpnie ilo[ u|ytych zmiennych w stosunku do maksymalnej dozwolonej ich ilo[ci. Rozmiar stosu mo|na zmieni za pomoc poleceniastacksize(nbmbts)gdzienbmbtsoznacza nowy po|dany roz- miar. Ta sama komenda bez argumentu pozwala uzyska rozmiar stosu mieszczcy maksymaln dopuszczaln liczb zmiennych. Tak wic chcc usun wektor v1 ze [rodowiska a wic zwolni miejsce w pamici nale|y u|y polecenia clear v1. Polecenie clear u|yte bez argumentu usunie wszystkie zmienne. Chcc natomiast usunc wicej ni| jedn zmienn mo|a poda ich nazwy jako kolejne argumenty poleceniaclear, na przykBad:clear v1 v2 v3.s 2.5 WywoBanie pomocy z linii poleceD Pomoc mo|na otrzyma wybierajc przyciskHelpz menu okna Scilab. Poczwszy od wersji 2.7 strony pomocy s w formacie HTML5. Mo|liwe jest u|ycie innego programu do przegldania dokumentacji w nastpujcy sposób -->global %browsehelp; %browsehelp=[]; -->help W wierszu pierwszym modyfikowana jest zmienna lokalna %browsehelpw taki sposób, |e zawiera ona macierz pust. Zatem wprowadzenie poleceniahelplub wybranie odpowiedniego przycisku, zostan zaproponowane inne alternatywne mo|liwo[ci. Mo|na t zmienn umie[ci w pliku.scilabdziki czemu uniknie si powy|szych manipulacji6. Wpisujc poleceniehelp (lub klikajc na przyciskHelp) pojawi si strona HTML zawierajca odno[niki do wszystkich stron pomocy (Scilab Programming, Graphic Library, Utilities and Elementary functions,...). Klikajc na wybran nazw otrzymuje si list szystkich funkcji zwizanych z danym tematem, do ka|dej z nich doBczono krótki opis. Wybierajc dan funkcj otrzymuje si stron z pytaniami dotyczcymi danej funkcji. Je|eli natomiast znana jest nazwa funkcji wystarczy wpisa komendhelpNazwaFunkcji w oknie Scilab-a aby wej[ na odpowiedn stron. Polecenieapropos sBowoKluczowe wy[wietli wszystkie strony na których pojawia si sBowoKluczowe w nazwach funkcji lub ich opisach. Przeszu- kiwanie opcji mo|e okaza sie niewygodne, gdy| czsto wiele funkcji zgrupowanych jest pod jedn nazw (na przykBad jako Elementary functions mojetutu wystpuje bardzo du|o ró|nych funkcji, które w jakimkolwiek stopniu na t nazw zasBuguj); nie wahaj sie zatem u|ywa poleceniaapropos 2.6 Generator prostych wykresów Przypu[my, |e chcemy narysowa wykres funkcji y = e-x sin(4x) dla x " [0, 2À]. Mo|na zdefiniowa podziaB przedziaBu poprzez funkcjlinespace -->x=linspace(0,2*%pi,101); a nastpnie policzy warto[ci funkcji dla ka|dego elementu podziaBu wykrzystujc w tym celu instrukcj -->y=exp(-x).*sin(4*x); 5 Format podstawowy to XML a z niego otrzymuje si HTML. 6 Je[li wasza przegldarka nie zostaBa uwzgldniona nale|y dokona modyfikacji plikuSCI/macros/util/browsehelp.sci 14 i ostatecznie -->plot(x,y, x , y , y=exp(-x)*sin(4x) ) gdzie trzy ostatnie BaDcuchy znaków (odpowiednio: opis osi rzdnych i odcitych oraz nazwa wykresu) s opcjonalne. Instrukcja ta pozwala narysowa krzyw przechodzc przez punkty, których wspóBrzdne okre[la wektorxna osi odcitych orazyna osi rzdnych. Je|eli punkty ukBadaj si cz[ciowo wzdBu| linii prostej wykres bdzie tym wierniejszy im wikszej ilo[ci punktów u|yjemy y=exp(-x)*sin(4x) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 0 1 2 3 4 5 6 7 x Rysunek 1: Prosty wykres. Uwaga: Ta instrukcja jest ograniczona, zatem nie nale|y sugerowa si jedynie tym rozwizaniem. W rozdziale dotyczcym grafiki wprowadza si instrukcjeplot2d, która jest o wiele lepsza. 2.7 Pisanie i wykonywanie skryptów Mo|liwe jest wpisanie cigu instrukcji do plikunazwaPlikui ich wywoBanie poprzez komend -->exec( nazwaPliku ) // lub exec("nazwaPliku") Prostsz metod jest wybranie pozycjiFile Operationznajdujc si w menuFileokna Scilab. Menu pozwala wybra insteresujcy nas plik (ewentualnie w celu jego modyfikazji); aby go wykona wystarczy wybraExec. Jako przykBad skryptu rozwa|my ponownie wykres funkcji e-x sin(4x) zadajc za ka|dym razem inny rozmiar przedziaBu [a, b] i ilo[ podprzedzaiBów. Napiszmy zatem w pliku nazywajcym si na przykBadscript1.scenastpujce instrukcje // pierwszy skrypt a = input(" Podaj lewy kraniec przyedzialu : "); b = input(" Podaj prawy kraniec przedzialu : "); n = input(" Podaj ilosc podprzedzialow : "); // obliczenie odcietych x = linspace(a,b,n+1); // obliczenie rzednych y = exp(-x).*sin(4*x); // maly rysunek plot(x,y, x , y , y=exp(-x)*sin(4x) ) 15 y Uwaga: " Aby odnalez si w naszych plikach Scilab-a zaleca si wykorzysta w nazwie pliku skryptowego rozszerzenie*.sce(w przypadku gdy plik zawiera funkcj, rozszerzeniem powinno by*.scu). " Pewne edytory mog by wyposa|one w narzdzia specyficzne dla edycji w Scilab-ie (zobacz strona Scialba). Dlaemacs istniej dwa, z krórych lepszy pochodzi od Alexandra Vigodnera, którego najnowsze wersje mo|na znalez pod adresem http://www.geocities.com/avezunchik/scilab.html " Skrypt jest zazwyczaj pierwotnym programem dla aplikacji napisanych w Scilab-ie. 2.8 Dodatkowe informacje 2.8.1 Skracanie instrukcji przy zapisie macierzowym Wcze[niej dowiedzieli[my si, |e mno|enie macierzy przez skalar (podobnie dzielenie macierzy przez skalar) jest zdefiniowane w Scilab-ie. Natomiast Scilab stosuje uproszczeD mniej oczywistych, takich jak dodawanie skalara i macierzy. Wyra|enieM + sgdzieMjest macierza asskalarem jest skrótem dla instrukcji M + s*ones(M) to znaczy skalar jest dodany do wszystkich elementów macierzy. Inny skrót: w wyra|eniu typu A./B (oznaczajcym dzielenie element po elemencie macierzy o tych samych wymiarach), je|eliAjest skalarem wówczas wyrarzenie to jest skróconym zapisem dla instukcji : A*ones(B)./B otrzymuje sie macierz [a/bij]. Takie uproszczenia pozwalaj na zapis bardziej syntetyczny w wielu przypadkach (por. wi- czenia). Na przykBad, je[li f jest funkcj zdefiniowan w [rodowisku, x jest wektorem a s jest zmienn, wówczas s./f(x) jest wektorem tych samych wymiarów co x, w którym i-ta wspóBrzdna równa jest s/f(xi). W ten sposób, aby wyznaczy wektor o wspóBrzdnych 1/f(xi), u|yjemy skBadni : 1./f(x) Taka skBadnia interpretuje 1. jako jedynk i wynik nie bdzie wBa[ciwy. Dobrym rozwizaniem aby otrzyma poszukiwany wynik jest ujcie liczby w nawias lub oddzielenie liczby spacj : (1)./f(x) // lub takze 1 ./f(x) lub 1.0./f(x) 2.8.2 PozostaBe uwagi dotyczce rozwizywania ukBadów równaD liniowych (*) 1. W przypadku gdy mamy wicej wyrazów wolnych, mo|na stosowa nastpujce polecenie : -->y1 = [1;0;0;0]; y2 = [1;2;3;4]; // przypadek gdy mamy 2 wyrazy wolne (mozna // ujac wiecej instrukcji w jednej lini) -->X=A\[y1,y2] // konkatenacja/laczenie y1 i y2 X = ! 4. - 0.8333333 ! ! - 3. 1.5 ! ! 1.3333333 - 0.5 ! ! - 0.25 0.0833333 ! Pierwsza kolumna macierzy X jest rozwizaniem ukBadu równaD Ax1 = y1, natomiast druga jest rozwizaniem ukBadu Ax2 = y2. 2. Wcze[niej zobaczyli[my, |e je|eli A jest macierz kwadratow (n,n) oraz b jest wektorem kolumnowym o n wspóBrzdnych (a wic macierz (n,1)) to : 16 x = A\b da nam rozwizanie ukBadu równaD liniowych postaci Ax = b. Je|eli natomiast macierz A oka|e sie macie| osobliw, Scilab zwróci bBd. Na przykBad : -->A=[1 2;1 2]; -->b=[1;1]; -->A\b !--error 19 singular matrix Podczas gdy macierz A jest zle uwarunkowana (lub ewentualnie zle zrównowa|ona) odpowiedz bdzie dostarczona ale wraz z komentarzem w postaci ostrze|enia zawierajcym oszacowanie odwrotno[ci uwarunkowania (cond(A) = ||A|| ||A-1||): -->A=[1 2;1 2+3*%eps]; -->A\b warning matrix is close to singular or badly scaled. results may be inaccurate. rcond = 7.4015D-17 ans = ! 1. ! ! 0. ! Z drugiej strony, je|eli macierz nie jest kwadratowa, ale posiada t sam liczb wirszy co wektor wyrazów wolnych, Sci- lab obliczy rozwizanie (w postaci wektora kolumnowego, którego liczba wspóBrzdnych bdzie równa liczbie kolumn macierzy A) nie informujc u|ytkownika o bBdzie. W efekcie, je|eli równanie Ax = b nie ma w ogólno[ci rozwi- 7 zania jednoznacznego , mo|na zawsze wybra wektor x który zweryfikuje pewne wBasno[ci (x o minimalne normie i jednocze[nie bdcy rozwizaniem min ||Ax - b||). W tym przypadku takie rozwizanie jest wynikiem dzieBanie innych algorytmów które umo|liwiaj (ewentualne) otrzymanie pseBdo-rozwizania8. Niedogodno[ jest taka, |e je[li zostanie popeBniony bBd podczas definiowania macierzy (na przykBad zdefiniowana zaostanie dodatkowa kolumna tak, |e macierz bdzie miaBa wymiar (n,n+1)) bBd ten nie zostanie wykryty natychmiast. Rozwa|my przykBad : -->A(2,3)=1 // zaburzenie A = ! 1. 2. 0. ! ! 1. 2. 1. ! -->A\b ans = ! 0. ! ! 0.5 ! ! - 3.140D-16 ! -->A*ans - b ans = 1.0D-15 * ! - 0.1110223 ! ! - 0.1110223 ! Poza zwróceniem naszej uwagi na konsekwencje tego typu zakBuceD, powy|szy przykBad ilustruje nastpujce fakty/aspekty : 7 niech (m, n) bdzie wymiarem macierzy A (gdzie n = m), mamy jedno rozwizanie wtedy i tylko wtedy, gdy m > n, KerA = {0} i w konsekwencji b " ImA ten ostatni warunek jest wyjtkiem je[li b jest zadane z góty w Km ; w ka|dym innym przypadku albo nie ma |adnego rozwizania albo isnieje nieskoDczenie wiele rozwizaD 8 w przypadkach trudnych, tzn. takich, gdzie macierz nie ma maksymalnego rzdu (rg(A) < min(n, m) gdzie n i m s dwuwymiarowe) nale|y wyznaczy rozwizanie korzystajac z pseBdo-odwrotnej macierzy do A (x = pinv(A)*b). Û 17 " x = A\y pozwala rozwiza problem mniej kwadratowy (kiedy macierz nie ma maksymalnego rzdu, nale|y le- piej u|y/wyko|ystax = pinv(A)*b, pseudo-inwersja zostanie obliczona poprzez dekompozycj pojedynczych warto[ci A (t dekompozycj mo|na otrzyma za pomoc instukcjisvd)); " instrukcjaA(2,3)=1(bBd zakBucenia. . . ) jest w rzeczywisto[ci skrótem od : A = [A,[0;1]] to znaczy Scilab domy[la si, |e macierz A ma by uzupeBniona (o 3 kolumn) ale brakuje mu jednego elementu. W takim przypadku uzupeBni je zerami. " element na pozycji (2,2) jest równy 2 + 3 µm. Epsilon maszynowy (µm) ma|e by zdefiniowany jako bardzo du|a liczba, dla której 1 •" µm = 1 w arytmetyce zmiennoprzecinkowej9. W konsekwencji musimy mie 2 •" 3µm > 2 aby na ekranie pojawiBa si 2. Wynika to z u|ytego domy[lnie formatu danych, mo|na go jednak zmodyfikawa u|ywajc instrukcjiformat: -->format( v ,19) -->A(2,2) ans = 2.0000000000000009 W ten sposób zmienimy format domy[lnyformat( v ,10)(patrzHelpaby pozna znaczenie argumentów). " Rozwizanie ukBadu Ax = b, które nie zostaBo przypisane |adnej konkretniej zmniennej, jest domy[lnie nazwane ans. Zmienn t mo|na nastpnie wyko|ysta do obliczenia warto[ci rezydualnej Ax - b. 3. Przy u|yciu Scilab-a mo|na równie| odpowiednio rozwiza ukBad równaD postaci xA = b gdzie x i b s wektorami wier- szowymi, a A jest macierz kwadratow (poprzez transpozycj wracamy do ukBadu klasycznego A¤"x¤" = b¤"), wystarczy pomno|y prawostronnie przez A-1 (odpowiednikiem tej operacji jest prawostronne podzielenie przez macierz A) : x = b/A I podobnie jak poprzednio, je[li A jest macierz prostoktn (gdzie liczba kolumn jest równa liczbie wierszy wektora b), Scilab wyznaczy rozwizanie. 2.8.3 Kilka prostych macierzy (*) Suma, iloczyn wspóBczynników macierzy, macierz pusta Aby zsymowa wspóBczynniki macierzy, u|ywamy funkcjisum: -->sum(1:6) // 1:6 = [1 2 3 4 5 6] : powinnismy zatem otrzymac 6*7/2 = 21 !!!!! ans = 21. Funkcja ta przyjmuje argument dodatkowy dla obliczenia tylko wierszy lub tylko kolumn : -->B = [1 2 3; 4 5 6] B = ! 1. 2. 3. ! ! 4. 5. 6. ! -->sum(B,"r") // oblicza sume wyrazow w poszczegolnych kolumnach, // otrzymujemy wektor wierszowy postaci ans = ! 5. 7. 9. ! -->sum(B,"c") // oblicza sume wyrazow w poszczegolnych wierszach, // otrzymujemy wektor kolumnowy postaci ans = ! 6. ! ! 15. ! 9 [ci[le mówic wszystkie liczby rzeczywiste x takie, |e m d" |x| d" M mog by zakodowane/zapisane za pomoc liczby zmiennoprzecinkowej fl(x) nastpujaco : |x-fl(x)| d" µm|x| gdzie m i M s odpowiednio najmniejszaÛi najwikszaÛliczb dodatni mo|liw do zakodowania za pomoc znormalizowanej Û liczby zmiennoprzecinkowej. 18 Jednym z wyjtkowo u|ytecznych obiektów jest macierz pusta, któr definiuje sie nastpujco : -->C =[] C = [] Macierz ta ma nastpujce wBasno[ci :[] + A = Ai[]*A = []. Stosujc teraz funkcjsumdo macierzy pustej otrzy- mujemy naturalnie rezultat : -->sum([]) ans = 0. wedBug matematycznej konwencji sumowania : n S = ui = ui gdy E = {1, 2, . . . , n} i"E i=1 je|eli E jest zbiorem pustym wówczas S = 0 W analogiczny sposób, aby otrzyma iloczyn elementów macierzy, stosuje sie funkcjprod : -->prod(1:5) // otrzymujemy 5! = 120 ans = 120. -->prod(B,"r") // napisz B aby zobaczyc macierz... ans = ! 4. 10. 18. ! -->prod(B,"c") ans = ! 6. ! ! 120. ! -->prod(B) ans = 720. -->prod([]) ans = 1. wynik otrzymuje si[ zgodnie z konwencj znan w matematyce : ui = 1, si E = ". i"E Suma i iloczyn skumulowany Funkcjecumsumicumprodobliczaj odpowiednio sum i iloczyn skumiulowany macierzy lub wektora : -->x = 1:6 x = ! 1. 2. 3. 4. 5. 6. ! -->cumsum(x) ans = ! 1. 3. 6. 10. 15. 21. ! -->cumprod(x) ans = ! 1. 2. 6. 24. 120. 720. ! 19 Dla macierzy kumulowanie dokonuje si[ zazwyczja jedynie w porzdku kolumnowym : -->x = [1 2 3;4 5 6] x = ! 1. 2. 3. ! ! 4. 5. 6. ! -->cumsum(x) ans = ! 1. 7. 15. ! ! 5. 12. 21. ! Podobnie jak w przypadku funkcjisumiprod, mo|na dokonywa kumulowania wedBug wierszy i kolumn : -->cumsum(x,"r") // suma skumulowana wedlug kolumn ! ans = ! 1. 2. 3. ! ! 5. 7. 9. ! -->cumsum(x,"c") // suma skumulowana wedlug wierszy ! ans = ! 1. 3. 6. ! ! 4. 9. 15. ! Pojawia si tutaj pozorna niezgodno[ pomidzy nazw a sposobem dokonywania obliczeD na wierszach i kolumnach10 : dla funkcjisumi prod argument informuje o ostatecznej formie otrzymanego wyniku (sum(x, r )pozwala otrzyma wektor wierszowy (r jak row, ang. wiersz) a wic sum ka|dej kolumny) lecz takie rozumowanie dziala jedynie w przypadku tych funkcji. minimum i maksimum wektora i macierzy Najmniejsz i najwiksz warto[ wektora lub macierzy mo|na wyznaczy posBugujc si funkcjamiminimax. Ich dzia- Banie jest takie jak w przypadku funkcjisumiprod gdzie dodatkowy argument okre[la czy minimum lub maksimum ma by wyznaczone dla wiersza czy dla kolumny. Dodatkowo mo|liwe jest okre[lenie wskaznika dla miejsca, w któtym znajduje sie owo minimum lub maksimum. PrzykBady : -->x = rand(1,5) x = ! 0.7738714 0.7888738 0.3247241 0.4342711 0.2505764 ! -->min(x) ans = 0.2505764 -->[xmin, imin] = min(x) imin = 5. // min jest osiagniete w x(5) xmin = 0.2505764 -->y = rand(2,3) y = ! 0.1493971 0.475374 0.8269121 ! ! 0.1849924 0.1413027 0.7530783 ! -->[ymin, imin] = min(y) imin = ! 2. 2. ! // min jest osiagniete dla y(2,2) ymin = 0.1413027 -->[ymin, imin] = min(y,"r") // minimum kazdej kolumny 10 sum,prod,cumsum,cumprod,mean,st_deviation,min,max 20 imin = ! 1. 2. 2. ! // => min to y(1,1) y(2,2) y(2,3) ymin = ! 0.1493971 0.1413027 0.7530783 ! -->[ymin, imin] = min(y,"c") // minimum kazdego wiersza imin = ! 1. ! // min to y(1,1) ! 2. ! // y(2,2) ymin = ! 0.1493971 ! ! 0.1413027 ! [rednia i odchylenie standardowe Funkcje mean i st_deviation pozwalaj wyznaczy [redni i odchylenie standardowe wspóBrzdnych wektora lub macierzy. FormuBa u|ywana dla odchylenia standardowego jest nastpujca : 1/2 n n 1 1 Ã(x) = (xi - x)2 , x = xi ¯ ¯ n - 1 n i=1 i=1 -->x = 1:6 x = ! 1. 2. 3. 4. 5. 6. ! -->mean(x) ans = 3.5 -->st_deviation(x) ans = 1.8708287 Podobnie mo|na wyznaczy [redni (a tak|e odchylenie standardowe) wiersza lub kolumny macierzy : -->x = [1 2 3;4 5 6] x = ! 1. 2. 3. ! ! 4. 5. 6. ! -->mean(x,"r") // srednia wyrazow kazdej z kolumn macierzy ans = ! 2.5 3.5 4.5 ! -->mean(x,"c") // srednia kazdego z wierszy ans = ! 2. ! ! 5. ! PrzeksztaBcenie macierzy. Funkcjamatrixumo|liwia takie przeksztaBcenie macierzy aby miaBa ona nowe wymiary (przy zachowaniu tych samych wspóBczynników). -->B = [1 2 3; 4 5 6] B = ! 1. 2. 3. ! ! 4. 5. 6. ! -->B_new = matrix(B,3,2) B_new = 21 ! 1. 5. ! ! 4. 3. ! ! 2. 6. ! Jest to mo|liwe dziki temu, |e wspóBczynniki macierzy zapamitywane s kolumnowo. Jednym z zastowowaD funkcjimatrix jest transpozycja wektora wierszowego na kolumnowy i odwrotnie. Nale|y zaznaczy, |e w ogólno[ci pozwala ona tak|e prze- ksztaBci macierzAzamieniajc (wektory wierszowe i kolumnowe) na wektor kolumnowyv: v = A(:), przykBad : -->A = rand(2,2) A = ! 0.8782165 0.5608486 ! ! 0.0683740 0.6623569 ! -->v=A(:) v = ! 0.8782165 ! ! 0.0683740 ! ! 0.5608486 ! ! 0.6623569 ! Wektory z postpem logarytmicznym. Czasami istnieje potrzeba operownia wektorem którego wspóBrzdne stanowi postp logarytmiczny (tzn. takie, |e stosunek dwuch ssiednichc jest staBy: xi+1/xi = const): mo|na u|y w tym przypadku funkcj logspace : logspace(a,b,n): pozwalajc otrzyma taki wektor o n wspóBrzdnych w którym pierwsza i ostatnia wsóBrzdna s odpowiednio 10a i 10b, np: -->logspace(-2,5,8) ans = ! 0.01 0.1 1. 10. 100. 1000. 10000. 100000. ! Warto[ci i wektory wBasne Funkcjaspecpozwala obliczy warto[ci wBasne macierzy (kwadratowej!) : -->A = rand(5,5) A = ! 0.2113249 0.6283918 0.5608486 0.2320748 0.3076091 ! ! 0.7560439 0.8497452 0.6623569 0.2312237 0.9329616 ! ! 0.0002211 0.6857310 0.7263507 0.2164633 0.2146008 ! ! 0.3303271 0.8782165 0.1985144 0.8833888 0.312642 ! ! 0.6653811 0.0683740 0.5442573 0.6525135 0.3616361 ! -->spec(A) ans = ! 2.4777836 ! ! - 0.0245759 + 0.5208514i ! ! - 0.0245759 - 0.5208514i ! ! 0.0696540 ! ! 0.5341598 ! oraz wy[wietla wyniki w postaci wektora kolumnowego (Scilab stosuje metod QR polegajc na interacyjnej dekompozycji) Schur macierzy). Wektory wBasne mog by otrzymane za pomocbdiag. W ogólno[ci aby wyznaczy warto[ci wBasne mo|na wyko|ysta funkcjgspec. 2.8.4 Funkcjesizeilength sizepozwala uzyska informacj o wymiarach macierzy (w kolejno[ci: liczba wierszy, liczba kolumn): -->[nl,nc]=size(B) // B jest macierza (2,3) z poprzedniego przykladu nc = 3. nl = 2. 22 -->x=5:-1:1 x = ! 5. 4. 3. 2. 1. ! -->size(x) ans = ! 1. 5. ! natomiastlengthmówi o liczbie elementów macierzy (rzeczywistych lub zespolonych). W ten sposób dla wektora wierszo- wego lub kolumnowego otrzymuje si dokBadnie liczb jego wspóBczynników : -->length(x) ans = 5. -->length(B) ans = 6. Obie te instrukcje u|ywane s po to aby pozna rozmiar macierzy lub wektora bdcego argumentem tych funkcji. Nale|y doda, |e size(A, r ) (lub size(A,1)) oraz size(A, c ) (lub size(A,2)) mówi o liczbie wierszy (rows) oraz kolumn (columns) macierzy A. 2.9 wiczenia 1. Zdefiniowa macierz wymiaru n nastpujco (zob. szczegóBy dotyczce funkcjidiagza pomocHelp): ëø öø 2 -1 0 · · · 0 ìø .. .. .. . ÷ø . ìø ÷ø . . . -1 . ìø ÷ø ìø .. .. .. ÷ø A = ìø ÷ø . . . 0 0 ìø ÷ø ìø . .. .. .. ÷ø . íø øø . . . -1 . 0 · · · 0 -1 2 2. Niech A bdzie macierz kwadratow; czym bdziediag(diag(A))? 3. Funkcjatrilpozwala wydoby tróktn doln cz[ macierzy (triudla cz[ci trójktnej górnej). Zdefiniowa dowoln macierz kwadratow A (np. za pomocrand) i skonstruowa pojedyncz instrukcj macierz trójktn górn T tak, |e tij = aij dla i > j (cz[ci trójktne górne macierzy A i T s jednakowe) i tak, |e tii = 1 (T jest diagonalnie jednostkowa). 4. NiechXbdzie macierz (lub wektorem) zdefiniowanym w [rodowisku. Napisa instrukcj pozwalajc obliczy macierz Y(tego zamego rozmiaru coX) w której element na pozycji (i, j) jest równy warto[ci f(Xij) w nastpujcych przypadkach : (a) f(x) = 2x2 - 3x + 1 (b) f(x) = |2x2 - 3x + 1| (c) f(x) = (x - 1)(x + 4) 1 (d) f(x) = 1+x2 sin x 5. Naszkicowa wykres funkcji f(x) = dla x " [0, 4À] (napisa skrypt). x 6. MaBa ilustracja prawa wielkich liczb : wyznaczy za pomoc funkcjirand n próbek wedlug prawa zero-jedynkowego w postaci wektora x, obliczy [redni skumulowan tego wektora, tzn. wektor x, którego ka|da z n wspóBrzdnych ¯ k 1 ma warto[ : xk = xi oraz naszkicowa wykres cigu otrzymanego w ten sposób. Zbada jego zachowanie, ¯ k i=1 zwikszajc warto[ n. 23 3 Programowanie w Scilabie Scilab, oprócz prostych instrukcji (pozwalajcych, z wyko|ystaniem instrukcji macierzowch, na bardzo syntetnyczne programo- wa bliskie jzykowi mamtematycznemu, co najmniej macierzowemu) dysponuje prostym aczkolwiek wystarczajco komplet- nym jzykiem programowania. Zasadnicz ró|nic w porównaniu z najpopularniejszymi dzi[ jzykami (C, C++, Pascal, . . . ) jest brak konieczno[ci deklarowania zmiennych: w przeciwieDstwie do operacji wcze[niej wykonywanych, nie musimy w |ad- nym momencie precyzowa rozmiaru macierzy ani jej typu (rzeczywista, zespolona, . . . ). To na interpreterze Scilaba spoczywa zadanie rozpoznawania ka|dego nowego obiektu. Takie podej[cie nie jest czsto brane pod uwag, gdy| z punktu widzenia programisty prowadzi mo|e do powstawania trudnych do wychwycenia bBdów. W przypadku jzyków takich jak Scilab (czy MATLAB) nie powoduje to problemów dziki u|ywaniu zapisu wektorowego/macierzowego oraz gotowych instrukcji prostych, pozwalajcych na znaczne zmniejszenie ilo[ci zmiennych oraz wielko[ci kodu (na przykBad dziki instrukcjom macierzowym unikamy znacznej ilo[ci ptli sBu|cych do definiowania macierzy). Inna zalet jzyków programowania tego typu jest dyspo- nowanie instrukcjami graficznymi (dziki czemu unikamy konieczno[ci ko|ystania z innego programu lub stosowania bibliotek graficznych) oraz bycie interpretowanym (co pozwala na Batwe wyszukiwanie bBdów bez pomocy debugera). Nie mniej jednak niedogodnym jest, |e programy napisane w Scilabie s wolniejsze od tych napisanych w C (ale program w C wymaga 50 razy wicej czasu na jego poprawne napisanie). Z innej strony, w aplikacjach gdzie ilo[ zmiennych jest szczególnie istotna (na przy- kBad macierz 1000 x 1000) lepiej powróci do jzyka kompilowanego. Istotnym punktem wpBywajcym na szybko[ wykonania jest konieczno[ stosowania maksymalnie czsto instrukcji prostych i instrukcji macierzowych (oznacza to ograniczenie do mi- nimum ilo[ ptli)11. W efekcie Scilab odwoBuje si tak|e do skompilowanych funkcji (fortran) dziki czekmu interpreter ma mniej pracy (pracuje mniej). Dla wykorzystania najlepszego z tych dwóch modeli (to znaczy szybko[ci jzyka kompilowanego i wygody [rodowiska takiego jak Scilab), u|ytkownik ma mo|liwo[ wywoBywania podprogramów w fortranie(77) lub C. 3.1 Ptle W Scilabie mamy mo|liwo[ u|ycia jednej z dwóch petli: fororazwhile. 3.1.1 Ptlafor PtlaforsBu|y do powtarzania pewnego cigu instrukcji; zmienna starujca ptli przyjmuje kolejne warto[ci wektora liniowego -->v=[1 -1 1 -1] -->y=0; for k=v, y=y+k, end Ilo[ iteracji (powtórzeD) okre[lona jest przez ilo[ skBadników wektora; w i-tej iteracji warto[k jest równav(i). Opisowo ptleforprzedstwi mo|na nastpujco dla i := istart a| do iend z krokiem istep powtarzaj : instrukcje koniec powtarzaj Jako wektora mo|na u|y wyra|eniaiStart:iStep:iEnd; je[li przyrostiStepma warto[ 1 to jak widzieli[my wcze[niej mo|e zosta pominity. Przedstawiona wcze[niej jako przykBad ptla for mo|e zosta tak|e napisana w bardziej naturalny sposób -->y=0; for i=1:4, y=y+v(i), end Kilka uwag: " Zmienna sterujca mo|e tak|e przyjmowa warto[ci pochodzce z macierzy. W takim przypadku ilo[ iteracji równa jest ilo[ci kolumn macierzy, a zmienna sterujca w i-tym przebiegu ptli jest równa i-tej kolumnie macierzy -->A=rand(3,3);y=zeros(3,1); for k=A, y = y + k, end " DokBadna skBadnia instrukcjiforwyglda nastpujco for zmienna = macierz, ciag instrukcji, end gdzie instrukcje oddzielone s od siebie przecinkiem (lub [rednikiem je[li nie chcemy widzie na ekranie rezultatów wy- konania instrukcji). Poniewa| lista zmiennych rozdzielonych przecinkiem równowa|na jest takiej samej li[cie rozdzielonej spacjami, dlatego w skryptach lub funkcjach cz[ciej stosuje si wygodniejszy (bardziej czytelny) zapis: 11 patrz równie|: sekcja "Kilka uwag dotyczcych szybko[ci" 24 for zmienna = macierz instrukcja1 instrukcja2 ........... instrukcjan end gdzie instrukcja mo|e koDczy si [rednikiem (bdzie tak zawsze jesli chcemy unikn pojawiania si rezultaów instrukcji na ekranie)12. 3.1.2 Ptlawhile Ptlawhilepozwala na powtarzanie cigu instrukcji tak dBugo dopóki prawdziwy jest pewnien warunek: -->x=1 ; while x<14,x=2*x, end Zadajc warunek mo|emy wykorzysta nastpujce operatory porównania: == równe < mniejsze > wiksze <= mniejsze lub równe >= wiksze lub równe ~= lub<> ró|ne Scilab posiada wbudowany typ logiczny, w którym dla oznaczenia prawdy stosuje si symbole%tlub%Tza[ faBszu%flub%F. Mo|na oczywi[cie zdefiniowa macierz lub wektor zBo|ony z warto[ci logicznych. Dostpne s nastpujce operatory logiczne : & i | lub ~ nie Formalna skBadnia ptliwhileprzedstawia si jak poni|ej while warunek, instrukcja_1, ... ,instrukcja_N , end lub (czsto w przypadku skryptu lub funkcji): while warunek instrukcja_1 ............. instrukcja_N end gdzie ka|da instrukcja_k mo|e zosta zakoDczona [rednikiem, warunek za[ jest wyra|eniem zwracajcym warto[ lo- giczn. 3.2 Instrukcja warunkowa S dwie instrukcje sterujce przebiegiem wykonania programu: if-then-else oraz select-case. 3.2.1 Instrukcjaif-then-else Spójrzmy na przykBad poni|ej : -->if x>0 then, y=-x,else,y=x,end // zmienna x powinna zostac wczesniej zdefiniowana Podobnie jak to miaBo miejsce dla instrukcji ptli przecinki nie s znakami obowizkowymi. Jak w wikszo[ci spotykaych jzyków programowania, w przypadku gdy nie podejmujemy |adnego dziaBania zwizanego z niespeBnieniem warunku, mo|emy cz[elsepomina. W przypadku gdy cz[elsezawiera by miaBa kolejn instrukcjif-then-elsezamiast ciaguelse orazifmo|na u|yelseifco prowadzi do ogólnego schematu tej instrukcji 12 jest tak tylko w przypadku, gdy mamy do czynienia ze skryptem, w przypadku funkcji, wynik instrukcji nie jest wy[wietlany nawet gdy nie jest ona zakoDczona [rednikiem; zasad t mo|na zmodyfikowa u|ywajac instrukcji mode Û 25 if warunek1 then ciag instrukcji wykonywanych gdy spelniony jest warunek1 elseif warunek2 then ciag instrukcji wykonywanych gdy spelniony jest warunek2 ... elseif warunekN then ciag instrukcji wykonywanych gdy spelniony jest warunekN else ciag instrukcji wykonywanych gdy nie jest spelniony zaden z warunk\ ow od 1 do N end Podobnie jak dla instrukcjiwhilewarunekjest wyra|eniem zwracajcym warto[ logiczn. 3.2.2 Instrukcjaselect-case(*) Spójrzmy na przykBad poni|ej (prosz sprawdzi jego dziaBanie dla ró|nych warto[ci zmiennejnum)13 -->num = 1, select num, case 1, y =  przypadek 1 , case 2, y =  przypadek 2 ,... -->else, y =  inne przypadki  , end który w skrypcie lub funkcji przyjmie raczej posta // zakladamy, ze zmienna num zostala dobrze okreslona select num case 1 y =  przypadek 1 case 2 y =  przypadek 2 else y =  inne przypadki end Scilab sukcesywnie porównuje warto[ zmiennej num z mo|liwymi warto[ciam przypadków (1, potem 2). Je[li zajdzie rów- no[, wykonywane s instrukcj poczwszy odcaserównego zmiennejnuma| do koDca instrukcjiselect-case(tutu jak dziaBa ten select-case??). Przekazanie sterowania do nieobowizkowej gaBzielsema miejsce w sytuacji gdy nie powiod si wszystkie wcze[niejsze testy. Formalnie instrukcja ta przedstawia si nastpujco select zmiennaTestjaca case wyrazenie1 ciag instrukcji wykonywany gdy zmiennaTestuaca rowna jest wyrazeniu1 ... case wyrazenieN ciag instrukcji wykonywany gdy zmiennaTestujaca rowna jest wyrazeniuN else ciag instrukcji wykonywany gdy zmiennaTestujaca nie jest rowna zadnemu z wyrazen od 1 do N end gdzie wyrazeniei jest wyra|eniem zwracajcym warto[, która mo|e zosta porównana ze zmiennaTestujaca. Kon- strukcjaselect-caserównowa|na jest nastpujcej konstrukcjiif-then-else if zmiennaTestujaca == wyrazenie1 then ciag instrukcji wykonywany gdy 13 Zmiennayjest BaDcuchem znaków  patrz kolejny rozdziaB. 26 zmiennaTestujaca rowna jest wyrazeniu1 ... elseif zmiennaTestujaca = wyrazenieN then ciag instrukcji wykonywany gdy zmiennaTestujaca rowna jest wyrazeniuN else ciag instrukcji wykonywany gdy nie jest speniony zaden z powyzszych warunkow end 3.3 Inne rodzaje zmiennych Do tej pory u|ywali[my zmiennych bdcych macierzami (wektorami, skalarami) zBo|onymi z warto[ci typu rzeczywistego, zespolonego lub logicznego. Oprócz nich mamy równie| do dyspozycji zmienne reprezentujce BaDcuchy znaków oraz listy. 3.3.1 AaDcuchy zanków W przykBadzie dotyczcym konstrukcji select-case zmienna y jest typu BaDcuch znaków. W Scilabie BaDcuchy znaków ograniczane s za pomoc znaku apostrofu lub cudzysBowu (je[li chcemy który[ z tych znaków wykorzysta w takim BaDcuchu, to nale|y napisa go dwukrotnie). Aby przypisa zmiennejczyToPrawdawarto[ Czy Scilab jest "cool" ? nale|y napisa -->czyToPrawda = "Czy Scilab jest ""cool"" ?" lub równowa|nie -->czyToPrawda =  Czy Scilab jest ""cool"" ? Mo|na tak|e zdefiniowa macierz zBo|on z BaDcuchów znaków -->M = ["a" "bc" "def"] M = !a bc def ! -->size(M) // w celu otrzymania wymiarow skladowych ans = ! 1. 3. ! -->length(M) ans = ! 1. 2. 3. ! Zauwa|my, |e w tym przypadkulengthnie zachowuje si tak samo jak dla macierzy liczbowej. Zwraca bowiem jako wynik macierz o takim samym wymiarze jak M, w której wspóBczynnik okre[lony przez wspóBrzdne (i, j) równy jest co do warto[ci ilo[ci znaków w BaDcuchu na pozycji (i, j). Do Bczenia (konkatenacji) BaDcuchów u|ywamy operatora+ -->s1 =  abc ; s2 =  def ; s = s1 + s2 s = abcdef podcig wydobywamy za[ (operacja extrakcja) przy pomocy funkcjipart -->part(s,3) ans = c -->part(s,3:4) ans = cd 27 Drugim argumentem funkcjipartjest wektor okre[lajcy indeksy znaków jakie mamy wydoby z BaDcucha. 3.3.2 Listy (*) Lista jest uporzdkowanym zbiorem obiektów Scilaba (macierzy i skalarów rzeczywistych, zespolonych czy te| zBo|onych z BDcuchów znaków, warto[ci logicznych, funkcji, list, . . . ; ka|demu elementowi przypisuje si numer). Dostpne s dwa rodzaje list: lista obiektow oraz lista typów. Oto przykBad listy obiektów : -->L=list(rand(2,2),["Obym skonczyl" " ten skrypt..."],[%t ; %f]) L = L(1) ! 0.2113249 0.0002211 ! ! 0.7560439 0.3303271 ! L(2) !Obym skonczyl ten skrypt... ! L(3) ! T ! ! F ! W ten oto sposób zdefiniowali[my list, w której pierwszy element jest macierz o wymiarze (2, 2), drugi wektorem BaDcuchów znakowych a trzeci wektorem o skBadnikach logicznych. Oto kilka podstawowych operacji na listach -->M = L(1) // extrakcja (wydobycie) pierwszego sk\l{}adnika M = ! 0.2113249 0.0002211 ! ! 0.7560439 0.3303271 ! -->L(1)(2,2) = 100; // modyfikowanie pierwszego skladnika -->L(1) ans = ! 0.2113249 0.0002211 ! ! 0.7560439 100. ! -->L(2)(4) = " przed wakacjami !"; // modyfikacja drugiego skladnika -->L(2) ans = !Obym skonczyl ten skrypt... przed wakacjami ! ! -->L(4)="Aby dodac czwarty skladnik" L = L(1) ! 0.2113249 0.0002211 ! ! 0.7560439 100. ! L(2) 28 !Obym skonczyl ten skrypt... przed wakacjami ! ! L(3) ! T ! ! F ! L(4) Aby dodac czwarty skladnik -->size(L) // ile elementow ma lista ans = 4. -->length(L) // podobnie jak wyzej ans = 4. -->L(2) = null() // usuniecie drugiego elementu L = L(1) ! 0.2113249 0.0002211 ! ! 0.7560439 100. ! L(2) ! T ! ! F ! L(3) Aby dodac czwarty skladnik -->Lbis=list(1,1:3) // definicja innej listy Lbis = Lbis(1) 1. Lbis(2) ! 1. 2. 3. ! -->L(3) = Lbis // 3 skladnik L jest teraz lista L = L(1) 29 ! 0.2113249 0.0002211 ! ! 0.7560439 100. ! L(2) ! T ! ! F ! L(3) L(3)(1) 1. L(3)(2) ! 1. 2. 3. ! Przejdzmy teraz do list typów. W tych listach pierwszy element jest BaDcuchem okre[lajcym typ elementów listy (co pozwala zdefiniowa nowy typ danych a nastpnie okre[li operatory na tym typie danych), kolejne elementy mog by dowolnymi obiek- tami Scilaba. Pierwszy element mo|e by tak|e wektorem zBo|onym z BaDcuchów znaków; pierwszy okre[la typ elementów listy, pozostaBe za[ sBu| jako indeksy elementów listy (w miejsce ich numerów na li[cie). Rozwa|my nastpujacy przykBad: chcemy reprezentowa wielo[cian (w którym wszystkie [ciany maj identyczn ilo[ krawdzi). W tym celu przechowujemy w macie- rzy wspóBrzdne wszystkich wierzchoBków (o wymiarze (3,ilo[WierzchoBków) do której odwoBywa si bdziemy u|ywajc ciagu coord. Nastpnie definiujemy macierz (o wymiarze (ilo[WierzchoBkówZciany, ilo[Zcian)) okre[lajc zwizek pomi- dzy [cian a wierzchoBkami: dla ka|dej sciany podajemy numery wierzchoBków j opisujcych w kolejno[ci zodnej z reguB korkocigu dla zewntrznego wektora normalnego. Do tej listy odwoBywa si bdziemy piszc fence. Dla sze[cianu (porównaj z 5 6 8 7 1 y 2 4 3 x Rysunek 2: Numéros des sommets du cube fig:cube rysunej 2) bdzie to wygldaBo jak poni|ej P=[ 0 0 1 1 0 0 1 1;... // wspolrzedne wierzcholkow 0 1 1 0 0 1 1 0;... 0 0 0 0 1 1 1 1]; connect = [1 5 3 2 1 1;... // zaleznosci pomiedzy wierzcholkami i scianami 2 8 7 6 5 4;... 30 3 7 8 7 6 8;... 4 6 4 3 2 5]; Pozostaje zdefiniowa list tutu typee, która zawiera bdzie wszystkie informacje o sze[cianie -->Cube = tlist(["polyedre","coord","fence"],P,connect) Cube = Cube(1) !polyedre coord fence ! Cube(2) ! 0. 0. 1. 1. 0. 0. 1. 1. ! ! 0. 1. 1. 0. 0. 1. 1. 0. ! ! 0. 0. 0. 0. 1. 1. 1. 1. ! Cube(3) ! 1. 5. 3. 2. 1. 1. ! ! 2. 8. 7. 6. 5. 4. ! ! 3. 7. 8. 7. 6. 8. ! ! 4. 6. 4. 3. 2. 5. ! W celu wskazania pewnego elementu z wykorzystniem jego numeru, mo|na wykorzysta odpowiadajcy jemy BaDcuch znaków: -->Cube.coord(:,2) // lub inaczej Cube("coord")(:,2) ans = ! 0. ! ! 1. ! ! 0. ! -->Cube.fence(:,1) ans = ! 1. ! ! 2. ! ! 3. ! ! 4. ! Pomijajc te szczególne mo|liwo[ci, listami typów manipuluje si tak samo jak innymi. Ich zalet jest mo|liwo[ samodzielnego zdefiniowania (rozszerzenia) operatorów+, -, / itd. na list danych tego typu (porównaj przyszB wersj tej dokumentacji, *, albo lepiej stron Scilaba z podobnymi informacjami. 3.3.3 Kilka przykBadów z wektorami i macierzami logicznymi Typ boolowski (u|ywany do reprezentacji warto[ci logicznych typu prawda i faBsz) okazuje si czsto przydatny ze wzgldów praktyznych podczas manipuowania macierzami. Gdy porównujemy dwie macierze tego samego wymiaru za pomoc jednego z operatorów porównania (<,>,<=,>=,==,~=) jako wynik otrzymujemy macierz o warto[ciach logicznych, równie| tego samego formatu, której zawarto[ jest wynikiem porównania element po elemencie odpowiadajcych sobie elementów porównywanych macierzy -->A = rand(2,3) A = ! 0.2113249 0.0002211 0.6653811 ! ! 0.7560439 0.3303271 0.6283918 ! -->B = rand(2,3) B = 31 ! 0.8497452 0.8782165 0.5608486 ! ! 0.6857310 0.0683740 0.6623569 ! -->A < B ans = ! T T F ! ! F F T ! Je[li natomiast porównujemy macierz z liczb, a wicA<sgdziesjest skalarem to faktycznie dokonujemy porównaniaA<s*one(A) -->A < 0.5 ans = ! T T F ! ! F T F ! Operatory logiczne równie| stosuje si dla macierzy na zasadzie element po elemencie -->b1 = [%t %f %t] b1 = ! T F T ! -->b2 = [%f %f %t] b2 = ! F F T ! -->b1 & b2 ans = ! F F T ! -->b1 | b2 ans = ! T F T ! -->~b1 ans = ! F T F ! Ponadto istniej dwie przydatne funkcje pozwalajce dokona wektoryzacji testu 1. Funkcja bool2s przeksztaBca macierz boolowsk w macierz o takich samych wymiarach przeksztaBcajc warto[ lo- giczna prawda na 1, faBsz za[ na 0 -->bool2s(b1) ans = ! 1. 0. 1. ! 2. Funkcjafind zwraca wspóBrzdne warto[ci prawda w macierzy boolowskiej -->v = rand(1,5) v = ! 0.5664249 0.4826472 0.3321719 0.5935095 0.5015342 ! -->find(v < 0.5) // v < 0.5 daje wektor boolowski ans = ! 2. 3. ! Stosowanie tych funkcji jest dosy kosztowne czasowo (porównaj "Kilka uwag zwizanych z szybko[ci wykonania"). Funkcje andorazoroznaczajce odpowiednio iloczyn logiczny (i) oraz sum logiczn (lub) wszystkich elementów macierzy boolow- skiej dajc w efekcie jedn warto[ logiczn. Funkcje to mog przyjmowa dodatkowy argument okre[lajcy czy odpowiednie dziaBanie logiczne ma zosta wykonane na wierszach czy kolumnach macierzy. 32 3.4 Funkcje W celu zdefiniowania funkcji w Scilabie, najbardziej odpowiedni metod jest zapisanie jej w utoworzonym pliku; bdzie mo|na do niego dopisywa kolejne funkcje grupujc na przykBad te zwizane z tym samym czy podobnym zagadnieniem czy aplikacj. Ka|da funkcja powinna rozpoczyna si w nastpujcy sposób function [y1,y2,y3,...,yn]=nazwaFunkcji(x1,....,xm) gdziexis argumentami funkcji, natomiastyiwarto[ciami przez ni zwracanymi. Dalej wystpuje ciaBo funkcji czyli wszyst- kie instrukcje niezbdne do wykonania przez funkcj okre[lonego zadania. Funkcja powinna koDczyc si sBowem kluczowym endfunction. Oto pierwszy przykBad function [y] = silnia1(n) // silnia; zakladamy, ze n jest liczba naturalna y = prod(1:n) endfunction ZaBó|my, |e powy|sz funkcj zapisali[my w pliku o nazwiesilnia1.sci14. Aby mo|na z niej korzysta w Scilabie nale|y j najpierw wczyta getf("silnia1.sci") // lub inaczej exec("silnia1.sci") Mo|na tego dokona tak|e z menuFile operationstak jak dla skryptów. Dopiero teraz mo|emy u|y zdefiniowanych w danym pliku funkcji, zarówno w  samodzielnych poleceniach jak i w skryptach czy innych funkcjach. -->m = silnia1(5) m = 120. -->n1=2; n2 =3; silnia1(n2) ans = 6. -->silnia1(n1*n2) ans = 720. Zanim przejdziemy do kolejnych przykBadów musimy ustali jeszcze pewn terminologi. Warto[ci zwracane przez funkcj (yi) oraz te przez ni zwracane (xi) nazywa bdziemy argumentami formalnymi. U|ywajc funkcji na przykBad w skrypcie czy innej funkcji argS = silnia1(argE) u|yte argumenty (argSorazargE) nazywa bdzimy argumentami efektywnymi. I tak w pierwszym przykBadzie argumentem wej[ciowym jest staBa o warto[ci 5, w drugim zmiennan2, natomiast w trzecim wyra|enien1*n2. Zwizek midzy argumentem efektywnym a formalnym mo|e obiawia si na ró|ne sposoby (por. nastpny paragraf, gdzie precyzuje si te zagadnienia przy u|yciu [rodowiska Scilaba). ce qui concerne Scilab). Drugi przykBad pokazuje funkcj znajdujc pierwiastki trójmianu kwadratowego trinom function [x1,x2] = trinom(a, b, c) // oblicza pierwiastki trojmianu kwadratowego postaci x^2 + b x + c = 0 // a, b oraz c musza byc liczbami rzeczywistymi lub zespolonymi roznymi od zera delta = b^2 - 4*a*c x1 = (-b - sqrt(delta))/(2*a) x2 = (-b + sqrt(delta))/(2*a) endfunction Oto wyniki trzech prób jej u|ycia -->[r1, r2] = trinom(1,2,1) r2 = - 1. 14 Tradycyjnie plik zawierajacy funkcj posiada rozszerznie .scinatomiast skrypt.sce Û 33 r1 = - 1. -->[r1, r2] = trinom(1,0,1) r2 = i r1 = - i -->trinom(1,0,1) // wywolanie bez przypisania do zmiennej ans = - i Zauwa|my, |e w trzecim przypadku nie widzimy drugiego z pierwiastków. Jest to normalne zachowanie w sytuacji gdy funkcja zwraca wicej ni| jedn warto[ i nie zostan one przypisane do zmiennych (tak jak ma to miejsce w drugim przypadku). W takiej sytuacji Scilab wykorzystuje domy[ln zmienn ans aby przechowywa wynik; ans bdzie przyjmowaBo kolejen zwracane warto[ci aby ostatecznie przyj warto[ równ tej zwróconej jako ostaniej. Teraz trzeci przykBad. Nale|y rozwin w punkcie t wielomian zapisany w bazie Newtona (Uwaga: Dla xi = 0 otrzymujemy baz kanoniczn.) p(t) = c1 + c2(t - x1) + c3(t - x1)(t - x2) + · · · + cn(t - x1) . . . (t - xn-1). U|ywajc rozkBadu na czynniki i wykonujc obliczenia od prawej do lewej (tutaj dla n = 4) p(t) = c1 + (t - x1)(c2 + (t - x2)(c3 + (t - x3)(c4))), otrzymujemy algorytm Hornera (1) p := c4 (2) p := c3 + (t - x3)p (3) p := c2 + (t - x2)p (4) p := c1 + (t - x1)p. Uogólniajc to na dowolne n otrzymujemy nastpujc funkcj Scilaba function [p]=myhorner(t,x,c) // rozwija wielomian c(1) + c(2)*(t-x(1)) + c(3)*(t-x(1))*(t-x(2)) + // ... + c(n)*(t-x(1))*...*(t-x(n-1)) // zgodnie z algorytmem Hornera n=length(c) p=c(n) for k=n-1:-1:1 p=c(k)+(t-x(k))*p end endfunction Je[li wektorycoef,xx,ttzostaBy dobrze okre[lone instrukcja val = myhorner(tt,xx,coef) spowoduje przypisanie dovalwarto[ci równej m-1 coef1 + coef2(tt - xx1) + · · · + coefm (tt - xxi) i=1 MaBe przypomnienie: instrukcjalengthzwraca iloczyn dwóch wymiarów macierzy dajc liczb elementów co w przypadku wektora (kolumnowego lub wierszowego) równowa|ne jest jego dBugo[ci. Instrukcja ta, wraz z instrukcjsizezwracajc ilo[ wierszy i ilo[ kolumn, pozwala zrezygnowa z przekazywania do funkcji informacji o wymiarze przekazywanych obiektów. 34 3.4.1 Przekazywanie parametrów (*) Przekazywanie parametrów odbywa si przez referencj je[li wewntrz funkcji nie s one modyfikowane oraz przez warto[15 w przeciwnym razie (to znaczy, parametry wej[ciowe funkcji nie mog by modyfikowane). Uwzgldnienie tego faktu w waszych programach mo|e spowodowa w znacznym stopniu przyspieszenie pewnych funkcji. Poni|ej przedstawiamy sztuczny, ale dobrze ilustrujcy zagadnienie, przykBad ukazujcy koszt zwizany z przekazywaniem parametrów przez warto[. function [w] = toto(v, k) w = 2*v(k) endfunction function [w] = titi(v, k) v(k) = 2*v(k) w = v(k) endfunction // skrypt testujacy m = 200000; ilPowtorzen = 2000; x = rand(m,1); timer(); for i=1:nb_rep; y = toto(x,300); end; t1=timer()/ilPowtorzen timer(); for i=1:nb_rep; y = titi(x,300); end; t2=timer()/ilPowtorzen Otrzymane wyniki bd ró|niBy si w zale|no[ci od maszyny na jakiej dziaBa Scilab; my otrzymali[my t1 = 0.00002 t2 = 0.00380 Wraz z zakoDczeniem funkcji destrukcji ulegaj wszystkie jej zmienne wewntrzne. 3.4.2 Wstrzymywanie funkcji Podczas debugowanie funkcji przydatna okazuje si funkcjadisp(v1, v2, ...) pozwalajca wy[wietli warto[ci zmien- nych v1, v2, . . . Nale|y mie na uwadze, |e funkcjadisp wy[wietla warto[ci swoich argumentów w odwrotnej kolejno[ci. W celu chwilowego wstrzymania wykonania programu mo|na u|y funkcjipause; po jej napotkaniu mamy mo|liwo[ spraw- dzenia warto[ci wszystkich zdefiniowanych do tej pory zmiennych (znak zachty-->Scilaba zmieni si na-1->). Polecenie resumepowoduje wznowienie przebiegu obliczeD. Istnieje tak|e inny sposób wskazywania miejsca wstrzymania programu (tak zwany break point, ang.) ni| dodawanie in- strukcjipause setbpt(NazwaFunkcji[, numerLinii]) Argumentem funkcjisetbptjestNazwaFunkcji, która ma by wstrzymywana oraz opcjonalny numer wiersza, w którym ma to nastpi (domysln warto[ci jest 1). Po napotkaniu przez Scilaba wskazanego miejsca zatrzymania, dalej wszystko odbywa si tak jak gdyby zastosowano instrukcjpausepo instrukcji wskazanej przeznumerLinii. Punkt zatrzymania mo|e zosta usnity przez delbpt(NazwaFunkcji[, numerLinii]) Je[li nie zostanie wskazana |adna linia, wszystkie punkty zatrzymania zostan usunite. Zdefiniowane do tej pory punkty zatrzy- mania mo|na zobaczyc korzystajc z dispbpt (). Oto przykBad zwizany z wcze[niej zdefiniowan funcjtrinom(patrz strona 33). ZaBó|my, |e wczytali[my uprzednio t funkcj czy to za pomoc funkcjigetfczy te|exec. -->setbpt("trinom") // pierwszy punkt zatrzymania -->[r1,r2] = trinom(1,2,7) // rozpoczynamy obliczenia => zatrzymanie Stop after row 1 in function trinom : -1->a // sprawdzamy wartosci zmiennych a = 1. -1->b 15 Przez warto[ czyli we wntrzu funkcji pracujemy na kopii przekazywanego argumentu. Przez referencj czyli we wntrzu funkcji pracujemy na oryginal- nym argumecie. (przyp. tBum.) 35 b = 2. -1->dispbpt() // sprawdzamy punkty zatrzymania breakpoints of function :trinom 1 -1->setbpt("trinom",5) // dodajemy kolejny punkt zatrzymania -1->x1 // x1 nie zostal jeszcze zdefiniowany x1 !--error 4 undefined variable : x1 -1->resume // wznawiamy wykonani az do nastepnego punktu zatrzymania Stop after row 5 in function trinom : -1->x1 // teraz mozemy sprawdzic jaka wartosc ma x1 x1 = - 1. - 2.4494897i -1->x2 // x2 nie zostal jeszcze zdefiniowany (nastapi to dopiero w linii 6) x2 !--error 4 undefined variable : x2 -1->dispbpt() // sprawdzamy punkty zatrzymania breakpoints of function :trinom 1 5 -1->resume // wznawiamy wykonanie funkcji r2 = - 1. + 2.4494897i r1 = - 1. - 2.4494897i -->delbpt("trinom") // usuwamy wszystkie punkty zatrzymania 4 Grafika W zakresie tworzenia grafiki Scilab posiada wiele mo|liwo[ci zarówno je[li bra pod uwag operacje niskopoziomowe jak te| bardziej zlo|one funkcje pozwalajce operowa skomplikowanymi obiektami. W dalszym cigu zostanie wyja[niona jedynie maBa cz[ dostpnych mo|liwo[ci. Uwaga: Dla osób znajcych instrukcje graficzne MATLABA Stephane Mottelet napisaBa (tutu czy napisal) bibliotek funkcji graficznych wywoBywanych w stylu MATLABA; dostpna jest ona pod adresem http://www.dma.utc.fr/~mottelet/scilab/ 4.1 Okna graficzne WywoBanie instrukcji graficznej takiej jak plot czy plot3d powoduje umieszczenie rysowanego obrazu w oknie raficznym o numerze 0. WywoBanie funkcji graficznej powoduje na ogóB powstanie nowego obrazu na ju| istniejcym, std potrzeba uprzedniego wymazania zawarto[ci okna przy pomocy funkcjixbasc(). Wymazania zawarto[ci okna mo|na tak|e dokona wybierajc z menuFileopcjclear. Manipulowanie oknami mo|liwe jest dzieki nastpujcym funkcjom 36 xset("window",num) okno o numerzenumstaje si oknem bie|cym; je[li |adne okno nie istniaBo to zostanie utworzone xselect() uaktywnia bie|ce okno; je[li |adne okno nie istniaBo to zostanie utworzone xbasc([num]) wymazuje zawarto[ okna o numerzenum; je[li zostanie pominity numer to wymazana zostanie zawarto[ bie|cego okna xdel([num]) powoduje zamknicie okna o numerzenum; je[li zostanie pominity numer to zamknite zostanie bie|ce okno Jako generaln zasad mo|na przyj, |e po wybraniu bie|cego okna (przezxset("window",num)) za pomoca instrukcji xset("nom",a1,a2,...)mo|emy dokonywa zmian w parametrach zwizanych z oknem. nom okre[la nazw parame- tru, który chcemy zmieni jak na przykBadthicknessgdy chcemy dostosowa grybo[ strzaBek czycolormapgdy chcemy zmieni u|ywan palet kolorów. Za nazw podajemy jeden lub wicej argumentów wmaganych do ustawienia nowej warto[ci parmetru. Zbiór tych parametrów nazywany jest kontekstem graficznym; ka|de okno posiada swój taki kontekst. W celu zdoby- cia wikszej ilo[ informacji na temat parametów (których jest caBkiem pokazny zbiór) odsyBsamy doHelp-u i tematuGraphic Library. W wikszo[ci sytuacji mog one by zmieniane interaktywnie przez menu graficzne, ukazujce si po wydaniu po- leceniaxset(). (Uwaga: W tym menu dostpne jest równie| podokno opisujce kolory; nie jest jenak mo|liwe wprowadzanie zmian do tego| opisu. Rodzina funkcji[a1,a2,...]=xget( nazwa )pozwala otrzyma parametry zwizane z pewnym kontekstem graficznym. 4.2 Wprowadzenie doplot2 Wcze[nie mieli[my ju| do czynienia z prosta instrukcj plot. Je[li jednak zamierzamy wykre[lic wicej krzywych lepiej stosowa funkcjplot2d. Najpro[ciej mo|na j u|y jak to pokazano poni|ej x=linspace(-1,1,61) ; // odci\k{e}te (wektor kolumnowy) y = x.^2; // rz\k{e}dne (r\ ownie\.z wektor kolumnowy) plot2d(x,y) // --> il lui faut des vecteurs colonnes ! tutu DoBó|my teraz inn krzyw . . . ybis = 1 - x.^2; plot2d(x,ybis) xtitle("Krzywe...") // dodajemy tytu\l{} . . . i jeszcze jedn, która tutu jest w innej skali16 ni| poprzednie yter = 2*y; plot2d(x,yter) Zauwa|my, |e Scilab przeskalowaB okno tak aby zmie[ciBa si w nim trzecia krzywa, ale tak|e odrysowaB dwie poprzednie krzywe w nowej skali (wydaje si to naturalne, ale mechanizm ten nie byB obecny a| do wersji 2.6 wprowadzajc czesto w bBad). Mo|liwe jest wykre[lenie tych trzech krzywych za jednym razem xbasc() // aby wyczy\ scic obszar rysowania plot2d(x,[y ybis yter]) // konkatenacja macierzy xtitle("Krzywe...","x","y") // tytu\l{} i nazwy dla obu osi Otrzymujemy w ten sposób wykres podobny do tego przedstawionego na rysunku 3. W celu jednoczesnego wykre[lenia kilku krzywych, instrukcja przyjmuje posta plot2d(Mx,My) gdzie Mx oraz My s macierzami o identycznym wymiarze, ilo[ kolumnncjest równa ilo[ci krzywych a i-ta krzywa jest otrzymywana na podstawie wektorówMx(:,i)(odcite) iMy(:,i) (rzdne). W sytuacji gdy wektor odcitych jest taki sam dla wszystkich krzywych (jak w naszym przypadku), mo|na poda go jeden raz zamiast powtarza nc razy (plot2d(x,My)zamiastplot2d([x x .. x],My)). 4.3 plot2dz argumentami opcjonalnymi Ogólna skBadnia przedstawia si nastpujco plot2d(Mx,My <,opt_arg>*) 16 Pod pojciem skali rozumiemy tutaj obszar w jakim zostanie wykre[lona krzywa oraz ewentualne dodatkowe cechy. 37 Courbes... y 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 x 0 -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 Rysunek 3: Les fonctions x2, 1 - x2 et 2x2 fig:1 gdzie<,opt_arg>*oznacza opcjonalny cig dodatkowych argumentów postaci17 sBowoKluczowe=warto[ podanych w dowolnej kolejno[ci. Poka|emy podstawowe epcje przy pomocy kilku przykBadów 1. Wybór koloru i definiowanie legendy W poprzednim przykBadzie mo|na byBo zaobserwowa, |e Scilab wykre[liB 3 krzywe przy u|yciu 3 ró|nych (pierwszych) kolorów okre[lonch w domy[lnej palecie kolorów 1 czarny 5 czerwony 23 fioletowy 2 niebieski 6 fioBkowo-ró|owy 26 brzowy 3 jasno-zielony 13 ciemno-zielony 29 ró|owy 4 cyan 16 jasno-niebieski tutu 32 jaune orangé 18 Instrukcjaxset() pozwala na ich zmian . Celem wybrania koloru u|ywamy zapisu styl=vect, gdzievect jest wektorem liniowym zawierajcym numery kolorów dla ka|dej krzywej. Legend otrzymujemy piszc leg=str, gdzie strjest BaDcuchem znaków postacileg1@leg2@... w którymlegijest podpisem dlai-tej krzywej. Oto przykBad (porówaj rysunek (4)) : x = linspace(0,15,200) ; y = besselj(1:5,x); xbasc() plot2d(x, y, style=[2 3 4 5 6], leg="J1@J2@J3@J4@J5") xtitle("Funkcje Bessela J1, J2,...","x","y") Poniewa| kolejno[ argumentów opcjonalnych nie jest istotna równie dobrze mo|na napisa plot2d(x, y, leg="J1@J2@J3@J4@J5", style=[2 3 4 5 6]) 2. Wkresy zawierajce symbole Czasami, gdy wa|ne jest wyrazne wskazanie punktów przez które wykres przechodzi, przydatne mog okaza si znaczniki tutu. Mo|emy wybra jeden z kilku rodzajów zancznika przedstawionych poni|ej styl 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 symbol · + × •" f& ³% ½% c& © 17 argumentFormalny = argumentEfektywny 18 tutu 38 Les fonctions de Bessel J1, J2,... y 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 x -0.4 0 2 4 6 8 10 12 14 16 J1 J4 J2 J5 J3 Rysunek 4: Wybrany styl i legenda fig:style Dla przykBadu spróbujmy (porównaj rysunek (5)) : x = linspace(0,2*%pi,40) ; y = sin(x); yp = sin(x) + 0.1*rand(x,"normal"); xbasc() plot2d(x, [yp y], style=[-2 2], leg="y=sin(x)+perturbation@y=sin(x)") 3. Okre[lanie skali Oto przykBad gdzie niezbdne jest narzucenie skali izometrycznej tutu bowiem chcemy narysowa okrg (patrz rysunek (6)). Dodatkowy parametr bdzie postaci frameflag=val, gdzieval powinna przyj warto[ 4 aby otrzyma skale izometryczn (dobran na podstawie warto[ci minimalnej i naksymalnej) t = linspace(0,2*%pi,60) ; x1 = 2*cos(t); y1 = sin(t); // elipsa x2 = cos(t); y2 = y1; // okr\k{a}g x3 = linspace(-2,2,60) ; y3 = erf(x3); // funkcja b\l{}\k{e}du legenda = "elipsa@okr\k{a}g@funkcja b\l{}\k{e}du" plot2d([x1 x2 x3],[y1 y2 y3],style=[1 2 3], frameflag=4, leg=legenda) xtitle("Znowu krzywe...","x","y") W pewnych przypadkach frameflag towarzyszy powinien parametr rect. Je[li na przykBad chcemy samodzielnie okre[li obszar rysowania (zwalniajc z tego punkcjplot2dtutu) musimy napisa rect = [xmin, ymin, xmax, ymax] w poBczeniu zframeflag=1jak w przykBadzie poni|ej x = linspace(-5,5,200) ; y = 1 ./(1 + x.^2); xbasc() plot2d(x, y, frameflag=1, rect=[-6,0,6,1.1]) xtitle("Funkcja Rungego tutu") Oto wszystkie mo|liwe warto[ci dlaframeflag 39 1.5 × 1.1 × × × × × × × × × 0.7 × × × × × × × 0.3 × × × × × -0.1 × × × × -0.5 × × × × × × × × × -0.9 × × × × × -1.3 0 1 2 3 4 5 6 7 × y=sin(x)+perturbation y=sin(x) Rysunek 5: dessin en trait plain et avec des symboles non reliés tutu fig:symbol frameflag=0 u|ycie wcze[niej zdefiniowanej skali (lub domy[lnej) frameflag=1 skala zadana przez kwadrat frameflag=2 skala obliczona jako maksimum i minimum zMxiMy frameflag=3 échelle isométrique calculée en fonction de rect tutu frameflag=4 skala izometryczna obliczona jako maksimum i minimum zMxetMy frameflag=5 identycznie 1 ale z ewntualnym dostosowaniem graduation tutu frameflag=6 identycznie 2 mais avec adaptation eventuelle pour la graduation frameflag=7 identycznie 1 ale wcze[niejsze krzywe s odrysowywane frameflag=8 identycznie 2 ale wcze[niejsze krzywe s odrysowywane Uwaga: W przypadkach 4 i 5 mo|e zaj[ (ewentualna) modyfikacja obszaru rysowania w taki sposób aby stopniowanie (tutu) (które jest zawsze tutu) tutu (tak zwane tutu). 4. Precyzowanie umieszczenis osi Rozmieszczenie osi mo|emy okre[la piszcaxesflag=val. W kolenym przykBadzie (patrz rysunek (7))valma warto[ 5 co powoduje, |e osie przetna si wpunkcie (0, 0)19 bez tutu boite englobante : x = linspace(-14,14,300) ; y = sinc(x); xbasc() plot2d(x, y, style=2, axesflag=5) xtitle("Funkcja sinc") Oto tabela podajca wszytkie mo|liwo[ci: axesflag=0 bez pudeBka, osi oraz jednostek axesflag=1 z pudeBkiem, osiami i jednostkami (x na dole, y po lewej) axesflag=2 z pudelkiem, ale bez osi i jednostek axesflag=3 z pudeBkiem, osiami i jednostkami (x na dole, y po lewej) axesflag=4 bez pudeBka, ale z osiami i jednostkami (punkt przecicia w [rodku) axesflag=5 bez pudeBka, ale z osiami i jednostkami (punkt przecicia dla x = 0 i y = 0) 19 Gdy punkt (0, 0) jest we wntrzu obszaru rysowania. 40 Encore des courbes ... y 1.413 1.010 0.606 0.202 -0.202 -0.606 -1.010 x -1.413 -2.000 -1.429 -0.857 -0.286 0.286 0.857 1.429 2.000 ellipse cercle fct erreur Rysunek 6: Elipsa, okrg i funkcja bBdu fig:echelle_iso 5. Skala logarytmiczna Odpowiedni parametr jest postacilogflag=str, gdziestr jest BaDcuchem zBo|onym z dówch znaków  pierwszy dotyczy osi OX, drugi OY i ka|dy przyjmuje warto[n(nie logarytmiczna) lubl(logarytmiczna). x = logspace(0,4,200) ; y = 1 ./x; xbasc() subplot(1,2,1) plot2d(x, y, style=2, logflag= "ln") xtitle("logflag=""ln""") subplot(1,2,2) plot2d(x, y, style=2, logflag= "ll") xtitle("logflag=""ll""") PrzykBad ten pokazuje tak|e jak przy pomocy funkcjisubplot(m,n,num)w jednym oknie umie[ci kilka (niezale|- nych) wykresów. Parametrminformuje o podziale okna w pionie namrównych cz[ci,ndecyduje o podziale w poziomie numjest natomiast kolejnym numerem okna spo[ród m × n okien. Okna numerowane s od lewej do prawej i od góry do doBu poczynajc od numeru 1. Std okno na pozycji (i, j) ma numer20 n × (i - 1) + j. Nic nie stoi na przeszkodzie aby modyfikowa siatk tutu celem dobranie najbardziej nam odpowiadajcego ukBadu wykresów. Na przykBad xbasc() subplot(1,2,1) titlepage("po lewej") subplot(3,2,2) titlepage(["po prawej";"powy\.zej"]) subplot(3,2,4) titlepage(["po prawej";"w centrum"]) subplot(3,2,6) titlepage(["po prawej";"poni\.zej"]) xselect() tutu okna w pionie na dwie cz[ci (lew i praw), okna po prawej podzielono w poziomie na trzy czsci. Funkcjsubplot mo|na rozumie jako dyrektyw pozwalajc wybra pewien podobszar okna graficznego. 20 A nie m × (i - 1) + j jak napisane jest w pomocy! 41 La fonction sinc 1.1 0.9 0.7 0.5 0.3 0.1 -14 -10 -6 -2 2 6 10 14 -0.1 -0.3 Rysunek 7: Umieszczenie osi otrzymane dlaaxesflag=5 fig:axesflag 6. SBowo kluczowe strf Pozwala ono zastpi zarazem frameflagi axesflag oraz, ze wzgldu na kompatybilno[ z wcze[niejszymi wersjamiplot2dzawiera tak|e flag okre[lajc czy ma ono zastosowanie do legendy, czy nie tutu. Podawana warto[ skBada si z trzech znakówxyz, gdzie x równe 0 (nie ma legendy) lub 1 (legenda tworzona jest przez podanie leg=val) ; y cyfra od 0 do 9, odpowiadajca warto[ci podawanej w frameflag ; z cyfra od 0 do 5, odpowiadajca warto[ci podawanej w axesflag. En fait il faut savoir l utiliser car les séquences optionnelles de nombreuses primitives de dessin ne disposent pas encore des mots clés frameflag i axesflag tutu. Dodatkowo jest to bardzo praktyczne je[li chcemy doda nowy wykres do ju| istniejcego bez zmieniania skali i ramki co mo|na osigna piszcstrf="000"(unikajc w ten sposób pisania frameflag=0, axesflag=0). 4.4 Inne wersjeplot2d:plot2d2,plot2d3 Zasadniczo u|ywa si ich jakplot2dtaka sama skBadnia, takie same argumenty opcjonalne. 1. plot2d2 tutu pozwala narysowa funkcj w oparciu o skalary: w miejsce wykresu prostego odcinka wprowadzamy punkty (xi, yi) et (xi+1, yi+1),plot2d2odcinek poziomy (od (xi, yi) do (xi+1, yi)) a nastpnie odcinek pionowy (od (xi+1, yi) do (xi+1, yi+1)). Oto przykBad (porównaj rysunek (9)): n = 10; x = (0:n) ; y = x; xbasc() plot2d2(x,y, style=2, frameflag=5, rect=[0,-1,n+1,n+1]) xtitle("plot2d2") 2. plot2d3 Rysuje diagram tutu batonow: dla ka|dego punktu (xi, yi) plot2d3rysuje jeden segment pionowy od (xi, 0) do (xi, yi) (porównaj rysunek (10)) : n = 6; x = (0:n) ; y = binomial(0.5,n) ; 42 logflag="ln" logflag="ll" 0 1.0 10 0.9 0.8 -1 10 0.7 0.6 -2 0.5 10 0.4 0.3 -3 10 0.2 0.1 -4 0 0 1 2 3 4 10 0 1 2 3 4 10 10 10 10 10 10 10 10 10 10 Rysunek 8: Wykorzystanie parametrulogflag fig:logflag xbasc() plot2d3(x,y, style=2, frameflag=5, rect=[-1,0,n+1,0.32]) xtitle("Prawdopodobienstwo dla prawa binomialnego tutu B(6,1/2)") 4.5 Rysowanie wikszej ilo[ci krzywych zBo|onych z ró|nej ilo[ci punktów Za pomocplot2di jego wariantów nie mo|na narysowa za jednym razem wikszej ilo[ci krzywych, które nie s dyskretyzo- wane na takiej samej ilo[ci przedziaBów (i chyba takich samych przedziaBów tutu) co pociga za soba konieczno[ kilkakrotnego u|ycia tej funkcji. Od wersji 2.6 mo|na oby si bez podawania skali bez obawy o przykre niespodzianki, poniewa| domy[l- nie (frameflag=8) wcze[niejsze krzywe s odrysowywane w przypadku zmiany skali. Dlatego je[li kto[ chce opanowa skal (chodzi tutaj chyba o to aby nie byla ona zmieniana automatycznie tutu), musi to zrobic przy pierwszym wywoBaniu a dla nastpnych u|ywaframeflag=021. Oto odpowiedni przykBad (porównaj rysunek (11)) x1 = linspace(0,1,61) ; x2 = linspace(0,1,31) ; x3 = linspace(0.1,0.9,12) ; y1 = x1.*(1-x1).*cos(2*%pi*x1); y2 = x2.*(1-x2); y3 = x3.*(1-x3) + 0.1*(rand(x3)-0.5); // identyczne jak y2, ale z zaburzeniem ymin = min([y1 ; y2 ; y3]); ymax = max([y1 ; y2 ; y3]); dy = (ymax - ymin)*0.05; // aby doda\ c margines tutu rect = [0,ymin - dy,1,ymax+dy]; // okno wykresu xbasc() // wyczyszczenie poprzednich wykres\ ow plot2d(x1, y1, style=1, frameflag=5, rect=rect) // 1-sze wywo\l{}anie, kt\ ore ustali skal\k{ plot2d(x2, y2, style=2, frameflag=0) // 2-ie i 3-ie wywo\l{}anie; plot2d(x3, y3, style=-1,frameflag=0) // u\.zywamy poprzedniej skali xtitle("Krzywe...","x","y") Uwaga: Wypróuj ten przykBad piszcframeflag=1w miejsceframeflag=5. Nie mo|na mie osobnej legendy dla ka|dej krzywej niemniej jest to mo|liwe do osignicia przy pomocy niewielkiej ilo[ci progamowania. 21 Metoda ta jest konieczna je[li u|ywamy skali iso-metrycznej. 43 plot2d2 11 9 7 5 3 1 -1 0 2 4 6 8 10 12 Rysunek 9: Wykorzystanieplot2d2 fig:plot2d2 4.6 Zabawa z kontekstem graficznym moze lepiej  praca ? :))) Je[li wyprobowali[my poprzednie przykBady bez wtpienia i z ochot bdziemy chcieli modyfikowa niekóre rzeczy jak rozmiar symboli, rozmiar czy styl u|tego fontu lub te| grubo[ strzaBek. 1. Fonty : Aby zmieni font nale|y u|y poni|szego wywoBania xset("font",font_id,fontsize_id) gdziefont_idifontsize_ids liczbami caBkowitymi odpowiadajcymi odpowiednio z rodzaj i wielko[ wybranego fontu. Bie|cy font mo|na otrzyma piszc f=xget("font") gdzie f jest wektorem, f(1) opisuje rodzaj fontu a f(2) jego wielko[. Mo|na prostozmienia/uzyska wielko[ za pomocxset("font size",size_id)ifontsize_id = xget("font size"). Oto te, które s teraz do- stpne Nazwa fontu Courier Symbol Times Times-Italic Times-Bold Times-Bold-Italic Identyfikator 0 1 2 3 4 5 Wielko[ 8 pts 10 pts 12 pts 14 pts 18 pts 24 pts Identyfikator 0 1 2 3 4 5 Uwagi: " Courier est U chasse fixe tutu; " font Symbol pozwala na u|ycie liter greckich (p odpowiada À, a  ±, etc...); " Times jest fontem domy[lnym a jego domy[lny rozmiar to 10 punktów. 2. rozmiar symboli: zmieniamy poleceniem xset("mark size",marksize_id) bie|cy otrzymujemy natomiast piszc 44 Probabilités pour la loi binomiale B(6,1/2) 0.32 0.28 0.24 0.20 0.16 0.12 0.08 0.04 0 -1 0 1 2 3 4 5 6 7 Rysunek 10: Wykorzystanie plot2d3 fig:plot2d3 marksize_id = xget("mark size") gdzie podobnie jak dla fontów rozmiar okre[lamy za pomoca cyfr od 0 do 5; 0 jest wielko[ci domy[ln. 3. grubo[ linii: zmieniamy / otrzymujemy piszc xset("thickness",thickness_id) thickness_id = xget("thickness") Grubo[ jest wielko[ci dodatni odpowiadajc liczbie pikseli (na grubo[) jak ma mie krzywa. Dokonujc zmiany grubo[ci kre[lonych linii wpBywamy tak|e, czy tego chcemy czy nie, na grubo[ kre[lonej ramki i skali na osiach. Wyj- [ciem w takiej sytuacji jest dwukrotne tworzenie rysunku. Za pierwszym razem pomijamy ramk i skal (axesflag=0). Nastpnie powracamy do normalnej grubo[ci i nie zmieniajc skali widocznego obszaru (frameflag=0) rysujemy co[ poza nim (na przykBad krzyw zredukowana do jednego punktu (-", -")): xset("thickness", 3) plot2d(x, y, axesflag=0, ...) xset("thickness", 1) plot2d(-%inf, -%inf, frameflag=0, ...) Zatem drugie wywoBanie sBu|y jedynie do narysowania ramki i skali. 4.7 Tworzenie histogramów Odpowiednia do tego celu funkcja scilaba nazywa sihistplota jej wywoBanie jest nastpujce histplot(n, X, <,opt_arg>*) gdzie " njest albo liczb caBkowit albo wektorem liniowym (w którym ni < ni+1): 1. w przypadku gdy n jest wektorem liniowym dane dzielone s na k klas Ci = [ni, ni+1[ (wektorn ma wic k + 1 skBadowych); 45 Des courbes... y 0.3 0.3 0.3 + + + + + + + + 0.2 0.2 0.2 + + 0.1 0.1 0.1 + + 0.0 0.0 0 -0.1 -0.1 -0.1 -0.2 -0.2 -0.2 x -0.3 -0.3 -0.3 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Rysunek 11: Jeszcze krzywe... fig:3 2. w przypadku gdynjest liczb caBkowit, dane dzielone s na n równoodlegBych klas ñø c1 = min(X), cn+1 = max(X) òø C1 = [c1, c2], Ci =]ci, ci+1], i = 2, ..., n, avec ci+1 = ci + "C óø "C = (cn+1 - c1)/n " Xwektor (liniowy lub kolumnowy) z danymi; " <,opt_arg>* cig opcjonalnych argumentów jak dla plot2d z dodatkow kombinacj normalization = val gdzie val jest staB (zmienn lub wyra|eniem) boolowskim (domy[lnie true). Kiedy histogram jest znormalizowany, tutu son intégrale vaut 1 et approche donc une densité (dans le cas contraire, la valeur d une plage correspond au nombre de composantes de X tombant dans cette plage). Plus précisemment, la plage de l histogramme correspondant U l intervalle Ci vaut donc (m étant le nombre de données, et "Ci = ni+1 - ni) : card {Xj "Ci} si normalization = vrai m"Ci card {Xj " Ci} si normalization = faux Poni|ej przedstawiamy maBy przykBad, tutu toujours avec la loi normale (patrz rysunek (12)) : X = rand(100000,1,"normal"); classes = linspace(-5,5,21); histplot(classes,X) // tutu on lui superpose le tracé de la densité de N(0,1) x = linspace(-5,5,60) ; y = exp(-x.^2/2)/sqrt(2*%pi); plot2d(x,y, style=2, frameflag=0, axesflag=0) 4.8 Zapisywanie grafiki w ró|nych formatach W celu zachowanie utworzonego rysunku wybieramy z menu File okna graficznego opcj Export a nastpnie wybieramy po|dany format. tutu - sprawdzi to Wikszo[ z nich opiera si na postscripcie... . Poczwszy od wersji 2.5 mamy tak|e mo|liwo[ eksportu grafiki do pliku typugif. 46 0.40 0.36 0.32 0.28 0.24 0.20 0.16 0.12 0.08 0.04 0 -5 -4 -3 -2 -1 0 1 2 3 4 5 Rysunek 12: Histogram próbki liczb losowych wybranych z rozkBadem N(0, 1) tutu fig:histo 4.9 Prosta animacja Realizacja animacji przy pomocy Scilaba jest do[ prosta, jako |e pozwala on na podwójne buforowanie dziki czemu unikamy efektu migotania, gdy| wykres tworzony jest  w ukryciu i dopiero potem pokazywany. Mamy do wyboru dwa  drivery majce wpByw na tworzenie wykresu na ekranie22 " Rec, który powoduje, |e wszystkie operacje graficzne zwizane s z oknem; jest to domy[lny drvier; " X11, tutu qui se contente simplement d afficher les graphiques (il est alors imposible de  zoomer ). Do celów tworzenia animacji najczsciej odpowiedni bdzie ten ostatni; wybieramy go piszc driver("X11")(w celu po- wrócenia do drivera domy[lnego piszemydriver({ec")). Przy podwójnym buforowaniu ka|dy kolejny obraz tworzony jest najpierw w pamici (powstaje wówczas tak zwana pixmapa) a dopiero pózniej przenoszony na ekran. Oto prosty schemat postpowania przy tworzeniu animacji driver("X11") // tutu pas d enregistrement des opérations graphiques xset("pixmap",1) // przej\ scie do trybu podw\ ojnego buforowania ........... // ewentualne instrukcje ustalaj\k{a}ce skal\k{e} for i=1:nb_dessins xset("wwpc") // wyczyszczenie pixmapy ....... ....... // tworzenie i-tego rysunku ....... xset("wshow") // przeniesienie rysunku na ekran end xset("pixmap",0) // powr\ ot do bezpo\ sredniego tworzenia rysunk\ ow na ekranie driver("Rec") // przej\ scie do domy\ slego drivera Uwaga: W powy|szym postpowaniu czyszczenie caBego obszaru rysowania przy pomocyxset("wwpc")nie czynno[ci wymagan. Zamiast tego mo|na u|y na przykBad funkcjixclea, co uchroni nas przed ka|dorazowym odrysowywaniem tych obszarów rysunku, które si nie zmieniaj. Aby tutu u|y technik maskowania nale|y za pomoc wywoBaniexset( alufunction ,nu zmieni funkcj sterujc wy[wietlaniem (gdzienumto liczba caBkowita zwizana z wybran funkcj); patrz Help i przykBady. Poni|ej przedstawiony zostanie przykBad animacji: poruszajcy si [rodek ci|ko[ci prostokta (o dBugo[ci L i szeroko[ci l) po okrgu o promieniu r i [rodku w punkcie (0, 0); prostokt dodatkowo obraca si wokóB swojego [rodka ci|ko[ci. tutu Il y a certain un nombre de détails qui se greffent autour du canevas général exposé ci-avant : " l ordreplot2dsert uniquement U régler l échelle (isométrique) pour les dessins ; 22 Oprócz  driverów pozwalajacych tworzy rysunki jakopostscript,figczygif. Û 47 " xset("background",1)impose la couleur 1 (du noir dans la carte des couleurs par défaut) comme couleur d arri re plan, mais il est recommandé d exécuter l instructionxbasr()pour mettre effectivement U jour la couleur du fond ; " le dessin consiste U appeler la fonctionxfpolysuivi dexpolypour dessiner le bord (avec ici 3 pixels ce qui est obtenu avecxset("thickness",3)) ; U chaque fois on change la couleur U utiliser avecxset(olor",num); " l instruction xset("default") repositionne le contexte graphique de la fentre U des valeurs par défaut ; ainsi la variablepixmapreprend la valeur 0,thicknessla valeur 1,backgroundsa valeur par défaut, etc... n = 4000; L = 0.6; l = 0.3; r = 0.7; nb_tours = 4; t = linspace(0,nb_tours*2*%pi,n) ; xg = r*cos(t); yg = r*sin(t); xy = [-L/2 L/2 L/2 -L/2;... // 4 punkty graniczne -l/2 -l/2 l/2 l/2]; xselect() driver("X11") xset("pixmap",1) plot2d(\%inf,\%inf, frameflag=3, rect=[-1,-1,1,1], axesflag=0) xset("background",1); // czarne t\l{}o xbasr() // tutu utile pour la mise a jour du background xset("thickness",3) // zwi\k{e}kszenie grubo\ sci kre\ slonych linii (3 pixele) xset("font",2,4) for i=1:n xset("wwpc") theta = 3*t(i); xyr = [cos(theta) -sin(theta);... sin(theta) cos(theta)]*xy; xset("color",2) xfpoly(xyr(1,:)+xg(i), xyr(2,:)+yg(i)) xset("color",5) xpoly(xyr(1,:)+xg(i), xyr(2,:)+yg(i),"lines",1) xset("color",32) xtitle("Animation simple") xset("wshow") // przeniesienie pixmapy na ekran end driver("Rec") // powrot do domy\ slnego drivera xset("default") // przywr\ ocenie domy\ slnego kontekstu graficznego 4.10 Powierzchnie Podstawow funkcj pozwalajc na tworzenie wykresów powierzchni jest plot3d23. Przy preprezentacji powierzchni za pomoca facettes tutu mamy mo|liwo[ okre[lenia innego koloru dla ka|dej z nich. Od wersji 2.6 mo|na tak|e, zarówno dla fscettes trójktnych jak i kwadratowych okre[li kolor dla ka|dego z wierzchoBków, przez co rendu tutu otrzymywane jest przez interpolacj kolorów okre[lonych w wierzchoBkach. 4.10.1 Wprowadzenie doplot3d Gdy powierzchnia opisane jest przez równanie typu z = f(x, y), szczególnie Batwo mo|na j narysowa je[li argumenty s z obszaru prostoktnego. Jako przykBad rozwa|my funkcj f(x, y) = cos(x)cos(y) dla (x, y) " [0, 2À] × [0, 2À] : x=linspace(0,2*\%pi,31); // dyskretyzacja zmiennej x (a tak\.ze y : b\k{e}dzie to samo) z=cos(x) *cos(x); // zbi\ or warto\ sci zmiennej z : macierz z(i,j) = f(x(i),y(j)) plot3d(x,x,z) // rysunek W efekcie otrzymamy co[ podobnego do rysunku (13)24. W najbardziej ogólnej formie funkcji plot3d lub plot3d1 u|ywamy piszc 23 plot3d1u|ywana analogicznie pozwala uzale|ni warto[ koloru od warto[ci przyjmowanej na osi OZ. 24 DokBadnie rzecz biorc rysunek ten przedstawia efekt u|ycia funkcjiplot3d1gdzie na potrzeby publikacji zamiast kolorów u|yto odcieni szaro[ci a tak|e ustawiono troch inny punkt widzenia. 48 Z 1 0 -1 0.0 0.0 3.1 3.1 Y X 6.3 6.3 Rysunek 13: Funkcja z = cos(x)cos(y) fig:4 plot3d(x,y,z <,opt_arg>*) plot3d1(x,y,z <,opt_arg>*) gdzie dlaplot2d,<,opt_arg>*ozancza cig argumentów opcjonalnych,opt_argprzyjmuje form sBowo_kluczowe=warto[. W najprostszym przypadku,xiys wektorami liniowymi ((1, nx) i (1, ny)) odpowiadajcymi dyskretyzacji zmiennej x oraz y, natomiastzjest macierz (nx, ny) tak, |e zi,j jest  wysoko[ci w punkcie (xi, yj). Opcjonalne argumenty to: 1. theta=val_theta i alpha=val_alpha to dwa kty (w stopniach) okre[lajce punkt widzenia we wspóBrzdnych sfe- rycznych (je[li O jest [rodkiem pudeBka englobante tutu, Oc kierunkiem patrzenia kamery, wówczas ± = kt(Oz, Oc) i ¸ = kt(0x, Oc2 ) gdzie Oc2 jest rzutem Oc na pBaszczyzn Oxy); 2. leg=val_leg pozwala okre[li nazw dla ka|ej z osi (na przykBad leg="x@y@z"), argument efektywny val_leg jest BaDcuchem znakowym, w którym@stanowi separator pomidzy nazwami; 3. flag=val_flag gdzie val_flag jest wektorem o trzech skBadowych [mode type box] pozwalajcym okre[li: (a) parametr mode zwizany jest z rysunkiem faces tutu i siatki: i. dla mode > 0, faces niewidoczne s usuwane25, siatka pozostaje widoczna; ii. dla mode = 0, otrzymujemy tutu un rendu (fr. fil de fer, ang. wireframe) de la surface; iii. dla mode < 0, faces niewidoczne s usuwane a siatka nie jest rysowana. Dodatnio okre[lona strona face tutu (patrz dalej) bdzie malowana z wykorzystaniem koloru numer mode za[ strona przeciwna przy u|yciu koloru, który mo|emy okre[li instrukcjxset("hidden3d",colorid)(domy[lnie jest to 4 kolor z palety). (b) parametr type pozwala okre[li skal: type otrzymana skala 0 u|ycie poprzedniej skali (tutu ou par défaut) 1 skala wraz zebox 2 skala otrzymana w oparciu o minimum i maximum zmiennych tutu 3 jak 1 ale skala jest izometryczna 4 jak 2 ale skala jest izometryczna 5 variante de 3 6 variante de 4 25 tutu actuellement c est l algorithme du peintre qui est utilisé, c-a-d qu un tri des facettes est effectué et les plus éloignées de l observateur sont dessinées en premier. 49 (c) parametr box kontroluje tutu le pourtour du graphe : box otrzymany efekt 0 juste le dessin de la surface 2 osie pod powierzchnia s rysowane 3 jak dla 2 z dodatkowym tutu la boite englobante 4 jak dla 3 z dodatkowym tutu la graduation des axes 4. ebox=val_ebox pozwala okre[li tutu la boite englobante, val_ebox jest wektorem o 6 skBadowych [xmin, xmax, ymin, ymax, zmin, zm Oto maBy przykBad, w którym u|ywa si prawie wszystkich parametrówplot3d. Jest to prosta animacja pozwalajca lepiej zrozumie zmiane punktu widzenia przy pomocy parametrówthetai alpha. W skrypcie tym u|ywamyflag=[2 4 4], co oznacza " mode = 2 powierzchnia bdzie rysowana (jej dodatno okre[lona strona) kolorem 2 oraz bdzie widoczna siatka; " type = 4 zostanie u|yta skala izometryczna obliczona na podstawie danych (jest to równowa|ne wyborowi type = 3 z parametrem ebox przyjmujcym warto[ci równe minimum i maksimum danych); " box = 4 rysunek bdzie zawieraB pudeBko tutu boite oraz stopniowanie et des graduations. x=linspace(-%pi,%pi,31); z=sin(x) *sin(x); n = 200; theta = linspace(30,390,n); // pe\l{}ny obr\ ot alpha = [linspace(60,0,n/2) linspace(0,80,n/2)]; // od g\ ory // do do\l{}u xselect() xset("pixmap",1) // aktywowanie podw\ ojnego buforowania driver("X11") // zmieniamy parametr theta for i=1:n xset("wwpc") // wyczyszczenie bie\.z\k{a}cego bufora plot3d(x,x,z,theta=theta(i),alpha=alpha(1),leg="x@y@z",flag=[2 4 4]) xtitle("zmiany punktu widzenia przy pomocy parametru theta") xset("wshow") end // zmieniamy parametr alpha for i=1:n xset("wwpc") // wyczyszczenie bie\.z\k{a}cego bufora plot3d(x,x,z,theta=theta(n),alpha=alpha(i),leg="x@y@z",flag=[2 4 4]) xtitle("zmiany punktu widzenia przy pomocy parametru alpha") xset("wshow") end xset("pixmap",0) driver("Rec") 4.10.2 Kolory Powrómy jeszcze na chwilk do ostatniego przykBadu i zamiastplot3dnapiszmyplot3d1uzale|niajc tym samym kolory od warto[ci zmiennej z. Narysowana powierzchnia powinna przypomina w tym momencie mozaike gdy| wykorzystywana plaeta kolorów domy[lnie nie jest  cigBa . Paleta kolorów jest macierz o wymiarze(nb_couleurs,3), gdzie i-ta linia definiuje intensywno[ci skBadowej czerwo- nej (warto[ zawarta w przedziale od 0 do 1), zielonej i niebieskiej dla i-tego koloru. Majc tak macierz, któr nazwiemy C, polecenie xset(olormap",C)pozwala j wczyta (zaBadowa) do kontekstu graficznego bie|cego okna graficznego. Funkcje, hotcolormaporaz greycolormapdostarczaj mapy ze stopniow zmian kolorów26. MaBa uwaga: je[li doko- nujemy zmiany palety kolorów po narysowaniu jakiego[ rysunku, zmiany nie bd widoczne natychmiast (jest to zachowania normalne); wystarczy zmieni rozmiar okna lub tutu d envoyer l ordrexbasr(numero_fenetre). Oto nowy przykBad x = linspace(0,2*%pi,31); z = cos(x) *cos(x); 26 Patrz tak|e sekcja Contributions na stronie Scilab. 50 C = hotcolormap(32); // hot colormap z 32 kolorami xset("colormap",C) xset("hidden3d",30) // wyb\ or koloru 30 do rysowania po ujemnie okre\ slonej stronie xbasc() plot3d1(x,x,z, flag=[1 4 4]) // spr\ obuj tak\.ze z flag=[-1 4 4] Uwaga : wplot3d1, wykorzystuje si jedynie znak parametru mode (je[li mode e" 0 siatka zostanie narysowana, nie zostanie za[ gdy mode < 0). 4.10.3 plot3dz des facettes Chcc u|y tej funkcji w jak najogólniejszej postaci nale|y poda opis powierzchni uwzgldniajc pojedyDczy tutu facettes. Okre[lany jest on przy pomocy 3 macierzyxf, yf, zfo wymiarze (nb_sommets_par_face, nb_faces), gdziexf(j,i),yf(j,i),zf( s wspóBrzdnymi j-tego wierzchoBka i-tego tutu facette. Poza tymi zmianami dalsze u|ycie jest takie samo jak w poprzednich przykBadach: plot3d(xf,yf,zf <,opt_arg>*) Orientacja tutu facettes jest odmienna od zwyczajowo przyjtej (patrz rysunek (14). P3 face négative pour Scilab P4 P2 face positive pour Scilab P1 Rysunek 14: orientacja tut facettes w Scilab-ie fig:face_orientation W celu zdefinowania kolorów dla ka|dego facette, trzeci argument musi by list : list(zf,colors)gdzie colors jest wektorem o rozmiarzenb_faces,colors(i)okre[lan numer (w palecie) koloru i-tego tutu facette. z P4 P3 P1 y P2 x Rysunek 15: Trój[cian fig:tetraedre Jak w pierwszym przykBadzie, narysujemy [ciany trój[cianu z rysunku (15), dla którego: îø ùø îø ùø îø ùø îø ùø 0 1 0 0 ðø ûø ðø ûø ðø ûø ðø ûø P1 = 0 , P2 = 0 , P3 = 1 , P4 = 0 , 0 0 0 1 51 gdzie definicja [cian jest nastpujaca (w ten sposób otrzymujemy [ciany zewntrzne z orientacj dodatni dla Scilab-a): f1 = (P1, P2, P3), f2 = (P2, P4, P3), f3 = (P1, P3, P4), f4 = (P1, P4, P2) Napiszmy wic: // f1 f2 f3 f4 xf = [ 0 1 0 0; 1 0 0 0; 0 0 0 1]; yf = [ 0 0 0 0; 0 0 1 0; 1 1 0 0]; zf = [ 0 0 0 0; 0 1 0 1; 0 0 1 0]; // tutu ouf ! xbasc() plot3d(xf,yf,list(zf,2:5), flag=[1 4 4], leg="x@y@z",alpha=30, theta=230) xselect() Otrzymany efekt powinien przypomina rysunek 16. Mo|na zauwa|y, |eplot3du|ywa prostej tutu projection orthographique et non une projection perspective plus réaliste. Rysunek 16: Trój[cian narysowany w Scilab-ie. fig:tetra Na bie|ce potrzeby, obliczenia tutu des facettes, mog by efektywniejsze dziki funkcjom: " eval3dpinf3ddla powierzchni zdefiniowanych przy pomocy x = x(u, v), y = y(u, v), z = z(u, v) (patrz 4.10.4) ; " genfac3ddla powierzchni definiowanych przy pomocy z = f(x, y) (przykBad pokazany jest troch dalej (4.10.5)). Je[li powierzchnia (wielo[cian) zdefiniowana jest w taki sam sposób jak trój[cian z przykBadu nie mo|na jego narysowa bezpo[rednio przy u|yciuplot3d. W celu otrzymania opisu oczekiwanego przez Scilab-a mo|na wykorzysta funkcj podobn do przedstawionej poni|ej: function [xf,yf,zf] = facettes_polyedre(P) // tutu transforme la structure Polyedre de l exemple // sur les tlist en description de facettes pour // visualisation avec plot3d [ns, nf] = size(P.face) // ns: ilo\ s\ c wierzcho\l{}k\ ow tworz\k{a}cych \ scian\k{e} // nf: ilo\ s\ c \ scian xf=zeros(P.face); yf=zeros(P.face); zf=zeros(P.face) 52 for j=1:ns num = connect(ns+1-j,:) // dla odwr\ ocenia orientacji xf(j,:) = P.coord(1, num) yf(j,:) = P.coord(2, num) zf(j,:) = P.coord(3, num) end endfunction Majc tak okre[lon funkcj, rysunek wielo[cianu otrzymamy piszc [xf,yf,zf] = facettes_polyedre(Cube); plot3d(xf,yf,list(zf,2:7), flag=[1 4 0],theta=50,alpha=60) 4.10.4 Rysowanie powierzchni opisanej przy pomocy x = x(u, v), y = y(u, v), z = z(u, v) doc:surfparam Odpowiedz: wzi dyskretyzacj dziedziny parametrów obliczy tutu les facettes przy pomocy funkcji (eval3dp). Ze wzgl- dów wydajno[ciowych, funkcja okre[lajca parametry powierzchni powinna by napisana  wektorowo . Je[li (u1, u2, . . . , um) i (v1, v2, . . . , vn) stanowi dyskretyzacj dziedziny parametrów, funkcja powinna by wywoBywana jednokrotnie z dwoma  du- |ymi wektorami rozmiaru m × n: U = (u1, u2, . . . , um, u1, u2, . . . , um, . . . . . . , u1, u2, . . . , um) 1 2 n V = (v1, v1, . . . , v1, v2, v2, . . . , v2, . . . . . . , vn, vn, . . . , vn) m fois v1 m fois v2 m fois vn W oparciu o te dwa wektory, funkcja powinna tworzy 3 wektory X, Y et Z wymiaru m × n takie, |e: Xk = x(Uk, Vk), Yk = y(Uk, Vk), Zk = z(Uk, Vk) Oto kilka przykBadów parametryzacji powierzchni, tutu écrite27 de façon U pouvoir tre utilisée aveceval3dp: function [x,y,z] = tore(theta, phi) // tutu paramétrisation classique d un tore de rayons R et r et d axe Oz R = 1; r = 0.2 x = (R + r*cos(phi)).*cos(theta) y = (R + r*cos(phi)).*sin(theta) z = r*sin(phi) endfunction function [x,y,z] = helice_torique(theta, phi) // tutu paramétrisation d une helice torique R = 1; r = 0.3 x = (R + r*cos(phi)).*cos(theta) y = (R + r*cos(phi)).*sin(theta) z = r*sin(phi) + 0.5*theta endfunction function [x,y,z] = moebius(theta, rho) // wst\k{e}ga Moëbius-a R = 1; x = (R + rho.*sin(theta/2)).*cos(theta) y = (R + rho.*sin(theta/2)).*sin(theta) z = rho.*cos(theta/2) endfunction function [x,y,z] = tore_bossele(theta, phi) // tutu paramétrisation d un tore dont le petit rayon r est variable avec theta R = 1; r = 0.2*(1+ 0.4*sin(8*theta)) x = (R + r.*cos(phi)).*cos(theta) 27 Piszemy je w naturalny sposób, pamitajac aby zamiast*i/napisa.*i./i wszystko powinno dziaBa! Û 53 y = (R + r.*cos(phi)).*sin(theta) z = r.*sin(phi) endfunction Oto przykBad wykorzystujcy ostatni powierzchni: // skrypt rysuj\k{a}cy powierzchni\k{e} opisan\k{a} za pomoc\k{a} r\ ownania parametrycznego theta = linspace(0, 2*%pi, 160); phi = linspace(0, -2*%pi, 20); [xf, yf, zf] = eval3dp(tore_bossele, theta, phi); // tutu calcul des facettes xbasc() plot3d1(xf,yf,zf) xselect() Je[li chcemy u|y kolorów a nie otrzymujemy ich na rysunku, spowodowane jest to niewBa[ciw orientacj; wystarczy odwróci kierunek jednego z wektorów stanowicego dyskretyzacj dziedziny. Rysunek 17: Un tore bosselé... tutu fig:6 tutu Funkcjanf3djest lekko podobna doeval3dp, ale majc dyskretyzacj u i v trzeba zdefiniowa soit-mme des matrices X, Y, Z tak, |e: Xi,j = x(ui, vj) Yi,j = y(ui, vj) Zi,j = z(ui, vj) et vos facettes s obtiennent alors avec[xf,yf,zf] = nf3d(X,Y,Z). Jako przykBad rozwa|my wstg Moëbius-a définit juste avant : nt = 120; nr = 10; rho = linspace(-0.5,0.5,nr); theta = linspace(0,2*%pi,nt); R = 1; X = (R + rho *sin(theta/2)).*(ones(nr,1)*cos(theta)); Y = (R + rho *sin(theta/2)).*(ones(nr,1)*sin(theta)); Z = rho *cos(theta/2); [xf,yf,zf] = nf3d(X,Y,Z); xbasc() plot3d(xf,yf,zf, flag=[2 4 6], alpha=60, theta=50) xselect() Uwaga: w celu otrzymania pawidBowej macierzy nale|aBo u|y funkcji ones, tutu ce qui ne rend pas le code tr s clair: funkcjaeval3dpjest Batwiejsza w u|yciu! 54 Rysunek 18: Wstga Moëbius-a fig:moebius 4.10.5 plot3dz interpolacj kolorów doc:interpcolor Od wersji 2.6 , mo|liwe jest okre[lenie koloru dla ka|dego wierzchoBka tutu d une facette. Wystarczy okte[li macierzcolors takiego samego wymiaru jakxf, yf, zfdajc opis ka|dego facette, to znaczy tak, |ecolors(i,j)jest kolorem zwi- zanym z i-tym wierzchoBkiem j-tego face, i poBczy z trzecim argumentem (zf) w list: plot3d(xf,yf,list(zf,colors) <,opt_arg>*) Oto przykBad pocztkowyplot3dbez rysowania siatki: " tutu une couleur par face pour le dessin de gauche, " tutu une couleur par sommet pour celui de droite. Aby obliczy wartosci kolorów, wykorzystano pomcnicz funkcj wi|c w sposób liniowy warto[ci na karcie graficznej tutu!!!. (u|yto funkcjidsearchdostpnej od wersji 2.7 ale mo|na Batwo tutu vous en passer). Zauwa|y tak|e bdzie mo|na wykorzy- tanie funkcjigenfac3dpozwalajcej na obliczenie tutu les facettes. // przyk\l{}ad wykorzystania plot3d z interpolacj\k{a} kolor\ ow function [col] = associe_couleur(val) // przypisanie koloru dla ka\.zdej z warto\ sci z val tutu przypsanie do warto\ sci?!!! n1 = 1 // numer 1-ego koloru n2 = xget("lastpattern") // numer ostatniego koloru nb_col = n2 - n1 + 1 classes = linspace(min(val),max(val),nb_col) col = dsearch(val, classes) endfunction x=linspace(0,2*\%pi,31); z=cos(x) *cos(x); [xf,yf,zf] = genfac3d(x,x,z); xset("colormap",graycolormap(64)) zmeanf = mean(zf,"r"); zcolf = associe_couleur(zmeanf); zcols = associe_couleur(zf); xbasc() xset("font",6,2) // font 6 (helvetica) tutu n est disponible // que dans la version cvs de scilab subplot(1,2,1) plot3d(xf,yf,list(zf,zcolf), flag=[-1 4 4]) 55 xtitle("Jeden kolor na tutu face") subplot(1,2,2) plot3d(xf,yf,list(zf,zcols), flag=[-1 4 4]) xtitle("Jeden kolor na wierzcho\l{}ek") xselect() Rysunek 19: Z i bez interpolacji kolorów. fig:shade 4.11 Krzywe w przestrzeni Podstawow funkcj pozwalajc rysowa krzywe w przestrzeni jestparam3d. Oto klasyczny przykBad tutu de l hélice: t = linspace(0,4*\%pi,100); x = cos(t); y = sin(t) ; z = t; param3d(x,y,z) // ewentualne wymazanie okna graficznego za pomoc\k{a} xbasc() tutu mais comme cette derni re ne permet que d affichier une seule courbe nous allons nous concentrer sur param3d1 qui permet de faire plus de choses. Oto jej skBadnia: param3d1(x,y,z <,opt_arg>*) param3d1(x,y,list(z,colors) <,opt_arg>*) Macierzex, yetzmusz by tego samego wymiaru (np,nc) a ilo[ krzywych (nc) jest okre[Bona przez ilo[ kolumn (jak dla plot2d). Parametry opcjonalne s takie same jak dla fnkcji plot3d, modulo le fait que flag tutu ne ne comporte pas de param tre mode. colorsjest wektorem okre[lajcym styl dla ka|dej krzywej (dokBadnie jakplot2d), to znaczy gdycolors(i)jest war- to[ci caBkowit dodatni, i-ta krzywa rysowana jest i-tym kolorem z bie|cej palety kolorów (tutu ou avec différents pointillés sur un terminal noir et blanc); je[li za[ zawarta jest pomidzy -9 i 0, otrzymujemy rysnek zBo|ony z punktów (tutu non reliés) zaznaczonych odpowiednim symbolem. Oto przykBad, który powinien doprowadzi nas do rysunku (20): t = linspace(0,4*\%pi,100) ; x1 = cos(t); y1 = sin(t) ; z1 = 0.1*t; // spirala x2 = x1 + 0.1*(1-rand(x1)); y2 = y1 + 0.1*(1-rand(y1)); z2 = z1 + 0.1*(1-rand(z1)); xbasc(); xset("font",2,3) param3d1([x1 x2],[y1 y2],list([z1 z2], [1,-9]), flag=[4 4]) xset("font",4,4) xtitle("Spirala z per\l{}ami") Jak dla plot2d funkcj t nale|y wywoBa kilkakrotnie je[li ró|ne krzywe nie maj takiej samej ilo[ci punktówon. Oto skrypt, wyja[niajcy jak utworzy dwie grupy punktów ze ró|nymi znakami i kolorami: n = 50; // ilo\ s\ c punkt\ ow P = rand(n,3); // losowanie punkt\ ow // tutu on impose les dimensions de la boite englobante ebox = [0 1 0 1 0 1]; 56 Rysunek 20: Krzywa i punkty w przestrzeni. fig:helice // rozdzielenie punkt\ ow na dwie grupy aby pokaza\ c jak przypisa\ c // r\ o\.zne symbole i kolory do punkt\ ow m = 30; P1 = P(1:m,:) ; P2 = P(m+1:n,:); // rysunek xbasc() // pierwsza grupa punkt\ ow xset("color",2) // tutu du bleu avec la carte par defaut niebieski w domy\ slnej palecie (?!) param3d1(P1(:,1),P1(:,2),list(P1(:,3), -9), alpha=60, theta=30,... leg="x@y@z", flag=[3 4], ebox=ebox) // tutu flag=[3 4] : 3 -> echelle iso se basant sur ebox // tutu 4 -> boite + graduation // druga grupa punkt\ ow xset("color",5) // tutu du rouge avec la carte par defaut param3d1(P2(:,1),P2(:,2),list(P2(:,3), -5), flag=[0 0]) // tutu -5 pour des triangles inverses // tutu [0 0] : echelle fixée et cadre dessiné avec l appel précédent xset("color",1) // aby ustawi\ c czarny jako bie\.z\k{a}cy kolor xtitle("Punkty...") xselect() 4.12 Ró|no[ci Istnieje jeszcze wiele prymitywów graficznych: 1. contour2dicontourpozwalaj rysowa linie poziomicowe dla funkcji z = f(x, y) okre[lonej na prostokcie; 2. grayplotiSgrayplotpozwalaj reprezentowa warto[ci funkcji tutu qui permettent de représenter les valeurs d une telle fonction en utilisant des couleurs ; 3. fecodgrywaja tak sam rol jak dwie poprzednie dla funkcji okre[lonej tutu est définie sur une triangulation plane ; 4. champpozwala okre[li pole wektorowe w 2D ; 5. tutu wiele funkcji wywoBywanych w tym rozdziale dopuszcza ró|ne parametry aby tworzy wykresy funkcji bardziej bezpo[rednio si on fournit une fonction scilab comme argument (le nom de ces fonctions commence par un ffplot2d, fcontour2d, fplot3d, fplot3d1, fchamp,...). 57 Aby zda sobie spraw z mo|liwo[ci28 wystarczy przejrze tutu rubrique (dziaB,sekcja (?!)) Graphic Library pomocy. Od wersji 2.7 istnieje nowy model graficzny  zorientowany obiektowo pozwalajcy modyfikowa wBa[ciwo[ci grafiki po ujrzeniu rysunku. Domy[lnie nie jest on aktywny, ale je[li kto[ chce poeksperymentowa nale|y wyda polecenie: set("figure_style","new") przed utworzeniem rysunku. Jako |e ten nowy tryb jest w trakcie rozwijania zalecane jest u|wanie wersji tutu cvs de scilab. Biblioteka Enrico Ségré, któr mo|na tutu  [cign z jego strony: http://www.weizmann.ac.il/~fesegre/ tutu compl te celle de Scilab et contient takze funkcje pozwalajce upro[ci tutu certaines taches. 5 Zastosowania i uzupeBnienia RozdziaB ten przedstawia metody rozwizywania w Scilabie pewnych typów problemów analizy numerycznej (aktualnie równaD ró|niczkowych...) oraz dostarcza uzupeBnieD/dodatkowych informacji niezbdnych podczas projektowania prostych symulacji stochastycznych. 5.1 Równania ró|niczkowe Scilab dysponuje pot|nym interfejsem dla rozwizywania numerycznego (w sposób przybli|ony) równaD ró|niczkowych przy zastosowaniu prostej instrucjiode. Rozwa|my zatem nastpujce równanie ró|niczkowe z warunkiem pocztkowym : u2 = f(t, u) u(t0) = u0 gdzie u(t) jest wektorem w Rn, f jest funkcj R × Rn -’! Rn, oraz u0 " Rn. ZakBadamy warunki brzegowe dla których rozwizanie istnieje i jest jednoznaczne do okresu T . 5.1.1 Podstawowe u|ycieode Zdefiniujmy funkcj f jako funkcj Scilaba w nastpujcyj sposób : function [f] = MojaFunkcja(t,u) // tutaj kod definiujacy f jako funkcje t i u. endfunction Rmq : Podobnie, je[li równanie jest ...autonome..., nale|y doprowadzi równanie do postaci w której mamy funkj zale|ne od zmiennej t Oto przykBad kodu dla równania Van der Pola : y2 2 = c(1 - y2)y2 - y kBadc u1(t) = y(t) oraz u2(t) = y2 (t) przeksztaBcamy je do ukBadu dwuch równaD ró|niczkowych pierwaszego rzdu : d u1(t) u2(t) = u2(t) c(1 - u2(t))u2(t) - u1(t) dt 1 function [f] = VanDerPol(t,u) // second membre pour Van der Pol (c = 0.4) f(1) = u(2) f(2) = 0.4*(1 - u(1)^2)*u(2) - u(1) endfunction Nastpnie wywoBujemyodeaby rozwiza równanie (ukBad równaD) wzgldem t0 w T , wychodzc od u0 (wektor kolumnowy) i chcc otrzyma rozwizanie w chwili t(1) = t0, t(2), ..., t(m) = T , wpiszemy instrukcj : t = linspace(t0,T,m); [U] = ode(u0,t0,t,MojaFunkcja) 28 tutu attention risque de noyade ! 58 Otrzymamy w ten sposób macierzUo wymiarach (n, m), w którejU(i,j)jest rozwizaniem cz[ciowym/czstkowym ui(t(j)) (i-ta skBadowa w chwili t(j)). Uwaga : Liczba skBadowych branych dlat(czyli momenty w których otrzymuje si rozwizanie) nie maj nic wspólnego z precyzj obliczeD. ............. Funkcjaodebazuje na wielu algorytmach pozwalajcych dostosowanie jej do wielu sytuacji... Aby wybra odpowiedni metod nale|y doda odpowiedni parametr w trakcie wywoBania (patrz Help). Domy[lnie (tzn. bez wybierania explicite jednej metody) stosowana jest metoda Adamsa predykcji, natomiast w przypadku gdy Scilab okre[li równanie jako sztywne29 algorytm zostaje zmieniony na metod Gear-a. Oto kompletny przykBad dla równania Van der Pola. W tym przypadku przestrzeD stanów jest pBaska, mo|na zatem otrzyma obraz dynamiki rysujc pole wektorowe w prostokcie [xmin, xmax] × [ymin, ymax] za pomoc instrukcji graficznejfchampw nastpujcy sposób: fchamp(MojaFunkcja,t,x,y) gdzie MojaFunkcja oznacza nazw funkcji Scilaba bdc skBadnikiem po prawej stronie równania ró|niczkowego, t jest momentem w którym chcemy naszkicowa pole (w przypadku cigBym równania mo|emy przyj dowoln warto[ np. 0) orazx iys wektorami wierszowymi o nx i ny wspóBrzdnych wyznaczacymi punkty siatki na której znajduj si kirunki wyznaczajce pole wektorowe. // 1/ kre\ slimy pole wektorowe odpowiadaj\k{a}ce r\ ownaniu Van der Pola n = 30; delta = 5 x = linspace(-delta,delta,n); // tutaj y = x xbasc() fchamp(VanDerPol,0,x,x) xselect() // 2/ rozwi\k{a}zujemy r\ ownanie r\ o\.zniczkowe m = 500 ; T = 30 ; t = linspace(0,T,m); // moment w kt\ orym znajduje sie rozwi\k{a}zanie u0 = [-2.5 ; 2.5]; // warunek pocz\k{a}tkowy [u] = ode(u0, 0, t, VanDerPol); plot2d(u(1,:) ,u(2,:) ,2,"000") 5.1.2 Van der Pol jeszcze raz W tej cz[ci zwrócimy nasz uwag na mo|liwo[ci graficzne Scilaba pozwalance otrzyma wiele |danych trajektorii bez ponow- nego uruchamiania skryptu z inn warto[ci argumentu u0 ......... Po wy[wietleniu pola wektorowego, ka|dy warunek pocztkowy bdzie zadany poprzez kliknicie lewym przyciskiem myszy30; pojawi si wówczas punkt na |danej pozycji pocztkowej. Jest to mo|liwe dziki funkcjixclick, której skBadnia jest nastpujca : [c_i,c_x,c_y]=xclick(); Jak wida, Scilab dyspinuje odpowiednimi procedurami, które umo|liwiaj miedzy innymi reagowanie na tzw zdarzenia graficzne jak kliknicie mysz . Kiedy takie zdarzenia ma miejsce, nastpuje automatczna lokalizacja pozycji punktu (w bie|cej skali) poprzez przypisanie jego wspóB|dnych do zmiennychc_xorazc_y. Trzeci argument, czylic_ioznacza odpowiedni klawisz myszy zgodnie z poni|sz tabelk : warto[ dlac_i klawisz 0 lewy 1 [rodkowy 2 prawy W skrypcie wyko|ystano kliknicie na lewy klawisz myszy jako ten, wyznaczajcy punkt pocztkowy dla pku zdarzeD. Na koniec nadano ka|dej z wyrysowanych trajektorii inny kolor (tablica couleur pozwala wybra jeden z dostpnych kolorów). Aby otrzyma skal izometryczn urzyjemy fchamp dodajc opcjonalny argument jak dla plot2d (nale|y u|y strf=wartosc_strf jako |e parametry frameflag i axesflag nie s ........). Ostatni aspekt : aby zaznaczy punkt pocztkowy obrysowano go naBum kóBkiem, natomiast aby pokaza wszystkie mo|liwo[ci okna graficznego wyko|ystano sze[cienny ukBad wspóBrzdnych. Ostatnia uwaga : podczas gdy jest wy[wietlone pole wektorowe, mo|na maksymalizowa okno graficzne! Klika- jc dwukrotnie otrzymano rysunek (21) wszystkie trajektorie zbie|ne na jednej sferze, bdcej ........ dla tego równania 29 równanie jest raide w przypadku gdy (mniej lub bardziej) Bczy si z metodami jednoznacznymi 30 jak sugerowano w jednym z artykuBów dotyczcych Scilaba, zamieszczonym w Linux Magazine 59 // 1/ wykres pola wektorowego dla rownania Van der Pola n = 30; delta_x = 6 delta_y = 4 x = linspace(-delta_x,delta_x,n); y = linspace(-delta_y,delta_y,n); xbasc() fchamp(VanDerPol,0,x,y, strf="041") xselect() // 2/ resolution de l equation differentielle m = 500 ; T = 30 ; t = linspace(0,T,m); couleurs = [21 2 3 4 5 6 19 28 32 9 13 22 18 21 12 30 27] // 17 couleurs num = -1 while %t [c_i,c_x,c_y]=xclick(); if c_i == 0 then plot2d(c_x, c_y, style=-9, strf="000") // un petit o pour marquer la C.I. u0 = [c_x;c_y]; [u] = ode(u0, 0, t, VanDerPol); num = modulo(num+1,length(couleurs)); plot2d(u(1,:) ,u(2,:) , style=couleurs(num+1), strf="000") elseif c_i == 2 then break end end 4.24 Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ 3.39 Ÿ Ÿ 2.54 1.70 Ÿ 0.85 Ÿ 0.00 -0.85 -1.70 -2.54 Ÿ Ÿ Ÿ -3.39 Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ -4.24 -6.0 -4.8 -3.6 -2.4 -1.2 0.0 1.2 2.4 3.6 4.8 6.0 Rysunek 21: Quelques trajectoires dans le plan de phase pour l équation de Van der Pol fig:vanderpol 5.1.3 Troche wicej oode 60 W tym drugim przykBadzie pokarzemy inne zastosowanie funkcjiodedla równaD z parametrem. ........ Oto nasze nowe równanie ró|niczkowe (Równanie Brusselator) : du1 = 2 - (6 + ë)u1 + u2x2 1 dt du2 = (5 + ë)u1 - u2u2 1 dt które przyjmuje jako punkt krytyczny Pstat = (2, (5 + ë)/2). Poniewa| warto[i parametru ë zmieniaj sie z ujemnej na dodat- ni, punkt stacjonarny zmienia swoj charakter (jest staBy lub zmienny, wchodzc, dla ë = 0 w zjawisko rozgaBzienia Hopf). Interesuj nas trajektorie z warunkiemi pocztkowymi z ssiedztwa tego punktu. Oto funkcja liczca to równanie : function [f] = Brusselator(t,u,eps) // f(1) = 2 - (6+eps)*u(1) + u(1)^2*u(2) f(2) = (5+eps)*u(1) - u(1)^2*u(2) endfunction Aby poda parametr, zastpimy w wywoBaniu ode nazw funkcji (tutaj Brusselator) przez uporzdkowan list zawierajc nazw funkcji oraz jej parametry : [x] = ode(x0,t0,t,list(MojaFunkcja, par1, par2, ...)) W naszym przypadku : [x] = ode(x0,t0,t,list(Brusselator, eps)) a nastpnie zastosujemyfchampaby narysowa pole. Aby ustali zakres tolerancji dla bBdu lokalnego rozwizania u|yjemy dodatkowych parametrówrtolorazatolzaraz po nazwie funckji (lub listy zawierajcej t funkcj wraz z jej parametrami. W ka|dym kroku czasowym, tk-1 ’! tk = tk-1 + "tk, obliczane jest oszacowanie bBdu lokalnego e (tzn. bBdu dla kroku czasowego wychodzc z warunku pocztkowego v(tk-1) = U(tk-1)) : tk e(tk) C" U(tk) - f(t, v(t))dt + U(tk-1) tk-1 (drugi skBadnik jest rozwizaniem dokBadnym zale|nym od rozwizania mumerycznego U(tk-1) otrzymanego w poprzednim kroku) a nastpnie sprawdza czy zawiera si on w zakresie tolerancji formuowane za pomoc dwóch parametrówrtoletatol : toli = rtoli " |Ui(tk)| + atoli, 1 d" i d" n w przypadku gdy dane s dwa wektory dBugo[i n dla tych praramertów, oraz: toli = rtol " |Ui(tk)| + atol, 1 d" i d" n gdy dane s skalary. Je|eli |ei(tk)| d" toli dla wszystkich tk, krok jest akceptowany a nastpnie obliczony zostaje nowy krok czasowy w taki sposób, |e kryterium naBo|one na przyszBy bBd bdzie pewnym sposobem jego realizacji< ?????. W przeciwnym razie, ponownie caBkuje sie wyra|enie dla t(k -1) przyjmujc nowy, nieco mniejszy krok(.........). Metody tego typu oprócz zmien- nych postpu czasowgo, manipuluj równie| kolejno[ci równaD w celu otrzymania dobrej skuteczno[ci/wydajno[ci informacji... Domy[lnie stosowanymi warto[ciami s rtol = 10-5 i atol = 10-7(za wyjtkiem typów wybierajcych metod Runge Kutta). Wa|na uwaga : caBkowanie ma|e czsto sie nie powiz... Poni|ej przedstawiamy przykBadowy skrypt, w którym punkty krytyczne oznaczono maBymi, czarnymi kwadratami (otrzymano je przy pomocy prostej funkcji graficznejxfrect) : // Brusselator eps = -4 P_stat = [2 ; (5+eps)/2]; // granice dla wykresu pola wektorowego delta_x = 6; delta_y = 4; x_min = P_stat(1) - delta_x; x_max = P_stat(1) + delta_x; y_min = P_stat(2) - delta_y; y_max = P_stat(2) + delta_y; n = 20; x = linspace(x_min, x_max, n); y = linspace(y_min, y_max, n); // 1/ wykres pola wektorowego xbasc() fchamp(list(Brusselator,eps),0,x,y, strf="041") xfrect(P_stat(1)-0.08,P_stat(2)+0.08,0.16,0.16) // dla zaznaczenia punktow krytycznych 61 xselect() // 2/ rozwiazanie rownania rozniczkowego m = 500 ; T = 5 ; rtol = 1.d-09; atol = 1.d-10; // zakres tolerancji bledu t = linspace(0,T,m); couleurs = [21 2 3 4 5 6 19 28 32 9 13 22 18 21 12 30 27] num = -1 while %t [c_i,c_x,c_y]=xclick(); if c_i == 0 then plot2d(c_x, c_y, style=-9, strf="000") // male o aby zanaczyc C.I. u0 = [c_x;c_y]; [u] = ode(u0, 0, t, rtol, atol, list(Brusselator,eps)); num = modulo(num+1,length(couleurs)); plot2d(u(1,:) ,u(2,:) , style=couleurs(num+1), strf="000") elseif c_i == 2 then break end end 4.74 Ÿ Ÿ Ÿ 3.89 Ÿ Ÿ Ÿ 3.04 Ÿ 2.20 Ÿ 1.35 0.50 Ÿ -0.35 -1.20 -2.04 Ÿ Ÿ Ÿ Ÿ -2.89 Ÿ Ÿ Ÿ Ÿ Ÿ Ÿ -3.74 -4.0 -2.8 -1.6 -0.4 0.8 2.0 3.2 4.4 5.6 6.8 8.0 Rysunek 22: Niektóre trajektorie w polu fazowym dla Brusselatora (ë = -4) fig:brusselator1 5.2 Generowniw liczb losowych 5.2.1 Funkcjarand Do tej pory funkcja ta sBu|yBa nam gBównie do wypeBniania liczbami losowymi naszych macierzy i wektorów... Funkcja ta u|ywa liniowego generatora nastpujco31 : ñø m = 231 òø Xn+1 = f(Xn) = (aXn + c) mod m, n e" 0, gdzie a = 843314861 óø c = 453816693 31 WedBug tego, co autor rozumiaB ogldanc kod 62 Jej okres jest z pewno[ci równy m (oznacza to, |e f jest permutacj cykliczn na przedziale [0, m - 1].) Zauwa|my, |e wszystkie generatory losowe w komputerze s .......... Aby przeksztaBci je do liczb rzeczywistych z przedziaBu [0, 1[, dzieli si otrzyman dan przez m (i otrzymuje si generator liczb rzeczywistych........). SkBadnik pocztkowy szeregu jest czsto zwany inicjatorem i przyjmuje domy[lnie warto[ X0 = 0. Zatem pierwsze wywoBanierand(pierwszy otrzymany wspóBczynnik je[li wypeBniamy macierz lub wektor) jest zawsze : u1 = 453816693/231 H" 0.2113249 Istnieje mo|liwo[ zmiany, w dowolnym momencie inicjatora za pomoc instrukcji : rand("seed",inicjator) gdzie inicjator jest zawarty w przedziale [0, m - 1]. Czsto istnieje potrzeba zainicjowania szeregu poprzez wybór inicjatora mniej lub bardziej przypadkowo (sytuacja taka, aby nie mi tej samen liczby za ka|dym razem), mo|emy chcie na przykBad otrzyma losowo dat lub godzine czy te| kombinaowa inicjatory<-??? Scilab dysponuje funkcj getdate która generuje wektor zBo|ony z 9 elementów. W[ród nich s : " drugi oznacza miesic (1-12), " szusty, dzieD miesica (1-31), " siudmy, godzin (0-23), " ósmy, minyty (0-59), " i dziewity, sekundy (0-61 ?). Aby otrzyma inicjator mo|na na przykBad dodawa powy|sze elementy midzy sob : v = getdate() rand("seed", sum(v([2 6 7 8 9]))) Zwrómy uwag równie| na mo|liwo[ zastpienia bierzcego inicjatora przez : germe = rand("seed") Poczwszy od rozkBadu jednostajnego na [0, 1[, mo|na otrzyma inne rozkBady a funkcjaranddostarcza dostateczny interface pozwalajcy otrzyma rozkBad normalny ([rednia 0 i wariancja 1). Aby przej[ od jednego do drugiego, stosuje sie nastpujc skBadnie : rand("normal") // aby otrzyma\ c rozk\l{}ad normalny rand("uniform") // aby powr\ oci\ c do rozk\l{}adu jednostajnego Domy[lnie generator przyjmuje rozkBad jednostajny ale jest rozsdnie w ka|dej symulacji upewni sie czyranddaje oczekiwany wynik u|ywajc jedn w tych instrukcji. Mo|na zraszt sprawdzi aktualny rozkad za pomoc : loi=rand("info") // rozk\l{}ad jest jeden z dw\ och mo\.zliwo\k{a}ci "jednostajny" lub "n Ponowne wywoBanierandmo|e przyj nastpujce formy : 1. A = rand(n,m)uzupeBnia macierzA(n,m) liczbami losowymi ; 2. je[li B jest macierz ju| zdefiniowan o wymiarach (n, m) to A = rand(B) daje ten sam efekt (pozwala unikn odzyskania<-???? wymiarów macierzyB) ; 3. wreszcie,u = rand()daje pojedyncz liczb losow. Do dwuch pierwszych metod mo|na doda argument aby wskaza dodatkowo rozkBad :A = rand(n,m,loi),A = rand(B,loi), gdzieloijest jednym z dwóch BaDcuchównormallubuniform. Wybrane zastosowania z u|yciem funkcji rand Poczwszy od rozkBadu jednostajnego, w prosty sposób mo|na otrzyma macierz (n,m) liczb wdBug : 1. rozkBad jednostajny na [a, b] : X = a + (b-a)*rand(n,m) 2. rozkBad jednostajny na liczbach caBkowitych z przedziaBu [n1, n2] : X = floor(n1 + (n2+1-n1)*rand(n,m)) 63 (losujemy nastpujce liczby rzeczywiste rozkBadu jednostajnego na przedziele rzczeywistym [n1, n2 + 1[ a nastpnie bie- rzemy ich cz[ caBkowit). Symulacja pojedynczej próby Bernulliego z prawdopodobieDstwem sukcesu p : succes = rand() < p otrzymujemy dziki prostej metodzie32 symulujcej rozkBad dwumianowy B(N, p) : X = sum(bool2s(rand(1,N) < p)) (bool2sprzekstaBca sykcesy w 1 i nie sumuje ich w funkcjisum). Z uwagi na fakt, |e wykonywanie iteracji przez Scilaba jest do[c wolne, mo|na otrzyma bezpo[rednio wektor (kolumnowy) skBadajcy si z m realizacji tego rozkBadu : X = sum(bool2s(rand(m,N) < p),"c") ale ko|y[tne bdzie u|ycie funkcjigrand, która u|ywa metod bardziej wyczynowej<-???. Z drugiej strony, je[li u|ywacie tego sposobu33, jest bardziej przej|y[cie aby kodowa go jako funkcj Scilaba. A oto prosta funkcja symulujce rozkBad geometryczny (ilo[ prób Bernulliego potrzebnych aby otrzyma sukces)34 : function [X] = G(p) // rozklad geometryczny X = 1 while rand() > p // pora\.zka X = X+1 end endfunction Wreszcie, w nawizaniu do rokBadu Normalnego N (0, 1), otrzymamy rozkBad Normalny N (µ, Ã2) ([rednia µ i odchylenie stan- dardowe Ã) za pomoc : rand("normal") X = mu + sigma*rand(n,m) // aby otrzymac macierz (n,m) tych liczb // lub przy uzyciu pojedynczej instrukcji : X = mu + sigma*rand(n,m,"normal") 5.2.2 Funkcjagrand W skomplikowanych symulacjach wyko|ystujcych wiele liczb losowych klasyczna funkcja rand z jej okresem rzdu 231(C" 2.147 109) mo|e okaza si troch za ciasna. Jest zatem wskazane u|yciegrandktóra równie| umo|liwia symulacj wszystkich klasycznych rozkBadów. grandu|ywa sie prawie w taki sam sposób corand, tzn. |e mo|na u|y jednej z dwóch nastpujcych skBadni (oczywi[cie dla drugiej macierzAmusi by zdefiniowana w momencie jaj wywoBania) : grand(n,m,loi, [p1, p2, ...]) grand(A,loi, [p1, p2, ...]) ooloioznacza BaDcuch znaków precyzujcy rozkBad po którym naspuj ewentualnie jego parametry. Kilka przykBadów (aby otrzyma prób n realizacji, pod postaci wektora kolumnowego) : 1. rozkBad jednostajny na liczbach caBkowitych z du|ego przedziaBu [0, m[ : X = grand(n,1,"lgi") gdzie m zale|y od generatora bazy (domy[lnie m = 232) ; 2. rozkBad jednostajny na liczbach caBkowitych z przedziaBu [k1, k2] : X = grand(n,1,"uin",k1,k2) (musi by speBniony warunek aby k2 - k1 d" 2147483561 bo w przeciwnym wypadku problem ten jest sygnalizowany jako bBd); 3. dla rozkBadu jednostajnego na przedziale [0, 1[ : 32 ale bardzo efektywnej dla duchych N ! 33 w ogólno[ci u|ywa bdziemy funkcjigrandpozwala otrzyma wiekszo[ klasycznych rozkBadów 34 ta funkcja jest bardziej efektywna dla maBych p. 64 X = grand(n,1,"def") 4. dla rozkBadu jednostajnego na przedziaBe [a, b[ : X = grand(n,1,"unf",a,b) 5. dla rozkBadu dwumianowego B(N, p) : X = grand(n,1,"bin",N,p) 6. dla rozkBadu geometrycznego G(p) : X = grand(n,1,"geom",p) 7. dla rozkBadu Poissona ze [redni µ : X = grand(n,1,"poi",mu) 8. dla rozkBadu wykBadniczego ze [redni » : X = grand(n,1,"exp",lambda) 9. dla rozkBadu normalnego ze [redni µ i odchylaniem standardowym à : X = grand(n,1,"nor",mu,sigma) Inne przykBady mo|na znalez na stronach pomocy. générateurs de base ??????.......... Pour opérer sur ces générateurs de base, vous pouvez utiliser les instructions suivantes : nom_gen = grand("getgen") permet de récupérer le (nom du) générateur courant grand([etgen",nom_gen) le générateurnom_gendevient le générateur courant etat = grand("getsd") permet de récupérer l état interne du générateur courant grand([etsd",e1,e2,..) impose l état interne du générateur courant La dimension de l état interne de chaque générateur dépend du type du générateur : de un entier pour urand U 624 entiers plus un index pour Mersenne Twister35. Si vous voulez refaire exactement la mme simulation il faut connaître l état initial (avant la simulation) du générateur utilisé et le sauvegarder d une façon ou d une autre. Exemple : grand("setgen","kiss") // kiss devient le générateur courant e = [1 2 3 4]; // etat que je vais imposer pour kiss // (il lui faut 4 entiers) grand("setsd",e(1),e(2),e(3),e(4)); // voila c est fait ! grand("getsd") // doit retourner le vecteur e X = grand(10,1,"def"); // 10 nombres s1 = sum(X); X = grand(10,1,"def"); // encore 10 nombres s2 = sum(X); s1 == s2 // en général s1 sera different de s2 grand("setsd",e(1),e(2),e(3),e(4)); // retour U l état initial X = grand(10,1,"def"); // de nouveau 10 nombres s3 = sum(X); s1 == s3 // s1 doit etre egal a s3 35 une procédure d initialisation avec un seul entier existe néanmoins pour ce générateur. 65 5.3 Dystrubuanty i ich odwrotno[ci Dystrybuanty s czsto u|ywane przy testach statystycznych (Ç2, ...) poniewa| pozwalaj na obliczenia, niech : r 1. dystrybuanta w 1 lub kilku punktach ; 2. jej odwrotno[ w 1 lub kilku punktach ; 3. jeden z parametrów rozkBadu jest dany przez pozostaBe i par (x, F (x)) ; W Helpie, dystrybuanty mo|na znalez pod hasBem Cumulative Distribution Functions... , wszystkie te funkcje poprzedzone s skrótemcdf. PrzykBadowo dla rozkBadu Normalnego N (µ, Ã2), funkcja która nas interesuje nosi nazwcdfnori jej skBadnia jest nastpujca : 1. [P,Q]=cdfnor("PQ",X,mu,sigma)aby otrzyma P = Fµ,Ã(X) i Q = 1 - P , X, mu oraz sigma mog by wektorami (o tych zamych wymiarach) i wówczas otrzymuje sie dlaPiQwektor za pomoc Pi = Fµ ,Ãi(Xi) ; i -1 2. [X]=cdfnor("X",mu,sigma,P,Q)aby otrzyma X = Fµ,Ã(P ) (podobnie jak poprzednio argumentami mog by -1 wektory o tych zamych wymiarach i wówczas otrzymuje si Xi = Fµ ,Ãi(Pi) ; i 3. [mu]=cdfnor("Mean",sigma,P,Q,X)aby otrzyma [redni ; 4. i ostatecznie[sigma]=cdfnor("Std",P,Q,X,mu)aby dosta odchylenie standardowe. Dwie ostatnie skBadnie funkcjonuj tak|e gdy argumentami s wektory o jednakowych wymiarach. Uwagi : " mo|liwo[ jednoczesnej pracy z p oraz q = 1 - p pozwala uzyska precyzj w obszarze gdzie p jest bliskie 0 lub 1. Dla p bliskiego 0 funkcja wyko|ystuje warto[ p, natomiast dla p bliskiego 1 funkcja wyko|ystuje warto[ q ; " laDcuch znaków pozwalajcy otrzyma funkcj odwrotn nie zawsze ma postaX... zobacznie na odpowiednich stronach pomocy. 5.4 Proste symulacje stochastyczne 5.4.1 Wprowadzenie i notacja Czsto symulacja polega przede wszystkim na otrzymaniu wektora : xm = (x1, ...., xm) którego skBadowe s rozumiane jako realizacje niezale|nych zmiennych losowych i podobnie rozkBad X1, X2, ...., Xm (bdziemy zapisywa przez X zmienn losow majc ten sam rozkBad <-???). W praktyce wektor xm otrzymuje sie bezpo[rednio lub po[rednio za pomoc funkcjirandlubgrand36. W przypadku próby xm poszukujemy zbli|onych charakterystyk jej rozkBadu takich jak warto[ oczekiwana, odchylenie stan- dardowe, dystrybuanta (lub funkcja gsto[ci) lub te| gdy stawiamy hipotez dotyczc rozkBadu, aby okre[li jej parametry, albo, gdy parametry s ju| znane, weryfikowa za pomoc testów statystycznych, |e nasza próba jest (mo|liwie) dobrze reprezenta- tywn zmienn losow która pod|a efektywnie za rozkBadem itd... <-??????????????????????????? Dla przypadków, które nas interesuj, znane s przewa|nie dokBadne wyniki teoretyczne, a symulacja sBu|y jedynie do zilustrowania wniosków, twierdzeD (lgn?, TCL =? CTG, itd...), dziaBania metody lub algorytmu, ... 5.4.2 PrzedziaBy ufno[ci Niekiedy o empirycznej warto[ci oczekiwanej otrzymanej z naszej próby (w Scilabie : x_bar_m = mean(xm)) chcieliby[my powiedzie, |e mie[ci si w pewnym przedziale. Ponadto chcieliby[my zna ów przedziaB aby móc zapisa, | : E[X] " Ic z prawdopodobieDstwem 1 - ± 36 w wikszo[ci przypadków próba statystyczna dotyczy miar fizycznych (temperatura, ci[nienie,...), biometrycznych (wzrost,waga), sonda|y, itd... otrzymane dane s zbierane w plikach (lub bazach danych); pewne programy (jak R) proponuj dla nich ukBad danych (dla których studenci bd mogli robi rcznie !) niestety Scilab nie dysponuje tak mo|liwo[cia; niemniej przypadków gdzie u|ywamy takich symulacji jest równie wiele, na przykBad aby przestudiowa Û zachowanie pewnych systemów wej[ciowych (lub zakBuceD) losowych lub podobnie aby rozwiza problemy czysto deterministyczne ale dla których metody analizy numerycznej s zbyt skomplikowane lub niemo|liwe do zaimplementowania<-??? . 66 gdzie najcz[ciej ± = 0.05 lub 0.01 (przedziaBy ufno[ci odpowiednio 95% i 99%). W celu wyznaczenia tych przedziaBów za punkt wyj[cia posBu|y na C.T.G.. Je[li przyjmiemy : m 1 ¯ Xm = Xi m i=1 ¯ za zmienn losow dla warto[ci [redniej (gdzie xm jest pojedyncz realizacj), wówczas prowo wielkich liczb mówi nam, |e Xm ¯ jest zbie|ne do E[X] i na mocy C.T.G (przy pewnych zaBo|eniach...) mamy : " b ¯ m(Xm - E[X]) 1 2 lim P (a < d" b) = " e-t /2dt m’!+" à 2À a (przyjmuje si, |e V ar[X] = Ã2). Dla dostatecznie du|ych m przyjmuje si, |e : " b ¯ m(Xm - E[X]) 1 2 P (a < d" b) H" " e-t /2dt à 2À a Je|eli chcemy aby przedziaB ufno[ci byB  symetryczny : a 1 2 " e-t /2dt = FN(0,1)(a) - FN(0,1)(-a) = 2FN(0,1)(a) - 1 = 1 - ± 2À -a wówczas : ± ± -1 -1 a± = FN(0,1)(1 - ) albo - FN(0,1)( ) 2 2 co zapisuje si w scilabie : a_alpha = cdfnor("X", 0, 1, 1-alpha/2, alpha/2) Otrzymujemy ostatecznie37 : a±Ã a±Ã E[X] " [xm - " , xm + " ] ¯ ¯ m m z prawdopodobieDstwem 1 - ± (je|eli przybli|enie granicy jest poprawne...). Problem pojawia si w momencie gdy nieznane jest odchylenie standardowe... Wyko|ystuje si wówczas etymator : m 1 ¯ Sm = (Xi - Xm)2 m - 1 i=1 Zastpujc odchylenie standardowe à przez sm (gdzie sm jest pojedyncz realizacj Sm), otrzymujemy przedziaB ufno[ci, który przyjmujemy równie| dla warto[ci empirycznej. Szczególne przypadki : 1. je[li Xi ma rozkBad normalny N(µ, Ã2) wówczas : ¯ " - µ Xm m <" t(m - 1) Sm gdzie t(k) ma rozkBad Studenta o k stopniach swobody. w tym przypadku wszelkie poprzednie przybli|enia trac moc (przybli|enie granicy i przybli|enie odchylenia standardowego) oraz otrzymuje sie : a±sm a±sm µ " [xm - " , xm + " ] avec probabilité 1 - ± ¯ ¯ m m gdzie sm jest empirycznym odchyleniem stndardowym z próby (sm = st_deviation(xm)w scilabie) gdzie a± jest warto[ci krytyczn rozkBadu Studenta : ± -1 a± = Ft(m-1)(1 - ) 2 otrzyman w scilacie38przez : a_alpha = cdft("T", m-1, 1-alpha/2, alpha/2) 2. w momencie gdy wariancja wyra|a si przez funkj warto[ci oczekiwanej (Bernouilli, Poisson, wykBadniczy,...) mo|na pzsBu|y die przybli|eniem odchylenia standardowego i otrzymuje si przedziaB ufno[ci poprzez rozwizanie odpowiedniej nierówno[ci. 37 dla przedziaBu z 95%, mamy a± C" 1.96 czsto zastpowane przez 2. 38 zobacz na odpowiednich stronach pomocy : wBa[ciwie funkcjecdfnie maj regularnej/uniwersalnej skBadni iXnie musi by zawsze oznaczeniem pozwa- lajcym otrzyma funkcj odwrotn do dystybuanty, tutaj jest toT! 67 5.4.3 Wykres dystrubuanty empirycznej Dystrybuant zmienne losowej X nazywamy funkcj : F (x) = Probabilité que X d" x Dystrybuanta empiryczna pochodzca z próby xm jest definiowana jako : m Fx (x) = card{xi <= x}/m m Jest to funkcja schodkowa, któr liczy si Batwo je[li posortujemy wektor xm w porzdku rosncym (mamy wic Fx (x) = i/m dla xi d" x < xi+1). Standardowym algorytmem sortujcym w scilabie jest funkcja sort która sortuje malejco39. Aby posortowa wektor xm w porzdku rosncym u|yjemy : xm = - sort(-xm) Aatwym sposobem narysowania wykresu jest funkcjaplot2d2. Oto przykBadowy kod : dystrybuanta_empiryczna(xm) // wykres dystrybuanty (empirycznej) // na pr\ obie xm m = length(xm) xm = - sort(-xm(:)) ym = (1:m) /m plot2d2(xm, ym, leg="dystrybuanta empiryczna") endfunction Zwrómy uwag na zapis : xm(:), który jest analogiczny do zapisywania wektora wierszowego. A teraz przykBad dla rozkBadu normalnego N (0, 1) : m = 100; xm = grand(m,1,"nor",0,1); xbasc() dystrybuanta_empiryczna(xm); // wykres fct dystrybuanty empirycznej // dane dla wykresu dystrybuanty "dokladnej" x = linspace(-4,4,100) ; y = cdfnor("PQ", x, zeros(x), ones(x)); plot2d(x, y, style=2) // nanosimy krzywa na pierwszy wykres xtitle("Dystrybyanty dokladna i empiryczna") który daje nastpujcy wykres (23). 5.4.4 Test Ç2 Niech xm = (x1, ..., xm) bdzie nasz prób dla dalszej analizy. Niech bdzie dana hipoteza H : zmienne losowe postaci (X1, ..., Xm) maj rozkBad L . Chcemy wiedzie czy hipoteza jest prawdziwa czy nie. Dla naszej próby mo|na bez wtpienia obliczy statystyki elementarne ([redni i odchylenie standardowe empiryczne) i je|eli s one dostatecznie bliskie warto[ci ocze- kiwanej oraz odchyleniu standardowemu rozkBadu L, mo|na wówczas zastosowa test statystyczny. Test Ç2 stosuje si dla roz- kBadów dyskretnych o skoDczonej liczbie warto[ci. Na przykBad zakBadajc, |e rozkBad L jest zadany przez {(vi, pi), 1 d" i d" n}. Test polega na obliczeniu wielko[ci : n (oi - mpi)2 i=1 y = mpi gdzie oi jest liczb realizacji xj równych vi i na przyrównaniu otrzymanej wielko[ci y do warto[ci progowej y±, test bdzie pozytywny40 je[li y d" y±. Je|eli do rówanania na y wstawiwy w miejsce próby x1, ..., xm, zmienne losowe X1, ..., Xm w ten sposób zdefiniujemy41 zmienn losow Y , któr mo|na przybli|a (dla dostatecznie du|ych m) rozkBadem Ç2 o n - 1 stopniach swobody. Warto[ progow otrzymujemy przez : -1 y± = F (1 - ±) Ç2 n-1 przy ± = 0.05 lub ± = 0.01. W Scilabie warto[ t otrzymamy poprzez : 39 zobacz tak|e funkcjgsort, która robi wicej rzeczy 40 z intuicyjnego punktu widzenia je|eli hipoteza jest dobra oczekujemy, |e oi nie odbiega zbytnio od mpi, zatem je|eli hipoteza jest faBszywa oczekujemy |e otrzymamy wysok warto[ y, skd odrzucimy hipotez dla y > y±... 41 poprzez zmienn losow wektorow O = (O1, ..., On) majc rozkBad wielomianowy<-????? 68 Fonctions de répartition exacte et empirique 1.0 1.0 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 repartition empirique Rysunek 23: Systrybuanta dokBadna i empiryczna rozkBadu normalnego fig:repart y_progowe = cdfchi("X", n-1, 1-alpha, alpha) Aby obliczy czsto[ oi za pomoc Scilaba mo|emy u|y nastpujcej metody : occ = zeros(n,1); for i=1:n occ(i) = sum(bool2s(xm == v(i))); // albo length(find(xm==v(i))) end if sum(occ) ~= m then, error("b\l{}\k{a}d oblicze\ n"), end I aby otrzyma wielko[ y mo|na napisa (stosujc zapis wktorowy42) : y = sum( (occ - m*p).^2 ./ (m*p) ) pod warunkiem, |ep(wektor przwdopodobieDstw rozkBadu L), bdzie tych samych wymiarów coocc(tutaj wektor kolumnowy ??????????). Uwagi : " przybli|anie rozkBadem Ç2 jest wa|ne gdy m jest dostatecznie du|e... |adamy czsto mpmin > 5 (pmin = mini pi) jako minimalny warunek zastosowania tego testu; w ten sposób mo|ecie weryfikowa to zaBo|enie i napisa wiadomo[ aby uprzedzi u|ytkownika je|eli nie jest speBnione. " mo|na Batwo pogrupowa te obliczenia w funkcji; " dla rozkBadu cigBego test mo|e si stosowa grupujc wielko[ci na przedziaBy, na przykBad dla rozkBadu U[0,1], stosuje si n równomiernie rozBo|onych przedziaBów; podobnie dla rozkBadu dyskretnego o nieskoDczonej liczbie wielko[ci, mo|na pogrupowa je w kolejki rozkBadu; podobn procedur mo|na zastosowa dla roskBadów skoDczonych dla których warynek zastosowania testu nie jest speBniony. " je|eli testujecie rozkBad w którym pewne parametry zostaBy wcze[niej obliczone na podstawie danych, nale|y okre[li liczb stopni swobody dla rozkBadu Ç2 ; na przykBad je[li spodziewamy si rozkBadu B(n - 1, p) i u|yjemy p = xm/(n - 1) ¯ -1 wówczas nieprzekraczaln warto[ci progow dla testowanej hipotezy zerowej bdzie y± = F (1 - ±) a nie y± = Ç2 n-2 -1 F (1 - ±). Ç2 n-1 42 wiczenie: przeksztaB t instrukcj wektorow aby zrozumie jak (i dlaczego) dziaBa. 69 5.4.5 Test KoBmogorowa-Smirnova Niech X bdzie rzeczywist zmienn losow dla której dystrybuanta rozkBadu jest funkcj cigB F oraz X1, X2,..., Xm, m niezale|nych zmiennych losowych. Dla realizacji Xi (mówimy wektor xm = (x1, ...., xm)), mo|na skonstruowa dystrybuant empiryczn która przy (m ’! +") dzy do dystrybuanty teoretycznej. Test KS polega na okre[leniu ró|nicy pomidzy dystrybu- ant teoretyczn i empiryczn (otrzyman na podstawie naszej próby (x1, ...., xm)) w nastpujcy sposób : " m km = m sup |F (x) - Fx (x)| -"<x<+" i na porównaniu jej z wielko[ci dopuszczaln . Je|eli zastpimy nasz realizac przez odpowiednie zmienne losowe, wówczas km staje si równie| zmienn losow (któr zapisujemy jako Km). Zgodnie z teori jej dystrybuanta jej rozkBadu jest nastpujca : +" 2 lim P (Km d" x) = H(x) = 1 - 2 (-1)j-1e-2j x2 m’!+" j=1 Jak przy te[cie Ç2, je|eli rozwa|ana hipoteza jest faBszywa, otrzymane warto[ci km bd du|e i odrzycimy t hipoteze gdy: km > H-1(1 - ±) 2 z ± = 0.05 lub 0.01 na przykBad. Je[li u|ywamy przybli|enia H(x) C" 1 - 2e-2x wówczas nieprzekraczalna warto[ progowa jest : 1 2 kseuil = ln( ) 2 ± Obliczenie km nie sprawia problemu je|eli posortuje si wektor (x1, x2, . . . , xm). Sortowaniw bdzie efektywne gdy zwrócimy uwag ne fakt aby : i i m m sup Fx (x) - F (x) = - F (xi), et sup F (x) - Fx (x) = F (xi+1) - m m x"[xi,xi+1[ x"[xi,xi+1[ dwie nastpne wielko[ci mo|na Batwo policzy : " " j + m km = m sup (Fx (x) - F (x)) = m max ( - F (xj)) 1d"jd"m m -"<x<+" " " j - 1 - m km = m sup (F (x) - Fx (x)) = m max (F (xj) - ) 1d"jd"m m -"<x<+" + - i otrzymujemy w ten sposób km = max(km, km). 5.4.6 wiczenia Czy to sprawka kostki ? Rzucamy 200 razy symetryczn kost do gry, do[wiadczenie jest nastpujce : rzucamy tyle razy a| wypadnie 1 (ale nie wicej ni| 10 rzutów). Otrzymali[my nastpujce wyniki : liczba rzutów 1 2 3 4 5 6 7 8 9 10 e" 11 ilo[ do[wiadczeD 36 25 26 27 12 12 8 7 8 9 30 na przykBad 36 razy jedynka pojawiBa sie w pierwszym rzycie, 25 razy jedynka pojawiBa sie w drugim |ucie, itd... Wykona test Ç2 aby udzieli odpowiedzi na postawione pytanie. Urna Polya Losujemy N razy z urny Polya zawierajcej pocztkowo r kul czerwonych i v kul zielonych. Ka|de losowanie polega na wycigniciu jednej kuli i jej odBo|eniu spowrotem do urny wraz z c kulami tego samego koloru. Przez Xk oznaczamy stosunek kul zielonych po k losowaniach do sumy kul, Vk oznacza ilo[ kul zielonych : v X0 = , V0 = v. v + r Je|eli v = r = c = 1, wynik bdzie nastpujcy : 1. E(XN ) = E(X0) = X0 = 1/2 ; 1 N+1 2. XN ma rozkBad normalny na zbiorze { , ....., } ; N+2 N+2 70 3. dla N ’! +", XN jest zbie|ne do rozkBadu jednostajnego na przedziale [0, 1). Wyko|ystacie symulacj aby zilustrowa dwa pierwsze wyniki. 1. Aby przeprowadzi ró|ne symulacje, mo|na napisa funkcj zale|n od parametru N i która wykonuje N kolejnych loso- waD. Funkcja ta zwraca wic XN i VN : function [XN, VN] = Urna_Polya(N) // symulacja N losowan z "Urny Polya" : VN = 1 ; V_plus_R = 2 ; XN = 0.5 for i=1:N u = rand() // losowanie jednej kuli V_plus_R = V_plus_R + 1 // zwieksza ilosc kul if (u <= XN) then // wylosowalislmy kule zielona : mamy XN prawdopodobienstwo // wylosowania kuli zielonej (1 - XN dla kuli czerwonej) VN = VN + 1 end XN = VN / V_plus_R // aktualizujemy stosunek kul zielonych end endfunction niestety skuteczne wykonanie statystyk wymaga wielokrotnego wywoBania tej funkcji, a zatem, z uwagi na stosunkowo powolne wykonywanie iteracji przez Scilab, nale|y naposa funkcj wykonujc m procesów równolegBych . Mo|na u|y funkcjifindaby naprawi[????? urny z których losujemy kule zielone. 2. Napisa skrypt znajdujcy poprzez symulac warto[ oczekiwan (wraz z jej przedzieBem ufno[ci). 3. Rozwin poprzedni skrypt testujc hipotez H dotyczc rozkBadu zmiennej losowej XN H : XN ma rozkBad jednostajny 1 N+1 na zbiorze { , ....., } za pomoc testu Ç2. N+2 N+2 4. Porówna wyniki graficznie, na przykBad rysujc na jednym rysunku dystrybuant empiryczn i teoretyczn, lub te| nary- sowa wykres funkcji gsto[ci rozkBadu Ç2 dla otrzymanych wielko[ci przez ten test, wielko[ci progowe zaznaczy liniami pionowymi. Most ????? Proces stochastyczny (w którym U oznacza rozkBad jednostajny na przedziale [0, 1]) : n 1 " Xn(t) = (1{U d"t} - t) i n i=1 jest taki, |e dla ustalonego t "]0, 1[ mamy : lim Xn(t) = Y (t) <" N (0, t(1 - t)) n’!+" Chcemy zilustrowa wyniki za pomoc symulacji. Schemat pracy : 1. Napisa funkcj Scilaba function [X] = pont_brownien(t,n) pozwalajc otrzyma realizacj Xn(t) ; w cigu wywoBamy t funkcj m razy z warto[ci n dostatecznie du|ym43 nale|y zatem napisa j bez u|ycia pentli. 2. Napisa sktypt scilaba wyko|ystujc t symulacj mostu Brownien : wykona m symulacji Xn(t) (z du|ym n) i przedsta- wi graficznie zachowanie rozkBadu (rysujc jego dystrybuant empiryczn i nakBadajc na ni dystrybuant teoretyczn Y (t)) a nastpnie zastosowa test KoBmogorowa-Smirnova. Mo|na umie[ci poprzedni funkcj na pocztku skryptu aby pracowa tylko z jednym plikiem. Uwaga : nie jest Batwo dobra odpowiednie wielko[ci dla m i n : kiedy m wierzy ...?????????????.......gdzie test si nie powiódB poniewa| uznaB, |e n nie jest dostatecznie du|e (poniewa| rozkBad normalny otrzymuje si z granicy gdy n ’! +"), nale|y zatem zwikszy n... Dla t = 0.3, mo|ecie na przykBad przetestowa z n = 1000 i m = 200, 500, 1000, i test daje dobre wyniki (|adko od|uca hipotez zerow) ale dla m = 10000 test jest prawie zawsze negatywny. Je|eli zwikszymy n (na przykBad do n = 10000 ale obliczenia zaczynaj by nieco dBu|sze na komputerze PC wykonanym w 2003) wówczas test ponownie staje si normalnie pozytywny. 43 aby przybli|y Y (t) ! 71 6 Zbiór drobiazgów W tej cz[ci przedstawiono/przytoczono przegld niektórych bBdów czsto spotykanych w Scilabie... 6.1 Definiowanie wektora i macierzy wspóBczynnik po wspóBczynniku Ten bBd jest jednym z najczstrzych. Rozwa|my nastpujcy skrypt : K = 100 // jedyny parametr w tym skrypcie for k=1:K x(k) = cos y(k) = cos innego end plot(x,y) Je|eli uruchomimy ten skrypt po raz pierwszy, zostan zdefiniowane dwa wektoryxiy. Pojawia si jeden maBy bBd poniewa| w ka|dej iteracji Scilab redefiniuje rozmiary tych wektorów (nie wie, |e ich wymiarem koDcowym bdzie (K,1)). Zauwa|my równie|, |e domy[lnie stworzy wektor kolumnowy. Podczas drugiero wywoBania (zmieniajc parametrK...) wektoryxiys znane i takie, |e k jest mniejsze od 100 (pocztkowa wielko[ parametru K) zatem zmieniane s rówenie| ich wspóBrzdne. W konsekwencji je|eli nowa warto[Kjest taka, |e : " K < 100wówczas nasze wektoryxiymaj zawsze 100 wspóBrzdnych (zmodyfikowane jest tylkoKpierwszych) a wykres nie pokazuje tego co oczekujemy ; " K > 100 nie pojawia si |aden problem (za wyjtkiem tego, |e rozmiar wektora jest za ka|dym razem aktualizowany poczwszy od 101 iteracji) Najlepszym sposobem jest zdefiniowanie wektorówxiyza pomoc inicjalizacji ich rodzaju : x = zeros(K,1) ; y = zeros(K,1) w ten sposób uniknie si wszelkich bBdów. Nasz skrypt zatem przyjmie posta : K = 100 // jedyny parametr w tym skrypcie x = zeros(K,1); y = zeros(K,1); for k=1:K x(k) = cos y(k) = cos innego end plot(x,y) 6.2 Na temat warto[ci zwracanych przez funkcj ZaBó|my, |e mamy zaimplementowan funkcj w Scilabie zwracajc dwa argumenty, na przykBad : function [x1,x2] = resol(a,b,c) // rozwi\k{a}zanie r\ ownania kwadratowego a x^2 + b x + c = 0 // formu\l{}ad poprawiona dla wi\k{e}kszo\ sci metod numerycznych // (w celu unikni\k{e}cia odejmowania dw\ och s\k{a}siednich liczb) if (a == 0) then error(" nie rozwa\.zamy przypadku gdy a=0 !") else delta = b^2 - 4*a*c if (delta < 0) then error(" nie rozwa\.zamy przypadku gdy delta < 0 ") else if (b < 0) then x1 = (-b + sqrt(delta))/(2*a) ; x2 = c/(a*x1) else x2 = (-b - sqrt(delta))/(2*a) ; x1 = c/(a*x2) end end end endfunction 72 W ogólnym przypadku je|eli wywoBamy funkcj w Scilab-ie w nastpujcy sposób : -->resol(1.e-08, 0.8, 1.e-08) ans = - 1.250D-08 jej wynik jest przypisany zmiennejans. Ale zmienna ta jest pojedyncz liczb, a funkcja zwraca dwie warto[ci z których tylko pierwsza jast przypisana zmienneans. Aby zobaczy drug warto[ u|yjemy skBadni : -->[x1,x2] = resol(1.e-08, 0.8, 1.e-08) x2 = - 80000000. x1 = - 1.250D-08 Inna, niebezpieczniejsza puBapka jest nastpujca. ZaBó|my, |e znamy precyzj z jak dziaBa ta formuBa (w stosunku do precyzji dziaBania formóB klasycznych). Aby dobra warto[ci a i c (na przykBad a = c = 10-k przyjmujc ró|ne warot[ci k) obli- czymy wszystkie pierwiastki dla kolejnych parametrów równania i ustawimy je jako dwa wektory sBu|ce do pózniejszej analizy. Najprostrzy schemat jest nastpujcy : b = 0.8; kmax = 20; k = 1:kmax; x1 = zeros(kmax,1); x2=zeros(kmax,1); for i = 1:kmax a = 10^(-k(i)); // c = a [x1(i), x2(i)] = resol(a, b, a); // BLAD ! end Taka skBadnia nie zadziaBa (jedynie w przypadku gdy funkcja zwraca jedn warto[). Nale|y wykona to w dwóch krokach : [rac1, rac2] = resol(a, b, a); x1(i) = rac1; x2(i) = rac2; Uwaga : Problem ten jest rozwizany w wercjach rozwijajcych (cvs) dla Scilaba. 6.3 Funkcja zostaBa zmidyfikowana ale... wszystko dziaBa tak jak przed modyfikacj ! By mo|e zapomnieli[cie zapisa zmiany w waszym edytorze albo, co jest bardziej prawdopodobne, zapomnieli[cie zaBadowa plik zawierajcy funkcj w Scilabie za pomoc instrukcjigetf(lubexec) ! Jedna maBa wskazówka : instrukcja getf (lub exec) jest z pewnoci zapisana niedaleko w historii komend, wystarczy zatem za pomoc strzaBki ‘! odszuka t instrukcj. 6.4 Problem zrand Domy[lnie funkcja rand dostarcza liczby losowe wedBug rozkBadu jednostajnego na przedziale [0, 1[, mo|na jednak otrzyma rozkBad normalny N (0, 1) za pomoc :rand(Dormal"). Chcc powróci do rozkBadu jednostajnego nale|y wpisa instrukcj rand("uniform"). Aby unikn tego problemu najlepiej pamita aby ka|de wywoBanie funkcjirand. precyzowaBo rozkBad który nas aktualnie interesuje (patrz poprzedni rozdziaB). 6.5 Wektory wierszowe, wektory kolumnowe... W kontek[cie macierzowym ich posta jest sprecyzowana, ale dla innych zastosowaD nie zawsze s one rozrózniane i stosuje sie funkcj która dziaBa w obydwu przypadkach. Aby przyspieszy obliczenia dobrze jest odwoBa si do skBakni macierzowej i wybra jedn lub drug form wektora. Mo|na wyko|ysta funkcjmatrixw nastpujcy sposób : x = matrix(x,1,length(x)) // aby otrzyma\ c wektor wierszowy x = matrix(x,length(x),1) // aby otrzyma\ c wektor kolumnowy Chcc otrzyma w prosty sposób wektor kolumnowy mo|na wyko[ysta instrukcj : x = x(:) 73 6.6 Operator porównania W pewnyc przypadkach Scilab akceptuje symbol=jako operator porównania : -->2 = 1 Warning: obsolete use of = instead of == ! ans = F ale lepiej jest zawsze u|ywa symbolu==. 6.7 Liczby zespolone a liczby rzeczywiste Scilab jest tak skonstruowany aby traktowa liczby rzeczywiste w taki sam sposób jak liczby zespolone ! Jest to calkiem prak- tyczne ale mo" niekiedy by przyczyn niespodzianek, na przykBad podczas szacowania funkcji rzeczywistej poza jej dziedzin |e (na przykBad x i log(x) dla x < 0, acos(x) i asin(x) dla x " [-1, 1], acosh(x) dla x < 1) ) poniewa| Scilab zwraca wówczas / oszacowanie cz[ci zespolonej tej funkcji. Aby dowiedzie si czy dziaBamy na zmiennej zespolonej czy rzeczywistej nale|y u|y funkcjiisreal: -->x = 1 x = 1. -->isreal(x) ans = T -->c = 1 + %i c = 1. + i -->isreal(c) ans = F -->c = 1 + 0*%i c = 1. -->isreal(c) ans = F 6.8 Proste instrukcje a funkcje Scilaba W niniejszym opracowaniu czstu stosuje sie zamiennie terminów prosta instrukcja i funkcja aby wskaza procedury dostpne w bie|cej werksji Scilaba. Istnieje jednak fundamentalna ró|nica pomidzy prost instrukcj która jest kodowana w fortranie 77 lub w C a funkcj (zwan równie| makro) która jest kodowana w jzyku Scilaba : funkcja jest wówczs rozumiana jako zmienna Scilaba, i dziki temu mo|na j u|y jako argumentu dla inne funkcji. Poczwszy od wersji 2.7, proste instrukcje bywaj zmiennymi w Scilabie (fptr). Niemniej jednak przysparzaj wiele problemów (aktualnie rozwizanych w rozwiniciach). Oto przykBad problemu jaki mo|na spotka : funkcja nastpujca pozwalajca przybli|a zbiór przy zastosowaniu metody Monte Carlo : function [im,sm]=MonteCarlo(a,b,f,m) // /b // przybli\.zenie | f(x) dx przy zastosowaniu metody Monte Carlo // /a // zwracane jest r\ ownie\.z empiryczne odchylenie standardowe xm = grand(m,1,"unf",a,b) ym = f(xm) 74 im = (b-a)*mean(ym) sm = (b-a)*st_deviation(ym) endfunction Jako argument f oczekuje si funkcji Scilaba ale ten kod powienien równie| dziaBa dla funkcji matematycznych które nale| do prostych instrukcji Scilaba44 poniewa| s one teraz rozwa|ane jako zmienne. Niemniej jednak test z u|yciem prostej instrukcji expnie powiedzie si : -->[I,sigma]=MonteCarlo(0,1,exp,10000) // bug ! WywoBujc funkcjMonteCarloza pomocgetfz opcj nie kompilowania : -->getf("MonteCarlo.sci","n") bBd ten [bug] (skorygowany w aktualnym cvs) bdzie wówczas ominity. 6.9 Obliczanie wyra|eD logicznych W przeciwieDstwie do jzyka C, obliczenie warto[ci logicznej wyra|eD : a lub b a i b poprzedzone jest wyliczeniem warto[ci logicznej a i b, dopiero potem zostanie wykonana na nich operacja sumy czy iloczynu logicznego (Priorytety operatorów zamieszczono w rozdziale Programowanie ). A Odpowiedzi do wiczeD z rozdziaBu 2 1. --> n = 5 // dla ustalenia warto\ sci n... --> A = 2*eye(n,n) - diag(ones(n-1,1),1) - diag(ones(n-1,1),-1) Szybszym sposobem jest u|ycie funkcjitoeplitz: -->n=5; // pour fixer une valeur a n... -->toeplitz([2 -1 zeros(1,n-2)]) 2. Je|eli A jest macierz (n, n),diag(A)zwróci wektor kolumnowy zawierajcy elementy diagonalne macierzy A (a zatem otrzymamy wektor kolumnowy o wymiarze n). tt diag(diag(A)) zwróci macierz kwadratow diagonaln o wymiarze n, w której elementy diagonalne bd takie jak w macierzy wyj[ciowej. 3. Oto jedna z mo|liwo[ci : --> A = rand(5,5) --> T = tril(A) - diag(diag(A)) + eye(A) 4. (a) --> Y = 2*X.^2 - 3*X + ones(X) --> Y = 2*X.^2 - 3*X + 1 // wyko\.zystuj\k{a}c zapis skr\ ocony --> Y = 1 + X.*(-3 + 2*X) // oraz schemat Hornera (b) --> Y = abs(1 + X.*(-3 + 2*X)) (c) --> Y = (X - 1).*(X + 4) // wyko\.zystuj\k{a}c zapis skr\ ocony (d) --> Y = ones(X)./(ones(X) + X.^2) --> Y = (1)./(1 + X.^2) // z skr\ otami 5. Oto skrypt : 44 na przykBadsin,expi wiele innych. 75 n = 101; // dyskretyzacja x = linspace(0,4*%pi,n); y = [1 , sin(x(2:n))./x(2:n)]; // aby zapobiec dzieleniu przez zero... plot(x,y,"x","y","y=sin(x)/x") 6. Mo|liwe rozwizanie (dla ró|nych warto[ci n) : n = 2000; x = rand(1,n); xbar = cumsum(x)./(1:n); plot(1:n, xbar, "n","xbar", ... "ilustracja lgn : xbar(n) -> 0.5 qd n -> + oo") B Rozwizania wiczeD z rozdziaBu 3 1. Klasyczny algorytm wyko|ystucy dwie ptle : function [x] = sol_tri_sup1(U,b) // // rozwi\k{a}zanie Ux = b gdzie U jest macierz\k{a} tr\ ojk\k{a}tn\k{a} g\ orn\k{a // [n,m] = size(U) // kilka wyj\k{a}tk\ ow .... if n ~= m then error( Macierz nie jest kwadratowa ) end [p,q] = size(b) if p ~= m then error( Wyraz wolny ma nieprawid\l{}owy wymiar ) end // pocz\k{a}tek algorytmu x = zeros(b) // rezerwuje miejsce dla x for i = n:-1:1 suma = b(i,:) for j = i+1:n suma = suma - U(i,j)*x(j,:) end if U(i,i) ~= 0 then x(i,:) = suma/U(i,i) else error( Macierz jest nieodwracalna ) end end endfunction A oto wersja wykorzystujca pojedyDcz ptl function [x] = sol_tri_sup2(U,b) // // [n,m] = size(U) // niekt\ ore wyj\k{a}tki .... if n ~= m then error( Macierz nie jest kwadratowa ) end 76 [p,q] = size(b) if p ~= m then error( Wektor wyraz\ ow wolnych ma z\l{}y wymiar) end // pocz\k{a}tek algorytmu x = zeros(b) // rezerwuje miejsce dla x for i = n:-1:1 suma = b(i,:) - U(i,i+1:n)*x(i+1:n,:) // zobacz komentarz ko\ ncowy if U(i,i) ~= 0 then x(i,:) = suma/U(i,i) else error( Nie mo\.zna odwr\ oci\ c macierzy ) end end endfunction Komentarz : w pierwszej iteracji (dla i = n) macierzeU(i,i+1:n)etx(i+1:n,:)s puste. Odpowiadaj one obiek- tom zdefiniowanym w Scilabie (macierz pusta), które oznaczone s przez[]. Dodawanie macierzy pustej jest zdefiniowane i daje : A = A + []. Zatem w pierwszej iteracji mamysuma = b(n,:) + []to znaczysuma = b(n,:). 2. // // skrypt rozwiazujacy x" + alpha*x + k*x = 0 // // Aby przeksztalcic rownanie do postaci ukladu rownan pierwszego rzedu // kladziemy : X(1,t) = x(t) i X(2,t) = x (t) // // Otrzymujemy wowczas : X (t) = A X(t) z : A = [0 1;-k -alpha] // k = 1; alpha = 0.1; T = 20; // moment koncowy n = 100; // dyskretyzacja czasowa : przedzial [0,T] zostanie // podzielony na n przedzialow t = linspace(0,T,n+1); // momenty czasowe : X(:,i) bedzie korespondowal z X(:,t(i)) dt = T/n; // zmiana czasu A = [0 1;-k -alpha]; X = zeros(2,n+1); X(:,1) = [1;1]; // warunki poczatkowe M = expm(A*dt); // oblicza exponent A dt // obliczenia for i=2:n+1 X(:,i) = M*X(:,i-1); end // wyswitlenie rezultat\ ow xset("window",0) xbasc() xselect() plot(t,X(1,:), czas , pozycja , Courbe x(t) ) xset("window",1) xbasc() xselect() plot(X(1,:),X(2,:), pozycja , predkosc , Trajektoria w przestrzeni stanow ) 3. function [i,info]=przedzial(t,x) // poszukuje dychotomii przedzialu i takiej, ze: x(i)<= t <= x(i+1) // jezeli t nie jest z przedzialu [x(1),x(n)] zwraca info = %f n=length(x) if t < x(1) | t > x(n) then 77 info = %f i = 0 // kladziemy wartosc domyslna else info = %t i_dolne=1 i_gorne=n while i_gorne - i_dolne > 1 itest = floor((i_gorne + i_dolne)/2 ) if ( t >= x(itest) ) then, i_dolne= itest, else, i_gorne=itest, end end i=i_dolne end endfunction 4. function [p]=myhorner(t,x,c) // rozwiniecie wielomianu c(1) + c(2)*(t-x(1)) + c(3)*(t-x(1))*(t-x(2)) + ... // algorytmem Hornera // t jest wektorem (lub macierza) moment\ ow czasu n=length(c) p=c(n)*ones(t) for k=n-1:-1:1 p=c(k)+(t-x(k)).*p end endfunction 5. Rozwinicie w szereg Fouriera : function [y]=signal_fourier(t,T,cs) // ta funkcja zwraca T-okresowy sygnal // t : wektor momentow dla kt\ orych obliczmy // sygnal y ( y(i) odpowiadaj\k{a}cy t(i) ) // T : okres sygnalu // cs : jest wektorem ktory wskazuje amplitude kazdej funkcji f(i,t,T) l=length(cs) y=zeros(t) for j=1:l y=y + cs(j)*f(j,t,T) end endfunction function [y]=f(i,t,T) // wielomiany trygonometryczne dla sygna\l{}u o okresie T : // jesli i jest parzyste : f(i)(t)=sin(2*pi*k*t/T) (z k=i/2) // jesli i jest nieparzyste : f(i)(t)=cos(2*pi*k*t/T) (z k=floor(i/2)) // stad w szczeg\ olnosci f(1)(t)=1 jest traktowany ponizej // jako przypadek szczegolny i niezbyt wa\.zny // t jest wektorem mament\ ow czasowych if i==1 then y=ones(t) else k=floor(i/2) if modulo(i,2)==0 then y=sin(2*%pi*k*t/T) else y=cos(2*%pi*k*t/T) 78 end end endfunction 6. Spotkanie : niech TA i TB bd momentami przybycia Pana A i Pani B. S to dwie zmienne losowe niezale|ne o rozkBadzie U([17, 18]) a spotkanie ma miejsce jezeli [TA, TA + 1/6] *" [TB, TB + 1/12] = ". Doswiadczenie spotkania odpowiada realizacji zmienne losowej wektorowej R = (TA, TB) która, zgodnie z hipotezami, ma rozkBad zero-jedynkowy na kwadra- cie [17, 18] × [17, 18]. PrawdopodobieDstwo spotkania wyznacza sie wic w oparciu o zale|no[ pomidzy powie|chni obszaru (odpowiadajc pojedynczemu spotkaniu) zdefiniowan jako: ñø TA d" TB + 1/12 òø TB d" TA + 1/6 óø TA, TB " [17, 18] oraz powie|chni kwadratow (1), przez co otrzymuje sie p = 67/288. Aby obliczy to prawdopodobieDstwo, mo|na przyj ze spotkanie powinno miec miejsce pomidzy poBódniem a godzin pierwsz. Poprzez symulacj otrzymuje sie m do[wiadczeD a prawdopodobieDstwo empiryczne jest liczb przypadków w których spotkanie ma miejce podzielon przez m. A oto rozwizanie : function [p] = rdv(m) tA = rand(m,1) tB = rand(m,1) spotkanie = tA+1/6 > tB & tB+1/12 > tA; p = sum(bool2s(spotkanie))/m endfunction mo|na tak|e wyko|ysta funkcjegrandopisan w rozdziale dotyczcym zastosowaD : function [p] = rdv(m) tA = grand(m,1,"unf",17,18) tB = grand(m,1,"unf",17,18) spotkanie = tA+1/6 > tB & tB+1/12 > tA; p = sum(bool2s(spotkanie))/m endfunction C Rozwizania wiczeD z rozdziaBu 4 Czy to sprawka kostki? Je|eli przez J oznaczymy zmienn losow bdc wynikiem do[wiadczenia, wówczas J ma rozkBad geometryczny G(p) gdzie p = 1/6 przy zaBo|eniu, |e kostka jest symetryczna. W ten sposób mo|emy okre[li prawdopodobieDstwo poszczególnych zmiennych losowych jako (przymujemy q = 1 - p = 5/6) : P (J = 1) = p, P (J = 2) = qp, P (J = 3) = q2p, ...., P (J = 10) = q9p, P (J > 10) = q10 W drugiej cz[ci tabeli podano czsto[ wystpowania poszczególnych zmiennych losowych, co w skrypcie zapiszemy nastpujco : occ= [36 ; 25 ; 26 ; 27 ; 12 ; 12 ; 8 ; 7 ; 8 ; 9 ; 30]; p = 1/6; q = 5/6; pr = [p*q.^(0:9) , q^10] ; y = sum( (occ - m*pr).^2 ./ (m*pr) ); y_progowy = cdfchi("X", 10, 0.95, 0.05); mprintf("\n\r Test dla :") mprintf("\n\r y = %g, et y_progowy (95 \%) = %g", y, y_progowy) mprintf("\n\r 200*min(pr) = %g", 200*min(pr)) Otrzmujemy y = 7.73915 gdy yseuil = 18.307, zatem kostka wydaje si by prawidBowa. Niemniej jednak mamy 200 × min(pr) C" 6.5 co jest znacznie poni|ej warunku pozwalajcego zastosowa test (nale|aBoby wykona jeszcze kilka dodatkowych do[wiadczeD). 79 Urna Polya *Funkcja scilaba function [XN, VN] = Urna_Polya_rownolegla(N,m) VN = ones(m,1) ; V_plus_R = 2 ; XN = 0.5*ones(m,1) for i=1:N u = grand(m,1,"de"f) // losowanie jednej kuli V_plus_R = V_plus_R + 1 // dodaje kule ind = find(u <= XN) // znajduje numer urny z kt\ orej // wylosowano kule zielona VN(ind) = VN(ind) + 1 // zwiekszamy liczbe kul zielonych w tej urnie XN = VN / V_plus_R // aktualizujemy proporcje kul zielonych end endfunction *Skrypt Scilaba // symulacja polya : N = 10; m = 5000; [XN, VN] = Urna_Polya_rownolegla(N,m); // 1/ obliczenie wartosci oczekiwanej i przedzialu ufnosci EN = mean(XN); // wartosc oczekiwana empiryczna sigma = st_deviation(XN); // odchylenie standardowe empiryczne delta = 2*sigma/sqrt(m); // zaokraglenie dla przedzialu empirycznego (2 dla 1.9599..) mprintf("\n\r E teoretyczne = 0.5"); mprintf("\n\r E oszacowane = %g",EN); mprintf("\n\r Przedzial (empiryczny) ufnosci przy 95\% : [%g,%g]",EN-delta,EN+delta); // 2/ test chi2 alpha = 0.05; p = 1/(N+1); // prawdopodobienstwo kazdego wyniku (rozklad jednostajny) occ = zeros(N+1,1); for i=1:N+1 occ(i) = sum(bool2s(XN == i/(N+2))); end if sum(occ) ~= m then, error(" Problem..."), end // mala poprawka Y = sum( (occ - m*p).^2 / (m*p) ); Y_seuil = cdfchi("X",N,1-alpha,alpha); mprintf("\n\r Test chi 2 : ") mprintf("\n\r ----------- ") mprintf("\n\r wielkosc otrzymana w tescie : %g", Y); mprintf("\n\r wielkosc progowa : %g", Y_progowe); if (Y > Y_progowe) then mprintf("\n\r Wniosek : Hipoteza odrzucona !") else mprintf("\n\r Wniosek : Hipoteza nie odrzucona !") end // 3/ ilustracja graficzna function [d] = gestosc_chi2(X,N) d = X.^(N/2 - 1).*exp(-X/2)/(2^(N/2)*gamma(N/2)) endfunction czestosc = occ/m; ymax = max([czestosc ; p]); rect = [0 0 1 ymax*1.05]; xbasc(); xset("font",2,2) subplot(2,1,1) plot2d3((1:N+1) /(N+2), czestosc, style=2 , ... frameflag=1, rect=rect, nax=[0 N+2 0 1]) plot2d((1:N+1) /(N+2),p*ones(N+1,1), style=-2, strf="000") xtitle("Czestosci empiryczne (kreski pionowe) i prawdopodobienstwa teoretyczne (krzy\.zyki)") subplot(2,1,2) // rysujemy g\k{e}sto\ s\ c chi2 ???a N ddl X = linspace(0,1.1*max([Y_seuil Y]),50) ; 80 D = densite_chi2(X,N); plot2d(X,D, style=1, leg="densite chi2", axesflag=2) // kreski pionowe dla Y plot2d3(Y, densite_chi2(Y,N), style=2, strf="000") xstring(Y,-0.01, "Y") // kreski pionowe dla Yseuil plot2d3(Y_seuil, densite_chi2(Y_seuil,N), style=5, strf="000") xstring(Y_seuil,-0.01, "Yseuil") xtitle("Pozycja Y w odniesieniu do warto\ sci Yseuil") Poni|ej przedstawiono rezultaty otrzymane dla N = 10 et m = 5000 (patrz wykres (24)) : E dok\l{}adne = 0.5 E oszacowane = 0.4959167 Przedzia\l{} ufno\ sci dla 95% : [0.4884901,0.5033433] Test du chi 2 : -------------- warto\ s\ c otrzymana w te\ scie : 8.3176 warto\ s\ c progowa : 18.307038 Wniosek : Hipoteza nie odrzucona ! Frequences empiriques (traits verticaux) et probabilites exactes (croix) 0.1 × × × × × × × × × × × 0.0 0.000 0.083 0.167 0.250 0.333 0.417 0.500 0.583 0.667 0.750 0.833 0.917 1.000 Positionnement de Y par rapport a la valeur Yseuil densite chi2 Y Yseuil Rysunek 24: Ilustracja dla testu Ç2 urny Polya... fig:polya Most ?brownien? *Funkcja function [X] = most_brownien(t,n) X = sum(bool2s(grand(n,1,"def") <= t) - t)/sqrt(n) endfunction 81 *Skrypt Skrypt ten wykonuje m symulacji zmiennej losowej Xn(t), przedstawia wykres funkcji rozkBadu empirycznego, na- kBadajc na niego wykres funkcji rozkBadu oczekiwanego (dla n = +") i wreszcie oblicza test KS : t = 0.3; sigma = sqrt(t*(1-t)); // oczekiwane odchylenie standardowe n = 4000; // n "du\.ze" m = 2000; // licczba symulacji X = zeros(m,1); // zainicjowanie wektora realizacji for k=1:m X(k) = most_brownien(t,n); // p\k{e}tla licz\k{a}ca kolejne realizacje end // wykres funkcji rozk\l{}adu empirycznego X = - sort(-X); // tri prcum = (1:m) /m; xbasc() plot2d2(X, prcum) x = linspace(min(X),max(X),60) ; // odci\k{e}te i [P,Q]=cdfnor("PQ",x,0*ones(x),sigma*ones(x)); // rz\k{e}dne dla funkcji dok\l{}adenj plot2d(x,P,style=2, strf="000") // rzutujemy je na wykres pocz\k{a}tkowy // zastosowanie testu KS alpha = 0.05; FX = cdfnor("PQ",X,0*ones(X),sigma*ones(X)); Dplus = max( (1:m) /m - FX ); Dminus = max( FX - (0:m-1) /m ); Km = sqrt(m)*max([Dplus ; Dminus]); K_progowe = sqrt(0.5*log(2/alpha)); // prezentacja wynik\ ow // mprintf("\n\r Test KS : ") mprintf("\n\r -------- ") mprintf("\n\r warto\ s\ c otrzymana w te\ scie : %g", Km); mprintf("\n\r warto\ s\ c progowa : %g",K_seuil); if (Km > K_progowe) then mprintf("\n\r Wniosek : Hipoteza odrzycona !") else mprintf("\n\r Wniosek : Hipoteza nie odrzucona !") end Oto rezultaty otrzymane dla n = 1000 i m = 4000 : Test KS : -------- warto\ s\ c otrzymana w te\ sci : 1.1204036 warto\ s\ c progowa : 1.2212382 Wniosek : Hipoteza nie odrzucona ! 82 Laboratorium Eksperyment¶ Matematycznych L E M ow EX1 Badanie szeregu harmonicznego _ Antoni Zochowski 1. Celem ¶ _ lasno¶ cwiczenia jest przyblizenie uczestnikom wà sci szeregu harmonicznego oraz szeregu Dirichleta, ktore wykorzystywane sa, czesto w kryterium por¶ ownawczym badania , zbiezno¶ Przypominamy, ze _ sci. _ n X 1 S(n) = k k=1 to ciag sum cze¶ sciowych szeregu harmonicznego, natomiast , , 1 X 1 D(x) = ; x > 1 kx k=1 to suma zbieznego szeregu Dirichleta. _ Zbadamy asymptotyczne zachowanie sie S(n) oraz przebieg funkcji D(x). , 2 . Wykonaj czynno¶ wstepne oraz komende Scilaba diary('nazwisko.ex1') . sci , , 3 . Utw¶ wektor (wykorzystujac poznane w ex0 komendy Scilaba) orz , s =[1; 2; :::::::; 10000] 4 . Utw¶ wektor zawierajacy odwrotno¶ orz sci , 1 1 sz =[1; ; :::::::; ] 2 10000 5 . Wprowadzimy teraz komende sum(wek) , ktora daje jako wynik sume skladowych ¶ à , , wektora lub element¶ macierzy. Wypr¶ ow obuj --> wek=[1,7,9,11] --> suma1=sum(wek) --> suma2=sum( wek(1:3) ) --> suma3=sum( wek(3:4) ) 6 . Utw¶ teraz skrypt harmonic.sci , do ktorego zapisz komendy wykonane w 3,4. orz ¶ Pamietaj o ¶ sli c sredniego. srednikach, je¶ nie chcesz oglada¶ wyniku po¶ , , 7 . Dopisz do skryptu komendy, kt¶ utworza, ciag sum cze¶ ore sciowych , , SM(k) =S(100 ¢ k); k =1 : 100 oraz ciag warto¶ argumentu sci x(k) =100 ¢ k +1; k = 1 : 100 (plus 1 wkr¶ sie wyja¶ Wykonaj skrypt komenda, otce sni). , 1 --> exec('harmonic.sci',1) 8 . Dopisz do skryptu komende rysujaca, wykres SM wzgledem x oraz na ko¶ ncu , , , pause xbasc() Wykonaj skrypt. Po obejrzeniu rysunku wydaj komende resume . , 9 . Co Ci ten wykres przypomina ? Moze logarytm ? Dopisz do skryptu tworzenie _ ciagu , z(k) =log (x(k)); k = 1 : 100 oraz wykresu SM wzgledem z. Czy otrzymaà s prosta, ? le¶ , 10 . Dopasuj metoda, najmniejszych kwadrat¶ wsp¶lczynniki w aproksymacji liniowej ow oà SM(k) =a ¢ z(k) +b Uwaga: Wyprowadzenie wzor¶ na a i b. Niech dane sa, ciagi yi; xi; i = 1; :::; n. ow , Minimalizujemy funkcje , n X F (a; b) = (axi + b ¡ yi)2: i=1 Stad, z warunku @F=@a = @F=@b = 0, przy oznaczeniach , X X X X X2 = x2; XY = xi ¢ yi; X = xi; Y = yi i dostajemy ukà r¶ n lad owna¶ a ¢ X2+b ¢ X = XY a ¢ X + b ¢ n = Y oraz, je¶ ¢ = (n ¢ X2 ¡ X ¢ X), sli a =(n ¢ XY ¡ X ¢ Y )=¢ b =(X2 ¢ Y ¡ XY ¢ X)=¢ 11 . Czy otrzymaà s a ¼ 1 ? Stala, b nazywamy staÃa, Eulera (obliczyà s jej przyblizenie) le¶ à l le¶ n X 1 ¼ log (n +1) +°: k k=1 Widzisz teraz skad sie wzieà (n +1). lo , , , 12 . Dopisz komendy tworzace wykres r¶znicy o_ , SM(k) ¡ a ¢ z(k) ¡ b Czy wida¶ zbiezno¶c ? c _ s¶ 13 . Utw¶ wzorujac sie na harmonic.sci , skrypt dirichlet.sci tablicujacy i orz, , , , rysujacy wykres , 10000 X 1 D(x) = x =2 : 0:2 : 5 kx k=1 2 Zadanie o robaczku (R.L. Graham, D.E.Knuth, O Patashnik, Matematyka konkretna, PWN 2002) Powolny, ale wytrwaBy robaczek R rozpoczyna wdrówk od pocztku metrowej gumowej ta[my i peBznie - z prdko[ci 1 centymetr na minut - w kierunku jej koDca. Pod koniec ka|dej minuty równie wytrwaBy wBa[ciciel ta[my W, którego wyBcznym celem w |yciu jest pokrzy|owanie planów R, rozciga ta[m równomiernie o 1 metr. Po jednej minucie R jest odlegBy o 1 cm od pocztku swojej drogi i o 99 cm od jej koDca. Wówczas W rozciga ta[m o jeden metr, ale tak, |e R zachowuje swoja wzgldn pozycj, czyli 1% od startu i 99% od finiszu. Zatem R znajduje si 2 cm od punktu wyj[cia i 198 od celu. Po nastpnej minucie R osiga 3 cm ta[my, ma przed sob 197 cm do przej[cia i w tym momencie W rozciga ta[m. Odpowiednie odlegBo[ci przyjmuj warto[ci 4,5 i 295,5 cm. I tak dalej.... Czy robaczek kiedykolwiek dobrnie do koDca ta[my?

Wyszukiwarka