Przeskocz do treści

Delta mi!

Loading

Jak to działa?

Maszyna różnicowa

Piotr Chrząstowski

o artykule ...

  • Publikacja w Delcie: luty 2017
  • Publikacja elektroniczna: 31 stycznia 2017
  • Wersja do druku [application/pdf]: (112 KB)
obrazek

Dlaczego w szkole tak dużo uczymy się o wielomianach? Są dwa podstawowe powody. Pierwszy z nich - całkiem zrozumiały - po prostu jest to niemal największa klasa funkcji, których wartości umiemy obliczać. Potrafimy jeszcze dzielić wartości wielomianów, ale z pozostałymi funkcjami, które występują w programie szkolnym, a później na studiach, w zasadzie mielibyśmy sporo problemów.

Co by było, gdybyśmy, porwani przez jakichś wrogich, a potężnych Marsjan, zostali postawieni w takiej sytuacji: ludzkość zostanie ocalona, jeśli wykażemy się odpowiednią inteligencją matematyczną i, dysponując jedynie kartką papieru i ołówkiem, obliczymy w ciągu jednego dnia wartość  ○ |sin 1 z dokładnością do 10 cyfr po przecinku? Kto z absolwentów liceów poradziłby sobie z takim wyzwaniem? Chyba niewielu Czytelników potrafiłoby aż z taką dokładnością wyznaczyć poprawną wartość: |0,0174524064. Na pewno nie za pomocą cyrkla, linijki i ołówka. Zresztą podobnie trudno by było, gdyby ci podli Marsjanie kazali nam obliczyć np. log210 czy |2π z tą samą dokładnością. Powstaje pytanie: jak sobie radzili z tymi problemami nasi poprzednicy, którzy wieki temu drukowali tabele różnych funkcji matematycznych: pierwiastków, sinusów i logarytmów z zaskakująco dobrymi dokładnościami, sięgającymi nieraz 8 cyfr znaczących?

Dochodzimy tu do drugiej niezwykłej własności wielomianów. Otóż, o czym dowiadują się studenci pierwszego roku kierunków ścisłych, wielomiany mają jeszcze jedną sympatyczną własność. Za ich pomocą możemy dowolnie dokładnie przybliżać wszystkie funkcje dostatecznie gładkie, czyli takie, których wykres jest linią ciągłą bez żadnych załamań (konkretnie funkcje klasy |C czyli nieskończenie wiele razy różniczkowalne), a do nich należą właśnie funkcje trygonometryczne, logarytmy, funkcje wykładnicze. Konkretnie: w każdym skończonym przedziale możemy tak dobrać współczynniki wielomianu, żeby wartości w tym przedziale różniły się od wartości interesującej nas funkcji o mniej, niż z góry zadana wartość. Dla funkcji sinus i dla przedziału [0,π /4] wystarczy wziąć wielomian

P(x) = x− 1x3+ -1-x5 −--1--x7+ ---1---x9− ---1----x11+ -----1-----x13 6 120 5040 362880 39916800 6227020800

(w tym wzorze w mianownikach mamy silnie wykładników, a x jest podane w mierze łukowej), aby otrzymać wartości sinx z dokładnością do 10 cyfr dziesiętnych po przecinku dla każdego x z tego przedziału. Znając ten wzór i dysponując paroma godzinami, na pewno byśmy sobie z marsjańskim problemem poradzili.

No to do roboty! Jak jednak wyznaczyć wiele wartości wielomianu tak, aby się jak najmniej napracować? Postawmy się w sytuacji autora tabel matematycznych dwa wieki temu - nie mamy komputerów, chcemy osiągnąć dużą dokładność i obliczyć wartości takiej funkcji, jak sinus, powiedzmy, co |-1-- 10000 radiana. Ile trzeba będzie wykonać mnożeń? One w końcu sprawiają więcej kłopotu niż dodawanie i odejmowanie. Powiedzmy, że zajmiemy się wielomianem stopnia |n, postaci P (x) = a0 +a1x + a2x2 + ...+ anxn. Na pierwszy rzut oka wygląda na to, że dla każdej wartości |x trzeba obliczyć wszystkie potęgi x od |2 do n, co kosztuje n − 1 mnożeń, potem otrzymane wartości  k x przemnożyć przez współczynniki ak dla |k = 1,...,n (kolejne |n mnożeń) i uzyskane liczby zsumować. Razem mamy więc 2n − 1 mnożeń i n dodawań.

Dość łatwo można ten wynik poprawić, stosując schemat Hornera. Jeśli przedstawimy nasz wielomian w postaci

P(x) = a + x(a + x(a + x(...+ x(a + xa )...))), 0 1 2 n− 1 n

