Występow. obl. numerycznych. 1) Błędy ich klasyfikacja i najprostsze spos. ich obliczania. Związki ilościowe wyst. w zjawiskach przyrodn. są opisywane na granice matematyczne za pom. równań wiążących pewne funkc., których część jest znana a część natomiast wymaga obliczenia. Równanie tego rodzaju można bez istotnego ograniczenia ogólności zapis. w postaci: U=Af (1) f-elem. pewnej przestrzeni funkcyjnej; A-operator czyli operacją określ. spos. przekszt. elem. A w elem. U, w szczególności A nazywa się funkcjonalem, J. U jest liczbą rzeczywistą. Zad., które trzeba rozwiązać, sprow. się do obl. z równ. (1) elem. U gdy dany jest elem. f. 1.Udowodnienie istn. i jednozn. rozw. równ. 2.Wyrażanie stabiln. zad. lub co na jedno wychodzi stabiln. rozw tego zad. Oznacza to że trzeba udowodnić iż dostatecznie małym zmianom danych wyjść f odpowiada dowolnie małej zmianie rozw. w rozpatrywanym obszarze. 3.Znalezienie rozw. w post. wzoru. Każde zagadnienie dla którego potrafimy pozytywnie rozw. problem 1 i 2 nazywa się zagadn. popr. postawionym. Poszukiw. rozw. U równ. (1) w postaci wzoru w wielu konkretnych przykładach związane jest z zasadniczymi trudnościami. Stąd prowadzi konieczność uciekania się do pomocy tzw. metod przybliżonych rozw. zad. Prowadzi to do otrzymania zamiast szukanej f. i liczb wartości przybliżonych. Wartości f. będziemy nazywać przybliżeniami f., a wartości przybliżone liczb liczbami przybliżonymi. Z przybliżeniami f. mamy do czynienia, np. wówczas gdy w celu uproszczenia obliczeń zastępujemy daną f. sumą skończoną ilości wyrazów jej szereg Fouriera lub szereg Taylora. Liczby przybliżone spotykamy wówczas, np. przy mierzeniu różnych wielkości fiz.: długości, powierzchni, objętości, masy, prędkości, itd. Liczby przybliżone zjawiają się również jako wynik działań artmet. na liczbach przybliżonych oraz na skutek odrzucania dost. małych wielkości podcz. obl. W związku z tym powst. wiele zagadnień: Jak wybrać przybliżoną metodę rozwiązywania równ typu (1)? Jaki będzie przy tym błąd otrzymanego rozw. przybliż. tzn. w jakim stopniu będzie się ono różnić od rozw. dokł.? Jak scharskt. błąd liczb przybl. miarę ich odchylenia od odpow. liczb dokł.? Jak oszacować błąd wyniku obl.? Jaki jest dop. bł. l. przybl., przy którym bł. wyniku działań artmet. na tych liczbach nie jest zbyt duży? W jaki sposób można wykonywać działania na l. przyb. żeby błąd wyniku był najmniejszy? Badanie wszystkich tych zagadnień jest przedmiotem metod num. Bł. bezwzględny i dokładna wartość bł. bezwzględnego. J. a jest wart. przybliżoną liczby a₀ to w skrócie zapiszemy a₀ ≈ a (a₀ jest przybliżalnie równa a). Dokładną wartością błędu liczby a lub przybliżonej równości (2) nazywamy wielkość ∆=a₀-a (3). Błędem bezwzględnym przybliżonej liczby a nazywamy każdą możliwie małą liczbę dodatnią α, taką że |∆|=|a₀-a|≤α (4) Np. niech będzie dana liczba a₀=0,5531. Gdy weźmiemy jako jej przybliżenie a=0,55 to dokładna wart. bł. wyn. ∆=0,0031. W pewnym przyp. dla wsk. bł. bezwzgl. równości przybliżonej (2) przyjmujemy zapis (a₀=a±α) (2'). J. spełn. jest nierówn. (4) to mówimy że przybliżona l a przedstawia l. a₀ z dokładn. do wielkości α. Bł. względny i dokł. wart. bł. bezwzględnego. Łatwo zauważyć że bł. bezwzgl. nie charakt. w pełni pomiarów. J. bł. bezwzgl. przy pomiarze dł. pudełka z zapałkami i przy pom. dł. biurka w obu przyp. jest bezporówn. bardziej dokł. niż pierwszy mimo że bł. bezwzgl. obu pomiarów jest taki sam. Naturalne jest wprowadzenie jeszcze jednej wielkości charakt. liczby przybliżone. Tą wielk. jest tzw. błędem przybliżonym liczby przybliżonej. Dokładną wielk. bł. wzgl. jest wielk. η=∆/a=(a₀-a)/a (5). Bł. wzgl. l. przybl. a lub przybl. równ. (2) nazywa się dowolną możliwie małą dodatnią liczbę δ taką że: |η|=∆/a=|(a₀ - a)/a|≤δ (6) Z określenia tego wynika że J. znany jest bł. bezwzgl. liczby przybliż. a to jako bł. wzgl. można przyjąć wielk.: δ=α/|a| (7). Bł. wzgl. zawsze wyraża się liczbą niemianowaną. W praktyce często bł. wzgl. wyrażamy w % wielkości zmiennej. Bł. wzgl. charakt. wyst. dokładnie l. przybl., a w szczególn. jakość pomiaru różnych wielkości. Przyjmiemy np. że przy średn. drutu i wys. drzwi możemy dopuścić bł. nie przekr. 1mm. otrzymując w rezult. pom. średn. drutu ≈5mm a wys. drzwi ≈2m łatwo stwierdzić że bł. wzgl. w 1 przyp. wyn. 20% a w 2 zaś 0,05%. Dochodzimy zatem do wn. że im mniejszy jest bł. wzgl., tym lepsza jest jakość pomiaru. Dokładanie cyfry i zapis przybliżony. Każda l. rzeczyw. a ma rozw. przy podst. q=2,3,8,10,16 tzn. można przedst. ją w spos.: a=±(a₁qⁿ+a₂q^n-1+...+a_mq^n-m+1+...) (8) gdzie: a₁-liczby całkowite (0≤a_i≤q). Daną liczbę a zapiszemy przy tym w postaci ciągu cyfr, im bardziej na prawo tym cyfra niższa i niższa jest jej pozycja tzn. mniejszy jest wykładnik potęgi, przy której cyfra wyst. we wzorze (8). Cyfry wyst. przy potęgach o nieujemnych wykładnikach są oddziel. „,” od cyfr niższych pozycji, np. l. 11,1 zapisana w postaci dwójkowego syst. liczb. jest równa 2¹+2⁰+2^-1=3½. Wszystkie cyfry liczby rozwinięcia przy pewnej podst., poczynając od pierwszej cyfry niezerowej nazywamy cyframi znaczącymi, np. l. 0,016 ma dwie cyfry znaczące 1 i 6. Ilość cyfr znacz. liczb można zwiększać lub zmn. przez dopisyw. zer na ostatnich pozycjach. Dana liczba nie zmienia się przy tym. dy/dx=lim_[∆x;0]∆y/∆x ; y(x₀+∆x)-y(x₀)=∆y ; dy/dx≈∆y/∆x=[y(x₀+∆x)-y(x₀)]/∆x ; ∫_[a;b]f(x)dx ; S_tr=½h[f(a)+f(b)]=½(b-a)[f(b)+f(a)]x=x₀ ; x_i-x_i-1=h_i ; ∆S_i=h/2[f(x_i)+f(x_i-1)] ; S_a0=∑∆S_i=h/2∑_[i=1,n][f(x_i)+f(x_i-1)] RYS. 1 W celu, żeby z ostatniej cyfry znaczącej można było sądzić o wartości błędu bezwzględnego obieramy parametr liczbowy Ω( ½≤Ω≤1) i umawiamy się, że cyfrę a_m liczby (8) nazywać będziemy dokładną, j. dla błędu bezwzględnego α tej liczby zachodzi nierówność α≤Ωq^n-m+1 (9) J. ten warunek nie jest spełniony to cyfrę a_m nazywamy wątpliwą. J. cyfra a_m jest dokładna, to również i wszystkie poprzednie cyfry są dokładne. W taki spos. pom. dokładnymi cyframi można zawsze znalęźć ostatnią (J. α>0) Ostatnia cyfra znacząca powinna być dokładna. Osiągamy to przez odrzucenie cyfr wątpliwych i dopisanie w razie konieczności po prawej stronie mnożnika q^k (kЄC). Np. w dziesiętnym syst. liczenia zapis a=3,14 przy α≤Ω*10^-2 jest poprawny zapis a=3,140 nie jest poprawny; zapis a=1200 przy α≤Ω≤1 - poprawny a przy α≤Ω*10 niepoprawny. Liczbę a=124237 z trzema cyframi dokładnymi trzeba zapisać w postaci a=124*10³, a liczbę a=0,00412 z dwiema dokładnymi cyframi 4 i 1 można zapisać poprawnie jak a=4,1*10^-3. Zaokrąglanie liczby przybliżonej (8) do ostatniej cyfry a_m polega na tym, że J. Θ=a_m+1q^n-m+a_m+2q^n-m+1+...<½ q^n-m+1 (10) to wartość Θ odrzuca się. J. jednak Θ ≥ ½ q^n-m+1, to odrzuca się Θ, aleprzy tym zwiększa się o jedność cyfrę a_m ostatniej z pozostawionych pozycji. W niektórych przyp. dziesiętnego syst. licz. przy Θ=½ q^n+m-1 zostawia się a_m zmiany gdy a_m jest parzysta i zwiększa się o jedność, gdy a_m jest nieparzysta. Błąd nieunikniony. Wielkości dane wyst. w każdym konkretnym zagadnieniu typu (1) z regóły nie są znane dokładnie, jedynie w przybliżeniu, z pewnym błędem. Ze wzgl. na ten błąd danych wyjściowych w końcowym wyniku rozwiązania zad. otrzym. nieuchronnie pewien błąd nazwany błędem nieuniknionym. Rozpatrzymy zagadnienia dot. wielk. bł. nieunikn. prostszych f. i działań artmet. Niech y⁰-f(x₁⁰,x₂⁰,x₃⁰,...,x_n⁰) (11) będzie różniczkowalna w spos. ciągły f. zmiennych niezal. x₁⁰,...x_n⁰. Wówczas przy podaniu niedokładnej ale przybliżonej wartości argumentów x₁=x₁⁰±α₁; x₂=x₂⁰±α₂ ;...; x_n=x_n⁰±α_n zamiast y⁰ otrzymujemy z f. (11) pewną wielkość y=y⁰±α_y, gdzie α_y - błąd bezwzgl. wielkości y. Powstaje pytanie, w jaki spos. znaleźć wielkość α_y(α_y), a także odpowiedni błąd względny δ_y, jako wielkości charakt. błąd nieunikniony f. y. Stosując wzór Lagrange'a dostaniemy y-y⁰=∑_[i=1;n][∂f(x)/∂x_i](x_i-x_i⁰) (12), gdzie: y(x)=f(x)=f(x₁,...,x_n), x - punkt o współrz. x₁,x₂,...,x_n. Nawiasy kwadratowe oznaczają, że ujęte w nich pochodne są brane w pewnym wewn. punkcjie x+Θ(x-x⁰), Θ-liczba zawarta między zerem a jednością. Rozważania będziemy prowadzić przy następujących uproszczonych założeniach: a)pochodne cząstkowe f. f zmieniają się dost. wolno; b)błędy wzgl. δ₁,δ₂,...,δ_n danych wyjściowych x₁,x₂,...,x_n są dost. małe. Ze wzoru (12) mamy: a)α_y≈∑_[i=1;n]|∂f(x) / ∂x_i|α_i (13) Jest to ogólny wzór przybliżony dla błędu bezwzględnego f.. Dzieląc obie str. rown. (13) przez y=f(x) oraz uwzględniając założenia a) i b) otrzymamy wzór na błąd względny δ_y≈∑_[i=1;n] |∂lnf(x)/∂x_i|α_i (14) lub δ_y≈∑ |∂lnf(x)/∂lnx_i|δ_i=∑{|x_i||∂lnf(x)/∂x_i|}δ_i (14') Z wzorów (13), (14) i (14') otrymujemy wielkości błędów działań artmet. Zasady obl. liczby dokładnych cyfr. Widzieliśmy już (pkt 4), że pozycja ostatniej dokładnej cyfry znaczącej liczby dokładnej i jej błąd bezwzględny α wyrażają się jedno przez drugie. Niech w liczbie przybliżonej a=±(a₁qⁿ+a₂q^n-1+...+a_mq^n-m+1) (8) wszystkie cyfry przy danym wyborze parametru Ω (½≤Ω≤1) będą dokładne. Znaczy to, że α≤Ωq^n-m+1 (9) Dzielac bie str. (9) przez |a| mamy δ≤Ωq^n-m+1/|a₁qⁿ+...+a_mq^n-m+1|≤Ωq^n-m+1/a₁qⁿ czyli δ≤Ω/a₁qⁿ (15) to dana liczba a z pierwszą cyfrą znacząca a₁ ma co najmniej m dokładnych cyfr znaczących. Istotnie mamy α=|a|δ≤δ(a₁+1)qⁿ≤Ωq^n-m+1 a to znaczy, że cyfra a_m jest dokładna. Z nierówności (15) i (15') wynikają stosunk. proste zasady obliczania liczb dokładnych. Błąd metody Dokładne rozw. wielu zad. typu (1) jest związ. z podst. trudnościami. Stosując jakąkolwiek przybliżoną metodę rozw. zad. U=Af (1) rozwiązujemy faktycznie inne „aproksymujące” zadanie U^--=A^—f^--(1') gdzie: f^---element zbliżony w pewnym sensie do elementu f, A^--- operator zbl. w p. sens. do op. A. Wielkość U-Ū (16) będziemy nazywać błędem metody. W ogólnym przyp. nie ma nie ma możliwości wskazania wartości lub oszacowania błędu metody. Zagadnienie to powinno być badane w zastosowaniu do każdej konkretnej metody rozwiązyw. dotyczącej szerokiej czasem zadań, a czasem i w zastosowaniu do szczególnego zadania. Obliczany np. metodą prostokątów całkę RYS. 2 I=∫_[0;1]f(x)dx Biorąc punkty x_i=ih, h=1/n(n-l) całkowita, i=0,1,2,...,n otrzymamy w przybliżeniu I_i=f(x_i)(x_i+1)=f_ih Dokładną wartość I≈∑_[i=0;n-1]I_i=∑_[i=0;n-1]f(x_i) bł. bezwzgl. metody określa równość ∆=∫_[0;1]f(x)dx-∑_[i=0;n-1]f(x_i)h=∑_[i=0;n-1]∫_{[x[i];x[i+1]]}f(x)dx-∑_[i=0;n-1]f(x_i)h=∑_[i=0;n-1]∫_{[x[i];x[i+1]]}[f(x)-f(x_i)]dx≤∑_[i=0;n-1]∫_{[x[i];x[i+1]]}|f(x)-f(x_i)|dx. U podstaw podanych tu części znanych metod przybliżeniowych leży jedna ogólna idea, która polega na tym że zamiast szukania rozw. zad. (1) w danej przestrzeni funkcyjnej poszukujemy rozwiązania tego zad. w węższej skończenie wymiarowej przestrzeni f. Bł. numeryczny. Pojęcie stabilności obliczeń. Przypuśćmy że jakąkolwiek metodą przybliżoną doprowadziliśmy rozwiązanie zagadnienia (1) do rozwiązania pomocniczego zagadnienia aproksymującego (1'). Żeby rozw. to zad. należy przede wszystkim podać algorytm rachunkowy. Oznacza to że trzeba przedstawić zbiór praw i reguł obliczeniowych określających przebieg obliczeń, prowadzący do szukanego wyniku. Np. pożądane jest tak dobierać metodę, żeby ilość potrzebnych działań rachunkowych była wedle możności najmniejsza. Czasem trzeba poświęcić prosty schemat obliczeniowy na korzyść takiej metody, której liczba działań rachunkowych jest najmniejsza. Np. rozwiązanie ukł. algebraicznych równań liniowych ∑_[j=1;n]a_ijx_j=b_i (17)(i=1,2,...,n) z wyznacznikiem ∆≠0 można zapisać za pomocą wzorów Cramera x_i=∆_i/∆ (18)(i=1,2,...,n) gdzie ∆_i - wyznacznik otrzymany z ∆ przez zastąp. i-tej kolumny elem. kolymną prawej str. ukł. (17). Należałoby wtedy obliczyć n+1 wyznaczników n-tego stopnia. Przy obliczaniu każdego z nich trzeba wykonać (n-1)n! mnożeń. Ogólna liczba koniecznych działań mnożenia i dzielenia wyniesie (n²-1)n!+n. Skorzystajmy teraz z metody eliminacji Gaussa. Pierwsze z równań (17) daje x_i=b_i/a_ii-a₁₂/a₁₁x₂-...-a_1n/a₁₁x_n (19). Wykomuje się przy tym n dzieleń. Podstawienie x_i w każde z n-1 następnych równań wymaga n mnożeń. Ogólnie rzecz biorąc, do otrzymania n-1 równań o niewiadomych x₂,x₃,...,x_n potrzeba n² mnożeń i dzieleń razem, należy zatem wykonać n/6(2n²+3n+1) mnożeń i dzieleń razem, żeby sprowadzić ukł. (17) do ukł. macierzy trójkątnej. Do rozwiązania ukł. trzeba wżiąć [n(n-1)/2 mnożeń. W taki spos. do rozw. ukł. (17) potrzeba n/3(n²+3n-1) mnożeń i dzieleń razem. Na pierwszy rzut oka mogłoby się wydawać że pierwsza metoda (Cramera) jest lepsza; przeprowadzone jednak obliczenia liczby potrzebnych mnożeń i dzieleń dowodzą jednak czegoś przeciwnego. Np. dla 10 równań metoda Cramera wymaga (10²-1)*10!+10 ≈ 36*10⁷, metoda Gaussa natomiast tylko 10/3(10²+3*10-1)=430. W taki spos. dzięki dobrze dobranej metodzie obliczeń można często znacznie skrócić pracę. Istnieją jednak jeszcze inne powody o charakt. bardziej zasadniczym, które zmuszają do tego, żeby dążyć do max zmn. liczby koniecznych działań rachunk. Podczas obl. z reguły zmuszeni jesteśmy zaokrąglać liczby lub po prostu odrzucać wielkości dost. małe. W wyn. tego zjawia się tzw. bł. numeryczny lub bł. obliczeń. Przy obliczeniach wyk. na maszynach cyfrowych wyst. 3 podst. rodzaje bł.: wejściowe, obcięcia i zaokrągleń. Bł. wejściowe (względnie bł. danych wejściowych) wyst. wówczas gdy dane liczbowe wprowadzone do pamięci lub rejestrów maszyny cyfrowej odbiegają od dokł. wartości tych danych. Błędy takie pojawiają się wtedy gdy dane wejściowe są wynikiem pomiarów wielk. fiz., mierzonych z pwen. bł. pomiaru. Bardziej jednak charakt. przyczyną wyst. bł. wejściowych jest skończona dł. słów binarnych, reprezentujących liczby w maszynie i nieuniknione w związku z tym wstępne zaokrąglanie. Zauważmy na koniec że wstępna zaokrąglanie wyst. w przypadku wszystkich liczb niewym. √2 π e itp. Bł. obcięcia powst. na skutek zmniejszania liczby działań. Postępujemy tak na ogół przy obl. sum nieskończ. (szeregów). Np. wartości wyrażenia e^x równego szeregowi 1+x+x²/2+...+x^N/N! o odpowiednio dobranej wartości N. Podobnie zastąpienie wartości pochodnej f. jej ilorazem różnicowym powoduje powstanie bł. obcięcia. Praktyka obliczeń stawia teorii metod num. jeszcze jedno ogólne żądanie o pierwszorzędnym znaczeniu. Jest to żądanie stabilności obliczeń względem błędu. Okazuje się często że dla otrzymania wyniku konieczne jest przeprowadzenie długotrwałego ciągu obl. związanych z daną metodą; ciąg ten zazwyczaj jest tym dłuższy, im większą dokł. chcemy uzyskać. J. popełnimy bł. w pewnym działaniu tych obliczeń to powinniśmy się spodziewać że bł. ten będzie maił wpływ na wynik nast. działań rachunku. J. bł. popełniony w pierwszych krokach rachunku nie zwiększa się w nast. krokach lub pozostaje tego samego rzędu to przebieg obl. nazywa się stabilnym względem bł. początkow. J. bł. ten wzrasta to przebieg obl. jest niestabilny. Metody numeryczne algebry liniowej. 1.Wiadomości podst. z algebry. 1.Macierze i wyznaczniki. Tablicę prostokątną zawierającą m wierszy i n kolumn nazywamy (m,n)-macierzą. Przez a_ij oznaczamy elem. m. stojący na przecięciu i-tego wiersza i j-tej kolumny. J. w m. A zmienimy role kolumn i wierszy to otrzymamy m. przestawioną A' m. A. J. A'=A to m. A naz. się symetryczną. M. której wszystkie elem. są rzeczyw. naz. się m. rzeczywistą. M. której wszystkie elem. są 0 naz. się m. zerową. Dwie m. A=[a_ij] i B=[b_ij] naz. równymi gdy a_ij=b_ij przy wszystkich i, j. M. A i E naz. sprzężonymi gdy a_ij=ē_ij przy wszystkich i, j (kreska u góry oznacza wielkość zespoloną sprzężoną do e_ij). Suma dwóch macierzy A+B=C o elem. c_ij=a_ij+b_ij. Iloczynem m. A przez liczbę α naz. m. o elem. c_ij=αa_ij. Iloczynem dwóch m. A i B naz. m. C o elem. c_ij=∑_[k=1;n]a_ikb_kj gdzie n - liczba kol. A=z zał. l. wierszy m. B. Dowilny wektor o składowych x₁,x₂,...,x_n można przedstawić jako wiersz (x₁,x₂,...,x_n) lub jako kol. o elem. (x₁,x₂,...,x_n). Dowolną (n,n)-macierz naz. kwadratową stopnia n. Elementy a_ii tej m. przedst. gł. przekątną. J. wyznacznik W=|A|=detA≠0 to m. kw. A naz. nieosobliwą. J. wszystkie elem. m. kw. leżące pod(nad) gł. przekątną są=0 to naz. ją m. górną(dolną) trójkątną. M. diagonalna - wszystkie elem. oprócz stojących na gł. przek.=0. M. jednostkowa - gł. przek. zawiera 1 a reszta elem.=0 i oznaczamy ją E - AE=EA=A. M. A^-1 naz. odwrotną j. AA^-1=E. Łatwo spr. że (AB)^-1=A^-1B^-1. M. dołączoną m. A naz. m. C=[A_ij*] gdzie A_ij* są dopełnieniami algebraicznymi elem. a_ij wyznacznika |A|. M. A naz. skośnosymetryczną, j. A=-A'. Gdy A^-1=A' to m. A naz. ortogonalną, gdy A=Ā' - hermitowską. Dwie m. A i B naz. podobnymi j. istn. m. nieosobliwa C taka że B=C^-1AC. Mówimy wtedy że m. B jest otrzymana z m. A przekształceniem przez podobieństwo. Przez wielomian względem m. A rozumiemy P(A)=a₀E+a₁A+a₂A²+...+a_nAⁿ gdzie a₁,a₂,...,a_n-liczby dane, A²=AA,...,Aⁿ=AA^n-1. Niektóre znane z algebry: 1)detAB=detAdetB 2)Na to żeby A miała A^-1 potrzeba i wystarcza żeby |A|=detA≠0. M. odwrotną zapis. w post. A^-1=C/|A| gdzie C -m. dołączona natomiast |A^-1|=1/|A| 3)J. detA≠0 to ukł. AX^=b^ (1) gdzie A=[a_ij]_nxn X^=[x_n]_nx1 b^=[b_n]_nx1 ma jedno rozw. które zapis w post. X^=A^-1b^. 4)Na to żeby ukł. równ. jednorodnych Ax^=0 miał rozw. ≠0 potrzeba i wystarcza żeby detA=C. Ukł. n równań jednorodnych Ax^=λx^ (2) (λ-wart. własna m.) może być zapis w post. Ax^-λx^=0 ⇒ (A-λE)x^=0 i dlatego ma niezerowe rozw. w tym i tylko w tym przypadku gdy wyznacznik D(λ)=det(A-λE)=|m.kw.n-tego st. i od każdego a na gł. przek. -λ| (3) jest=0 (D(λ)=0). Równanie D(λ)=0 naz. równ. charakterystycznym lub wiekowym. D(λ) - wielomianem charakterystycznym m. A. Pierwiastki równ. D(λ)=0 naz. wartościami własnymi lub liczbami własnymi m. A, a odpowiadające im rozwiąz. niezerowe ukł. równ. (2) wektorami własnymi m. A. M. podobne mają jednakowe wielomiany charakt.: |C^-1AC-λE|=|C^-1AC-λC^-1EC|=|C^-1||A-λE||C|=|A-λE|. Łatwo zauważyć że m. przestawiona A' ma taki sam wielomian charakt. jak m. A ponieważ wartość wyznacznika przy zamianie wierszy przez kolumny i przeciwnie nie zmienia się. TW1. Niech λ₁,λ₂,...,λ_n będą wartościami własnymi m. A, a φ^₁,φ^₂,..., φ^_n odpowiadającymi im wektorami własnymi i Ψ^₁, Ψ^₂,...,Ψ^_n odpow. im wekt. wł. m. podst. A. Wówczas przy λ_i≠λ_j wektory φ^_i i φ^_j są ortogonalne tzn. φ_iΨ_j=∑_[k=1;n]ξ_ikΨ_jk=0 gdzie ξ_i1, ξ_i2,..., ξ_in -składowe wektorów φ^_i, Ψ_j1, Ψ_j2,...,Ψ_jn -składowe wektorów Ψ_j TW2. Wekt. własne odpow. różnym wartościom wlasnym są liniowo niezależne. Wekt. X^₁, X^₂,..., X^_n, naz. wekt. liniowymi niezależnymi, j. tożsamość ∑_[i=1;n]c_iX^_i=0 gdzie c₁, c₂,..., c_n -liczby, zachodzi tylko przy c₁=c₂=...=c_n=0. Rzeczywista (n,n) m. symetryczna ma n rzeczywistych liniowo niezależnych ortogonalnych wektorów własnych; wszystkie jej wartości własne są rzeczywiste. RYS3 TW3. Każdy n-wymiarowy wektor X^ może być przedst. w post. kombinacji liniowej n wektorów własnych rzeczywistej m. symetrycznej X^=∑_[i=1;n]c_iφ^_i gdzie c_i (i=1,...,n) -stałe. Często jest wygodnie normować wektory ortogonalne tj. dzielić każdy wekt. przez pierwiastek z sumy kwadratów jego składowych. TW4. J. n wektorów X^_i=(x_i1,x_i2,...,x_in) przedst. ukł. ortonormalny to ukł. n wektorów Y^_i=(x_1i,x_2i,...,x_ni) (i=1,2,...,n) jest też ukł. ortonormalnym. TW5. Dowolną (n,n) m. A można przedst. w post. iloczynu 2 m. trójkątnych - górnej i dolnej j. a₁₁≠0, |a₁₁ a₁₂; a₂₁ a₂₂|≠0, |A|≠0 TW6. J. A ma n-liniowo niezależnych wektorów własnych to przekształceniem przez podobieństwo można ją doprowadzić do postaci diagonalnej. Przy tym elementami gł. przek. w tej m. będą jej wart. własne. TW7. (Cayley'a - Hamiltona) J. wartości własne m. A spełn. równ. algebraiczne λⁿ+g_n-1λ^n-1+...+g₁λ+g₀=0 to Aⁿ+g_n-1A^n-1+...+g₁A+g₀E=0 2.Przestrzenie liniowe unormowane. Niech X^₁,X^₂,...,X^_n będzie ukł. n wektorów liniowo niezależnych o składowych (x_1i,x_2i,...,x_ni) (i=1,...,n). Wówczas zbiór wszystkich możliwych kombinacji liniowych wektorów X^=C₁X^₁+C₂X^₂,+...+C_nX^_n (4) gdzie C₁,C₂,...,C_n-dowolne stałe; przedstawia przestrzeń liniową n-wymiarową R. TW8. J. każdy z k wektorów Y^₁,Y^₂,...,Y^_k jest kombinacją innych m<k wektorów X^₁,X^₂,...,X^_m to wektory Y^₁,Y^₂,...,Y^_k są liniowo zależne. Każdy zbiór n liniowo niezależnych wektorów przestrzeni R nazywa się jej bazą. Np. jako bazę można wziąć wektory X^₁,X^₂,...,X^_n lub wektory e^₁=[100...0]^T,e^₂=[010...0]^T,...,e^_n=[00...1]^T Wypisanie powyżej n wektorów nazywamy bazą wyjściową przestrzeni R. Na mocy tw. 8 każdy wekt. jest kombinacją liniową wektorów bazy. Współczynniki tej kombinacji liniowej są określone jednoznacznie i naz. się składowymi (współrzędnymi) wektora w danej bazie. Np. X^=[x₁x₂...x_n]^T=x₁e^₁+...+x_ne^_n (4') gdzie x₁,x₂,...,x_n składowe wekt. X^ w bazie wyjściowej e^₁,e^₂,...,e^_n. Zobaczmy jak przekształcają się składowe wekt. X^ przy przejściu od bazy wyjściowej e^₁,e^₂,...,e^_n do dowolnej bazy X^₁,X^₂,...,X^_n. Niech C będzie m. której kolumny przedst. składowe wekt. X^₁,X^₂,...,X^_n w bazie wyjściowej C=[c_ij]_nxn. M. ta jest oczywiście nieosobliwa pon wekt. X^₁,X^₂,...,X^_n są liniowo niezależne; naz. się ona m. przekształcenia. Podstawiając X_i=∑_[j=1;n]c_jie^_j w równ. (4) otrzym. [x₁x₂...x_n]^T=[c₁₁c₁+c₁₂c₂+...+c_1nc_n;... ;c_n1c₁+c_n2c₂+...+c_nnc_n]=[c_ij]_nxn*[c₁...c_n]^T. W taki spos. składowe wyjściowe wekt. X^ otrzymujemy z jego skład. przekszt. przekształceniem linowym z m. C; X^=CX^', gdzie przez X^' oznaczono kolumnę utworzoną ze składowych wekt. X^ w bazie przekszt. Niech będzie dane przekszt. liniowe przestrzeni R w siebie Y^=AX^ gdzie A-(n,n) m. nieosobliwa nazwana m. przekszt. liniowego X^, Y^-wekt. wyznaczone przez swoje składowe w bazie wyjściowej. Zobaczmy jak zmienia się m. A przekszt. linowego przy przejściu od bazy wyjściowej e^₁,e^₂,...,e^_n do bazy przekształconej X^₁,X^₂,...,X^_n. Niech X^`,Y^` oznacz. wektory X^,Y^ w bazie. Wówczas X^=CX^';Y^=CY^' skąd wynika Y^`=C<sup>-1</sup>Y^`=C^-1AX^=C^-1ACX^` W ten spos. jednemu i temu samemu przekszt. w różnych bazach odpowiadają m. podobne ze współcz. podobieństwa C. Wartości własne i wektory własne m. A naz. także wart. własnymi i odpowiednio wektorami własnymi przekształcenia Y^=AX^. J. przekszt. liniowe ma w jakiejś bazie macierz diagonalną to baza skł. się z wekt. wł., elem. zaś stojące na przek. tej m. przedst. wartości wł. I przeciwnie j. w przestrzeni istn. baza złożona z wekt. wł. przekszt. to w tej bazie m. przekszt. A jest diagonalną. Wszystko to wynika bezpośrednio z równości AX^=λX^ określającej wartości własne i wektory własne przekształcenia Y^=AX^. Normą (analogia do długości) wektora X lub -ogólnie-dowolnego elem. przestrzeni liniowej nazywamy przyporządk. temu wekt. nieujemną liczbę rzeczyw. ||X|| spełn. nast. aksjomaty normy: 1)||X||>0 przy X≠0, ||0||=0 2)||αX||=|α| ||X|| α -liczba 3)||X+Y||≤||X||+||Y|| Każda przestrzeń liniowa w której wprowadzono pojęcie normy naz. przest. lin. unormowaną. Iloczynem skalarnym wektorów naz. przyporządkow. każdej parze wekt. X^ i Y^ liczby rzeczyw. (zespolonej) nazywanej iloczynem skalarowym XY=(X,Y) Aksjomaty skalarowego mnożenia: 1)XX>0 przy X≠0 i 00=0 2)XY=YX (XY=YX -w przestrzeni zespolonej) 3)(αX)Y=α(X,Y) α -liczba 4)(X+X₁)Y=XY+X₁Y Przestrzeń e której wprowadzono pojęcie iloczynu skalarowego i norma ||X||=√(XX) naz. przest. euklidesową. J. przestrzeń jest zespoloną a norma ||X||=√(XX) to naz. przest. unitarną. Np. w przestrz. R przyjmujemy XY=∑_[i=1;n]x_iŷ_i (też w zespolonej). Przez x₁,...,x_n y₁,...,y_n oznaczono składowe wekt. X i Y. W tejże samej przestrzeni można wziąć na równi z normą wektora określoną równością ||X||=√(XX) w charakterze ||X|| jedną z wielkośi max_[i]|x_i| (i-1,...,n) lub ∑_[1;n]|x_i|. Dla odróżnienia jednej od drugiej tych różnych postaci norm oznaczymy je wskaźnikami u dołu ||X||₁=max_[i]|x_i| (5), ||X||₂=∑_[i=1;n]|x_i| (6), ||X||₃=√(XX)=√(∑_[1;n]|x_i|) (7). Dwa wektory X i Y naz ortogonalnymi, j. XY=0. Układ wektorów Y₁,...,Y_K naz. ukł. ortogonalnym j. wszystkie te wekt. są parami ortogonalne; j. przy tym norma każdego z nich=1, to ukł. naz. ortogonalnym. Dla m. A tak samo jak dla wekt. może być wprowadzone pojęcie normy ||A|| spełniającej aksjomaty 1,2,3 normy a warunek ||AX||≤M||X|| gdzie X -dowilny wektor n-wymiarowej przestrzeni euklidesowej R. Określenie to jest równoważne równości ||A||=max|| A X/||X|| || (8) Konkretna postać normy macierzy zależy oczywiście od postaci normy przyjętej dla wektorów. Np. można wykazać że w przyp. normy (5)(6)(7) dla wekt. odpowiednie normy macierzy A są ||A||₁=max_[i]∑_[j=1;n]|a_ij| (9), ||A||₂=max_[j]∑_[j=1;n]|a_ij| (10), ||A||₃=√λ₁ (11) gdzie λ₁ -najw. wart. własna m. Istotnie w przyp (5) ||X||=1, wówczas mamy ||AX||=max|∑_[k=1;n]|a_ikx_k|≤max∑_[k=1;n]|a_ik||x_k|≤max_[i]∑_[k=1;n]|a_ik| zatem max||AX||=max_[i]∑_[k=1;n]|a_ik|. W przyp. (7) zauważymy że słuszne jest prawo przestawienia m. w iloczynie skalarowym (AX,Y)=(X,A*Y) gdzie A*-m. której elem. są zespolone sprzężone z A' przestawiona macierzy A. Niech λ₁≥λ₂≥...≥λ_n będą wartościami własnymi X₁,X₂,...,X_n ortonormalnym ukł. wektorów własnych. Obierzmy dowolny wekt. X z normą=1 i zapiszmy go za pomocą wekt. własnych X=c₁X₁+...+c_nX_n wtedy (X,X)=c₁²+c₂²+...+c_n²=1. Następnie ||AX||²=(X,A'AX)=(c₁X₁+...+c_nX_n, c₁λ₁X1+...+c_nλ_nX_n)=λ₁c₁²+λ₂c₂²+...+λ_nc_n²=λ₁(c₁²+...+c_n²)=λ₁ a poza tym, dla X=X₁ mamy ||AX₁||²=(X₁,A'AX₁)=(X₁,λ₁X₁)=λ₁ W taki spos. max||AX||=√λ₁ (11). Mówimy że ciąg wektorów X_k (ciąg m. A_k) przy k→∞ dąży do wektora X (do m. A) lim_[k→∞]X_k=X (lim_[k→∞]A_k=A), jeśli każda składowa wekt. X_k (każdy elem. m. A_k) przy k→∞ dąży do odpowiedniej składowej wekt. X (do odp. elem. m. A). Żeby lim_[k→∞]X_k=X potrzeba i wystarcza by lim_[k→∞]||X_k-X||=0 (lim_[k→∞]A_k=A ; lim_[k→∞]||A_k-A||=0) TW12. Na to by A^m→0 przy m→∞ potrz. i wyst. by wszystkie wartości własne m. A były co do wart. bezwzgl.<1. |λ₁|<1 (i=1,...,n) TW13. Wartość bezwzgl. każdej wart. wł. m. A jest większa od jej normy (dowolniej) |λ_i|≤||A|| (i=1,...,n) TW14. Na to by szereg E+A+A²+...+A^m+... był zbieżny potrz. i wyst. by A^m→0, gdy m→∞. Przy spełnieniu tego warunku zachodzi równość E+A+A²+...+A^m+...=(E-A)^-1 (12). Istotnie konieczność powyższego warunku jest oczywista. Pokażemy, że jest on war. dost. Na mocy twierdzenia (12) mamy |λ_i|<1 (i=1,...,n) co oznacza że równanie AX=1X (λ=1) nie ma rozwiąz. niezerowego i wyznacznik |E-A|≠0. Dlatego istnieje m. (E-A)^-1 odwrotna do (E-A). Mnożąc prawostronnie przez (E-A)^-1 tożsamość (E+A+A²+...+A^m)(E-A)=E-A^m+1 otrzymujemy E+A+A²+...+A^m=(E-A)^-1-A^m+1(E-A)^-1 Stąd wynika że S_m=E+A+...+A^m⇒(E-A)^-1 ponieważ A^m+1→0 gdy m→∞. Zauważmy jeszcze że (E-A)^-1-S_m=A^m+1+A^m+2+... ||(E-A)^-1-S_m||≤||A||^m+1+||A||^m+2+...=||A||^m+1/(1-||A||) Metody numeryczne rozw. ukł. równ. algebraicznych liniowych. 1.Wstęp. Rozwiązanie układu a₁₁x₁+a₁₂x₂+...+a_1nx_n=a_{1 n+1} ;...; a_n1x₁+a_n2x₂+...+a_nnx_n=a_{n n+1} (14) lub w postaci macierzowej AX=f (14') gdzie A=[a₁₁ ... a_1n ;...; a_n1 ... a_nn] X=[x₁ ... x_n]^T f=[a_{1 n+1} ... a_{n n+1}]^T nie przedstawia trudności z teoretycznego punktu widzenia i w przypadku gdy wyznacznik |A|=W≠0 wyraża się wzorami Cramera x_i=|A_i|/|A| (15) gdzie |A_i|-wyznacznik otrzymany z |A| drogą zamiany jegi i-tej kol. przez kol. wyr. wolnych. Ale rozwiązanie ukł. (14) przez obliczenie wart. liczbowych wg wzorów (15) przedst. znaczne trudności i do chwili obecnej jest przedm. zainteresowania matematyków. Metody rozw. można podzielić na 2 duże grupy: metody bezpośrednie i iteracyjne. M. bezpośrednie charakt. to że prowadzą zawsze do rozwiąz. ukł. j. ono istn. Ale duż liczba działań arytmet. które należy przy tym przeprow. pociąga za sobą niebezpieczeństwo popełnienia bł. co przy znacznej liczbie równań może dać bezwartościowy wynik. M. iteracyjne nie posiadają ostatniej wady, ale zastosow. do odkeśl. klasy ukł. konkretnej z tych metod nie zawsze daje proces zbieżny. Z metod bezpośr. przede wszystkim należy wymienić metodę eliminacji Gaussa i metodę rzutu ortogonalnego. Najbardziej rozpowsz. z met. iteracyjnych jest met. iteracji prostej, met. Seidla i met. relaksacji (met. poprawek). Podst. idea met. elimin. Gaussa polega na tym że sprowadzamy ukł. (14) do ukł. równoważnego z macierzą trójk. górną (postępowanie proste eliminacji). Z otrzymanego w taki spos. ukł. znajdujemy niewiadome przez kolejne podstawienia (postępow. odwrotne eliminacji). Przyjmując że nazywany elem. kierującym współczynnik a₁₁≠0 dzielimy pierwsze z równań ukł. przez a₁₁. Następnie mnożąc otrzymane w taki spos. równ. przez a₂₁, a₃₁,...,a_n1 odejmujemy od odpowiednich równań ukł. (14). W rezultacie dochodzimy do układu x₁+b₁₂x₂+...+b_1n=b_{1 n+1} ; c₂₂x₂+a_23*1x₃+...a_2n*1x_n=a_{2 n+1} ;...; c_n2x₂+a_n3*1, x₃+...a_nn*1x_n=a_{n n+1*1} (14₁) Postępując z tym ukł. (bez uwzględnienia pierwszego równ.) tak samo jak z ukł. (14) otrzymamy ukł. x₂+b₂₃x₃+b₂₄x₄+...+b_2nx_n=b_{2 n+1} ; c₃₃x₃+a_34*2x₄+...+a_3n*2x_n=a_{3 n+1*2} ;...; c_n3x₃+a_n4*2x₄+...a_nn*2x_n=a_{n n+1*2} (14₂) Powtarzając to postępowanie otrzymamy w końcu ukł. x_n-1+b_{n-1 n}x_n=b_{n-1 n+1} ; c_nnx_n=a_{n n+1*n-1} ; x_n=b_{n n+1} Biorąc w każdym ukł. (14₁)(14₂)...(14_n) pierwsze równania otrzymujemy szukany układ równań z m. trójkątną x₁+b₁₂x₂+b₁₃x₃+...+b_1nx_n=b_{1 n+1} ; x₂+b₂₃x₃+...b_2nx_n=b_{2 n+1} ;...; x_n-1+b_n-1x_n=b_{n-1 n+1} ; x_n=b_{n n+1} (16). Podany schemat eliminacji znany jest pod nazwą schematu jedynego dzielenia. Ten schemat zawsze istn. w tych przypadkach gdy elem. kierujące każdego z ukł. (14₁)(14₂)...(14_n) są ≠0. J. jednak którykolwiek z nich zeruje się to w odpowiednim ukł. wystarczy zmienić numeracje równań tak by można było powtórzyć proces eliminowania. W pewnych przypadkach dla otrzymania większej dokładności obieramy w charakt. współczynnika kierującego ukł. największy co do wart. bezwzgl. elem. A. Drogą eliminow. odpowiadającej niewiadomej otrzymujemy ukł. równ. analogiczny do ukł. (14₁). Z tego ukł. w taki sam spos. wybieramy w charakt. kierującego najw. co do wart. bezwzgl. współczynnik tego ukł. i otrzym. ukł. analogiczny do (14₂). Powtarzając to postępow. otrzymujemy szuk. ukł. równ. z m. trójkątną. Taki schemat eliminacji jest znany pod nazwą schematu z wyborem elementów największych. 2.Metoda eliminacji Gaussa. Rozkład LU. Etap pierwszy (zwany etapem elimin. zmiennych). Dany jest ukł. równań. A⁽¹⁾x^=b^⁽¹⁾, postaci a₁₁⁽¹⁾x₁+...+a_1n⁽¹⁾x_n=b₁⁽¹⁾ ; a₂₁⁽¹⁾x₁+...+a_2n⁽¹⁾x_n=b₂⁽¹⁾ ;...; a_n1⁽¹⁾x₁+...+a_nn⁽¹⁾x_n=b_n⁽¹⁾ (1). Odejmując od i-tego wiersza tegi ukł. (i=2,3,...,n) wiersz pierwszy pomnożony przez a_i1⁽¹⁾/a₁₁⁽¹⁾ otrzymujemy ukł. A⁽²⁾x^=b^⁽²⁾ postaci a₁₁⁽²⁾x₁+...+a_1n⁽²⁾x_n=b₁⁽²⁾ ; a₂₂⁽²⁾x₂+...+a_2n⁽²⁾x_n=b₂⁽²⁾ ;...; a_n2⁽²⁾x₂+...a_nn⁽²⁾x_n=b_n⁽²⁾ Wyeliminowaliśmy w ten spos. pierwszą niewiadomą x₁ z równań leżących w wierszach o nr i=2,3,...,n. Podobnie eliminujemy x₂ z równań leżących w wierszach 3,4,...,n odejmując od i-tego wiersza (i-3,4,...,n) wierszdrugi pomnożony przez a_i2⁽²⁾/a₂₂⁽²⁾. Postępując w ten spos. otrzymamy kolejne przekszt. ukł. równ. A⁽³⁾x^=b^⁽³⁾, A⁽⁴⁾x^=b^⁽⁴⁾, ..., i po (n-1) eliminacjach uzyskamy trójkątny ukł. równań A⁽ⁿ⁾ℋ ^=b^⁽ⁿ⁾, postaci a₁₁⁽ⁿ⁾x₁+...+a_1n⁽ⁿ⁾x_n=b₁⁽ⁿ⁾ ; a₂₂⁽ⁿ⁾x₂+...+a_2n⁽ⁿ⁾x_n=b₂⁽ⁿ⁾ ;...; a_nn⁽ⁿ⁾x_n=b_n⁽ⁿ⁾ (3). Etap drugi (zwany postępow. odwrotnym). Otrzymany po przekształceniu ukł. rozwiązujemy zgodnie ze wzorami (4) x_n=b_n⁽ⁿ⁾/a_nn⁽ⁿ⁾, x_i=(b_i⁽ⁿ⁾-a_in⁽ⁿ⁾x_n-...-a_{i i+1}⁽ⁿ⁾x_i+1)/a_ii⁽ⁿ⁾ (i=n-1,n-2,...,1). Można obliczyć że w celu wyznaczenia tą metodą rozwiązania x^ należy wykonać łącznie M=n³/3+n²-n/3 mnożeń i dzieleń oraz D=n³/3+n²/2-5n/6 dodawań. Jeśli nie musimy zachować danych A⁽¹⁾ i b⁽¹⁾ do dalszych obliczeń lub sprawdzenia to min. liczba zajętych komórek pamięci P=n²+n, w przeciwnym razie P'=2P. Liczba wykonyw. w met. Gaussa jest bez porównania mniejsza niż w metodzie Cramera. Tam bowiem należało obliczyć (n+1) wyznaczników, które w przyp. stosow. rozwijania wyznacznika wg kolumny lub wiersza wymagają co najmniej (n+1)! mnożeń. Zatem w przypadku 15 równań należałoby wyk. 4922789888000 ≈ 5*10¹² mnożeń i maszyna cyfrowa wykonująca 10⁶ mnożeń na sek. musiałaby pracować nieprzerwanie przez ponad rok, natomiast w met. Gaussa M=1345 mnożeń, a więc komputer wykona pracę w niecałe 0,01 sek. W wielu zagadnieniach num. dotyczących m. kwadratowych A, poszukujemy takich m. trójkątnych L i U by A=LU. Met. eimin. Gaussa umożliwia wyznaczenie takiego rozkładu. Zauważmy, że przekształcenie układu A⁽¹⁾x^=b^⁽¹⁾ do postaci A⁽²⁾x^=b^⁽²⁾ jest równoważne pomnożeniu obu str. ukł. (1) A⁽¹⁾x^=b^⁽¹⁾ przez m. [1 0 0 ... 0 ;-l₂₁ 1 0 0 ... 0 ;-l₃₁ 0 1 0 0 ... 0 ;-l₃₁ 0 0 1 0 0 ... 0 ;...;-l_n1 0 0 ... 1]=L⁽¹⁾, l_i1=a_i1⁽¹⁾/a₁₁⁽¹⁾, i=2,3,...,n. Mamy zatem 2 równości: L⁽¹⁾A⁽¹⁾=A⁽²⁾, L⁽¹⁾b^⁽¹⁾=b^⁽²⁾. Analogicznie otrzym. 2 kolejne równości: L⁽²⁾A⁽²⁾=A⁽³⁾, L⁽²⁾b^⁽²⁾=b^⁽³⁾ gdziel_i2=a_i2⁽²⁾/a₂₂⁽²⁾, L⁽²⁾=[1 0 ... 0 ; 0 1 0 ... 0 ; 0 -l₃₂ 1 0 ... 0 ;...; 0 -l_n2 ... 1] i=3,4,...,n. Postępując w ten spos. otrzymamy L^(n-1)L^(n-2)...L⁽¹⁾A⁽¹⁾=A⁽ⁿ⁾ L^(n-1)L^(n-2)...L⁽¹⁾b^⁽¹⁾=b^⁽ⁿ⁾. Macierze L⁽ⁱ⁾, i=1,2,...,n-1 są nieosobliwe, możemy zatem napisać A⁽¹⁾=(L⁽¹⁾)^-1(L⁽²⁾)^-1...(L^(n-1))^-1A⁽ⁿ⁾ (4) Zauważmy że (L⁽¹⁾)^-1=[1 0 ... 0 ; l₂₁ 1 0 ... 0 ; l₃₁ 0 1 0 ... 0 ;...; l_n1 ... 1], (L⁽²⁾)^-1=[1 0 ... 0 ; 0 1 0 ... 0 ; 0 l₃₂ 1 ... 0 ;...; 0 l_n2 0 ... 1] a ponad to L=(L⁽¹⁾)^-1(L⁽²⁾)^-1(L^(n-1))^-1=[1 0 ... 0 ; l₂₁ 1 0 ... 0 ; l₃₁ l₃₂ 1 ... ;...; l_n1 l_n2 l_n3 ... 1]. J. oznaczymy m. A⁽ⁿ⁾ symbolem U to rozkład (4) m. A⁽¹⁾ na iloczyn możemy zapisać następująco A⁽¹⁾=LU, gdzie L -m. trójkątna dolna z jedynkami na diagonalnej, U -m. trójkątna górna. J. etap elimin. zmiennych met. Gaussa można wykonać do końca to ukł. równań Ax^=b^ można zapisać jako Lux^=b^. Etap ten jest zatem równoważny przekształceniu wyjściowego ukł. do postaci Ux^=L^-1b^=b^⁽ⁿ⁾, a wyznaczenie rozwiąz. postaci x^ polega na rozwiąz. 2 ukł. z macierzami trójkątnymi Ly^=b^, Ux^=y^. Warianty stabilne. Met. Gaussa nie jest dozwolona. Np. [0 2 2 ; 3 3 0 ; 1 0 1]x^=[b₁ b₂ b₃]^T M. tego ukł. jest nieosobliwa, istn. zatem dokł. jedno rozw. Jednak ukł. tego nie można sprow. do postaci trójkątnej pon. a₁₁=0=u₄. Łatwo stwierdzić że wersja algorytmu Gaussa z częściowym wyborem elelm. podst., zwana czasem met. Gaussa-Crouta. (GCW) jest met. niezawodną przy założeniu dokładnych działań arytmet. nie nastąpi zatrzymanie procesu obliczeń z powodu dzielenia przez zero. Podamy ważne przykł. m., w stos. do których nie trzeba stos. strategii częściowego wyboru elem. podst. M. kwadratowa A naz. diagonalną dominującą gdy moduły elem. diagonali są niemniejsze od sumy modułów pozostałych elem. stojących w tym samym wierszu |a_ii|≥∑_[k=1;n]|a_ik|, i=1,..,n Np. A=[3 1 0 ; 0 2 2 ; 4 1 5] jest diagonalnie dominująca. M. A jest diag. dom. kolumnowo, j. A'=A^T jest diagonalnie dominująca |a_ii|≥ (>) ∑_[k=1;n]|a_ik|, i=1,...,n 3.Metoda eliminacji Jordana. Ukł. równ. postaci (1) A⁽¹⁾x^=b^⁽¹⁾ przekształc. następująco. Pierwsze równanie dzielimy obustronnie przez a₁₁⁽¹⁾, a następnie do i-tego wiersza, i=2,3,...,n odejmujemy wiersz pierwszy pomnożony przez a_i1⁽¹⁾, otrzymując układ A⁽²⁾x^=b^⁽²⁾: x₁+a₁₂⁽²⁾x₂+...+a_1n⁽²⁾x_n=b₁⁽²⁾; a₂₂⁽²⁾x₂+...+a_2n⁽²⁾x_n=b₂⁽²⁾; ...; a_n2⁽²⁾x₂+...+a_nn⁽²⁾x_n=b_n⁽²⁾. Następnie drugie równanie dzielimy obustronnie przez a₂₂⁽²⁾ i od i-tego wiersza, i=1,3,4,...,n odejmujemy wiersz drugi pomnożony przez a_i2⁽²⁾, otrzymując ukł. A⁽³⁾x^=b^⁽³⁾: x₁+a₁₃⁽³⁾x₃+...+a_1n⁽³⁾x_n=b₁⁽²⁾; x₂+a₂₃⁽³⁾x₃+...+a_2n⁽³⁾x_n=b₂⁽²⁾; ...; a_n3⁽³⁾x₃+...+a_nn⁽³⁾x_n=b_n⁽²⁾. Po (n-1) eliminacjach otrzymujemy ukł. x₁=b₁⁽ⁿ⁾; x₂=b₂⁽ⁿ⁾; ...; x_n=b_n⁽ⁿ⁾ -a wiec gotowe rozwiązanie. Liczba działań jest tu około półtora raza większa niż dla met. Gaussa. Nie można bowiem dla met. Jordana znaleźć odpowiednika rozkładu A=LU jak dla met. Gaussa. 4.Układy z m.trójdiagonalną. Tx^=d^, T=[b₁ c₁ 0 ... 0 ; a₂ b₂ c₂ 0 ... 0 ; 0 a₃ b₃ c₃ 0 … 0 ; 0 0 a₄ b₄ c₄ 0 … 0 ;…; 0 … 0 a_n b_n]; b₁x₁+c₁x₂=d₁; a_ix_i-1+b_ix_i+c_ix_i+1=d_i, i=2,3,…,n-1; a_nx_n-1+b_nx_n=d_n. Ukł. równ. Tx^=d^ można rozwiązać nadzwyczaj szybko, wykon. małą liczbę działań arytm. Zauważmy na wstępie że jeśli można dokonać rozkładu T na iloczyn LU, to macierze L i U mają specyficzną postać: L=[1 0 ... 0 ; l₂ 1 0 ... 0 ; 0 l₃ 1 0 ... 0 ;...; 0 ... 0 l_n 1]; U=[u₁ c₁ 0 ... 0 ; 0 u₂ c₂ 0 ... 0 ; 0 0 u₃ c₃ 0 ... 0 ; 0 ... 0 u_n c_n] Elementy l₂,…,l_n oraz u₁,...,u_n można obl. rekurencyjnie z wzorów u₁=b₁, l_i=a_i/u_i-1, u_i=b_i-l_ic_i-1, dla i=2,3,...,n. Liczba działań wynosi M=2n-2, D=n-1, a liczba zajętych komórek pamięci P=3n-2. Aby znaleźć rozw. ukł. Tx^=d^, czyli LUx^=d^, należy rozw. kolejno 2 ukł. Lx^=d^ oraz Ux^=d^ według wzorów y₁=d₁, y_i=d_i-l_iy_i-1 (i=2,3,...,n) oraz x_n=y_n/u_n, x_i=(y_i-c_ix_i+1)/u_i dla i=n-1, n-2, ...,1. Łącznie należy wykonać dodatkowo M=3n-2, D=2n-2 działań. J. m. T jest diagonalnie dominująca kolumnowo |b₁|≥|a₂|, |b_i|≥|c_i-1|+|a_i+1| dla i=2,3,...,n-1 |b_n|≥|c_n-1|, to rozkład T=LU jest równoważny wersji rozkładu z częściowym wyborem elementu podst., która jest niezawodna. Ponadto elementy m. L spełniają warunek |l_i|≤1 i=2,...,n co poprawia oszacowanie błędów zaokrągleń. J. T jest diagonalnie dominująca |b_i|≥|c_i|, |b_i|≥|a_i|+|c_i|, i=2,...,n |b_n|≥|a_n|, to met elimin. Gaussa też jest niezawodna, ale stosujemy inną wersję rozkładu T=LU. L=[l₁ 0 ... 0 ; a₂ l₂ 0 … 0 ; 0 … 0 a_n l_n], U=[1 u₁ 0 ... 0 ; 0 1 u₂ 0 … 0 ; 0 … 0 1]; l₁=b₁; y₁=d₁/l; u_i=c_i/l_i; l_i=b_i+1-u_ia_i+1; y_i+1=(d_i+1-a_i+1y_i)/l_i; i=1,...,n-1 oraz x_n=y_n, x_i=y_i-u_ix_i+1 dla i=n-i,n-2,...,1. 5.Układy z m.symetryczną. Rozkład Banachiewicza (Cholesky'ego) Dla symetrycznych m. dodatnio określonych istn. jeszcze jeden typ rozkł. na iloczyn m. trójkątnych zwany rozkładem Banachiewicza (w wielu publikacjach anglosaskich -Cholesky'ego) A=L L^T (L^T=L'). M. L jest m. trójkątną dolną przy czym na digonali mogą być elem. nieskończenie=1. Pon. m. L*(sprzężone)=-L też spełnia warunek A=(L L^T)* dla A=L L^T, więc rozkład taki nie jest jednoznaczny. J. jednak założyć że liczby na diagonali m. L są dodatnie to otrzymujemy jednoznaczne wzory obliczeniowe w sposób następujący. Obliczamy kolejno: l_ii=√(a_ii-∑_[k-1;l-1]l_ik²), i=1,...,n (5.1) l_ji=(a_ji-∑_[k=1;i-1](l_jk / l_ik))/l_ii, j=i+1,...,n Np. Dla m. A=[4 2 2 ; 2 5 3 ; 2 3 6] znajdujemy i=1: l₁₁=√4=2, l₂₁=2/2=1, l₃₁=2/2=1; i=2: l₂₂=√(5-1*1)=2, l₃₂=(3-1)/2=1; i=3: l₃₃=√(6-1*1-1*1)=2, zatem A=[4 2 2 ; 2 5 3 ; 2 3 6]=[2 0 0 ; 1 2 0 ; 1 1 2]*[2 1 1 ; 0 2 1 ; 0 0 2]. Dodatnia określoność A gwarantuje że elem. l_ii na diagonali są niezerowe a bład wyznaczenia rozkładu A=LL^T jest niewielki, zwłaszcza jeśli działanie wykonujemy w zwiększonej precyzji. 6.Obliczanie wyznacznika i odwracanie m. J. zachodzi konieczność obliczenia wyznacznika najlepiej dokonać rozkładu m. A=LU (stosując częściowy wybór elem. podstawowego), a wówczas detA=detLU=detL*detU. Ponieważ elementy na diagonali m. L są=1 (l_ii=1) więc detL=1 wtedy detA=detLU=detU (6.1). Wyznacznik m. trójk. U=iloczynowi elem. na diagonali i można wyznaczyć go z dużą dokł. wyk. (n-1) mnożeń. detA=detU=u₁₁u₂₂...u_nn (6.2) Wyznaczanie m. odwrotnej A^-1 met. GCW. Rozpatrzymy ukł. równ. {a₁₁x₁+a₁₂x₂+...+a_1nx_n=1; a₂₁x₁+a₂₂x₂+...+a_2nx_n=0 ;...; a_n1x₁+a_n2x₂+...+a_nnx_n=0} (6.3) Oznaczmy rozw. rego ukł. przez x₁⁽¹⁾,...,x_n⁽¹⁾. Wówczas na podst. wzorów Cramera mamy |A|x₁⁽¹⁾=A_11, ..., |A|x_n⁽¹⁾=A_1n (6.4) gdzie A_1j (j=1,...,n) są dopełnieniami agbraicznymi odpowiednich elem. a_1j. J. oznaczymy przez x₁⁽²⁾,...,x_n⁽²⁾ rozw. ukł. (6.3), którego prawe strony są równe 0,1,0...0 to w podobny spos. otrzymamy |A|x₁⁽²⁾=A_21, ..., |A|x_n⁽²⁾=A_2n (6.5) itd. Oznaczając wreszcie przez x₁⁽ⁿ⁾,...,x_n⁽ⁿ⁾ rozw. ukł., którego prawe strony są równe 0,0,...,0,1 otrzym. |A|x₁⁽ⁿ⁾=A_n1, ..., |A|x_n⁽ⁿ⁾=A_nn (6.6) Utwórzmy obecnie tablicę której kol. są rozw. wskazanych rówm.: [x₁⁽¹⁾,...,x₁⁽ⁿ⁾ ;x₂⁽¹⁾,...,x₂⁽ⁿ⁾ ;...;x_n⁽¹⁾,...,x_n⁽ⁿ⁾] Ta tablicz przedstawia m. odwrotną A^-1 pon A^-1=[x₁⁽¹⁾,...,x₁⁽ⁿ⁾ ;x₂⁽¹⁾,...,x₂⁽ⁿ⁾ ;...;x_n⁽¹⁾,...,x_n⁽ⁿ⁾]=[A₁₁/|A| A₂₁/|A| ... A_n1/|A| ;...; A_1n/|A| A_2n/|A| ... A_nn/|A|]. Każdy ukł. równ. typu (6.3) rozw. za pom. met. GCW. 7.Met.iteracyjne. Met.iteracji prostej. Niech ukł. równ. lin. będzie dany w post. x₁=b₁₁x₁+...+b_1nx_n+g₁ ;...; x_n=b_n1x₁+...+b_nnx_n+g_n lub w post. macierzowej (7.1) X=BX+G albo (E-B)X=G. W charakt. pierwszego przybliżenia obieramy dowolny wekt X₀ (np. X₀=G) i podst. go w prawą str. równości (7.1). Z lewej str. równości otrzymamy pewien wektor X₁=BX₀+G. Postępując z X₁ tak samo jak z X₀ otrzym. wekt. X₂. Powtarzając to postępow. otrzym. ciąg wektorów X_m+1=BX_m+G (7.2). Oczywiście że jeśli przy (m=0,1,...) m→∞ mamy że X_m→X to wekt. X będzie rozw. równ. (7.1). Dla zbieżności tego procesu konieczne jest żeby m. B była mała w tym lub innym sensie, tj. żeby w m. (E-B) elem. wyst. na gł. diagonali przewyższały pozostałe elem. i były zbliżone do jedności. TW: (o zbieżności met. iteracji prostej). Dla zbieżności met. iteracji prostej przy dowolnym wekt. X₀ i przy dowolnej wart. wekt. G wyrazów wolnych potrzeba i wystarcza by wszystkie wart. wł. m. B były co do wart. bezwzgl.<1. Dowód Wychodząc ze wzoru (7.2) możemy wyrazić wekt. X_m przez początkowy wekt. X₀. Otrzym.: X_m=BX_m-1+G=B(BX_m-2+G)+G=B²X_m-2+(B+E)G=B²(BX_m-3+G)+(B+E)G=B³X_m-3+(B²+B+E)G=... B^mX₀+(B^m-1+B^m-2+…+B+E)G. W taki spos. mamy X_m=B^mX₀+(E+B+B²+...+B^m-1)G (7.3). Stąd od razu wynika że dla istn. lim_[m→∞]X_m przy dowolnym wekt. G i dowolnym X₀ wystarczy żeby szereg E+B+B²+...+B^m-1+... (7.4) był zbieżny. Warunkiem koniecznym i dostat. zbieżności szeregu (7.4) jest to by wszystkie wart. własne λ_i m. B były co do wart. |λ_i|<1. Tym samym tw. zostało udowodnione. Sprawdzenie spełnienia warunków tw. z reguły jest kłopotliwe dla tego w praktyce bywają bardzo przydatnespecjalne kryteria dostateczne zbieżności, wynikające z następującegi twierdzenia: TW: Dla zbieżności procesu iteracji prostej wystarcza by którakolwiek z norm m. B była<1, ||B||<1. Np. Należy rozw. ukł. {10x₁+2x₂+x₃=10; x₁+10x₂+2x₃=12; x₁+x₂+10x₃=8} Pon. elem. na gł. diagonali |a_ii|>∑_{[j=1,i≠j;3]}|a_ij|, i=1,2,3 więc met. iteracji prostej jest zbieżna. Rozwiązując ukł. względem niewiadomych stojących na gł. przek. otrzym.: {x₁=-0,2x₂-0,1x₃+1; x₂=-0,1x₁-0,2x₃+1,2; x₃=-0,1x₁-0,1x₂+0,8} ⇒ B=[0 -0,2 -0,1 ;-0,1 0 -0,2 ;-0,1 -0,1 0] ⇒ ||B||<1 Obliczenia kolejnych przybliżeń X_m wykonujemy na podst. wz. (7.2) X_m=BX_m-1+G. Obierając jako pierwsze przybliżenie X₀=0 otrzym. nast. wyniki RYS4 Widzimy że met. iteracji prostej w danym przyp. jest dość szybko zbieżna. 8.Met.Seidala różni się od met. iteracji prostej tylko tym że przy obliczaniu k-tej składowej wekt. X_m+1 (m+1)-ego przybliżenia wykorzystujemy obliczone już za pom. ukł. (7.1) k-1 pierwszych składowych tego wekt. tzn. obliczenie kolejnych przybliżeń prowadzimy wg wzorów x_i^(m+1)=∑_[j=1;i-1]b_ijx_j^(m+1)+∑_[j=1;m]b_ijx_j^(m)+G_i (8.1) TW: J. pierwsza norma m. B ||B||₁=max_[i]∑_[j=1;n]|b_ij|=ν<1 to met. Seidla jest zbieżna. 9.Met.Jacobiego. Niech m. A ukł. równ. Ax=b będzie taką m. że można ją rozłożyć na składniki A=L+D+U gdzie L jest m. poddiagomalną, D -m. diagonalną, U -m. naddiagonalną: A=[1 2 3 ; 4 5 6 ; 7 8 9]=(L)[0 0 0 ; 4 0 0 ; 7 8 0]+(D)[1 0 0 ; 0 5 0 ; 0 0 9]+(U)[0 2 3 ; 0 0 6 ; 0 0 0] wówczas Ax=b ⇒ Lx+Dx+Ux=b ⇒ Dx=-(L+U)x+b. Otrzym. wz. iteracyjny Jacobiego Dx^(k+1)=-(L+U)x^(k)+b (9.1) k=0,1,... Aby móc zast. wz. (9.1) należy tak zmienić kolejność równ. ukł. Ax=b by na diagonali były elem. niezerowe. W przyp. gdy wymiar m. jest duży wprowadzenie takiej zmiany kolejn. równ. jest trudne. Istn. jednak metody automatycznego przestawienia wierszy m. A. Jedna z tych met.: a) spośród kol., w któeych na diagonali znajduje się elem. zerowy, wybieramy tę, w której jest największa liczba 0; b) w kol. tej wybieramy elem. o maks. module i tak przestawiamy wiersze by wybrany elem. znalazł się na diagonali: wiersz ten ustalamy i omijamy go w dalszym badaniu wartości elem. c) spośród pozostałych kol., w których na diagonali znajduje się elem. zerowy wybieramy tem, której jest największa l. 0; d) dokonując kolejnych przestawień wierszy doprowadzamy ostatecznie wyjściową m. do takiej postaci m., w której na diagonali są elem. ≠0. Met. Jacobiego jest niezawodna. Przybliżone rozw. równań nieliniowych i ich układów. 1.Jedno równanie z jedną niewidomą. Jak wiadomo nie każde równanie można rozwiązać dokładnie. Nie można podać metody na rozwiązanie dowolnego równania algebraicznego stopnia wyższego niż czwarty. Dokładne rozw. równ. nie zawsze jest też konieczne pon. w równaniach często mamy jedynie przybliżone wart. współczyn. i zadanie dokładnego wyznaczenia pierwiastków traci sens. Rozpatrzymy zadanie przybliżonego obliczania pierwiastków rzeczywistych równania f(x)=0. Gdy f(x) jest ciągła na pewnym przedziale (skończonym lub nieskończonym). Większość met. przybliżonego rozw. równań polega na poprawianiu kolejnych przybliżonego pierwiastków. Ponadto większość met. można stos. jedynie wtedy, gdy znamy przedział, w którym znajduje się pojedynczy pierwiastek. O takim przedziale mówi się, że jest przedziałem izolacji pierwiastka. Dla wyznaczenia przedziałów izolacji można wykorzystać wykres f. y=f(x), albo j. jest to trudne, równanie wyjściowe f(x)=0 należy przekształcić do równoważnej mu postaci φ(x)=ψ(x) gdzie wykresy y₁=φ(x) i y₂=ψ(x) są prostsze do wykonania. Odcięte punktów przecięcia tych wykresów będą pierwiastkami równoważnego równania. J. f. f(x) jest wielomianem algebraicznym to istn. metody analitycznego izolowania pierwiastków. Będziemy zakładali że f(x) jest klasy C² (f(x)ЄC²) a także będziemy przyjmować dodatkowe warunki na jej pierwszą f'(x) i jej drugą f''(x) pochodną. Mamy do czynienia z 3 źródłami błędów uzyskiwanych przybliżeń: 1)błąd przybliżenia wynikający ze stosowanej met. 2)wpływ błędów współczyn. równania 3)dokładność obliczeń przybliżonych 1.1 Metoda połowienia. Mamy dane równania f(x)=0 (1) przy czym f(x) jest ciągła na przedziale domkniętym [a;b] wewnątrz którego znajduje się dokł. jeden pierwiastek i na którego końcach wartości f(a) i f(b) mają przeciwne znaki ⇒ f(a)f(b)<0. Aby znaleźć pierwiastek dzielimy przedział [a;b] na połowy punktem x₁=(a+b)/2. J. f(x₁)=0, to x₁ jest szukanym pierwiastkiem a jak f(x₁)≠0 to z otrzymanych dwóch przedziałów [a;x₁] i [x₁;b] wybieramy ten na końcach którego f(x) ma przeciwne znaki. Z kolei ten przedział dzielimy na 2 połowy, ponownie badamy wartość f(x₂) i znaki f(x) na końcach otrzymanych przedziałów itd. W wyniku takiego postępow. po pewnej ilości kroków albo otrzymamy pierwiastek dokładny f(x_n)=0 albo ciąg przedziałów takich że f(x_i)f(x_i+1)<0 (2)przy czym x_i oraz x_i+1 są odpowiednio początkiem i końcem i-tego przedziału a jego długość |x_i=x_i+1|=(b-a)/2ⁱ (3). Pon. lewe końce ciągu przedziałów tworzą ciąg niemalejący i ograniczony z góry, a prawe końce ciąg nierosnący i ograniczony z dołu więc z (3) wynika że istnieje ich wspólna granica lim_[x→∞]|x_i-x_i+1|=lim_[x→∞](b-a)/2ⁱ=0 ⇒ lim_[x→∞]x_i=lim_[x→∞]x_i+1=ξ pierwiastek α. Ten spos. znajdowania pierw. naz. met. połowienia także met. równego podziału albo met. bisekcji. Np. Znaleźć pierwiastek x³+x²-3x-3=0 położony w przedziale [1,2]. Mamy f(x)=x³+x²-3x-3 zatem f(1)=-4 i f(2)=3 i f(1)f(2)<0. Pon. pochodna f'(x)=3x²+2x-3>0 na [1,2] więc jest to przedział izolacji pierwiastka. Postępując zgodnie z met. połow. obliczamy kolejne przybliżenia są zestawione w tablicy RYS.5 Schemat blokowy: RYS.6 Zmienne A i B określają końce przedziału w każdej z kolejnych iteracji: YA i YB są to wart. f. y=f(x) odpowiednio w tych punktach a zmienna EPS to przyjęta dokładność z jaką chcemy obliczyć pierwiastek. Gdy długość przedziału w danej iteracji (B-A) jest mniejsza od EPS następuje przerwanie procedury. Oczywiście to założone EPS powinno być większe od maksymalnej granicznej dokładności obliczeń, w przeciwnym razie gdy EPS będzie mniejsze możemy otrzymać wynik jedynie z maksymalną graniczną dokładnością. 2.Reguła falsi. Nazwa metody pochodzi od łacińskich słów: regula - linia i falsus - fałszywy; jest to zatem met. fałszywego założenia liniowości f.. Przyjmijmy że w przedziale [a;b] ma dokładnie jeden pierwiastek i że f(a)f(b)<0. Ponadto niech f(x)ЄC² i f'(x), f''(x) mają stały znak na [a;b]. Wobec tych założeń wykres f. y=f(x) może być tylko jednym z czterech rodzajów: RYS.7 Rozpatrzymy przypadek gdy na [a;b] pochodne f'(x)>0 i f''(x)>0. Dla pozostałych przypadków rozumowanie jest analogiczne. Przez punkty A(a,f(a)) i B(b,f(b)) prowadzimy cięciwę o równaniu y-a=[(f(a)-f(b))/(b-a)](x-a) (4). Odciętą X, punktu w którym cięciwa AB przecina oś OX przyjmujemy za pierwsze przybliżenie szukanego pierw. Stąd pon. y(x₁)=0 mamy x₁=a-[f(a)/(f(b)-f(a))](b-a) (5). J. f(x₁)=0 to oczywiście x₁ jest pierw. α i zadanie zostało zakończone. Załóżmy że f(x₁)≠0. Przez punkt C(x₁,f(x₁)) oraz ten z punktów A, B którego rzędna ma przeciwny znak niż f(x₁) prowadzimy następną cięciwę (rys.7(5)). Odcięta x₂ punktu w którym ta cięciwa przetnie oś OX da nam drugie przybliżenie pierw. α. Dla f'(x)>0, f''(x)>0 w kolejnych przybliżeniach punktu B pozostaje nieruchomy. (Gdy f'(x)>0 i f''(x)<0 (rys.7(4)) wówczas nieruchomy byłby punkt A). J. przybliżenie x₂ jest nadal niewystarczające to przez punkty B i D(x₂,f(x₂)) prowadzimy trzecią cięciwę co daje nam trzecie przybliżenie x₃ itd. Otrzymujemy kolejne wyrazy ciągu x₁,x₂,...,x_n,... (6) określonego wzorem rekurencyjnym: x₀=a, x_k+1=x_k-[f(x_k)/(f(b)-f(x_k))](b-x_k) (7) k=1,2,...,n. Można wykazać że przy przyjętych założeniach ciąg (6) {x_i} jest rosnący i ograniczony a więc zbieżny i jego granica jest szukanym pierw. α. Przy nieruchomym punkcie B, kolejne wyrazy ciągu {x_i} są mniejsze od α oraz f(x_i)<0 dla każdego i. Przechodząc do granicy dla n→∞ z równości (7) otrzymujemy g=g-[f(g)/(f(b)-f(g))](b-g) przy czym g=lim_[g→∞]x_n (a<g<b). Stąd f(g)=0 i g≡α. 3.Metoda Newtona. Załóżmy że αЄ[a;b], f(a)f(b)<0 oraz f'(x), f''(x) mają stały znak na [a;b]. W met. Newtona zwanej też met. stycznych z tego końca przedziału [a;b] w którym f(x) ma ten sam znak co f''(x) prowadzimy styczną do wykresy f. y=f(x) i punkt x₁, w którym ta styczna przecina oś OX przyjmujemy za pierwsze przybliżenie pierwiastka α. RYS.8 Oczywiście f(x₁) ma ten sam znak co rzędna punktu z którego końca wyprowadzona była styczna. J. otrzymane przybliżenie jest za mało dokładne to z punktu o współrzędnych (x₁,f(x₁)) prowadzimy następną styczną. Punkty x₂ przecięcia tej stycznej z osią OX jest drugim przybliżeniem. W ten sposób otrzymujemy kolejne wyrazy ciągu przybliżeń x₁, x₂, ..., x_n... Zbieżnego monotonicznie do α. Na rys.8(6) pokazano interpretację met. Newtona dla przyp. gdy zarówno f'(x)>0 i f''(x)>0. Styczną prowadzimy z punktu B₀. Jej równanie ma postać: y-f(b)=f'(b)(x-b). Przyjmując y=0 otrzymamy x₁=b-f(b)/f'(b) (8). Z tw. Lagrange'a mamy f(x₁)-f(α)=f'(c₁)(x₁-α), przy czym c₁Є(α,x₁). Pon. f(α)=0 i f'(c₁)>0, więc f(x₁)>0, a stąd wynika że styczna poprowadzona w punkcie B₁(x₁,f(x₁)) będzie miała identyczne własności co styczna w punkcie B₀. Z równości drugiej stycznej y-f(x₁)=f'(x₁)(x-x₁) otrzymamy przy x=x₂, y(x₂)=0 x₂=x₁-f(x₁)/f'(x₁) (9). Ogólnie równanie stycznej w punkcie B_n(x_n,f(x_n)), n=0,1,2,... ma postać y-f(x_n)=f'(x_n)(x-x_n) co pozwala na podanie wzoru rekurencyjnego kolejne wyrazy ciągu przybliżeń: x_n+1=x_n=f(x_n)/f'(x_n) (10). Ciąg {x_n} jest ciągiem malejącym (x_n+1<x_n) ograniczonym z dołu (α<x_n) a więc zbieżnym. Przechodząc do granicy w (10) dla n→∞ mamy g=g-f(g)/f'(g) czyli f(g)=0, a więc g=α. Przy rozwiązyw. równ. met. stycznych (Newtona) wygodnie jest skorzystać z TW: J. mamy przedział [a;b] taki że: 1)f(a) i f(b) mają przeciwne znaki; 2)f''(x) jest ciągła i nie zmienia znaku na [a;b]; 3)styczne do krzywej y=f(x) poprowadzono w punktach o odciętych a i b przecinają oś OX wewn. przedziału [a;b], wówczas równanie f(x)=0 ma dokładnie jeden pierw. α w przedziale [a;b] i met. Newtona jest zbieżna do α dla dowolnego punktu startowego x₀Є[a;b]. Znanym przykł. zast. met. Newtona jest algorytm obliczania pierw. kwadratowego x²-c=0. Mamy tutaj f(x)=x²-c i f'(x)=2x, zatem na podst. wzoru (10) x_n+1=x_n-(x_n²-c)/2x_n (x_n≠0) i po uproszczeniu otrzymamy x_n+1=½(x_n+c/x_n) (11). Wszystkie założenia powyższego tw. są spełnione na przedziale 0<a<c^½ i b>½(a+c/a). Dotychczasowe rozważania prowadziliśmy przy założeniu że szukany pierwiastek α jest pojedynczy. Omówimy teraz pokrótce metodę w przypadku gdy α jest pierw. wielokrotnym. DEF. Liczbę α naz. r-krotnym (r≥2) pierw. równania f(x)=0 jest (r-1)-krotnym pierw. równania f'(x)=0. J. α jest pierw. o parzystej krotności to met. połowienia traci sens, natomiast dla krotności nieparzystej (r≥3) gdy f(a)b(b)<0 pozostaje słuszna. Podobnie przedstawia się sprawa dla met. regula falsi oraz met. siecznych, przy czym dla tej ostatniej obniża się rząd metody. Met. Newtona pozwala na znalezienie pierwiastka zarówno o nieparzystej jak i o parzystej krotności, jeżeli tylko istnieje pewne lewo lub prawo stronne sąsiedztwo szukanego pierw. o stałym znaku f'(x) i f''(x). W przyp. wielokrotności jednym ze spos. modyfikacji metody Newtona jest przekształcenie zależności (10) do postaci x_n+1=x_n-rf(x_n)/f'(x_n) (12) gdzie r jest krotnością pierw. (znana a priori). Ponieważ na ogół krotność nie jest znana, często we wzorach (7), (10) stosuje się podstawienie: zamiast f. f(x) wprowadza się f. U(x)=f(x)/f'(x). Równanie U(x)=0 ma identyczne pierwiastki jak równanie f(x)=0, ale wszystkie te pierwiastki są pojedyncze. Wzór (10) przyjmuje wtedy postać x_n+1=x_n-U(x_n)/U'(x_n), gdzie U'(x_n)=1-f''(x_n)U(x_n)/f(x_n) (13). 4.Metoda siecznych. Metodę regiła falsi można znacznie ulepszyć, j. zrezygnujemy z żądania, aby f(x) miała w punktach wytyczających następna cięciwę różne znaki, natomiast do wyznaczenia (n+1)-ego przybliżenia będziemy korzystali zawsze z punktów x_n i x_n-1. Mamy wówczas x_n+1=x_n-f(x_n)(x_n-x_n-1)/f(x_n)-f(x_n-1) (14) Ta metoda (14) naz. się met. siecznych. Jej zbieżność jest znacznie szybsza; niestety zdarzają się przypadki, gdy może ona nie być zbieżna, gdy początkowe przybliżenia nie leżą dostatecznie blisko pierwiastka. Przy stosowaniu met. siecznych jest niezbędne aby pierwsza iteracja zaczynała się od punktów, w których f. ma różne znaki, w przeciwnym razie możemy wykryć nieistniejący pierwiastek RYS.16, co jest szczególnie niebezpieczne przy obliczeniach na maszynach cyfrowych. W met. siecznych istotne znaczenie ma maksymalna graniczna dokładność wynikająca z przyjętej arytmetyki. Gdy bowiem różnica x_n+1-x_n jest tego samego rzędu co oszacowanie błędu następne przybliżenie może już być całkowicie błędne. Dlatego też za dodatkowe kryterium przerwania iteracji należy przyjmować wartości |f(x_n)|, tak aby tworzyły one ciąg malejący (w końcowej fazie obliczeń). Program obliczeniowy powinien także zapewniać możliwość przerwania iteracji, j. różnica między kolejnymi przybliżeniami zamiast maleć zaczyna wzrastać. Wtedy powtórna lokalizacja pierwiastka

