1 Introduction to Numerical AnalysisMarek Kręglewski
2 Course content Week 1 Solutions of nonlinear equations in one variable: the bisection algorithm. Week 2 The Newton-Raphson method, the secant method. Fixed point iteration. Week 3 Numerical integration: trapezoidal rule and Simpson’s rule. Week 4 Numerical differentiation: forward and backward-difference formula. Three-point formula of numerical differentiation. Week 5 Initial value-problem for differential equations: Euler’s method, the Runge-Kutta methods. Week 6 Taylor expansion – error of a numerical method. The Richardson’s extrapolation. Week 7 Week 8 Polynomial interpolation: Newton and Lagrange polynomials. Week 9 Methods for solving linear systems: linear systems of equations, Cramer’s rule, Gaussian elimination. Week 10 Approximation theory: least-squares approximation. Week 11 Linear algebra, matrix inversion and the determinant of a matrix. Week 12 The similarity transformations. Eigenvalues and eigenvectors. Week 13 Iterative techniques in matrix algebra: Jacobi iterative method. Week 14 Optimization. Week 15 Round-off errors: absolute error, relative error, significant digits.
3 Course content 2 LABORATORY CLASSES MS Excel – general introductionApplication of the MS Excel in solving numerical problems MANUALS: E. Steiner, Mathematics for chemists, Oxford. A. Ralston, Introduction to numerical analysis.
4 Solution of equation in one variable x=f(x)START READ x , ε, A y=x x=y y=½(x+A/x) |x-y|< ε NO YES WRITE y STOP Trace of operations
5 Algorithm notation START and STOP of a sequential algorithmINPUT and OUTPUT operations SUBSTITUTION operations CONDITIONAL operation = LOOP ? SUBSTITUTION variable = expression Calculate the value of the expression and save it under the name of the variable
6 Convergent process: x=½(x+4/x)y 4 2.5 2.05 2 Iteration process
7 Divergent process: x=6-x*xy 2.1 1.59 3.4719 E+11 E+23 E+47 E+95
8 Solution of equation in one variable Bisection methodSolution of an equation f(x)=0, i.e. search for zero points of the function f(x). Search for the a zero point in the range
9 Bisection Algorithm START READ a, b, ε NO f(a)*f(b)<0WRITE: incorrect range YES p=(a+b)/2 NO f(a)*f(p)<0 YES b=p a=p |a-b|<ε NO YES WRITE a,b STOP Trace of operations
10 Differential CalculusDerivative of a function – a measure how rapidly the dependent variable changes with changes of the independent variable y=y(x) y2 Tangent line tan(α) (slope) y = y2-y1 α y1 x = x2-x1 x1 x2 derivative
11 Differential CalculusFind the derivative of the function y = a x2 Let x = x2-x1 and y = y(x2)-y(x1) y = a(x2)2-a(x1)2 = a(x1+x)2-a(x1)2 = a[(x1)2+2x1x+(x)2]-a(x1)2 = = a[2x1x+(x)2] After dividing by x In the limit as x2 → x1 (i.e. x → 0) The derivative of the function y=ax2 is dy/dx=2ax
12 Differential CalculusDerivatives of some elementary functions (a is a constant): Function y=y(x) Derivative dy/dx=y’(x) xn n xn-1 ax ax ln(a) ln(x) 1/x sin(x) cos(x) -sin(x) a Let y(x) and z(x) are differentiable functions of x: Composite function f(u(x))
13 Solution of equation in one variable Newton-Raphson methodThe search of a zero point begins at any point x0, if: 1) the function f(x) and its first derivative are continuous 2) the first derivative is different from zero y zero point x x1 x3 x2 x0 The expansion inTaylor series:
14 Newton-Raphson algorithmSTART READ x0 , ε x0=x1 x1=x0 - f(x0) / f ’(x0) |x0-x1|< ε NO YES WRITE x1 STOP Trace of operations
15 Solution of equation in one variable Secant MethodThe search for the zero point begins from a pair of points(x0, x1), if: 1) the function f(x) is continuous 2) f(x0) f(x1), when x0x1 y zero point x x2 x3 x1 x0 The first derivative from the Newton-Raphson method approximated with an expression:
16 Secant method algorithmSTART READ x0 , x1 , ε q0=f(x0) q1=f(x1) x0=x1; x1=x2 q0=q1 ; q1=f(x2) x2=x1 – q1(x1-x0) /(q1-q0) |x2-x1|< ε NO YES WRITE x2 Trace of operations STOP
17 Integral Calculus – principal factsThe antiderivative F(x) of f(x) is the function such that dF(x)/dx=f(x) The indefinite integral is the same thing as the antiderivative function A definite integral is the limit of a sum of terms f(x)x
18 Integral Calculus - examplesA car moves with constant velocity v(t)=50 km/h. Calculate the distance it covers in 2 hours. A stone is falling with the acceleration g(t) = 10 m/s2. At the begining its velocity is 0 m/s. Calculate the distance the stone covers between 2nd and 4th second of the fall.
19 Numerical integrationTrapezoidal rule h T2 Tm a b
20 Numerical integrationSimpson’s rule m must be even Sm/2 a b
21 Analytical integration – an examplef(x)=x3 f(x)=x4
22 Numerical integration – an exampleCalculation results f(x) x x3 x4 10 1000 10000 11 1331 14641 12 1728 20736 x3 x4 T(h=2) 2728 30736 T(h=1) 2695 30009 S(h=1) 2684 29766,67 I (accurate) 29766,4 f(x)=x3 f(x)=x4 h T(h) T(h)-I 2 2728 44 1 2695 11 Errors of the trapezoidal rule error ~ h2 h T(h) T(h)-I 2 30736 969,6 1 30009 242,6
23 Geometrical series When a=1 i) The sum is equal toii) is a series expansion of the function
24 Taylor series expansion at x=0constants Thus
25 Taylor series expansionconstants Thus
26 Series expansion of a functionCalculate the value f(6) using the Taylor series expansion Call the Taylor series
27 Różniczkowanie numerycznePrzybliżenia jednostronne: Średnia P i L (różnica centralna):
28 Różniczkowanie – błąd metody_ Pochodna jednostronna Pochodna centralna pochodna błąd ~ h1 pochodna błąd ~ h2
29 Przykład – obliczenie pochodnejOblicz pochodną ln(x) w punkcie x=3 metodą pochodnej centralnej oraz jednostronnej dla różnych długości kroków: f(x)=ln(x) ln'(3)=1/3 ln(3)= f'(x)=[f(x+h)-f(x-h)]/(2*h) h xh f(xh) f'(3) błąd h^2 błąd/h^2 1 4 2 0.5 3.5 0.25 2.5 0.1 3.1 0.01 2.9 f'(x)=[f(x+h)-f(x)]/h x+h f(x+h) f’(3) błąd/h Zmniejszenie kroku zmniejsza błąd, przy czym szybciej błąd maleje w metodzie różnic centralnych
30 Równanie różniczkowe I rzęduRównanie różniczkowe opisujące rozpad promieniotwórczy Propozycja rozwiązania: Sprawdzanie poprawności: Podstawienie do równania: Lewa strona równa prawej, gdy: Wartość a wyznaczana z warunku początkowego: Ostateczne rozwiązanie analityczne: k – stała szybkości rozpadu promieniotwórczego
31 Rozpad promieniotwórczyRównanie różniczkowe opisujące rozpad promieniotwórczy Rozwiązanie analityczne: Okres połowicznego rozpadu :
32 Równanie różniczkowe – metoda EuleraRównanie (f jest znaną funkcją): Wzór przybliżony na pochodną: Po przekształceniu: Uproszczony zapis: Ostatni wzór pozwala na obliczanie wartości funkcji y punkt po punkcie. Wartość funkcji w punkcie zerowym y0 określają warunki początkowe.
33 Równanie różniczkowe I rzędut N dN/dt Nanalit 1000 -5000 1 0.1 500 -2500 2 0.2 250 -1250 3 0.3 125 -625 4 0.4 62.5 -312.5 5 0.5 31.25 82.085 6 0.6 15.625 7 0.7 7.8125 8 0.8 9 0.9 11.109 10 11 1.1 12 1.2 13 1.3 14 1.4 15 1.5 16 1.6 17 1.7 18 1.8 19 1.9 20 0.0454
34 Równanie różniczkowe II rzęduF Drgania harmoniczne Fp = ma a - przyspieszenie Fw = -kx x - wychylenie Przyjmijmy: m=1 k=1 Równowaga sił Fp = Fw a=-x x Rozwiązania szczególne równania: Rozwiązanie ogólne równania: Stałe c1 i c2 wyznaczane z warunków początkowych
35 Równanie różniczkowe II rzędux -1 1 Warunki początkowe: Rozwiązanie ogólne z uwzględnieniem warunków początkowych:
36 Rozwiązanie numeryczne IKorzystamy z przybliżonych wzorów na pochodne: Oznaczamy: Z postaci równania wynika:
37 Rozwiązanie numeryczne I c.d.Gdy t=0 : k t(k) x(k) v(k) a(k) 1 -1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
38 Rozwiązanie numeryczne IIKorzystamy z przybliżonych wzorów na pochodne centralne: Oznaczamy: Z postaci równania wynika:
39 Rozwiązanie numeryczne II c.d.Gdy t=0 : k t(k) x(k) v(k+1/2) a(k) 1 -1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
40 Ekstrapolacja RichardsonaCzy wykonując obliczenia ze skończona długością kroku h można oszacować wynika graniczny dla h 0 ? F(h) – wartość obliczona dla długości kroku h a0 = F(0) hipotetyczna wartość dla zerowej długości kroku p – rząd błędu metody numerycznej Obliczamy wynik numeryczny F dla dwóch różnych kroków h i (qh)
41 Ekstrapolacja Richardsona c.d.odejmujemy stronami a0 też jest obarczone błędem i postępowanie można prowadzić dalej. Najczęściej ekstrapolację stosujemy dla q=2, a wtedy:
42 Ekstrapolacja Richardsona przykład 1Wyniki numeryczne metodą trapezów: h T(h) 2 2728 1 2695
43 Ekstrapolacja Richardsona przykład 2f(x)=ln(x) f'(x)=[f(x+h)-f(x-h)]/(2*h) ln'(3)=1/3 h P(h) D/3 a0 0.8 3.8 2.2 0.4 3.4 2.6 0.2 3.2 2.8 0.1 3.1 2.9 błąd metody różnic centralnych h2, czyli p=2. = P(h)-P(2h)
44 Interpolacja wielomianemDana jest funkcja f(x) w postaci tablicy, tzn. znamy jej wartości w (n+1) punktach (węzłach) f(x0), f(x1), f(x2), ..., f(xn). Zadanie: znaleźć wielomian n-tego stopnia taki, że: w(x0)= f(x0) w(x1)= f(x1) ... w(xn)= f(xn) wn(x) nazywamy wielomianem interpolacyjnym. Cele interpolacji: łatwe zapamiętanie postaci funkcji (współczynniki) wykonywanie operacji matematycznych na wielomianie wyznaczanie pośrednich wartości funkcji
45 Obliczanie wartości wielomianuPostać naturalna wielomianu Obliczanie wartości wielomianu wg schematu Hornera
46 Obliczanie wartości wielomianuSTART Algorytm Wczytaj n , {ai}, x w=an i=n-1 w=w*x+ai i=i-1 TAK i≥0 NIE Wypisz w STOP
47 Ślad działań w3(x)=1+3x-2x2+4x3 n=3 a0=1 a1=3 a2=-2 a3=4Oblicz wartość wielomianu w punkcie x=3. n w i 3 4 2 4*3-2=10 1 10*3+3=33 33*3+1=100 -1 Wartość wielomianu w punkcie x=3 wynosi 100.
48 Postać Newtona wielomianuNiech x0, x1, x2,..., xn-1 są danymi liczbami, dla których wartości wielomianu są określone (dane). Tworzymy wielomiany pomocnicze pk (k=0,1,2,...,n) takie, że p0(x) = 1 p1(x) = x-x0 p2(x) = (x-x0)(x-x1) ... pk(x)= (x-x0)(x-x1)... (x-xk-1) Wielomian wn(x) przedstawiamy jako Jak wyznaczyć współczynniki bk?
49 Wyznaczanie współczynników bkx f(x) f[xl,xl+1] f[xl,xl+1,xl+2] x0 f(x0) x1 f(x1) x2 f(x2) ... xn f(xn)
50 Przykład b0 = 100 b1 = 183 b2 = 58 b3 = 4 x f(x) f[x0,x1] f[x0,..,x2]466 183 7 1296 415 58 9 2782 743 82 4 b0 = 100 b1 = 183 b2 = 58 b3 = 4
51 Interpolacja liniowa y f1 f0 x x0 x1 Prosta: w1(x)=a0+a1xf(x0) = f0 = a0+a1x0 (/ x1) f(x1) = f1 = a0+a1x1 (/ x0) Wyznacz a0, a1 f1 -f0 = a1x1 – a1x a1=(f1-f0)/(x1-x0) f0x1-f1x0 = a0x1 – a0x a0=(f0x1 –f1x0 )/(x1-x0) w1(x)= [(f0x1 –f1x0 )/(x1-x0)] + [(f1-f0)/(x1-x0)] x w1(x)= [(f0x1 –f0x0 +f0x0 –f1x0 )/(x1-x0)] + [(f1-f0)/(x1-x0)] x w1(x)=f0 + [(f1-f0)/(x1-x0)] (x-x0) to postać Newtona dla w1(x) = b0 p0(x) + b1 p1(x) , gdzie p0(x) = b0 = f0 p1(x) = x-x b1 = (f1-f0)/(x1-x0) f0 x x0 x1
52 Zjawisko Rungego Przy interpolacji wielomianem wysokiego stopnia, np. 10-tego dla funkcji w przedziale [-1,1] dla węzłów równoodległych xi = -1 + i *0, i = 0,1,2,...,10 x f(x) w(x) -1 -0.8 -0.6 0.1 -0.4 0.2 0.5 -0.2 1.5 2.5 1 1.48E-15 -2.5 -12.5 -25 -31.25 0.4 -1.5 25 62.5 93.75 0.6 -0.5 -1.5E-15 -93.75 0.8
53 Zjawisko Rungego Porównanie wykresu funkcji i wielomianu:
54 Rozwiązywanie układu równańPrzykład: x1+2x2+3x3=1 2x1+3x2+4x3=1 3x1+4x2+ x3=1 x1=-1 x2=1 x3=0
55 Macierze a przekształcenia geometrycznex1 x2 y1 inwersja P y2 x1 x2 y1 odbicie w płaszczyźnie P y2 x1 x2 y1 obrót P φ Macierze transformacji geometrycznych są macierzami ortogonalnymi
56 Przekształcenie macierzy przez podobieństwoIstnieje odwzorowanie A, które przekształca x y: y’ x x’ y Jeżeli wektory x’ i y’ przekształcane są do wektorów x i y poprzez transformację Q, jak wygląda odwzorowanie wektora x’ w wektor y’ ? Jeżeli oraz , to Jeżeli macierz Q jest nieosobliwa, to Macierze A i B są swoimi transformatami przekształconymi przez podobieństwo
57 Przekształcenie - przykład2 1’ 2’ φ=-45° 1
58 Przekształcenie - przykład2 1’ 2’ x -45° 1 y (x1,x2)=(1,2) (y1,y2)=(3,-1) 1+2=3 1-2=-1 (x’1,x’2)= (y’1,y’2)=
59 Równanie charakterystyczne macierzyλ – skalar , A(nn) I(nn) K(nn) K = A – λI macierz charakterystyczna macierzy A detK = K(λ) = det(A - λI)=A - λI= równanie charakterystyczne macierzy K(λ) = λn + an-1 λn-1 + an-2 λn a1 λ + a0 = 0 Pierwiastki wielomianu K(λ): λ1 , λ2 , ... , λn-1 , λn nazywamy pierwiastkami charakterystycznymi (wartościami własnymi) macierzy A. Jeżeli B = Q-1AQ, to macierz charakterystyczna macierzy B K = B – λI = Q-1AQ - Q-1IQ = Q-1(A - I)Q , a wyznacznik detK =B - λI= Q-1 A - λI Q = A - λI= 0 Dwie macierze związane przekształceniem przez podobieństwo mają te same pierwiastki charakterystyczne.
60 Pierwiastki charakterystyczne
61 Macierz diagonalna Jeżeli istnieje takie przekształcenie przez podobieństwo, które macierz A sprowadzi do macierzy diagonalnej D , to elementy na przekątnej macierzy diagonalnej są zarazem pierwiastkami charakterystycznymi (wartościami własnymi) macierzy A.
62 Przykład diagonalizacji macierzyAby wyzerować elementy niediagonalne: Po przekształceniu otrzymujemy macierz:
63 Pierwiastki i wektory charakterystyczneC-1AC jest przekształceniem diagonalizującym macierz A. Kolumny macierzy C są wektorami charakterystycznymi. Jeżeli macierz C jest ortogonalna, to C-1=CT , a C-1AC = CTAC. Obustronne pomnożenie macierzy A przez wektor charakterystyczny daje wartość charakterystyczną: Ogólnie:
64 Regresja liniowa Regresja liniowa: y=a*x+b(x1,y1) (x2,y2) Regresja liniowa: y=a*x+b Zadanie: Wyznaczyć optymalne wartości a oraz b.
65 Regresja liniowa Podstawowe założenia:Rozkład yi wokół linii prostej jest losowy Wariancja σy2 jest niezależna od x Metoda najmniejszych kwadratów: Wyznaczamy min Φ(a,b) względem a oraz b:
66 Regresja liniowa Rozwiązanie układu równań ze względu na a, b:
67 Regresja liniowa Estymata wariancji dla wartości yi:Estymata wariancji dla parametrów a oraz b: Współczynnik korelacji liniowej dla próbki r Wartość r zawiera się między -1 i +1. r>0 wskazuje na zależność dodatnią, a r<0 na zależność ujemną między x oraz y. r=0 wskazuje na brak zależności liniowej między x oraz y.
68 Regresja liniowa - przykład
69 x [km] y [kg] x*x x*y y-a*x-b (y-a*x-b)^2 x-xsr y-ysr 1 -2 -0.4 0.16 -4 18 3 -10 9 -30 0.8 0.64 10 5 -20 25 -100 7 49 -210 -0.8 2 -38 81 -342 0.4 4 -18 Sum: 165 -684 0.00 1.6 a= -4.6 kg/km b= 3 kg s^2= 0.5333 s= kg sa^2= 0.0133 sa= 0.1155 sb^2= 0.44 sb= 0.6633 xsr= cov(x,y)= ysr= var(x)= 8.0000 var(y)= r(x,y)=
70 Więcej o korelacji - kwadrantyIV III μy I II μx Kwadranty: I x-μx<0 y-μy<0 (x-μx)(y-μy)>0 II x-μx>0 y-μy<0 (x-μx)(y-μy)<0 III x-μx>0 y-μy>0 (x-μx)(y-μy)>0 IV x-μx<0 y-μy>0 (x-μx)(y-μy)<0
71 Współczynnik korelacji liniowejx r=-1 -1
72 Regresja liniowa jako układ równańNiewiadome: a , b Szukamy rozwiązania takiego, aby uzyskać Zapis macierzowy:
73 Układ równań nadmiarowyPoszukujemy rozwiązania a, dla którego T jest minimalne. Tak obliczone wartości parametrów a zapewniają minimalizację sumy kwadratów odchyleń od prostej
74 Przykład przedstawienia macierzowego
75 Wariancje Wariancja dla zmiennej yWariancje i kowariancja dla parametrów Współczynnik korelacji liniowej
76 Jakobian W regresji liniowej funkcja modelu to prosta y= a*x + b.Jakobian to macierz pochodnych po parametrach a, b we wszystkich punktach danych i = 1,2,...,n Jeżeli do danych chcielibyśmy dopasować wielomian 2-go stopnia y= a0 + a1*x + a2*x2 , to jakobian miałby postać:
77 Rozkład złożonego pasmaPasmo doświadczalne Należy dopasować do pasma krzywe Gaussa w postaci a - wysokość b - położenie c - szerokość
78 Metoda najmniejszych kwadratów{ak}, k=1:M , M dopasowanych parametrów Funkcja błędu (suma po n punktach): Φ{ak} = j [yj(dośw) - yj({ak}]2 Zadanie Minimalizować Φ modyfikując zbiór {ak} startując z wartości początkowych {ak}0
79 Funkcja błędu i jakobianElementy jakobianu Rozkład na N pasm
80 Algorytm Poprawiona wartość {ak}
81 Metoda najmniejszych kwadratów
82 Calculation precisionSources of errors: input data errors round-off errors cut-off errors the model errors accidental errors Absolute and relative errors: approximate value accurate value absolute error relative error
83 Round-off and cut off errorsround-off cut-off round-off to t digit after the decimal point the resulting absolute error ½·10-t Example above: ½·10-3 = 0.,240 0,0005 How to round-off numbers ending with 5? in addition the errors cancel
84 Przenoszenie się błędówDodawanie i odejmowanie Jaki jest błąd sumy? Jaki jest błąd różnicy?
85 Przenoszenie się błędówDodawanie i odejmowanie Podobnie: Błąd bezwzględny sumy lub różnicy równa się sumie błędów bezwzględnych składników.
86 Znoszenie się składników przy odejmowaniubłąd bezwzględny błąd względny
87 Przenoszenie się błędówMnożenie i dzielenie Podobnie: Błąd względny iloczynu lub ilorazu równa się sumie błędów względnych czynników.
88 Wykorzystanie zasad przenoszenia błędówOblicz pierwiastki równania kwadratowego wykonując obliczenia z dokładnością do 5 cyfr znaczących. tylko 2 cyfry znaczące 5 cyfr znaczących
89 Wykorzystanie zasad przenoszenia błędówWykorzystanie wzorów Viete’a
90 Błędy maksymalne złożonych wyrażeńDana zależność funkcyjna Parametry xi obarczone błędami. Jaki jest błąd maksymalny wielkości złożonej y?
91 Przykład szacowania błędu maksymalnego
92 Błędy standardowe złożonych wyrażeńDana zależność funkcyjna sx to błędy standardowe zmiennych x. Jaki jest błąd standardowy wielkości złożonej y?
93 Przykład szacowania błędu standardowego