wówczas - o dziwo! - liczba mnożeń spadnie nam do n, a liczba dodawań się nie zmieni. W tej metodzie obliczenia zaczynamy od wewnątrz: mnożymy x przez |an, dodajemy an−1, mnożymy przez x, dodajemy |an−2 itd., aż dodamy na końcu a0. Wydaje się, że tego wyniku poprawić już nie można.

A jednak! Przedstawimy teraz sposób, za pomocą którego będziemy obliczali wartości wielomianu w kolejnych liczbach, nie wykonując, poza kilkoma początkowymi wartościami, żadnego mnożenia! Metodę zilustrujemy na przykładzie konkretnego wielomianu  3 2 V (x) = 2x − x − 4x + 7. Będziemy wyznaczali wartości tego wielomianu dla kolejnych liczb naturalnych. Pierwsze cztery wartości to, dla x = 0,1,2,3, odpowiednio 7,4,11,40, co obliczamy w pamięci lub schematem Hornera.

Teraz pora na magię. Obliczmy różnice między tymi wartościami: będą to kolejno liczby |−3,7,29. Powtórzmy to jeszcze raz, tym razem dla tych trzech wartości: |10,22. Na koniec jeszcze raz obliczmy różnicę między tymi dwiema wartościami: 12. Zapamiętajmy tę wartość, będziemy z niej wielokrotnie korzystać. Jesteśmy już przygotowani do wyznaczenia wartości V (4) za pomocą 3 dodawań. Skupiamy się na ostatnich wartościach z przedstawionych tutaj ciągów, zapisanych w odwrotnej kolejności, czyli |12,22,29,40. Dalsze trzy wartości otrzymamy, sumując je w następujący sposób: 12+ 22 = 34, 34+ 29 = 63,63+ 40 = 103. Ostatnia otrzymana wartość to właśnie V (4). Teraz, żeby otrzymać |V(5), powtarzamy nasz algorytm dla właśnie wygenerowanych liczb, zaczynając od tej samej dwunastki, co poprzednio i dostajemy kolejno |12 + 34 = 46, 46+ 63 = 109, 109+ 103 = 212, więc V (5) = 212. Biorąc znowu dwunastkę i sumując ostatnio otrzymane wyniki, otrzymamy kolejno wartości 58, 167,379, z których ostatnia liczba jest, jak się już domyślamy, równa |V(6). Następnie tą samą metodą wyznaczamy po kolei |V(7),V (8),V (9)..., za każdym razem wykonując tylko trzy dodawania i żadnego mnożenia!

Co się tu dzieje? Wszystko dzięki niepozornemu operatorowi różnicy skończonej, który oznaczymy przez |∆. Ten operator działa na wielomianach i tworzy z nich nowe wielomiany według następującego wzoru: |∆(P(x)) = P(x + 1)− P(x) dla dowolnego wielomianu P (x). Czyli dla każdego argumentu x wartość ∆ (P(x)) jest równa przyrostowi wartości wielomianu, gdy argument zwiększymy o 1. Obliczmy zatem dla naszego wielomianu

∆(V (x)) = 2(x+1)3 −(x+1)2−4(x+1)+7 − 2x3+x2+4x −7 = 6x2+4x −3.

Okazuje się, że otrzymaliśmy wielomian stopnia o jeden niższego niż oryginał. Chwila zastanowienia i widzimy, że tak będzie dla każdego wielomianu trzeciego stopnia. Operator ∆, działając na taki wielomian, zmniejszy jego stopień do kwadratowego. I ogólnie, jeśli wielomian P(x) jest stopnia |n, to ∆ (P(x)) będzie wielomianem stopnia n − 1. Można to sprawdzić, rozwijając P (x + 1) ze wzoru dwumianowego Newtona; współczynnik przy najwyższej potędze to an, więc wyraz z xn skróci się przy odjęciu |P(x). Jeszcze krótki namysł i zobaczymy, że współczynnikiem przy najwyższej potędze (czyli przy  n−1 x ) wielomianu |∆(P (x)) będzie |n⋅an.