Układy równań nieliniowych Metoda Newtona - Raphsona Rozważmy układ równań (X^^(P+1)=X^^(P)-[f^']^-1f^) {f₁(x₁, x₂, ...,x_n)=0; f₂(x₁, x₂, ...,x_n)=0; ...; f_n(x₁, x₂, ...,x_n)=0} (1) albo f^(x^)=0 (2) gdzie x^=[x₁ x₂ ... x_n]^T, f^=[f₁ f₂ ... f_n]^T Niech mamy P-ξ przybliżenie x^^(P)=(x₁^(P), x₂^(P), ..., x_n^(P))^T w pewnym otoczeniu K(α^,r^) punktu α w którym f(α^)=0 wtedy dokładnie pierwiastek można przedstawić w postaci x^=x^^(P)+ξ^^(P) (2) gdzie ξ^^(P)=(ξ₁^(P), ξ₂^(P), ...,ξ_n^(P))^T -błąd pierwiastka f(x^^(P)+ξ^^(P))=0 (3) Niech f. f^(x^) będzie różniczkowalna w otoczeniu K(α^,r^). Załóżmy, że pochodna f^'(x^) jest ciągła w punkcie α^, a pochodna f^'(α) jest nieosobliwa. Więc punkt α jest punktem przyciągania metody iteracyjnej. Interoplacja 1.1 Sformułowanie zagadnienia interpolacji. Na przedziale [a;b] danych jest n+1 różnych punktów x₀, x₁, ...,x_n -które naz. węzłami interpolacji, oraz wartości pewnej f. y=f(x) w tych punktach f(x₀)=y₀, f(x₁)=y₁, ..., f(x_n)=y_n RYS.9 Zadaniem interpolacji (1.1) jest wyznaczenie przybliżonych wartości f. w punktach nie będących węzłami oraz oszacowanie błędu tych przybliżonych wartości. W tym celu należy znaleźć f. F(x) zwaną f. interpolującą, która w węzłach interpolacji przyjmuje takie same wart. co f. y=f(x). Można powiedzieć, że interpolacja jest w pewnym sensie zadaniem odwrotnym do tablicowania f.. Przy tablicowaniu f. budujemy tablicę wart., przy interpolacji natomiast na podstawę tablicy wartości f. określamy jej postać analityczną. F. interpolującej poszukuje się najczęściej jako f. pewnej określonej postaci. Mogą to być wielomiany algebraiczne, wielom. trygonometryczne, funkcje sklejane (splajny) itp. Zależnie od postaci f. interpolującej sformułowane powyżej zagadnienie interpolacyjne może mieć rozwiązania lub mieć nieskończenie wiele rozwiązań. Jest bardzo ważne aby istn. jedna f. interpolująca. Wzory interpolacyjne są punktem wyjścia do wyprowadzenia wielu metod stosowanych w innych działach metod numerycznych np. do różniczkowania numerycznego, kwadratur czy numerycznego rozwiązywania równań różniczkowych. 1.2 Interpolacja za pomocą wielomianów. Poszukamy wielomianu W_n(x) stopnia co najwyżej n, spełniającego warunki (1.1) f(x_i)=y_i (i=0,1,...,n) TW: Istn. dokł. jeden wielom. stopnia co najwyżej n≥0 który w punktach x₀, x₁, ..., x_n przyjmuje wartości y₀, y₁, ..., y_n. Dowód: Przyjmijmy że węzły interpolacji są rozmieszczone w zupełnie dowolny sposób w przedziale [a;b]. Mamy dane n+1 węzłów, w których są znane wartości pewnej f. y=f(x). Szukany wielomian zapisujemy w postaci W_n(x)=a₀+a₁x+...+a_nxⁿ korzystając z warunków (1.1) f(x_i)=y_i otrzymujemy układ (n+1) równań z (n+1) niewiadomymi współczynnikami a₀, a₁, ..., a_n. ∑_[k=0;n]a_kx_i^k=y_i, i=0,1,...,n (1.3) M. współcz. tego ukł. ma postać A=[1 x₀ x₀² ... x₀ⁿ; 1 x₁ x₁² ... x₁ⁿ; ...; 1 x_n x_n² ... x_nⁿ] (1.4) Wyznacznik D=|A|=detA jest wyznacznikiem Vandermonde'a więc przy x_i≠x_j (i≠j) zawsze D=II(x_i-x_j)≠0 (1.5) 0≤j<i≤n. Zatem ukł. (1.3) ma dokł. jedno rozw., a wartości a_k wg tw. Cramera są określone wzorem a_k=1/D∑_[j=0;n]y_kD_kj (1.6) gdzie D_kj (j=0,1,...,n) są kolejnymi dopełnieniami algebraicznymi elementów k-tej kol. m. A. 1.2.1 Wzór interpolacyjny Lagrange'a Podstawiając (1.6) do (1.2) i grupując razem te składniki, w których wyst. jednakowe y_k otrzymujemy W_n(x)=y₀Φ₀(x)+y₁Φ₁(x)+...+y_nΦ_n(x) (1.8) gdzie funkcje Φ_i(x) są wielomianami stopnia co najwyżej n. Jak poprzednio zakładamy że węzły są rozmieszczone w dowolnych sposób. Poszukujemy wielomianu W_n(x) stopnia co najwyżej n określonego równaniem (1.8). Dla każdego i=0,1,2,...,n zachodzi zależność W_n(x_i)=y₀Φ₀(x_i)+y₁Φ₁(x_i)+...+y_nΦ_n(x_i) stąd wynika że Φ_n(x_i)={0, gdy i≠j ; 1, gdy j=i} (1.9). Aby określić f. Φ_j(x) należy znaleźć wielomian stopnia n, który w punktach x₀, ..., x_j-1, x_j+1, ..., x_n jest tożsamościowo=0, a w punkcie x_j=1. Stąd Φ_j(x)=λ(x-x₀) ...(x-x_j-1)(x-x_j+1) ...(x-x_n) (1.11) a pon. Φ_j(x_j)=1 to 1=λ(x_j-x₀) ...(x_j-x_j-1)(x_j-x_j+1) ...(x_j-x_n) (1.12). Po wygórowaniu stałej λ ze wzorów (1.11) - (1.12) otrzymujemy Φ_j(x)=(x-x₀) ...(x-x_j-1)(x-x_j+1) ...(x-x_n) / (x_j-x₀) ...(x_j-x_j-1)(x_j-x_j+1) ...(x_j-x_n) (1.13) Zatem W_n(x)=y₀ (x-x₁)...(x-x_n) / (x₀-x₁)...(x₀-x_n)+y₁ (x-x₀)(x-x₂)...(x-x_n) / (x₁-x₀)(x₁-x₂)...(x₁-x_n)+...+y_n (x-x₀)...(x-x_n-1) / (x_n-x₀)...(x₀-x_n-1)=∑_[j=0;n]y_j (x-x₀) ...(x-x_j-1)(x-x_j+1) ...(x-x_n) / (x_j-x₀) ...(x_j-x_j-1)(x_j-x_j+1) ...(xj-x_n) (1.14) Przyjmując oznaczenie W_n(x)=(x-x₀)...(x-x_n) możemy (1.14) zapisać w postaci W_n(x)=∑_[j=0;n]y_jW_n(x)/(x-x_j)W_n'(x_j) (1.1.6) gdzie y_j=y(x_j), a W_n'(x_j) jest wartością pochodnej wielomianu W_n(x) w punkcie x_j (będącym zerem tego wielomianu) co można łatwo pokazać rozpatrując wzór Taylora w otoczeniu punktu x_j. Otrzymamy wielomian spełniający wszystkie warunki wymagane przy interpolacji. Wzór (1.14) nazywamy wzorem nterpolacyjnym Lagrange'a. Np.: Znaleźć wielomian interpolacyjny który w punktach -2, 1, 2, 4 przyjmuje wartości 3, 1,-3, 8. Stosujemy wzór interpolacyjny Lagrange'a (1.14) dla n=3 Mamy W₃(x)=3(x-1)(x-2)(x-4) / (-2-1)(-2-2)(-2-4)+1(x+2)(x-2)(x-4) / (1+2)(1-2)(1-4) - 3(x+2)(x-1)(x-4) / (2+2)(2-1)(2-4)+8(x+2)(x-2)(x-2) / (4+2)(4-1)(4-2)=2x³/3 - 3x²/2 - 25x/6+6. 1.2.2. Oszacowanie błędu wzoru interpolacyjnego. Powstaje pytanie z jaką dokładnością wielomian Lagrange'a przybliża f. f(x) w pozostałych punktach x≠x_i, xЄ[a;b] czyli określenie wielkości błędu ξ(x)=f(x)-W_n(x)=? (1.17). Niech f(x) ma pochodne rzędu n+1 włącznie. Wprowadzamy f. pomocniczą Ψ(u)=f(u)-W_n(u)+K(u-x₀)(u-x₁) ...(u-x_n) (1.18) w której K jest pewną stałą .Oczywiście Ψ(x₀)=Ψ(x₁)=...=Ψ(x_n)=0. Współczynnik K dobieramy w taki spos. aby pierwiastkiem f. Ψ(u) był również punkt x^--, różny od węzłów interpolacji. Stąd K=f(x^--)-W_n(x^--) / (x^---x₀)...(x^---x_n)=f(x^--)-W_n(x^--)/W_n(x^--) (1.19) Mianownik (1.19) jest ≠0 dla każdego x^--≠x_i (x₀<x₁<...<x_i<x_n). F. Ψ(u) ma więc n+2 miejsc zerowych w przedziale [a;b], w punktach x_o, x₁, ..., x_n, x^--. Na podstawie tw. Rolle'a możemy stwierdzić, że pochodna Ψ'(u) ma co najmniej jedno miejsce zerowe w każdym z przedziałów których końcami są miejsca zerowe f. Ψ(u), a zatem wewnątrz przedziału o końcach min(x;x₀), max(x;x_n) ma co najmniej n+1miejsc zerowych: ζ₀⁽¹⁾, ζ₁⁽¹⁾, ..., ζ_n⁽¹⁾. Stosując ponownie tw. Rolle'a do pochodnej Ψ'(u), stwierdzamy że istn. con. n punktów ζ₀⁽²⁾, ζ₁⁽²⁾, ..., ζ_n-1⁽²⁾ takich że Ψ''(ζ₀⁽²⁾)=...=Ψ''(ζ_n-1⁽²⁾)=0 (1.20). Kontynuując to rozumowanie dochodzimy do wniosku że istn. con. jeden taki punkt w przedziale (min(x;x₀); max(x;x_n)) że Ψ⁽ⁿ⁺¹⁾(ζ)=0 (1.21) Ponieważ W_nⁿ⁺¹(x)=0, (1.22) W_nⁿ⁺¹(x)=(n+1)!. Więc Ψ⁽ⁿ⁺¹⁾(u)=f⁽ⁿ⁺¹⁾(u)-K(n+1)! (1.23). Z (1.2.1) dla u=ζ otrzymujemy K=f⁽ⁿ⁺¹⁾(ζ)/(n+1)! Stąd f(x)-W_n(x)=f⁽ⁿ⁺¹⁾(ζ)W_n(x) / (n+1)! (1.24) J. kres górny M_n+1=max_[xЄ[a;b]]|f⁽ⁿ⁺¹⁾(x)|, to (1.25) |f(x)-W_n(x)|≤M_n+1|W_n(x)| / (n+1)! (1.26). Wyrażenie (1.26) może służyć do oszacowania błędu bezwzględnego wzoru interpolacyjnego. W₄(x)=1+1(x-0)-2/3(x-0)+2/3(x-0)(x-2)(x-3)-2/9(x-0)(x-2)(x-3)(x-4)=-2x⁴/9+8x³/3-88x²/9+35x/3+1. J. tylko możma oszacować wartość fⁿ⁺¹(ζ). Np.: Z jaką dokładnością można obliczyć wartość sin(π/36) stosując wzór interpolacyjny Lagrange'a, j. dane są wartości sin0, sin(π/6), sin(π/4) i sin(π/3). Mamy f(x)=sinx, n=3, a=0, b=π/3, f⁽⁴⁾(x)=sinx, M₄=√3/2|sin(π/36)-W₄(π/36)|≤√3/(2*4!)|(π/36-0)(π/36-π/4)(π/36-π3)|≈9*10^-4 1.2.3.Wzór interpolacyjny Newtona dla nierównych odstępów argumentu. Przyjmijmy, że f. f(x) jest określona za pom. wielom. algebraicznych. Do tego celu niezbędne jest wprowadzenie pojęcia ilorazów różnicowych. Wyrażenia f(x₀,x₁)=[f(x₁)-f(x₀)]/(x₁-x₀); f(x₁,x₂)=[f(x₂)-f(x₁)]/(x₂-x₁); f(x_n-1,x_n)=[f(x_n)-f(x_n-1)]/(x_n-x_n-1) nazywamy ilorazami różnicowymi pierwszego rzędu. Analogicznie definiujemy ilorazy różnicowe rzędu drugiego f(x₀,x₁,x₂)=[f(x₁,x₂)-f(x₀,x₁)]/(x₂-x₀); f(x_n-2,x_n-1,x_n)=[f(x_n-1,x_n)-f(x_n-2,x_n-1)]/(x_n-x_n-2). Ogólnie iloraz różnicowy rzędu n tworzymy z ilorazu różnicowego rzędu n-1 za pomocą wzoru rekurencyjnego f(x_i, x_i+1, ..., x_i+n)=[f(x_i+1, ..., x_i+n)-f(x_i, ..., x_1+n-1)]/(x_i+n-x_i) dla n=1,2,... oraz i=0,1,... Z ilorazów różnicowych tworzy się zazwyczaj tablicę

Ilorazy różnicowe
xi	f(xi)	rzędu 1	rzędu 2	rzędu 3	rzędu 4	rzędu 5
x0 x1 x2 x3 x4 x5	f(x0) f(x1) f(x2) f(x3) f(x4) f(x5)	f(x0,x1) f(x1,x2) f(x2,x3) f(x3,x4) f(x4,x5)	f(x0,x1,x2) f(x1,x2,x3) f(x2,x3,x4) f(x3,x4,x5)	f(x0,x1,x2,x3) f(x1,x2,x3,x4) f(x2,x3,x4,x5)	f(x0,x1,x2,x3,x4) f(x1,x2,x3,x4,x5)	f(x0,x1,x2,x3,x4,x5)

Załóżmy, że f. f(x) jest okr. tablicą wartości (x₀,f(x₀)), (x₁,f(x₁)), ..., (x_n,f(x_n)), gdzie x₀,...,x_n są węzłami interpolacji. Niech W_n(x) będzie wielom. interpolac. W_n(x_i)=f(x_i), I=0,1,…,n W_n(x)=W₀(x)+[W₁(x)-W₀(x)]+[W₂(x)-W₁(x)]+...+[W_n(x)-W_n-1(x)] (1.27). Oczywiście różnica W_k(x)-W_k-1(x) jest wielom., który można przedstawić w postaci W_k(x)-W_k-1(x)=A(x_k-x₀)(x_k-x₁)...(x_k-x_k-1) a ponieważ W_k(x_k)=f(x_k), więc A=[f(x_k)/(x_k-x₀)...(x_k-x_k-1)]-[(∑_[i=0;k-1]f(x_i)W_k-1(x_k)/(x_k-x_i)W_k-1'(x_i))/(x_k-x₀)...(x_k-x_k-1)] (1.28) ⇒ A=f(x_k)/W_k-1(x_k)+∑_[i=0;k-1]f(x_i)/(x_k-x_i)W_k-1'(x_i)=∑_[i=0;k-1]f(x_i)/W_k(x_i) (1.29). Z (1.27) i (1.29) mamy W_k(x)-W_k-1(x)=(∑_[i=0;k-1]f(x_i)/W_k(x_i))W_k-1(x) (1.30). Można udowodnić zależność dotyczącą ilorazu rzędu n (met. indukcji) f(x_i, x_i+1, ..., x_i+n)=f(x_i)/(x_i-x_i+1)...(x_i-x_n+i)+f(x_i+1)/(x_i+1-x_i)...(x_i+1-x_i+n)+...+f(x_i+n)/(x_i+n-x_i)...(x_i+n-x_i+n) (1.31). Przy i=0 z (1.31) mamy f(x₀,x₁,...,x_n)=∑_[j=0;j=n]f(x_j)/(x_j-x₀)...(x_j-x_j-1)(x_j-x_j+1)...(x_j-x_n) dla n=1,2,... i stąd W_n(x)=f(x₀)+f(x₀,x₁)W₀(x)+f(x₀,x₁,x₂)W₁(x)+...+f(x₀,x₁,...,x_n)W_n-1(x) (1.32). Wzór interpolacyjny (1.32) nosi nazwę wzoru interpolacyjnego Newtona dla nierównych odstępów argumentów albo wzoru interpolacyjnego Newtona z ilorazami różnicowymi. Np.: Napisać wzór interpolacyjny dla f(x) okereślony tablicą wartości f(0)=1, f(2)=3, f(3)=2, f(4)=5, f(6)=7

xi	f(xi)	1	2	3	4
0 2 3 4 6	1 3 2 5 7	1 -1 3 1	-2/3 2 -2/3	2/3 -2/3	-2/9

1.2.5 Różnice progresywne i wsteczne. Obecnie przyjmujemy że wartości argumentów x₀, x₁, ..., x_n (węzły interpolacji) tworzą ciąg arytmetyczny i że przyrost argumentu jest wielkością stałą x₁=x₀+h, ..., x_n=x₀+nh. Takie założenie bardzo upraszcza wzory. Różnice progresywne. Niech będą dane wartości f. f(x) w punktach x₀, x₁, ..., x_n. Wyrażenie Δf=f(x+h)=f(x) (1.33) nazywamy różnicą progresywną rzędu pierwszego f. f(x). Podobnie definiujemy różnice progresywne wyższych rzędów Δⁿf=Δ(Δ^n-1f), n=2,3,... (1.34). Np.: Δ²f=Δ(Δf)=Δ[f(x+h)-f(x)]=[f(x+2h)-f(x+h)]-[f(x+h)-f(x)]=f(x+2h)-2f(x+h)+f(x). Wzór (1.34) można także zapisać w postaci Δ⁰f(x)=f(x), Δⁿf(x)=Δ^n-1f(x+h)-Δ^n-1f(x) (1.3.5). Różnicę progresywną dowolnego rzędu f. f można więc przedstawić jako kombinację liniową wartości tej f.. W przypadku wielomianu rzędu n różnica rzędu pierwszego jest wielomianem stopnia n-1. Np.: Obliczyć różnice progresywne wielomianu W₄(x)=x⁴-x-1 przyjmując h=1. Mamy Δw=[(x+1)⁴-(x+1)-1]-[x⁴-x-1]=4x³+6x²+4x; Δ²w=[4(x+1)³+6(x+1)²+4(x+1)]=[4x³+6x²+4x]=12x²+24x+14; Δ³w=[12(x+1)²+24(x+1)+14]-[12x²+24x+14]=24x+36; Δ⁴w=[24(x+1)+36]-[24x+36]=24; Δ⁵w=0. Jak łatwo wykazać, w przypadku wielomianu n-tego stopnia n-ta różnica progresywna jest stałą i ma postać ΔⁿW_n(x)=n!a₀hⁿ=const. Stosując interpolację f. f(x) mamy określoną za pomocą stablicowanych wartości y_if(x_i) a punkty x_i (i=0,1,...,n) są rozmieszczone w równych odstępach. Wówczas Δy_i=y_i+1-y_i; Δ²y_i=Δ(Δy_i)=Δy_i+1-Δy_i; Δⁿy_i=Δ(Δ^n-1y_i)=Δ^n-1y_i+1-Δ^n-1y_i (1.36). Można udowodnić (met. Indukcji) że Δⁿy_i=(ⁿ₀)y_i+n-(ⁿ₁)y_i+n-1+(ⁿ₂)y_i+n-2-...-(-1)ⁿ(ⁿ_n-1)y_i+1+(-1)ⁿ(ⁿ_n)y_i albo Δⁿy_i=∑_[j=0;n](-1)^j(ⁿ_j)y_i+n-j (1.37). Analogicznie dowodzi się że każdą wart. y_k (k=1,2,...,n) można obliczyć za pom. różnic Δy₀, Δ²y₀, ..., Δ^ky₀ oraz wartości y₀: y_k=y₀+(^k₁)Δy₀+(^k₂)Δ²y₀+...+(^k_k)Δ^ky₀ (1.38). Operator Δ ma wiele własności podobnych do właściwości operatora różniczkowego: Δ(g_k±f_k)=Δg_k±Δf_k; Δ(af_k)=aΔf_kΔ^m(Δⁿf_k)=Δ^m+nf_k; Δ(f_kg_k)=f_k+1Δg_k+g_k+1Δf_kΔ f_k/g_k=(Δf_kg_k-f_kΔg_k)/g_k+1g_k dla g_k+1g_k≠0 Δc=0 (c-const). Obliczone różnice kolejnych rzędów najwygodniej jest układać w post. trójkątnej tablicy różnic. Np.: y=x³-1, x₀=0 z krokiem 1 dla n=5 RYS.10 Różnice wsteczne. Korzystając z przyjętych oznaczeń, pierwszą różnicę wsteczną f. f(x) w punkcie x_i określamy jako pierwszą różnicę progresywną tej f. w punkcie x_i-1 i oznaczamy ją symbolicznie ∇y_i, ∇¹y_i, ∇f(x_i) albo ∇¹f(x_i). Zatem ∇y_i=∇y_i-1=f(x_i)-f(x_i-1) (1.38) i=1,2,...,n. Ogólnie, k-tą różnicę wsteczną f. f(x) w punkcie x_i określamy jako różnicę różnic (k-1)-szych w punktach x_i oraz x_i-1, a więc ∇^ky_i=∇^k-1y_i-∇^k-1y_i-1 (1.39) k=1,2,...,n, i=k,k+1,...,n przy czym ∇⁰y_i=y_i. Można bez trudności wykazać, że ∇^ky_i=∇^ky_i-k (k=0,1,...,n). Tablica różnic wstecznych jest identyczna z tablicą różnic progresywnych - zmianie ulegają jedynie symbole różnic. Problem, czy korzystać z różnic progresywnych czy z różnic wstecznych jest uzależniony najczęściej od tego, w jakiej części przedziału (x₀;x_n) chcemy badać f. f(x). J. interesuje nas początek przedziału (otoczenie punktu x₀) korzystamy z różnic progresywnych. J. natomiast interesuje nas zachowanie się f. f(x) na otoczeniu końca przedziału to korzystamy z różnic wstecznych i tablicę rozbudowujemy w kierunku malejących wskaźników x_i (x_n,x_n-1,...). W tym drugim przypadku punkty w których jest określona f. f(x) wygodnie jest oznaczać symbolami: x_-n, x_-n+1, ..., x_-1, x₀. sposób sporządzanie tablicy różnic wstecznych: RYS.11 1.2.6 Wzory interpolacyjne Newtona dla równoległych wartości argumentu. Niech będą dane wartości f. y_i=f(x_i), i=0,1,...,n; przy czym punkty x_isą rozmieszczone w jednakowych odstępach x_i=x₀+ih, i=0,1,...,n; h=const. W punkcie (1.2.4) podaliśmy wzór interpolacyjny Newtona przy założeniu dowolnego rozmieszczenia węzłów W_n(x)=f(x₀)+f(x₀;x)ω₀(x)+...+f(x₀;...;x_n)ω_n-1(x) (1.40) gdzie ω_k(x)=(x-x_o)(x-x₁)...(x-x_k). W przypadku węzłów równoległych ilorazy różnicowe kolejnych rzędów wynoszą odpowiednio: f(x₀)=a₀=y₀; f(x₀;x₁)=(y₁-y₀)/h; f(x₀;x₁;x₂)=[(y₂-y₁)/h - (y₁-y₀)/h]/2h=(y₂-2y₁+y₀)/2h²; ...; f(x₀;x₁;...;x_n)=(y_n-ny_n-1+(ⁿ₂)y_n-2-...+(-1)ⁿy₀)/(n!hⁿ) (1.41). korzystając ze wzoru (1.36) na różnice progresywne rzędu n, otrzymujemy a₀=y₀; a₁=Δy₀/h; a₂=Δ²y₀/2h²; ...; a_k=Δ^ky₀/k!h^k (1.42) przy czym 0!=1, Δ⁰y₀=y₀. Podstawiając (1.42) do (1.40) otrzymujemy W_n(x)=y₀+Δy₀(x-x₀)/h+Δ²y₀(x-x₀)(x-x₁)/2!h²+...+Δⁿy₀(x-x₀)(x-x₁)...(x-x_n-1)/n!hⁿ (1.43). Wzór nosi nazwę pierwszego wzoru interpolacyjnego Newtona. Jest to wielomian stopnia nie większego niż n, a ponad to W_n(x_k)=y_k. W praktycznych zastosowaniach znacznie łatwiej jest stos. wielomian Newtona, wprowadzając nową zmienna g=(x-x₀)/h. Ponieważ (x-x₁)/h=[x-(x₀+h)]/h=g-1; x-x₂/h=[x-(x₀+2h)]/h=g-2; (x-x_n-1)/h=g-(n-1) więc wzór (1.34) można sprowadzić do postaci W_n(x)=W_n(x₀+gh)=y₀+gΔy₀/1!+g(g-1)Δ²y₀/2!+...+g(g-1)...(g-n+1)Δⁿy₀/n! (1.44). Wielkość g=(x-x₀/h określa liczbę kroków potrzebnych do przejścia od punktu x₀ do punktu x. Wzór (1.44) wygodnie jest stos. do interpolacji f. na otoczeniu wartości początkowej x₀ (dla g<1). Przechodząc do przedziału x₁<x<x₂ w którym g>1 wygodniej jest przyjąć za x₀ następny węzeł x₁. J. posługujemy się trójkątną tablicą różnic, to do wzoru (1.44) wchodzą różnice „schodzące” w dół po odpowiedniej przekątnej. Dlatego też stosownie tego wzoru jest wygodne dla górnej, początkowej części tablicy, gdzie mamy dostateczną liczbę różnic. Wzór (1.44) nazywany jest tu także wzorem Newtona na interpolację w przód. Pierwszy wzór Newtona jest niewygodny do interpolacji w pobliżu końca tablicy. Stosuje się tu inny wzór. Szukamy wielomianu zapis. w postaci W_n(x)=a₀+a₁(x-x_n)+a₂(x-x_n)(x-x_n-1)+...+a_n(x-x_n)(x-x_n-1)...(x-x₁) (1.45). Współczynniki a_i, jak poprzednio, wyznaczamy ze wzorów na ilorazy różnicowe a_k=Δ^ky_n-k/k!h^k (1.46) i korzystając ze wprowadzonego pojęci różnic wstecznych mamy a_k=∇^ky_n/k!h^k (1.47). Podstawiając (1.47) do (1.45) otrzymujemy W_n(x)=y_n+Δy_n-1(x-x_n)/1!h+Δ²y_n-2(x-x_n)(x-x_n-1)/2!h²+...+Δⁿy₀(x-x_n)...(x-x₁)/n!hⁿ (1.48). Wzór (1.48) naz. drugim wzorem Newtona. Wprowadzamy podstawienie x=x_n+gh. Wówczas W_n(x)=y_n+gΔy_n-1+g(g+1)Δ²y_n-2/2!+...+g(g+1)(g+2)...(g+n-1)Δⁿy₀/n! (1.49). Skorzystajmy teraz ze wzorów (1.47) określających współczynniki a_k za pom. różnic wstecznych a ponadto pon. z różnic wstecznych korzystamy wtedy, gdy interesuje nas zachowanie się f. na końcowej części przedziału, zmienimy numerację punktów x_i pisząc x_i-n zamiast x_i (i=0,1,...,n) x_i=x₀+ih (i=0,-1,...,-n). Wówczas W_n(x)=y₀+∇y₀(x-x₀)/h+∇²y₀(x-x₀)(x-x_-1)/2!h²+...+∇ⁿy₀(x-x₀)(x-x_-1)...(x-x_-n+1)/n!hⁿ (1.50). Punktom x_-n, x_-n+1, ..., x_-1, x₀ odpowiadają punkty x₀, x₁, ..., x_n. Po podstawieniu do (1.50) g=(x₀-x)/h mamy W_n(x)=y₀-g∇y₀+g(g-1)∇²y₀/2-...+(-1)ⁿg(g-1)(g-2)...(g-n+1)∇ⁿy₀/h! (1.51). Wzór naz. się często wzorem interpolacyjnym Newtona na interpolację wstecz. 1.2.7 Zbieżność procesów interpolacyjnych. powstaje pytanie, czy przybliżenie polepszy się gdy zwiększymy liczbę węzłów? Intuicyjnie wydaje się oczywiste że j. będziemy mieli większą liczbę zmierzonych wartości f., to można będzie dokładniej opisać matematycznie tę zależność funkcyjną. Nie zawsze jednak tak jest, zwłaszcza wtedy, gdy do opisu matematycznego stos. interpolację. Pogorszenie wyników interpolacji przy zwiększaniu liczby węzłów pokażemy na przykładzie f. y=|x| na przydziale [-1;1]. Znajdziemy wielomiany interpolacyjne i narysujemy ich wykresy dla n=2,4,10. Mamy dla n=2, x₀=-1, h=1, W₂(x)=x²; dla n=4, x₀=-1, h=0,5, W₄(x)=-4x⁴/3+7x²/3; dla n=10, x₀=-1, h=0,2, W₁₀(x)=390625x¹⁰/5184-1015625x⁸/6048+221875x⁶/1728+6835x⁴/162+11527x²/1792 RYS.12 Na rys. są podane wykresy wielomianów interpolacyjnych. Widać wyraźnie, że początkowo ze wzrostem liczby węzłów przybliżenie polepsza się, lecz przy dalszym wzroście n przybliżenie zaczyna się pogarszać, zwłaszcza przy krańcach przedziału. Wadą f. y=|x| jest fakt że w punkcie x=0 nie ma ona pochodnej. Można jednakże łatwo znaleźć f. nie mającą tej wady, np.: y=1/(1+25x²) w przypadku której mamy do czynienia z tym samym zjawiskiem. Takie zachowanie się wielomianu interpolacyjnego (zwłaszcza przy końcach przedziału) jest zjawiskiem typowym dla interpolacji za pomocą wielomianu wysokich stopni przy stałych odległościach węzłów. Jest to tak zwane zjawisko Kungego - jest ono przykładem źle uwarunkowanego zadania. Interpolacja tego typu na środkowych częściach przedziału [x₀;x_n] jest natomiast dobra i użyteczna. Również w przypadku interpolowania f., której wykres znacznie różni się od wykresu wielomianu (np. f. x=1/x) interpolacja przy dużej liczbie węzłów daje nie najlepsze wyniki (pojawiają się ekstrema w f. interpolującej). W takich przypadkach lepsze wyniki daje interpolacja przedziałowa (np. wielomianami stopnia drugiego). Problem wyboru odpowiedniego wzoru interpolacyjnego ma szczególne znaczenie przy obliczeniach „ręcznych”, gdy nie możemy zbytnio zagęszczać węzłów. Stos. się wówczas również inne wzory interpolacyjne zawierające różnice położone nad i pod poziomym wierszem, w którym powinna znajdować się podana wartość zmiennej x. Są to tzw. wzory z różnicami centralnymi; zaliczamy do nich wzory: Gaussa, Stirlinga i Bessela, których nie będziemy omawiać. 1.3. Interpolacja za pom. f. sklejanych. 1.3.1 Określenie f. sklejanych. Niech w przedziale [a;b] danych będzie (n+1) punktów x₀, x₁, ..., x_n przy czym a=x₀<x₁<...<x_n=b (1). Punkty x_i określają pewien podział przedziału [a;b] na n podprzedziałów. Podział ten oznaczymy Δ_n. F. S(x)=S(x,Δ_n) określoną na przedziale [a;b] naz. f. sklejaną stopnia co najwyżej m (m≥1); j.: 1)S(x) jest wiel. stopnia co najwyżej m na każdym podprzedziale (x_i;x_i+1), i=0,1,...,n-1; 2)S(x)ЄC^(m-1)([a;b]). Punkty x_i naz. węzłami f. sklejanej. Zbiór wszystkich f. sklejanych stopnia m o węzłach x_i oznaczymy S_m(Δ_n). Niech S(x)ЄS_m(Δ_n). Na każdym przedziale (x_i;x_i+1), i=0,1,...,n-1 f. S(x) jest wielom. stop. co najw. m S(x)-C_imx^m+C_im-1x^m-1+...+C_i1x+C_i0 (2) xЄ(x_i;x_i+1). Mamy więc n(m+1) dowolnych stałych C_ij. Żądanie ciągłości pochodnych rzędu 0,1,...,n-1 w każdym wewnętrznym x_i, i=1,...,n-1 daje m(n-1) warunków. Tak więc f. S(x) zależy od n(m+1)-m(n-1) parametrów. W obliczeniach na maszynach cyfrowych funkcje bardzo często przybliża się f.mi sklejanymi. Związane jest to z tym, że wartości tych f. są łatwe do wyznaczenia oraz że są one zbieżne dla licznych klas f. (jest to jedna z zalet f. sklejanych w porównaniu z wielomianami). W praktyce często stos. się funkcje sklejane stopnia trzeciego. Dla wielu zagadnień są one wystarczająco gładkie, szybkość ich zbieżności jest wystarczająca a ponadto wyznacza się je stosunkowo łatwo. Postać interpolującej f. sklejanej w dużym stopniu zależy od zagadnienia. Często wygodnie jest przedstawić poszukiwaną f. S(x)ЄS_m(Δ_n) w postaci kombinacji liniowej elementów bazy przestrzeni. S_m(Δ_n). Wyznaczymy teraz bazę przestrzeni f. S(x) stopnia trzeciego z węzłami równoodległymi x_i=x₀+ih, h=(b-a)/n, i=0,1,...,n. Oznaczmy dodatkowo przez x_-3, x_-2, x_-1, x_n+1, x_n+2, x_n+3 punkty x_i=x₀+ih dla i=-3,-2,-1, n+1, n+2, n+3. i określmy f. ϕ_i³(x), i=-1, 0, 1, ..., n, n+1: ϕ_i³(x)=1/h³ {(x-x_i-2)³ dla xЄ[x_i-2;x_i-1]; h³+3h²(x-x_i-1)+3h(x-x_i-1)²-3(x-x_i-1)³ dla xЄ[x_i-1;x_i]; h³+3h²(x_i+1-x)+3h(x_i+1-x)²-3(x_i+1-x)³ dla xЄ[x_i;x_i+1]; (x_i+2-x)³ dla xЄ[x_i+1;x_i+2]; 0 dla pozostałych xЄR} (3). Funkcje (3) są klasy C²((-∞;+∞)). W tablicy

	xj-2	xj-1	xj	xj+1	xj+2
ϕj^3(x)	0	1	4	1	0
(ϕj^3(x))'	0	3/h	0	-3/h	0
(ϕj^3(x))''	0	6/h^2	-12/h^2	6/h^2	0

podano wartości f. ϕ_j³(x) oraz jej pierwszej i drugiej pochodnej w punktach x_k dla k=j-2, j-1, j, j+1, j+2 (poza przedziałem [x_j-2;x_j+2] f. ta jest tożsamościowo=0), a na RYS.13 pokazano jej wykres. TW: Funkcje ϕ_i^—(x), i=-1,0,1,...,n+1, określone na przedziale [a;b] w nast. spos. ϕ_i^—(x)=ϕ_i³(x) (4), a≤x≤b stanowią bazę przestrzeni f. sklejanych trzeciego stopnia S₃(Δ_n). Każdą f. S(x)ЄS₃(Δ_n) można więc przedstawić w postaci kombinacji liniowej S(x)=∑_[i=-1;n+1]C_iϕ_i³(x), a≤x≤b (5) gdzie C_i są liczbami rzeczywistymi. 1.3.2 Interpolacyjne f. sklejane stopnia 3. Niech f(x)ЄC([a;b]). F. S(x)ЄS₃(Δ_n) naz. interpolacyjną f. sklejaną stopnia 3-go dla f. f(x), j. S(x_i)=f(x_i)=y_i (6), i=0,1,...,n; n≥2. Jak wiemy, f. S(x) stopnia 3-go zależy od (n+3) parametrów. Tak więc f. sklejana stopnia 3-go mają dwa stopnie swobody, wobec czego nakładamy na nie dwa dodatkowe warunki. Wybór tych warunków zależy zarówno od własności f(x), jak i od posiadanej info. o tej f.. Najczęściej rozpatrujemy następujące dodatkowe warunki: S'(a+0)=α, S'(b-0)=β (7) lub S''(a+0)=α₂, S''(b-0)=β₂ (8) gdzie α₁, β₁, α₂, β₂ są ustalonymi liczbami. J. f(x) ma pochodne w punktach a, b i są one znane, to możemy je przyjąć jako α₁, β₁, α₂, β₂. Na f. S(x) możemy nakładać inne warunki. Np.: j. f(x) jest okresowa, o okresie (b-a), to najczęściej żądamy, aby f. S(x) można było przedłużyć na przedziale (-∞;+∞) w taki spos., że będzie ona f. okresową o okresie (b-a), różniczkowalną w spos. ciągły. W takiej sytuacji na f. sklejaną nakładamy warunki S⁽ⁱ⁾(a+0)=S⁽ⁱ⁾(b-0), i=1,2 (9). Omówimy spos. wyznaczania oraz problem istnienia i jednoznaczności interpolacyjnej f. sklejanej dla węzłów dowolnie rozłożonych. Oznaczymy M_j=S''(x_j), j=0,1,...,n (10). Zgodnie z określeniem f. sklejanej 3-go stopnia, pochodna S''(x) jest f. ciągłą na przedziale [a;b] i liniowa na każdym podprzedziale [x_j-1;x_j], można więc przedstawić ją w postaci S''(x)=M_j-1(x_j-x)/h_j+M_j(x-x_j-1)/h_j, xЄ[x_j-1;x_j] (11) h_j=x_j-x_j-1. Całkując stronami (11), otrzymujemy (dla xЄ[x_j-1;x_j]) S'(x)=-M_j-1(x_j-x)²/2h_j+M_j(x-x_j-1)²/2h_j+A_j (12) S(x)=M_j-1(x_j-x)³/6h_j+M_j(x-x_j-1)³/6h_j+a_j(x-x_j-1)+B_j (13) gdzie A_j, B_j są stałymi. Nakładając na (12) i (13) warunki interpolacji S(x_j-1)=M_j-1h_j²/6+B_j=y_j-1; S(x_j)=M_jh_j²/6+A_jh_j+B_j=y_j wyznaczamy stałe A_j, B_j. Mamy B_j=y_j-1-M_j-1h_j²/6, A_j=(y_j-y_j-1)/h_j-h_j(M_j-M_j-1)/6 (14). Nałożymy teraz warunek na S(x), aby pochodna S'(x) była f. ciągłą na [a;b]. W tym celu obliczymy granice jednostronne S'(x_j-0)=h_jM_j-1/6+h_jM_j/3+(y_j-y_j-1)/h_j (15) S'(x_j+0)=-h_j+1M_j/3-h_j+1M_j+1/6+(y_j+1-y_j)/h_j+1 i zażądamy, aby S'(x_j-0)=S'(x_j+0), j=1,...,n-1. Po podst. (15) do (16) otrzymujemy ukł. (n-1) równań: h_jM_j-1/6+(h_j+h_j+1)M_j/3+h_j+1M_j+1/6=(y_j+1-y_j)/h_j+1-(y_j-y_j-1)/h_j (17), j=1,...,n-1. Równania te wygodnie jest zapisać w postaci μ_jM_j-1+2M_j+λ_jM_j+1=d_j, j=1,...,n-1 (18) gdzie: λ_j=h_j+1/(h_j+h_j+1), μ_j=1-λ_j; d_j=[6/(h_j+h_j+1)][(y_j+1-y_j)/h_j+1-(y_j-y_j-1)/h_j]=6f(x_j-1;x_j;x_j+1). J. M_j, j=0,1,...,n spełniają powyższy ukł. to f. S(x) określona wzorami (13), (14) na każdym podprzedziale [x_j-1;x_j], j=1,...,n jest interpolacyjną f. sklejaną stopnia 3-go. Do ukł. (18) należy dołączyć jeszcze dwa równania wynikające ze spełnienia przez f. S(x) jednego z dodatkowych warunków (7-9). Równania te mają postać: dla warunków (7) 2M₀+M₁=d₀, M_n-1+2M_n=d_n gdzie d₀=6/h₁[(y₁-y₀)/h₁-α₁), d_n=6/h_n[β₁-(y_n-y_n-1)/h_n) (19) dla warunków (8) M₀=α₂, M_n=β₂ (20) dla warunków (9) M₀=M_n, λ_nM₁+μ_nM_n-1+2M_n=d^~~_n gdzie λ_n=h₁/(h₁+h_n), μ_n=1-λ_n, d^~~_n=[(y₁-y₀)/h₁-(y_n-y_n-1)/h_n]6/(h₁+h_n) (21). Zapiszemy teraz ukł. równ. (18) i (19) w postaci macierzowej [2 1 0 ... 0; μ₁ 2 λ₁ ... 0; 0 μ₂ 2 ... 0; ...; 0 ... 2 λ_n-1; 0 ... 1 2][M₀ M₁ ... M_n-1 M_n]^T=[d₁ d₂ ... d_n-1 d_n] (22). We wszystkich trzech przypadkach (19) - (21) m. współczynników układów, które należy rozwiązać są m. silnie dominującymi. stąd wynika, że układy te mają jednoznaczne rozwiązanie, istn. zatem dokł. jedna interpolacyjna f. sklejana stopnia trzeciego, spełniająca jeden z dodatkowych warunków (19), (20), (21). Dla układów z macierzami trójdiagonalnymi silnie, diagonalnie dominującymi są opracowane niezawodne algorytmy ich rozwiązywania, a więc wyznaczanie interpolacyjnych f. sklejanych w przedstawiony powyżej spos. nie nastręcza numerycznych trudności. Zauważmy że j. węzły są równoodległe to λ_j=μ_j=½. Aproksymacja 1.Wstęp. W poprzednim rozdziale omówiliśmy zagadnienie interpolacji, czyli poszukiwania pewnej f. interpolującej F(x), która na danym, dyskretnym zbiorze argumentów pokrywa się z wartościami f. interpolowanej f(x). Takie postępow. jest często niedogodne. Zajmiemy się teraz zagadnieniem bardziej ogólnym, w którym warunek aby f. dana i f. szukana przyjmowały dokładnie te same wartości na zb. danych z góry punktów węzłowych, nie musi być spełniony. F. f(x), znaną lub określoną tablicą wartości, będziemy aproksymować (zastępować) inna f. F(x), zwaną f. aproksymującą lub przybliżeniem f. f(x). Oczywiście przybliżenie takie powoduje pojawienie się bł. zwanych bł. aproksymacji (przybliżenia) i problem oszacowania tych błędów oraz ich wielkości mają istotny wpływ na wybór metody aproksymacyjnej. Gdy zbiór, na którym mierzymy błąd aproksymacji, jest zbiorem dyskretnym, wówczas jest to apr. punktowa, gdy zaś jest to przedział - apr. integralna. Niech f(x) będzie f. którą chcemy aproksymować, X - prwną przestrzenią liniową unormowaną (skończenie lub nieskończenie unormowaną; f(x)ЄX) a X_m - m-wymiarową podprzestrzenią liniową p. X. Apr. f. f(x) polega na wyznaczeniu takich współczynników a₀, a₁, ..., a_m f. F(x)=a₀ϕ₀(x)+a₁ϕ₁(x)+...+a_mϕ_m(x) (1) gdzie ϕ₀,ϕ₁,...,ϕ_m są funkcjami bazowymi m+1 wymiarowej podprz. liniowej X_m+1 aby f. F(x) spełniała pewne warunki, np. minimalizowała normę różnicy ||f(x)-F(x)||. Apr. określona wzorem (1) nosi nazwę apr. liniowej (wzór (1) naz. się też wilom. uogólnionym). W obliczeniach na maszynach cyfrowych duże znaczenie ma też apr. wymierna, określona wzorem: F(x)=[a₀ϕ₀(x)+a₁ϕ₁(x)+...+a_mϕ_m(x)]/[b₀ψ₀(x)+b₁ψ₁(x)+...+b_mψ_m(x)] (2) gdzie ϕ_i(x), ψ_j(x) (i=0,1,...,n) (j=0,1,...,m) są elementami tej samej bazy k-wymiarowej podprz. linowej (k=max(n,m)), natomiast a_i i b_i są stałymi współczynnikami, które należy wyznaczyć. Zajmiemy się teraz zagadnieniem wyboru odpow. podprz. X_m i związanej z nią bazy (f. bazowych ϕ_k(x)). J. f. aproksymowana f(x) jest f. cigłą na przedziale [a;b] f(x)ЄC[a;b], to f. ϕ_k(x) będą elementami pewnej m+1 wymiarowej podprz. C[a;b]. Jedną z często obieranych podprz. X_m jest podprz. f. trygonometrycznych z bazą 1, sinx, cosx, sin2x, cos2x, ..., sinkx, coskx szczególnie przydatna, gdy apr. f. f(x) jest f. okresową. Inną, też często obieraną jest podprz. wilomianów (4) 1, x, x², ..., x^m - z bazą jednomianów. Interpolacja jest oczywiście także jedną z met. apr. f. Przy interp. wzorami Newtona używaliśmy jako bazy wielomianów ϕ₀(x)=1, ϕ₁(x)=(x-x₀), ..., ϕ_m(x)=(x-x₀)(x-x₁)...(x-x_m-1). Przyjmijmy, że w p. X jest określona norma ||•||. Mówimy, że f. F)x) dobrze przybliża f. f(x), j. norma ||F9x)-f(x)|| jest mała. Normą taką może być np.: norma Czebyszewa ||f||=max_[a;b]|f(x)|; norma L₂ ||f||₂=(∫_[a;b]|f(x)|²dx)^½; norma L₂ z wagą ||f||_2,w=(∫_[a;b]w(x)|f(x)|²dx)^½ gdzie x(x) jest ciągłą, nieujemną f. wagową, dodatnią poza zb. miary zero. Gdy f. f(x) jest określona na dyskretnym zb. punktów wówczas można dane wartości f(x_i) przyjąć za współrzędne wekt. f i wówczas rozpatrujemy normę ||f||=(∑_[i=0;n][f(x_i)]²)^½ Zagadnienie najlepszej apr. przy wybranych f. bazowych ϕ_k(x) sprowadza się do znalezienia wartości współczynników a_k takich, aby otrzymać wyrażenia ||f(x)-(a₀ϕ₀(x)+...+a_mϕ_m(x))|| (5) i aby istn. jedyne możliwe rozw. tego zagadnienia ze wzgl. na a_k. Apr. średniokwadratowa. Dla f. f(x) okr. na przedziale [a;b] poszukujemy minimum całki ||F(x)-f(x)||=∫_[a;b]W(x)[F(x)-f(x)]²dx (6) a dla f. f(x) danej na dyskretnym zb. argumentów poszukujemy min. sumy (met. najmn. kwadratów) ||F(x)-f(x)||=∑_[i=0;n]W(x_i)[F(x_i)-f(x_i)]² (7) W(x_i)≥0 dla i=0,1,...,n. Apr. jednostajna. Dla f. f(x) poszukujemy f. F(x) dającej najmn. maksimum różnicy między F(x) a f(x) na całym przedziale ||F(x)-f(x)||=max_[xЄ[a;b]]|F(x)-f(x)| (8) TW1: J. f. f(x) jest ciągła na [a;b] (skończonym) to dla każdego ξ>0 można dobrać takie n, że jest możliwe utworzenie wielom. P_n(x) stopnia n (n=n(ξ)), który spełnia nierówność |f(x)-P_n(x)|<ξ (9) na całym [a;b]. TW2: J. f. f(x) jest ciągła na R i okresową o okresie 2π, to dla każdego ξ dodatniego istn. wielom. trygonometryczny S_n(x)=a₀+∑_[k=1;n](a_kcoskx+b_ksinkx) (10), n=n(ξ) spełniający dla wszystkich x nierówność |f(x)-S_n(x)|<ξ. 2.Apr. średniokwadratowa. Niech będzie dana f. y=f(x), która na pewnym zb. X punktów x₀, x₁, ..., x_n przyjmuje wart. y₀, y₁, ..., y_n. Wartości te możemy znać tylko w przybliżeniu z pewnymi błędami (np.: jako wyniki pomiarów obarczone błędami obserwacji) przy czym błędy te mogą być niejednokrotnie dosyć duże, co w istotny spos. wpływa na jakość apr. Będziemy poszukiwać takiej f. F(x) przybliżającej daną f. f(x), która umożliwi wygładzenie f. f(x) i pozwoli z zakłóconych błędami danych wart. f. przybliżanej otrzymać gładką f. przybliżającą z dużym prawdopodobieństwem mało odchylającą się od f. przybliżanej zarówno między węzłami, jak i w węzłach x₀, x₁, ..., x_n, j. tylko przyjmiemy, że f. przybliżana ma dość gładki przebieg. Niech ϕ_i(x), j=0,1,...,m. będzie układem f. bazowych podprz. X. Poszukujemy wielom. uogólnionego F(x), RYS.14 będącego najlepszym przybliżeniem średniokwadratowym f. f(x) na zb. X=(x_j), tj. f. postaci F(x)=∑_[i=0;n]a_iϕ_i(x) (11) przy czym współczynniki a_i są tak określone aby wyrażenie (7) było min. Oznaczmy H(a₀, a₁, ..., a_m)=∑_[j=0;n]W(x_j)[f(x_j)-∑a_iϕ_i(x_j)]²=∑_[j=0;n]W_jR_j² (12) gdzie f. wagowa W(x) jest ustalona z góry i taka że W(x_i)>0 dla i=0,1,...,n, a R_j jest odchyleniem w punkcie x_j. Aby znaleźć a_i obliczamy pochodne cząstkowe f. H względem zmiennych a_i. Z warunku ekstremów f. wielu zmiennych ∂H/∂a_k=0 (13), k=0,1,...,m otrzymamy wówczas układ (m+1) równań liniowych z m+1 niewiadomymi a_i ∂H/∂a_k=-2∑_[j=0;n]W(x_j)[f(x_j)-∑_[i=0;m]a_iϕ_i(x)]ϕ_k(x_j)=0 (14) zwany ukł. normalnym. Ponieważ f. ϕ_j(x) tworzą bazę p. X_m, więc układ (14) ma wyznacznik różny od zera i rozwiązanie tego ukł. daje min. f. H. W zapisie macierzowym ukł. (14) przyjmuje postać D^TDA=D^Tf (15), gdzie D=[ϕ₀(x₀) ... ϕ_m(x₀); ϕ₀(x₁) ... ϕ_m(x₁); ...; ϕ₀(x_n) ... ϕ_m(x_n)], A^=[a₀ a₁ … a_m]^T, f^=[f(x₀) f(x₁) … f(x_n)]^T. M. współczynników ukł. (15) jest m. symetryczną i dodatnio określona co zapewnia jednoznaczność rozwiązania. Ukł. normalny (14) lub (15) można otrzymać z równania (11) po podst. wart. węzłowych x_i (i=0,1,...,n). Otrzymamy wówczas nieokreślony ukł. (n+1) równań z (m+1) niewiadomymi DA=f (16) z którego po pomnożeniu przez D^T otrzymujemy (15). 2.1 Apr. wielomianowa. J. jako f. bazowe przyjmiemy ciąg jednomianów (xⁱ), i=0,1,...,m to wzór (15) przyjmie postać ∑_[j=0;n][f(x_j)-∑a_ix_jⁱ]x_j^k=0 (17) k=0,1,...,m. Zmieniając kolejność sumowania mamy ∑_[j=0;n]f(x_j)x_j^k=∑_[i=0;m]a_i(∑_[j=0;n]x_j^i+k), (18) k=0,1,...,m i oznaczając g_ik=∑_[j=0;n]x_j^i+k, ρ_k=∑_[j=0;n]f(x_j)x_j^k (13) otrzymamy ukł. normalny postaci ∑_[i=0;m]a_ig_ik=ρ_k (20) k=0,1,...,m. Można wykazać, że j. punkty x₀, x₁, ..., x_n są różne i m≤n to wyznacznik ukł. (20) jest różny od zera, a więc ukł. tem ma jednoznaczne rozwiązanie. J. m=n, to wielom. aproksymac. f(x) pokrywa się z wielom. interpolac. dla punktów x₀, x₁, ..., x_n i wówczas H=0. W praktyce stopień wieom. jest i powinien być znacznie mniejszy od liczby punktów x_k (k=0,1,...,n), wtedy bowiem korzystamy z dużej ilości inf., uzyskując równocześnie prostrze (niskiego stopnia) f. aproksymujące. Np.: Niech dane będą wyniki pomiarów zestawione w tablicy