Wprowadźmy teraz potęgi operatora |∆ ze względu na składanie. Ustalmy, że |∆0 jest operatorem identycznościowym, czyli z definicji |∆0(P(x)) = P(x). Przyjmijmy, że |∆n+1(P(x)) = ∆(∆ n(P(x))) dla |n⩾ 0. Co się stanie, jeśli będziemy iterować operator ∆ dla wielomianu |P(x)? Pierwsza iteracja - to już wiemy - będzie wielomianem stopnia |n− 1 o współczynniku przy najwyższej potędze równym nan. Druga - wielomianem stopnia |n− 2 z najwyższym współczynnikiem równym |n(n −1)a , n trzecia da nam wielomian stopnia n −3 o najwyższym współczynniku równym n(n − 1)(n −2)an i tak dalej, aż w końcu |n-ta iteracja, czyli |∆n(P (x)) będzie wielomianem stopnia zerowego o współczynniku przy najwyższej potędze (czyli wyrazie stałym) równym |n(n −1) ⋅... ⋅2 ⋅1⋅a = n!a . n n Zatem postać wielomianu ∆ n(P (x)) zależy jedynie od stopnia tego wielomianu i od współczynnika an. Jest więc dość łatwa do wyliczenia. Przykładowo, dla naszego wielomianu V (x) mamy

  •  0 3 2 ∆ (V(x)) = 2x −x − 4x +7,
  • ∆ 1(V (x)) = 6x2 + 4x − 3,
  •  2 ∆ (V(x)) = 12x +10,
  • ∆ 3(V(x)) = 12 (co było do przewidzenia od samego początku, bo 12 = 3!⋅2 )

i kolejne iteracje ∆ dadzą już tylko wielomiany zerowe. Teraz pora na najważniejsze. Widzimy, że bezpośrednio z definicji operatora |∆ wynika wzór

P (x + 1) = P (x)+ ∆ (P(x)),

co zapiszemy w dogodniejszej dla dalszych rozważań postaci

 1 P (x) = P (x −1) +∆ (P (x −1)).

Ale |∆1(P(x − 1)) to też wielomian zmiennej x, więc proces możemy iterować, bo w szczególności  1 1 2 ∆ (P (x− 1)) =∆ (P (x− 2)) + ∆ (P(x − 2)). Idąc tak dalej, otrzymujemy

pict

a ostatni wyraz, jak wiemy, równy jest |n!⋅an. Widać więc, że do wyznaczenia wartości P (x) wystarczy do wartości P (x −1) dodać wartości wszystkich |i-tych potęg |∆ dla argumentów x −i, gdzie |i przebiega zbiór {1,...,n}. A te wartości możemy wyznaczać kolejno, wykonując n dodawań.

Stwórzmy więc dla naszego wielomianu V (x) tabelę takich wartości. Zaczniemy od obliczenia "na piechotę" wartości V (0),V (1),V(2),V (3).

----------------------------------------------- |x |V (x) |∆1(V (x)) |∆ 2(V(x)) |∆ 3(V (x)) | |--|-------|----------|-----------|----------|- |0-|--7----|----------|-----------|----------|- |1-|--4----|----------|-----------|----------|- |2 | 11 | | | | |--|-------|----------|-----------|----------|- -3----40---------------------------------------

Obliczmy teraz wartości ∆ (V (x)) dla |x = 0,1,2.

|--|-------|-1--------|--2--------|--3-------|- |x-|V-(x)--|∆-(V-(x))-|∆--(V(x))--|∆-(V-(x))-|- |0-|--7----|---−-3----|-----------|----------|- |1 | 4 | 7 | | | |--|-------|----------|-----------|----------|- |2-|--11---|---29-----|-----------|----------|- -3----40---------------------------------------

Następnie wyznaczmy wartości ∆2(x) dla x = 0,1 ; więcej nam nie potrzeba. Od razu wstawmy dla porządku 12 w miejsce ∆ 3(V(x)) dla x = 0 - wiemy, że będzie to | 12 = 2 ⋅3!, ale gdybyśmy nie dowierzali, zawsze możemy tę wartość uzyskać z ostatniej wypełnionej kolumny, odejmując 10 od 22.

|--|-------|----------|-----------|----------|- |x-|V-(x)--|∆1(V-(x))-|∆-2(V(x))--|∆-3(V-(x))-|- |0 | 7 | − 3 | 10 | 12 | |--|-------|----------|-----------|----------|- |1-|--4----|----7-----|----22-----|----------|- |2-|--11---|---29-----|-----------|----------|- |3 | 40 | | | | -----------------------------------------------

W powyższej tabeli została wytłuszczona przekątna, za pomocą której będziemy generowali kolejną przekątną. Zaczniemy od domyślnej dwunastki w ostatniej kolumnie (cała kolumna składa się z samych dwunastek), którą wstawimy dla |x = 1. Teraz, poczynając od góry, będziemy do każdego z wytłuszczonych elementów przekątnej dodawali jego prawego sąsiada i wynik zapisywali pod spodem, wytłuszczając go:

|--|-------|----------|-----------|----------|- -x--V-(x)---∆1(V-(x))--∆-2(V(x))---∆-3(V-(x))--- |0 | 7 | − 3 | 10 | 12 | |--|-------|----------|-----------|----------|- |1-|--4----|----7-----|----22-----|----12----|- |2-|--11---|---29-----|----34-----|----------|- |3 | 40 | 63 | | | |--|-------|----------|-----------|----------|- -4---103---------------------------------------

Proszę! Wykonaliśmy trzy dodawania i mamy wartość |V(4). Jedziemy dalej:

|x-|V-(x)--|∆1(V-(x))-|∆-2(V(x))--|∆-3(V-(x))-|- |--|-------|----------|-----------|----------|- |0-|--7----|---−-3----|----10-----|----12----|- |1-|--4----|----7-----|----22-----|----12----|- |2 | 11 | 29 | 34 | 12 | |--|-------|----------|-----------|----------|- |3-|--40---|---63-----|----46-----|----------|- |4-|-103---|---109----|-----------|----------|- |5 | 212 | | | | -----------------------------------------------

I jeszcze V (6), skoro nam tak dobrze idzie:

|--|-------|-1--------|--2--------|--3-------|- |x-|V-(x)--|∆-(V-(x))-|∆--(V(x))--|∆-(V-(x))-|- |0-|--7----|---−-3----|----10-----|----12----|- -1----4---------7----------22----------12------ |2 | 11 | 29 | 34 | 12 | |--|-------|----------|-----------|----------|- |3-|--40---|---63-----|----46-----|----12----|- |4-|-103---|---109----|----58-----|----------|- |5 | 212 | 167 | | | |--|-------|----------|-----------|----------|- -6---379---------------------------------------

Działa! Trzy dodawania i mamy obliczoną kolejną wartość wielomianu trzeciego stopnia.

No dobra. To było proste, ale mieliśmy sytuację komfortową, kiedy to wartości wielomianu obliczaliśmy dla kolejnych liczb naturalnych. Co by jednak było, gdyby argumenty były wymierne? Otóż nic by się nie zmieniło. Okazuje się, że jeśli operator |∆ zastąpimy operatorem ∆ ε dla dowolnego wymiernego dodatniego ε i zdefiniujemy ∆ (P(x)) = P(x + ε)− P (x), ε to zachowamy podstawową cechę operatora |∆: zastosowanie tego operatora zmniejsza stopień wielomianu o 1. To nam wystarcza do tego, żeby powtórzyć naszą procedurę dla dowolnego początkowego x i dla dowolnego |ε> 0. Obliczenie kolejnych |M wartości dla argumentów różniących się o z góry ustaloną stałą dla dowolnego wielomianu stopnia n wymaga zatem obliczenia jego |n kolejnych wartości "na piechotę", a następnie zastosowania −n M razy naszego schematu, za każdym razem wymagającego jedynie |n dodawań, bez wykonywania jakiegokolwiek mnożenia! Niewiarygodne!

Na to, żeby metodę opisaną powyżej wykorzystać do budowy maszyny, przypuszczalnie wpadł jako pierwszy inżynier heskiej armii Johann Helfrich von Müller (1742-1830). Pomysł ten rozwinął i częściowo zrealizował angielski matematyk Charles Babbage (1791-1871). Jest to bohater niezwykły, wizjoner, twórca modelu komputera w postaci, której praktycznie dziś powszechnie używamy. Żył długo, choć jego życie nie było usłane różami. Dość młodo został profesorem matematyki w Cambridge, głównie za prace w dziedzinie kryptografii - ciągnęło go zawsze w kierunku zastosowań matematyki w informatyce, choć informatyki wtedy jeszcze nie było. W latach 20. XIX wieku wpadł na pomysł skonstruowania maszyny, która wykonywałaby automatycznie działania prowadzące do wyznaczania kolejnych wartości wielomianów według podanej wyżej metody. Sam mechanizm arytmometru był znany od prawie 200 lat dzięki konstrukcjom Schickarda, Pascala i Leibniza. Chodziło o to, żeby ustawić kilka arytmometrów tak, aby mogły sobie nawzajem przekazywać obliczone sumy w odpowiedniej kolejności. Efektem prac Babbage'a była maszyna różnicowa, której model przedstawiony jest na zdjęciu powyżej. Wymagała ona ręcznego wykonania początkowych czynności. Babbage na tym nie poprzestał. Zaczął budowę kolejnej wersji maszyny różnicowej, która byłaby znacznie większa i... nigdy jej nie ukończył. Nie dlatego, żeby napotkał jakieś nieoczekiwane przeszkody. Po prostu w trakcie prac nad nią wpadł na pomysł, który tak go pochłonął, że zarzucił zaczęte prace. Ale jest to temat na osobną historię.