x	1	3	4	6	8	9	11	14
y	1	2	4	4	5	7	8	9

Poszukujemy zależności między x i y postaci ax+by=1 przy czym zadanie polega na znalezieniu najlepszych a i b. RYS.15 Przy tak sformowanym zad. pojawia się problem: czy rozpatrywać odchylenia x czy y. Przedyskutujmy obydwa przypadki pokazując jednocześnie że prowadzą one do różnych wniosków. Gdy przyjmiemy, że wartości y są obarczone pewnymi bł. wówczas minimalizujemy wyrażenie ∑_[i=1;n](y_i-y_i^--)² (21) gdzie y_i^--=-ax_i/b+1/b=a₁x_i+a₀. J. natomiast przyjmiemy że wartości x są obarczone bł. to powinniśmy minimalizować wyrażenie ∑_[i=1;n](x_i-x_i^--)² (22) a₁x_i^--+by_i=1. W przypadku minimaliz. wyrażenia (21) otrzymamy układ normalny postaci {na₀+(∑_[i=1;n]x_i)a₁=∑_[i=1;n]y_i (n=8); (∑_[i=1;n]x_i)a₀+(∑_[i=1;n]x_i²)a₁=∑_[i=1;n]x_iy_i} w którym a₁=-a/b, a₀=1/b i stąd {a₀=[(∑_[1;n]y_i)(∑_[1;n]x_i²)-(∑_[1;n]x_i)(∑_[1;n]x_iy_i)]/[n(∑_[1;n]x_i²)-(∑_[1;n]x_i)²; a₁=[n(∑_[1;n]x_iy_i)-(∑_[1;n]x_i)(∑_[1;n]y_i)]/n∑_[1;n]x_i²-(∑_[1;n]x_i)²} (24). Układ normalny (24) dla danych z tablicy ma rozwiązanie a₀=6/11 i a₁=7/11 i stąd równanie ma postać 11y-7x=6 (25). Minimalizując w taki spos. wyrażenie (22) otrzymamy równanie 2x-3y=-1 (26). Jak widać, równania (25) i (26) są różne, chociaż ich wykresy prawie się pokrywają. 3.Apr. za pomocą wielomianów ortogonalnych. Przy apr. wielom. z f. bazowymi 1, x, x², ..., x^k ze wzrostem stopnia wielom. k obliczenia stają się coraz bardziej pracochłonne, a ponadto ich wyniki są niepewne. Dodatkową trudnością jest fakt że zmiana stopnia wielom. wymaga ponownego rozwiązania ukł. normalnego, co też przemawia przeciwko stosowaniu ciągu funkcji ϕ_i(x)=xⁱ do apr. Obydwie trudności można usunąć używając do apr. wielomianów ortogonalnych. Zdefiniujmy zatem ortogonalność f. na dyskretnym zb. punktów. DEF. F. f(x) i g(x) naz. ortogonalnymi na dyskr. zb. punktów x₀, x₁, ..., x_n j. ∑_[i=0;n]f(x_i)g(x_j)=0 (1) przy czym ∑_[i=0;n][f(x_i)]²>0, ∑_[i=0;n][g(x_i)]²>0 (2). Analogicznie ciąg funkcyjny {ϕ_m(x)}≡ϕ₀(x), ϕ₁(x), ..., ϕ_m(x), ... (3) naz. ortogonalnym na zb. punktów x₀, x₁, ..., x_n (4), j. ∑_[i=0;n]ϕ_j(x_i)ϕ_k(x_i)=0 (5) dla j≠k. Zakładamy, że nie wszystkie punkty x_i są miejscami zerowymi f. ϕ_j(x). Zastosowanie wielom. ortogonalnych usuwa jedną z podst. trudności obliczeniowych. Przy apr. wielom. ortogonalnymi m. ukł. normalnego jest m. diagonalna, a jej elementy położone na głównej przekątnej wynoszą a_jj=∑_[i=0;n]ϕ_j²(x_i). Załóżmy że mamy n+1 równoodległych punktów x_i (x_i=x₀+ih), i=0,1,...,n). Niech q=(x-x₀)/h przeprowadzimy te punkty odpowiednio w kolejne liczby całkowite od 0 do n. x₀→0; x₁→q=(x₁-x₀)/h=1; x_n→(x_n-x₀)/h=n=q Poszukujemy ciągu wielomianów {F_i⁽ⁿ⁾(q)}: F₀⁽ⁿ⁾(q), F₁⁽ⁿ⁾(q), ..., F_m⁽ⁿ⁾(q) (6) (dolny indeks oznacza stopień, m≤n) spełniających na zb. punktów q=0,1,...,n warunek ortogonalności ∑_[i=0;n]F_j⁽ⁿ⁾(i)F_k⁽ⁿ⁾(i)=0 (7) dla j≠k przy czym F_kⁿ(q)=a₀+a₁q+a₂q(q-1)+...+a_kq(q-1)...(q-k+1) (8). Wówczas F_k⁽ⁿ⁾(q)=a₀+a₁q⁽¹⁾+a₂q⁽²⁾+...+a_kq^(k). Często używa się unormowanego ciągu wielom. spełniających warunek F^—_k⁽ⁿ⁾(q)=∑_[s=0;k](-1)^s(^k_s)(^k+s_s) q^[s]/n^[n] (10) k=0,1,...,m gdzie (^k_s) są symbolami Newtona. Wzór (10) określa wielom. Grama stopni k=0,1,...,m dla (n+1) równoodległych węzłów. Wzór apr. Grama ma postać: y_m(x)=∑_[k=0;m]C_kF^—_k⁽ⁿ⁾(q)/S_k=∑_[k=0;m]C_kF^—_k⁽ⁿ⁾((x-x₀)/h)/S_k (11) S_k=∑_[q=0;n][F^—_k⁽ⁿ⁾(q)]², C_k=∑_[i=0;n]y_iF^—_k⁽ⁿ⁾(x_i) dla m≤n. Np.:Dana jest tabl. wartości f(x). Znaleźć najlepszy, w sensie aproksymacji średniokwadratowej, wielom. przybliżający f. f(x). Tabl.1

x	1	1,5	2	2,5	3
y	3	4,75	7	9,75	13

Ponieważ węzły są równoodległe, posłużymy się wielomianami Grama. Konstruujemy wielom. apr. dla n=4. Wartości F_k(q), q=0,1,2,3,4 obliczamy ze wzoru (11) F₀(q)=1, F₁(q)=1-2q/n, n≥1, F₂(q)=1-6q/n+6q(q-1)/n(n-1), n≥2, F₃(q)=1-12q/n+30q(q-1)/n(n-1)-20q(q-1)(q-2)/n(n-1)(n-2), n≥3, F₄(q)=1-20q/n+90q(q-1)/n(n-1)-140q(q-1)(q-2)/n(n-1)(n-2)+70q(q-1)(q-2)(q-3)/n(n-1)(n-2)(n-3), n≥4 i zestawiamy w postaci tabl. 2

q	xi	yi	F0(q)	F1(q)	F2(q)	F3(q)	F4(q)
0	1	3	1	1	1	1	1
1	1,5	4,75	1	0,5	-0,5	-2	-4
2	2	7	1	0	-1	0	6
3	2,5	9,75	1	-0,5	-0,5	2	-4
4	3	13	1	-1	1	-1	1

Obliczamy wartości C_i, S_i oraz b_i=C_i/S_i (tabl. 3)

Ci	Si	bi
37,5	5	7,5
-12,5	2,5	-5
1,75	3,5	0,5
0	10	0
0	70	0

Dla q=(x-1)/0,5 i ze wzoru (11) otrzymujemy y(x)=∑_[i=0;4]b_iF_i(q)=x²+x+1. Apr. średniokwadratowa f. ciągłych. Niech f(x) ciągła na przedziale [a;b]. Wielom. apr. P(x)=a₀ϕ₀(x)+a₁ϕ₁(x)+…+a_nϕ_n(x) gdzie {ϕ_I(x)} są elementami bazy pewnej podprzestrzeni funkcji całkowych z kwadratem na [a;b]. Podobnie jak w przyp., gdy f(x) była określona na dyskretnym zb. punktów, apr. średniokw. f. ciągłych polega na znalezieniu takiego ciągu współczynników a_i (i=0,1,...,n) aby otrzymać minimum normy ||F(x)-f(x)||=∫_[a;b]W(x)[F(x)-f(x)]²dx. Przyjmijmy zatem, tak jak w przyp. apr. dyskretnej, że W(x)≡1 i H_n=∫_[a;b]{P(x)-f(x)]²dx=∫_[a;b][∑_[i=0;n]a_iϕ_i(x)-f(x)]²dx (12). Aby rozwiązać zadanie obl. pochodne cząstkowe ∂H_n/∂a_i, i=0,1,...,n, przyrównujemy je do zera i z otrzymanego ukł. równ. wyznaczamy współczynniki a_i. Rachunki będą znacznie prostrze, j. zamiast ϕ_i(x)-xⁱ użyjemy ciągu f. ortogonalnych z wagą. Może to być np. ciąg wielomianów Legendre'a P_n(x)=(1/2ⁿn!)(dⁿ/dxⁿ)(x²-1)ⁿ (13), n=0,1,... Przykładem ortogonalnego układu f. z wagą na przedziale [-1;1], W(x)=1/√(1-x²) są wielom. Czebyszewa. 4.Apr. trygonometryczna. W zagadnieniach apr. często spotykamy się z przyp. gdy f. f(x) jest okresowa. Taką f. wygodniej i lepiej jest apr. nie wielom. algebraicznymi, a wielom. tryg. Jako przykład ukł. ortogonalnego f. można podać układ (baza) 1, cosx, sinx, cos2x, sin2x, ..., cosmx, sinmx ortogonalny na odcinku [0;2π] ponieważ przez całkowanie można łatwo sprawdzić słuszność następujących równości dla dowolnych całkowitych k i n. ∫_[0;2π]sinnx sinkx dx={0, k≠n, k=n=0; π, k=n} (13) ∫_[0;2π]coskx sinnx dx=0, ∫_[0;2π]cosnx coskx dx={0, k≠n; π, k=n≠0; 2π, k=n=0} Niech dana f. ciągła f(x) ma dla ustalenia wagi okres równy 2π. J. okres=L to wyst. wykonać podstawienie x₁=Lx/2π, żeby sprowadzić zagadnienie do przypadku okresu=2π. Obieramy wielom. apr. T_m(x) w postaci T_m(x)=∑_[k=0;m](a_kcoskx+b_ksinkx) (14). Współczynniki a_k i b_k zapiszemy to zgodnie ze wzorami: H_m=∫_[a;b][f(x)-∑a_iϕ_i(x)]dx, ∂H_m/∂a_k=0, k=1,2,...,n; ∑a_i(ϕ_i,ϕ_k)=(f,ϕ_k) (15), k=1,2,...,m; (ϕ_i ϕ_k)=∫_[a;b]ϕ_i(x)ϕ_k(x)dx=0, i≠k; (f,ϕ_k)=∫_[a;b]f(x)ϕ_k(x)dx, a_k=(f,ϕ_k)/(ϕ_k,ϕ_k)=(∫_[a;b]fϕ_kdx)/(∫_[a;b]ϕ_k²dx) i (13) a₀=(1/2π)∫_[0;2π]f(x)dx, a_k=(1/π)∫_[0;2π]f(x)coskx dx, b_k=(1/π)∫_[0;2π]f(x)sinkx dx (16). Są to znane z analizy współczynniki Fouriera f. f(x). W szczególności a₀=(1/π)∫_[0;π]f(x)dx, a_k=(2/π)∫_[0;π]f(x)coskx dx (17), (b_k=0, k=1,2,...,m) j. f(x) jest parzysta oraz a₀=0, a_k=0, b_k=(2/π)∫_[0;π]f(x)sinkx dx (18) j. f(x) jest nie parzysta. Np.: Za pom. wielom. tryg. dziewiątego stopnia należy otrzymać przybliżenie następującej f. okresowej f(x)={x(π-x), 0≤x≤π ⇒ -x²+πx; x(π+x),-π≤x≤0 ⇒ x²+πx} RYS.16 Na mocy nieparzystości danej f., ze wzorów (18) b_k=(2/π)∫_[0;π]x(π-x)sinkx dx=(2/π)[-x(π-x)coskx/k|₀^π+(1/k)∫_[0;π](π-2x)coskx dx]=(2/πk)[(π-2x)sinkx/k|₀^π+(2/k)∫_[0;π]sinkxdx]=4/πk²∫_[0;π]sinkx dx=(-4/πk³)coskx|₀^π=(4/πk³)[1-(-1)^k] zatem b_2n=0, b_2n-1=8/π(2n-1)³. Szukany wielom. tryg. przyjmuje postać: T₉(x)=8/π[sinx+sin3x/27+sin5x/125+sin7x/343+sin9x/729]. Za pom. wielom. tryg. można apr. też f. nieokresowe. W tym przypadku przybliżone przedstawienie f. f(x) nie na całej prostej (-∞;∞), lecz tylko na odcinku dł. 2π albo okresu. Przejdziemy obecnie do rozpatrzenia przypadku punktów dyskretnych. W przyp. f. okresowej o okresie 2π znamy wartości w równoodległych punktach (węzach) odcinka [0;2π] x_j=2πj/n (j=0,1,...,n-1). Dobierzemy tak współczynniki a_k, b_k wielom. T_m(x) przy n>2m, żeby wyrażenie S=∑_[j=0;n-1][f(x_j)-T_m(x_j)]² miało wartość najmniejszą. Przyrównując do zera pochodne cząstkowe S względem a_q i b_q otrzymujemy następujący układ równań: {∑_[k=0;m][a_k∑_[j=0;n-1]coskx_jcosqx_j+b_k∑_[j=0;n-1]sinkx_jcosqx_j]=∑_[j=0;n-1]f(x_j)cosqx_j; ∑_[k=0;m][a_k∑_[j=0;n-1]coskx_jsinqx_j+b_k∑_[j=0;n-1]sinkx_jsinqx_j]=∑_[j=0;n-1]f(x_j)sinqx_j} (20) (q=0,1,...,m). Łatwo sprawdzić następujące równości dla k,q=0,1,...,m: ∑_[j=0;n-1]coskx_j={0, k≠0; n, k=0}; ∑_[j=0;n-1]sinkx_j=0 (21); ∑_[j=0;n-1]coskx_jsinqx_j=0 (22); ∑_[j=0;n-1]coskx_jcosqx_j={0, k≠q; n/2, k=q≠0; n, k=q=0} (23); ∑_[j=0;n-1]sinkx_jsinqx_j={0, k≠q; n/2, k=q=0; 0, k=q=0} (24). Istotnie, słuszność równości (21) jest oczywista. Dla udowodnienia ich przy k≠0 rozpatrzymy sumę ∑_[j=0;n-1]coskx_j+I∑_[j=0;n-1]sinkx_j=∑_[j=0;n]e^{Ikx (z indeksem j)}=∑e^Ik2πj/n=(1-e^i2kπ)/(1-e^i2kπ/n)=0. Przyrównując tu do zera część równości (21). W celu udowodnienia równości (22) zapiszemy lewą jej stronę w postaci ∑_[j=0;n-1]coskx_jsinqx_j=½∑_[j=0;n-1][sin(k+q)x_j+sin(q-k)x_j]=½∑_[j=0;n-1]sin(k+q)x_j+½∑_[j=0;n-1]sin(k-q)x_j. Stąd na mocy równości (21) wynika, że strona prawa tej równości=0. Analogicznie dowodzimy słuszności równości (23,24). Z równania (20) wykorzystując równości (21-24) znajdujemy a₀=1/n∑_[j=0;n-1]f(x_j), a_k=2/n∑_[j=0;n-1]f(x_j)coskx_j, b_k=2/n∑_[j=0;n-1]f(x_j)sinkx_j (25) (k=1,2,...,m). J. w (25) przejdziemy do granicy przy n→∞, to przekształcają się one we wzory (16) określające współczynniki Fouriera f. f(x). J. rozpatrywana f. f(x) jest parzysta, to wszystkie b_k=0 (k=1,...,m). W przyp. zaś jej nieparzystości a_k=0 (k=0,1,...,m) f(0)=f(π). W tych przyp. szczególnych przy obliczaniu ≠0 współczynników na podstawie wzorów (25) można wykonywać sumowanie tylko w połowie przedziału [a;b]=[0;2π] i następnie mnożyć wynik przez dwa. Np.: Za pom. wielom. tryg. stopnia 3 znaleźć przybliżenie nieparzystej funkcji f(x) danej na odcinku [0;π] następującym ciągiem wartości:

x	0	π/4	π/2	3π/4	π
y	0	5	7	6	0

Tu n=8 na [-π;π] Na mocy nieparzystości a_k=0, a dla określenia b_k na podstawie wzorów (25) mamy b_k=¼∑_[j=0;7]f(x_j)sinkx_j=½∑_[j=0;3]f(πj/4)sin(kπj/4). Obliczone współczynniki b_k podano w tablicy

x	f(x)	sinx	f(x)sinx	sin2x	fsin2x	sin3x	fsin3x
0	0	0	0	0	0	0	0
π/4	5	√2/2	5√2/2	1	5	√2/2	5√2/2
π/2	7	1	7	0	0	-1	-7
3π/4	6	√2/2	3√2	-1	-6	√2/2	3√2

2b₁=11√2/2+7, 2b₂=-1, 2b₃=11√2/2-7, b₁=11√2/4+7/2, b₂=-˝, b₃=11√2/4-7/2. Szukany wielom. tryg. przyjmuje postać T₃(x)=(11√2/4+7/2)sinx-½sin2x+(11√2/4-7/2)sin3x 4.Uwagi końcowe Często zdarza się że dla f. określonej tablicą wartości albo wykresem (np. wyniki pomiarów) chcemy dobrać wyrażenie analityczne, które będzie w przybliżeniu opisywać daną f. Proces dobierania takiego wzoru empirycznego składa się z 2 części: wyboru postaci wzoru i określenia wartości liczbowych dla występujących w nim parametrów, tak aby przybliżenie to było jak najlepsze. Czasami mamy pewne inf. dotyczące postaci poszukiwanego wzoru, ale najczęściej postać wzoru dobieramy „na oko” na zasadzie podobieństwa wykresu analizowanej f. do znanych f. Staramy się przy tym, aby f. te miały stosunkowo prostą postać. Takimi często używanymi wzorcami są m. in.: y=ax^b+c; y=exp(ax²+bx+c); y=ae^bx+c; y²=ax²+bx+c; y=ax²+bx+c; y=(ax+b)/(cx+d); y=1/(ax²+bx+c); y=x/(ax²+bx+c); y=a+b/x+c/x²; y=ax^be^cx; y=ae^bx+ce^dx. Ponieważ przyjęta „metoda” określania podobieństwa wykresów jest bardzo niedokładna po wybraniu postaci wzoru, ale przed doborem parametrów, sprawdza się jego przydatność za pom. met. wyrównywania. Wielkość parametrów można określić albo w spos. dokładny omówioną met. najmniejszych kwadratów albo w spos. przybliżony met. przeciętnych. Całkowanie numeryczne. 1.Wstęp. Uwagi ogólne o całkowaniu numerycznym. Ponieważ wyznaczenie f. pierwotnej w wielu przypadkach jest bardzo trudne lub wręcz niemożliwe stosowanie metod przybliżonych jest więc często konieczne. Ponadto, j. f. podcałkowa jest określona za pom. tablicy to pojęcie f. pierwotnej traci sens i wówczas możemy obliczyć tylko przybliżoną wartość całki. Gdy przedział całkowania jest skończony, wówczas wiele sposobów przybliżonego obliczania całek polega na tym, że f. podcałkową F(x) zastępujemy interpolującą f. ϕ(x), którą łatwo można całkować. Niech ϕ(x) będzie wielom. interpolacyjnym Lagrange'a dla f. F(x) z węzłami interpolacji x₀,x₁,...,x_N ϕ(x)=L_N(x)=∑_[k=0;N]ψ_k(x)F(x_k); ψ_k(x)=(x-x₀)...(x-x_k-1)(x-x_k+1)...(x-x_N)/(x_k-x₀)...(x_k-x_k-1)(x_k-x_k+1)...(x_k-x_N). Podstawiając w miejsce f. podcałkowej F(x) wielom. ϕ(x), otrzymamy ∫_[a;b]F(x)dx ≈ ∫_[a;b]ϕ(x)dx=∑_[K=0;N]A_kF(x_k), A_K=∫_[a;b]ψ_k(x)dx. J. F(x)-ϕ(x)|<ξ i xЄ[a,b], to |∫_[a;b](F(x)-ϕ(x))dx|≤ξ(b-a). Przedstawiony sposób pozwala zatem obl. całkę z dowolną dokł. j. tylko f. F(x) można przybliżyć wielom. z dowolną dokładnością; oprócz wielom. można stosować też inne f. interpolacyjne ϕ(x), np.: f. wymierne, trygonom. Podobny wzór przybliżonego całkowania otrzymamy gdy przedział całkowania podzielimy na podprzedziały i w każdym podprzedziale zastosujemy powyższy sposób obliczania całki. J. f. F(x) ma osobliwości utrudniające dobre przybliżenie (np. gdy F(x) jest nieograniczona lub jej pochodne niskiego rzędu nie istn. w przedziale całkowania) to postępowanie można zmodyfikować. F. podcałkową przedstawiamy w postaci iloczynu F(x)=p(x)f(x), gdzie f(x) jest już f., którą można dobrze przybliżyć, a p(x) ma wszelkie osobliwości F(x) utrudniające to przybliżenie. Niech teraz ϕ(x) będzie wielom. interpolacyjnym Lagrange'a dla f. f(x) z węzłami interpolacji x₀,x₁,...,x_N. Podstawiając w miejsce f(x) wielom. ϕ(x) otrzymamy ∫_[a;b]F(x)dx=∫_[a;b]p(x)f(x)dx ≈ ∫p(x)(x)dx=∑_[K=0;N]A_k'f(x_k); A_k'=∫_[a;b]p(x)ψ_k(x). Ten zmodyfikowany sposób można zastosować i do obliczania całek na nieskończonym przedziale całkowania. F. p(x) określa wówczas szybkość malenia F(x), gdy x→±∞. Rozbicie f. podcałkowej F(x) na iloczyn f(x) i p(x) niekoniecznie musi być związane z wydzieleniem osobliwości F(x). W met. Gaussa f. p(x) dobiera się często ze wzgl. na wybór odpowiedniego typu wielom. ortogonalnych. Do przybliżonego obliczania całek I(f)=∫_[a;b]p(x)f(x)dx (1) będziemy stosować wzory S(f)=∑A_kf(x_k), x_kЄ[a;b] (2) przy czym współczynniki A_k nie zależą od f(x). Wzór (2) naz. kwadraturą, a x_k węzłami kwadratury. F. p(x) naz. f. wagową lub wagą. Błąd przybliżenia całki (1) sumą (2) czyli błąd metody, oznaczymy jako E(f)=I(f)-S(f) (3). Za krytrium dokładności kwadratury przyjmujemy zgodność S(w) i I(w), gdy w jest wielom. Mówimy, że wzór (2) jest rzędu r, j.: 1) I(w)=S(x) dla wszystkich wielom. W(x) stopnia mniejszego niż r. 2) istnieje wielom. W(x) stopnia r taki, że I(w)≠S(w). Zakładać będziemy, że r≥1. 2.Kwadratury Newtona-Cotesa. Rozpatrzymy następujący sposób obliczania całki (1) S(f)=I(L_n)=∑A_kf(x_k), A_k=∫_[a;b]p(x)ψ_k(x)dx (5) gdzie L_n jest wielom. interpolacyjnym Lagrange'a dla f. f(x) z węzłami interpolacji x₀,x₁,...,x_N L_n(x)=∑ψ_k(x)f(x_k) (6), wzór na błąd przybliżenia f(x) wielom. L_n(x) jest R_n+1(x)=f(x)-L_w(x)=1/(N+1)!w_n+1(x)f⁽ⁿ⁺¹⁾(ζ), ζЄ[a;b] (7) gdzie w_N+1(x)=(x-x₀)(x-x₁)...(x-x_N) TW Rząd kwadratury (2) wynosi co najmniej (N+1) ⇔ kwadratura ta jest postaci (5). Kwadratury (5) z węzłami równoodległymi naz. kwadraturami Newtona-Cotesa. Niech końce przedziału są węzłami i waga p(x)≡1. S(f)=∑_[k=0;N]A_kf_k, A_k=∫_[a;b]ψ_k(x)dx=[h(-1)^N-k/k!(N-k)!]∫_[0;N]t(t-1)...(t-N)/(t-k) dt (8) h=(b-a)/N, f_k=f(a+kh) Podst. własności kwadratury: 1. J. N jest liczbą nieparzystą, to kwadratura (8) nie jest dokładna dla wszystkich wielom. stopnia (N+1) i (N+1) jest jej rzędem dokładności. Gdy jest dokładna dla wielom. stopnia (N+1) i jej rząd wynosi (N+2). Dla dowolnego N rząd kwadratury wynosi zatem 2[N/2]+2 (9) gdzie [N/2] oznacza całkowitą część liczby N/2. 2. J. fЄC^P([a;b]), gdzie r=2[N/2]+2, to błąd kwadratury można przedstawić w postaci F(f)=C_Pf_x^(P)(ζ), ζЄ[a;c], współczynnik C_P nie zależy od f. 3. Błąd kwadratury zależy od f^(P)=d^Pf/dx^P, a jak wiadomo pochodne wysokiego rzędu trudno jest oszacować. Inną trudnością jest to że pochodne te często rosną, nieograniczenie; ten szybki wzrost może już wystąpić dla pochodnych niskiego rzędu. Ponadto dla dużych N kwadratura jest mało odporna na niedokładność obliczeń na maszynach cyfrowych 4. Zaznaczamy, że A_k zależą od N (A_k=A_k^N). Dla pewnych k mamy lim_[N→∞]|A_k^N|=∞, a więc metoda nie jest zbieżna w klasie f. ciągłych. Podamy teraz wzory (8) gdy N=1 i N=2. Dla N=1, mamy A₀=-h∫_[0;1](t-1)dt=h/2, A₁=h∫_[0;1]tdt=h/2 a więc S(f)=h/2(f₀+f₁) (9). Błąd tego wzoru na podstawie (7) E(f)=1/2!∫_[a;b](x-a)(x-b)f_x''(ζ)dx (8) gdy fЄC²([a;b]), ζЄ[a,b]. F. W₂(x)=(x-a)(x-b) nie zmienia znaku na [a;b], możemy zatem do (8) zast. uogólniona tw. o wartości średniej całek. E(f)=f''(ζ₁)/2∫_[a;b](x-a)(x-b)dx=-h³f''(ζ₁)/12 (9) ζ₁Є[a;b], h=b-a RYS.17 Dla N=2 otrzymujemy A₀=h/3, A₁=4h/3, A₂=h/3 S(f)=h/3(f₀+f₁+f₂) (10). Jest to znany wzór parabol, zwany też wzorem Simpsona. RYS.18 Wzór (10) jest też dokładny dla wielom. st. 3-go. Niech f(x)ЄC⁴([a;b]) i niech L₃(x) będzie wilom. interpolacyjnym dla f(x) z węzłami a, (a+b)/2, b i x₄, gdzie x₄ jest dowolnym ale różnym od pozostałych węzłów punktem z przedziału [a;b]. Mamy wówczas f(x)-L₃(x)=(x-a)(a-(a+b)/2)(x-b)(x-x₄)f_x''''(ζ)/4! (11) ζЄ[a;b]. Przechodząc w (11) do granicy z x₄→(a+b)/2 otrzymamy f(x)-L₃^~~(x)=(x-a)(x-(a+b)/2)(x-b)f_x⁽⁴⁾( ζ^~~)/4!, (12) ζ^~~Є[a;b]. Błąd kwadratury E(f)=1/4!∫_[a;b](x-a)(x-(a+b)/2)²(x-b)f⁽⁴⁾(ζ^~~)dx (13). Po zastosowaniu do (13) uogólnionego tw. o wartości średniej dla całek otrzymamy E(f)=f⁽⁴⁾(ζ₁)/4!∫_[a;b](x-a)(x-(a+b)/2)²(x-b)dx=-h⁵f⁽⁴⁾(ζ₁) (14) gdzie ζ₁Є[a;b], h=(b-a)/2. 3.Kwadratury złożone Newtona-Cotesa. Z powodów, które przedstawiliśmy w poprzednim punkcie, kwadratur Newtona-Cotesa wysokich rzędów używa się rzadko. Dla N równych 1 i 2 współczynniki C_P we wzorze E(f)=C_Pf^(P)(ζ), określającym błąd wynoszą odpowiednio (b-a)³/12, (b-a)⁵/2880. Błąd tych kwadratur jest zatem proporcjonalny do pewnej potęgi długości przedziału całkowania i przy dużych przedziałach stosow. kwadratur nawet niskiego rzędu może nie zapewnić żądanej dokładności. W praktyce kwadratury te stosuje się w nast. spos.: 1. przedział całkowania [a;b] dzielimy na pewną liczbę podprzedziałów 2. w każdym podprzedziale stos. kwadraturę niskiego rzędu i wyniki sumujemy. Kwadraturę która jest sumą kwadratur na podprzedziałach naz. kwadraturą złożoną. Omówimy wzór złożony trapezów RYS.19. Przedział całkowania [a;b] dzielimy na m części o dł. h=(b-a)/m. W każdym podprzedziale stosujemy wzór trapezów i wyniki sumujemy. Otrzymamy S(f)=∑_[k=0;m-1]h(f_k+f_k+1)/2=h(f₀/2+f₁+...+f_m-1+f_m/2) (15) gdzie f_i=f(a+ih). Niech fЄC²([a;b]). Błąd złożonego wzoru trapezów wynosi E(f)=-h³/12∑_[k=0;m-1]f''(ζ_k)=(-(b-a)³/12m²)(1/m)∑_[k=0;m-1]f''(ζ_k) (16) przy czym ζ_kЄ[a+kh;a+(k+1)h]. Czynnik 1/m∑_[k=0;m-1]f''(ζ_k) jest średnią arytmetyczną wartości drugiej pochodnej w punktach ζ_k. Jest oczywiste że wartość ta znajduje się między najmniejszą i największą wartością f''(x) na [a;b]. Ponieważ pochodna f''(x) jest ciągła na [a;b] więc istn. punkt ζЄ[a;b] taki że 1/m∑_[k=0;m-1]f''(ζ_k)=f''(ζ) Zatem E(f)=(-(b-a)/12m²)f''(ζ) (17). Wzór złożony parabol. W tym przypadku dzielimy [a;b] na m równych części o dł. h=(b-1)/m, przy czym m jest parzyste. W przedziałach [a;a+2h]=[x₀;x₂], ..., [a+(m-2)h;b]=[x_m-2;x_m] o dł. 2h stosujemy wzór parabol, wyniki sumujemy S(f)=h/3∑_[k=1;m/2](f_2k=2+4f_2k-1+f_2k)=h/3[f₀+f_m+2(f₂+f₄+f_m-2)+4(f₁+f₃+...+f_m-1)]. Niech fЄC⁴([a;b]). Błąd E(f)=-h⁵/90∑_[k=1;m/2]f⁽⁴⁾(ζ_k), ζЄ[a+2(k-1)k;a+2kh]. Istn. punkt ζЄ[a;b] taki, że 2/m∑_[k=1;m/2]f⁽⁴⁾(ζ_k)=f⁽⁴⁾(ζ). Stąd mamy E(f)=(-(b-a)/180m⁴)f⁽⁴⁾(ζ) (18). Rozwiązywanie przybliżone równań różniczkowych zwyczajnych. W przypadku nawet całkiem prostych równań różniczkowych zwyczajnych nie zawsze można wyrazić szukane rozw. w postaci wzoru. Dlatego duże znaczenie mają met. przybliżone rozwiązywania. Met. przybliżone można podzielić na: analityczne i numeryczne. Do met. analitycznych zaliczamy met. kolejnych przybliżeń Picarda i met. szeregów potęgowych, a do met numerycznych met. Eulera i jej odmiany oraz met. Rungego-Kutty, Adamsa i Störmera. W przypadku zagadnień brzegowych rozpatrzymy numeryczne met. różnicowe. 1. Met. kolejnych przybliżeń. Jedną z prostych co do pomysłu (jednak nie w praktycznym realizowaniu) przybliżonego rozwiązania zagadnienia Cauchy'ego dla równania dy/dx=f(x,y) (1) z warunkiem początkowym y(x₀)=y₀ (2) jest met. kolejnych przybliżeń. Rozwiązanie zagadnienia (1-2) jest równoważne rozwiązaniu równania całkowego y(x)=y₀+∫_[x0;x]f(t,y)dt (3) niewiadomej f. y dowolną f. - przybliżenie zerowe, np. y=y₀ - po wykonaniu całkowania otrzymujemy pierwsze przybliżenie y₁(x)=y₀+∫_[x0;x]f(t,y₀)dt. Postępując z y₁(x) tak samo jak z y₀ dostajemy drugie przybliżenie y₂(x)=y₀+∫_[x0;x]f(t,y₁(t))dt Powtarzając ten proces otrzymamy przybliżenie y_n+1(x)=y₀+∫_[x0;x]f(t,y_n(t))dt (4) Przy dość ogólnych założeniach co do f. f(x,y) ciąg {y_n} jest zbieżne jednostajnie do rozwiązania zagadnienia (1), (2) -gdy n→∞. Np.: Rozwiązać zagadnienie y'=x-y, y(x₀)=1. Jako przybliżenie zerowe obieramy y₀=1, wówczas: y₁(x)=1+∫_[0;x](t-1)dt=1-x+x²/2; y₂(x)=1+∫_[x0=0;x](t-1+t-t²/2)dt=1-x+x²-x³/6; ...; y_n=1-x+2x²/2!+2x³/3!+...+(-1)ⁿ2x⁴/n!-(-1)ⁿxⁿ⁺¹/(n+1)! Przy n→∞ y_n→2e^-x-(1-x), co pokrywa się z dokł. rozw. zad. 2.Zwyczajna Metoda Eulera. Podczas numerycznego rozw. równania (1) przy warunku początkowym y(x₀)=y₀ należy znaleźć w punktach x₀, x₁, x₂, ..., x_n,... przybliżenia y_n wartości rozwiązania dokładnego y(x_n). Różnicę Δx_n=x_n+1-x_n=h naz. krokiem siatki. Gdy f jest stałą wówczas x_n=x₀+nh (n=0,1,...). Można w przybliżeniu przyjąć, prawa str. równ. y'=f(x,y) na każdym z odcinków zawartych między punktami podziału jest stała. Na tej podst. szukane rozw. na pierwszym odcinku w przybliżeniu przyjmuje postać f. liniowej y=y₀+(x-x₀)f(x₀,y₀) (5), w szczególności w punkcie x₁ otrzymujemy y₁=y₀+hf(x₀,y₀). Równość (5) oznacza że na odcinku [x₀,x₀+h] szukaną krzywą całkową y(x) zastępujemy w przybliżeniu odcinkiem prostoliniowym wychodzącym z punktu (x₀,y₀) i o współczynniku kątowym równym f(x₀,y₀). RYS.20 Analogicznie wychodząc z punktu (x₁,y₁) znajdujemy wart. przybliżoną y₂ wg wzoru y₂=y₁+hf(x₁,y₁) itd. Dla punktu x_n=x₀+nh mamy y_n=y_n-1+hf(x_n-1,y_n-1) (6). W ten sposób jako przybliżenie szukanej krzywej całkowej otrzymamy linię łamaną o wierzchołkach w węzłach (x₀,y₀), (x₁,y₁), ..., (x_n,y_n). Metoda (wzór (6)) zwyczajna Eulera daje dosyć grube przybliżenia rozwiązania zagadnienia Cauchy'ego (1), (2) i wykorzystywana jest zazwyczaj w tych przypadkach kiedy trzeba otrzymać poglądowe przedstawienie rozwiązania na niedużym odcinku. 3.Ulepszuna met. Eulera. Istota tej met. jest następująca: Obliczamy najpierw wartości pomocnicze szukanej f. y_n+½ w punktach x_n+½=x_n+h/2 za pomocą wzoru y_n+½= y_n+f(x_n,y_n)h/2. Następnie znajdujemy wartości prawej strony f w punkcie x_n+½, f(x_n+½,y_n+½) i potem przyjmujemy y_n+1=y_n+hf_n+½ (7) RYS.21 Metoda ulepszona Eulera (7) jest nieco dokładniejsza w porównaniu ze zwyczajną met. Eulera. 4.Ulepszona met. Eulera-Cauchy'ego. Znajdujemy przybliżenie szukanego rozw. na postawie wzoru y_n+1=y_n+h/2[f(x_n,y_n|+f(x_n+1,y_n+1)] (8) tgα_n+1=f(x_n+1,y_n+1); tgα_n=f(x_n,y_n). RYS.22 Jak łatwo widać prawa strona (8) zależy od niewiadomej y_n+1. Dlatego wykorzystujemy wzór iteracyjny y_n+1^(k+1)=y_n+h/2[f(x_n,y_n)+f(x_n+1,y_n+1^(k)] (9) (k=0,1,...). Na pierwszym kroku zerowe przybliżenie można znaleźć na podstawie met. zwyczajnej Eulera y_n+1⁽⁰⁾=y_n+hf_n. Następnie można znaleźć y_n+1⁽¹⁾=y_n+h/2[f(x,y_n)+f(x_n+1,y_n+1⁽⁰⁾)]; y_n+1⁽²⁾=y_n+h/2[f(x_n,y_n)+f(x_n+1,y_n+1⁽¹⁾)] itd. Zakończymy proces iteracyjny gdy |y_n+1^(k+1)-y_n+1^(k)|≤ξ, gdzie ξ-błąd. 5.Metoda Rungego-Kutty (RK). Metoda RK jest jedną z najczęściej stosowanych met. o większej dokładności. Idea tej met. ma dużo wspólnego z ideą met. Eulera. Przedstawmy y(x) w otoczeniu każdego punktu x=x_n wzorem Taylora i obliczmy współczynniki rozwinięcia bezpośrednio na podst. f(x,y) wykorzystując warunki y(x₀)-y₀ (2) y(x)=y_n+hdy/dx|_Xn+(h²/2!)(d²y/dx²)|_Xn+(h³/3!)(d³y/dx³)|_Xn+... (10) xЄ[x_n,x_n+1]. W zależności od tego do jakich wyrazów ograniczamy się we wzorze (10), otrzymujemy taką lub inną dokładność szukanego rozwiązania przybliżonego. Np.: zachowując wyrazy do rzędu h włącznie dojdziemy do met. Eulera. W met. RK ograniczamy się do 4 lub 5 wyrazów rozwinięcia (zachowując wyrazy z potęgami do h³ lub h⁴ włącznie). Rozpatrzymy met. RK trzeciego rzędu, zachowujemy w (10) wyrazy do h³ y(x)≈y_n+h(dy/dx)|_Xn+(h²/2!)(d²y/dx²)|_Xn+(h³/3!)(d³y/dx³) (11) przyjmujemy y_n+1=y_n+λ_n (12), gdzie λ_n=h(dy/dx)+(h²/2)(d²y/dx²)+(h³/6)(d³y/dx³). Wielkość λ_n przedstawiamy za pom. kombinacji liniowych postaci λ_n=αk₁+βk₂+γk₃+δk₄ (14), gdzie α,β,γ,δ -współczynniki nieokreślone, k₁,k₂,k₃,k₄ -liczby okr. równościami: {k₁=hf(x_n,y_n); k₂=hf(x_n+h/2,y_n+k₁/2); k₃=hf(x_n+h/2,y_n+k₂/2); k₄=hf(x_n+h,y_n+k₃)} (15). Wyrazimy pochodne λ_n przez prawą str. f(x,y). dy/dx=f(x,y); d²y/dx²=∂f/∂x+(∂f/∂y)(dy/dx) itd. Dla prostoty wprowadzimy operator D=∂/∂x+f(∂/∂y) wówczas d²y/dx²=Df, d³y/dx³=D(Df)=∂²f/∂x²+2f(∂²f/∂x∂y)+f²(∂²f/∂y²)+∂f/∂y(∂f/∂x+f(∂f/∂y))=D²f+(∂f/∂y)Df. Podstawiając znalezione wartości pochodnych w równości (13) dla λ_n mamy λ_n=hf+(h²/2)Df+h³/3(D²f+(∂f/∂y)Df) (16). Przedstawimy obecnie k₂,k₃,k₄ jako f. dwóch zmiennych za pom. szeregu Taylora k₂=f(x_n+h/2,y_n+k₁/2)h=[f+((h/2)∂/∂x+(k₁/2)∂/∂y)f+½((h/2)∂/∂x+(k₁/2)∂/∂y)f+...]h. Ograniczając się do trzeciej potęgi h i uwzględniając że k₁=hf, mamy k₂=hf+(h²/2)Df+(h³/8)D²f, k₃=hf+(h²/2)Df+(h³/4)(∂f/∂y)Df +(h³/8)D²f, k₄=hf+h²Df+(h³/2)(∂f/∂y)Df +(h³/2)D²f. Utworzymy obecnie sumę λ_n=αk₁+βk₂+γk₃+δk₄. Porównując w ostatniej równości i w wyrażeniach (16) dla λ_n współczynniki przy jednakowych potęgach h otrzymamy ukł. równań: {α+β+γ+δ=1 (przy hf); β/2+γ/2+δ=½ (przy h²Df); β/8+γ/8+δ/2=1/6 (przy h³D²f); γ/4+δ/2=1/6 (przy h³(∂f/∂y)Df)} Ukł. ten ma rozw. α=δ=1/6, β=γ=1/3. W ten spos. λ_n=1/6(k₁+2k₂+2k₃+k₄) (17) a do obliczania y_n+1 mamy wzór y_n+1=y_n+1/6(k₁+2k₂+2k₃+k₄) (18), gdzie k₁,k₂,k₃,k₄ -liczby określone przez (15). Dla równania y'=f(x,y) Bieberbach otrzymał za pom. rozwinięcia Taylora następujące oszacowanie błędu met. RK |y_n+1-y(x_n+1)|<6MN|x_n+1-x_n|⁵|N⁵-1|/|N-1|, gdzie M i N -stałe takie że w obszarze |x-x_n|<a, |y-y_n|<b spełnione są nierówności |f(x,y)|<M, |∂^i+kf/∂xⁱ∂y^k|<N/M^k-1 (i+k≤3), |x-x_n|N<1, aM<b, h=x_n+1-x_n≤a. W ogólnym przypadku trudno jest wskazać dokładne oszacowanie bł. met. RK. Dlatego w praktyce korzysta się często z przybliżonego oszacowania bł. tej met. wychodząc z różnych rozwiązań pośrednich. Starając się zwiększyć dokł. rozw. kosztem zmniejszenia kroku kierujemy się w praktyce następującą zasadą: tworzymy różnicę |k₂-k₃| i |k₁-k₂| oraz żądamy żeby pierwsza z nich nie przekraczała kilku procent wartości drugiej. J. ten warunek nie jest spełniony to krok podziału należy zmniejszyć. Odpowiednia liczba Θ=|(k₂-k₃)/(k₁-k₂)| jest swego rodzaju „miarą niepewności” lub „wskaźnikiem kroku”. W celu grubego oszacowania bł. metody wykorzystujemy często zas. Rungego: przypuśćmy dla ustalenia uwagi, że mamy do czynienia z met. RK o dokł. h^k i niech x będzie punktem w którym interesuje nas rozwiązanie przybliżone. Znajdujemy to rozw. wykorzystując krok h a następnie krok 2h. Przy kroku h dla punktu x₁=x₀+h mamy y(x₁)=y₁+Ah^k+1, A=f^k+1(ζ)/(k+1)! (x₀≤ζ≤x₁). Bł. wynosi Ah^k+1. W celu orientacyjnego obliczenia bł. w p. x=x₀+2nh można przyjąć że na każdym kroku o dł. h czynimy bł. proporcjonalny do h^k+1 i dla p. x bł. sumaryczny jest równy Ah^k+12n czyli y(x)=y_2n+A2nh^k+1 (19). W obliczeniach z krokiem 2h również zakładamy że na każdym kroku dopuszczamy bł. proporcjonalny do (2h)^k+1 z tym samym współczynnikiem proporcjonalności A i w p. x bł. sumaryczny jest równy An(2h)^k+1 to jest y(x)=y_n+An(2h)^k+1 (20). Z równości (19) i (20) |y_2n-y(x_2n)|≈|y_2n-y_n|/(2^k-1). W ten spos. dochodzimy do zas. Rungego. Szukany bł. met. RK jest w przybliżeniu równy modułowi różnicy rozwiązań o kroku h i kroku 2h dzielonej przez 2^k-1 przy k=4 mamy 2⁴-1=15. (błąd kroku ~ h⁵). 6.Rozwiązywanie ukł. równań różniczkowych. Zagadnienia początkowe: Znaleźć f. y_i=y_i(x) i=1,2,...,m spełniające na przedziale [x₀=a;b] ukł. równ. różniczkowych {dy₁/dx=f₁(x,y₁,...,y_m); dy₂/dx=f₂(x,y₁,...,y_m); ...; dy_m/dx=f_m(x,y₁,...,y_m)} (21) i warunek początkowy y(x₀)=y_i,0 (22) (i=1,...,m) albo w postaci wektorowej dy^/dx=f^(x,y^), y^(x₀)=y₀^ (23) y^=[y₁...y_m]^T, f^=[f₁...f_m]^T. Jak i dla jednego równania skalarnego y'=f(x,y) mamy metodą Eulera wg której y_n+1^=y_n^+hf^(x_n,y_n^), n≥1. y^(x₀)=y₀^ (24). Ulepszona met. Eulera-Cauchy'ego y_n+1^^(k+1)=y_n^+h/2[f^(x_n,y_n^)+f(x_n+1,y_n+1^^(k))] (25) y_n+1^⁽⁰⁾=y_n^+f^(x_n,y_n^), (k=1,2,...). Metody typu Rungego-Kutty y_n+1^=y_n^+∑_[i=1;s]w_ik_i^, k₁^=hf^(x_n,y_n^) (26), k_i^=hf^(x_n+a_ih,y_n^+∑_[j=1;i-1]b_ijk_j^), w_i,a_i,b_ij -stałe. 7.Met. różnicowe rozwiązywania zagadnień brzegowych. Przez liniowe zagadnienie brzegowe rzędu n rozumiemy równanie różniczkowe rzędu n≥2 względem f. niewiadomej y(x) L(y)=r(x) (27), gdzie L(y)=∑_[j=0;n]f_j(x)y⁽ⁱ⁾ przy n warunkach brzegowych s_j=γ_j, (j=1,2,...,n), s_j(y)=∑_[k=0;n-1][α_jky^(k)(a)+β_jky^(k)(b)] (28), gdzie: f_j(x), r(x) -dane f. ciągłe na obcinku [a;b], α_jk,β_jk,γ_j -dane stałe. J. r(x)=0 to równanie (27) naz. jednorodnym, j. γ_j=0 to odpowiednie warunki brzegowe naz. jednorodnymi. Zagadnienie brzegowe naz. jednorodnym, gdy r(x)=0 i γ_j=0. Rozwiązanie zagadnienia brzegowego polega na znalezieniu takiej f. y(x) która spełnia równanie (27) i warunki (28). Niech będzie dane równanie różniczkowe rzędu drugiego L(y)=f₀(x)y+f₁(x)y'+f₂(x)y''(x)=r(x) z warunkami brzegowymi: γ₁=s₁=α₁₀y(a)+β₁₀y(b)+α₁₁y'(a)+β₁₁y'(b); γ₂=s₂=α₂₀y(a)+β₂₀y(b)+α₂₁y'(a)+β₂₁y'(b) (29). Z najczęściej stosowanych zagadnień brzegowych jest zagadnieniem trzech typów na końcach przedziału [a;b]. Pierwsze zagadnienie na wartości y(x) na końcach [a;b] y(a)=y_a; y(b)=y_b (30) α₁₀=1, β₂₀=1, α_jk=β_jk=0 dla k=1, β₁₀=0, α₂₀=0. Drugie zagadnienie y'(a)= γ_a; y'(b)= γ_b (31). Trzecim zagadnieniem nazywamy zadanie kombinację liniową: {α₁₀y(a)+α₁₁y'(a)=γ₁; β₂₀y(b)+β₂₁y'(b)=γ₂} (32) U podstaw met. różnicowej leży zastąpienie zależności różniczkowych związkami różnic skończonych. Oznaczmy punkty przedziału [a;b] przez x_i=x₀+ih, h=(b-a)/n, (i=0,1,...,n). Niech y(x_i) jest wartością szukanej f., a jej wartości przybliżone to y_i. Zastępujemy pochodne ilorazami różnicowymi y'(x_i)≈[y(x_i+1)-y(x_i-1)]/2h=(y_i+1-y_i-1)/2h (33); y''(x_i)≈[y(x_i+1)-2y(x_i)+y(x_i-1)]/h²=(y_i+1-2y_i+y_i-1)/h² (34). W przypadku ogólnym dla k parzystego mamy y^(k)(x_i)≈(1/h^k)Δ^ky_i-k/2 i nieparzystym y^(k)(x_i)=1/2h^k(Δ^ky_i-[(k+1)/2]+Δ^ky_i-[(k+1)/2]). Podstawiając do równania (27) zamiast pochodnych iloraz różnicowy otrzymamy tzw. równanie różnicowe. Postępując dokładnie tak samo z warunkami brzegowymi sprowadzamy rozwiązanie zagadnienia brzegowego do rozwiązania ukł. równań algebraicznych względem y_i. Np.: y''+x²y+2=0; y(-1)=y(1)=0; a=-1, b=1. Dzielimy przedział [a;b]=[-1;1] na części o dł. h=½. RYS.23 Na podstawie warunków brzegowych i symetrii zadania mamy y_-2=y₂=0; y_-1=y₁. Dlatego zagadnienie sprowadza się do znalezienia 2 niewiadomych y₀ i y₁. Równanie różnicowe zapisujemy w postaci: {2+(y₁-2y₀+y₁)/h²=0, dla x=0; ¼y₁+2+(0-2y₁+y₀)/h⁴=0, dla x=½} ⇒ {y₀=1,05; y₁=0,8} 8.Met. „przeganiania”. Przy rozwiązywaniu wielu specjalnych zagadnień, np. przy obliczaniu reaktorów jądrowych, mamy często do czynienia z równaniami różnicowymi postaci A_iy_i-1-C_iy_i+B_iy_i+1=-F_i (35), przy warunkach brzegowych -C₀y₀+B₀y₁=-F₀ (36); A_Ny_N-1-C_Ny_N=-F_n (37), gdzie i=1,2,...,N-1; A_i,B_i,C_i,F_i -dane liczby dodatnie, y_i -f. szukana. Do rozwiązywania ukł. równań liniowych algebraicznych z macierzą trój-diagonalną, bardzo wygodne jest szukanie rozw. ukł. (35) w postaci y_i=α_i+1y_i+1+β_i+1 (38). Dla i-1 mamy y_i-1=α_iy_i+β_i=α_i(α_i+1y_i+1+β_i+1)+β_i (39). podstawiając (38)(39) w (35) mamy (A_iα_iα_i+1-C_iα_i+1+β_i)y_i+1+(A_iα_iβ_i+1+A_iβ_i-C_iβ_i+1+F_i)=0. Żeby były spełnione równości =0 A_iα_iα_i+1-C_iα_i+1+β_i=0 ⇒ α_i+1=B_i/(C_i-α_iA_i) (40); A_iα_iβ_i+1+A_iβ_i-C_iβ_i+1+F_i=0 ⇒ β_i+1=(A_iβ_i+F_i)/(C_i-α_iA_i) (41). Współczynniki „przeginania” α_i+1, β_i+1 można obliczać rekurencyjnie ze wzorów (40)(41) j. dane są α₁ i β_i. Dla znalezienia ich wykorzystujemy warunek brzegowy (36) i wzór (38) dla i=0 y₀=y₁B₀/C₀+F₁/C₀ (36') y₀=α₁y₁+β₁ (38'). Z porównania (36') i (38') mamy α₁=B₀/C₀, β₁=F₁/C₀. Następnie wykorzystujemy warunek brzegowy (36) i wzór (38) dla i=N-1, mamy A_Ny_N-1-C_Ny_N==F_N (36'') ma rozw. względem y_N=(A_Nβ_N+F_N)/(C_N-α_NA_N) (42). Teraz można obliczyć rekurencyjnie niewiadomą y_i na podstawie wzoru (38) dla i=N-1, N-2, ..., y₀.

---