METODA ELIMINACJI GAUSSA

Układ zapisujemy w postaci macierzy C

{p=1,2,…,n

{{i=p+1, p+2, …,n

{{c_ij^(p)=c_ij^(p-1)

{a₁₁x₁+a₁₂x₂+…+a_1nx_n=b₁

{a₂₁x₁+a₂₂x₂+…+a_2nx_n=b₂

{…

{a_n1x₁+a_n2x₂+…+a_nnx_n=b_n

[c₁₁ c₁₂ … c_1n c_1,n+1]

[c₂₁ c₂₂ … c_2n c_2,n+1]

C= [… ]=C⁽⁰⁾

[c_n1 c_n2 … c_nn c_n,n+1]

Przekształcamy układ do układu trójkątnego (przez odejmowanie od równań 2,3…n równanie 1 podzielone przez a₁₁ i pomnożone przez a_i1, gdzie i-numer równania).

{ c₁₁x₁+c₁₂x₂+…+c_1nx_n=c_1,n+1

{ c⁽¹⁾₂₂x₂+…+c⁽¹⁾_2nx_n=c⁽¹⁾_2,n+1

{ …=…

{ c^(n-1)_nnx­_n=c^(n-1)_n,n+1

Indeksy górne oznaczają numer operacji (0), (1)

{p=1,2,…,n

{{i=p+1,p+2…,n

{{c^(p)_ij=c^(p-1)_ij – (c^(p-1)_ip/c^(p-1)_pp)* c^(p-1)_pj

j=p+1,p+2,…,n+1

Najlepiej jest dzielić przez możliwie duże liczby, gdyż dzielenie przez liczby bliskie 0 może prowadzić do błędu. Aby to zrobić sortujemy równania w taki sposób, aby dzielić przez największą dostępną; wprowadzamy sortowanie na każdym etapie. Jest to wybór elementu głównego.

UWARUNKOWANIE ZADAŃ OBLICZENIOWYCH

Badanie uwarunkowania zadania obliczeniowego związane jest z analizą wpływu zaburzeń dawnych na zaburzenie rozwiązania.

Weźmy pod uwagę następujący układ równań wyrażony w postaci macierzowej: AX=Y (1)

Przeanalizujemy wpływ prawej strony na rozwiązanie, tj. wpływ ΔY na względny błąd ΔX.

Mamy więc A(X+ΔX)= (Y+ΔY) (2) Z równań (1) i (2) wynika: ΔX=A^-1ΔY , stąd i z równia (2) wynika nierówność: ||ΔX||<||A^-1|| ||ΔY|| oraz: ||ΔX||/||A||<=||X|| i ||X||=||A^-1Y|| z tych zależności uzyskujemy ||ΔX||/||X||<=(||A^-1|| ||Y||)/||A^-1Y|| ||ΔY||/||Y||<=||A^-1||||A|| ||ΔY||/||Y|| (3) Liczba będąca mnożnikiem po prawej stronie nierówności (3) jest uwarunkowaniem układu równań (condition number) cond(A)=||A||||A^-1||

Cond(A)=λmax(A)/λmin(A)

Cond(A)>=1 wiec im mniejsze wartości uwarunkowania macierzy cond(A) (im uwarunkowanie bliższe 1) tym lepiej, jesli jest ono dużą liczbą np 10¹⁰ należyw wybrać metodę regularyzacji.

METODY DOKŁADNE – pobierając dane wejściowe dążą do podania rozwiązania finalnego (pętle występują wewnątrz).

M. ITERACYJNE-najpierw proponujemy szereg rozwiązań a następnie otrzymujemy przybliżenia (kolejno). Otrzymujemy ciąg quasi rozwiązań, każdy kolejny wektor będzie lepszym rozwiązaniem, zaprzestajemy rozwiązań, zaprzestajemy rozwiązań, kiedy otrzymany wynik spełnia dane warunki.

W przypadku bardzo dużych układów równań opłaca się wykorzystywać metody iteracyjne prowadzące do ciągu wektorów X⁰,X¹,…Xⁿ określanych z równania macierzowego: Xⁱ⁺¹=WXⁱ+Z

Jeśli ciąg przybliżeń (Xⁿ) zmierza do granicy X to proces przybliżeń nazywamy zbieżnym z metodą efektywną. Np. z układu równań:

{3x₁+2x₂+1x₃=6

{x₁+2x₂+x₃=4

{2x₁+1x₂+x₃=1

Możemy uzyskać następujące równanie macierzowe.

{x₁} [0 -2/3 -1/3] {x₁} {2}

{x₂}=[-1/2 0 ½ ] {x₂} +{2}

{x₃} [-2 -1 0 ] {x₃} {1}

Metody obliczania macierzy

[3 2 1] [3 0 0] [0 2 1]

A=[0 2 0] = [0 2 0]+[1 0 -1]= D+R

[2 1 1] [0 0 1] [2 1 0]

Zatem AX=B => DX+RX=B => DX=-RX+B

Zakładając, że det(D)różny od0 możemy pomnożyć ostatnie równanie przez D^-1otrzymując X=-D^-1 RX+D^-1B czyli X=WX+Z gdzie W=-D^-1R12+D^-1B Poprawne rozwiązanie równania X=WX+Z otrzymujemy, gdy λmax[W]<1. Aby spełnić ten warunek wystarczy aby ||W||<1.

NORMY MACIERZY

-maksimum sumy wartości bezwzględnych z elementów kolumn analizowanej macierzy ||W||=max_iϵ<1,n>sumaⁿ_i=1|W_ij|

-maksimum sumy wartości bezwzględnych z elementów wierszy analizowanej macierzy

||W||=max_jϵ<1,n>sumaⁿ_i=1|W_ij|

-pierwiastek kwadratowy z sumy kwadratów elementów analizowanej macierzy ||W||=pierwiastek(sumaⁿ_i=1 sumaⁿ_j=1 W²_ij)

Metody ITERACYJNE:

Iteracja prosta xⁱ⁺¹=wxⁱ+z, i=0,1… xⁱ⁺¹=u oraz xⁱ=x

Np. {u_i=0,2x₂+0,3x₃+1 czyli {u₁} [0 0,2 0,3] {x₁} {1}

{ u_i=0,4x₁+0,5x₃+2 {u₂} = [0,4 0 0,5] {x₂}+{2}

{ u_i=0,1x₁+0,6x₂+3 {u₃} [0,1 0,6 0] {x₃} {3}

Metoda Gaussa-Siedla istota logarytmu polega na wykorzystaniu doliczonych k pierwszych składowych wektora niewiadomych xⁱ⁺¹ do oblicznia k+1 składowej itd. Xⁱ⁺¹=w_u Xⁱ + W_l Xⁱ⁺¹+z i=0,1…

W_u-macierz trójkątna górna W_l-macierz trójkątna dolna; gdzie W=W_u+W_l , W_ij=0 dla i=1(1)n

Metoda nadrelaksacji jest ulepszeniem metody Gaussa-Siedla mającym na celu poprawić zbieżność iteracji. Polega na sukcesywnym wprowadzaniu w miejsce składowych n_i po prawej stronie układu wartości x_i+α(u_i-x_i) gdzie α to wsp. Nadrelaksacji Xⁱ⁺¹=W_uXⁱ+W_i[Xⁱ+α(Xⁱ⁺¹-Xⁱ)]+z gdzie αϵ<1,2> oraz α_opt=1,8, dla α=1 mamy metodę G-S

INTERPOLACJA funkcji – jest to metoda obliczania przybliżonej wartości funkcji w danym punkcie przedziału <x₀,x_n>, gdy znana jest wartość tej funkcji w przedziałach x₀, x₁, x₂,…,x_n zwanych węzłami interpolacji y₀=f(x_0­), y₁=f(x₁), …. y_n=f(x_n) zadaniem interpolacji jest oblicznie przybliżonej wartości funkcji y=f(x) w punktach wewnątrz przedziału <x_0, x_n> nie będących węzłami interpolacji.

INTERPOLACJA WIELOMIANOWa

Daną funkcję przybliża się przy pomocy wielomianu uogólnionego w postaci:w(x)=sumaⁿ_i=1a_i ϕ_i(x) wektor funkcji: φ(x)=[ϕ₀(x), ϕ₁(x), …ϕ_n(x)] Wektor współczynników funkcji interpolującej: A^T=[a₀,a₁,…,a_n]^T wtedy: w(x)=φ(x) A - (wielomian uogólniony,>> XA=Y, W_(x)= ɸ_(x)A=ɸ_(x)X^-1Y (<<wielo interpolacyjny) gdzie: ɸ_(x)-macierz bazowa, X^-1-m. interpolacyjna, Y- wektor wartości funkcji w węzłach

INTERPOLACJA Lagrange’a

W tej interpolacji dla n+1 węzłów interplacji (x_i­,y_i), i=0(1)n przyjmuje się następujące funkcje bazowe ϕ_i(x)=Πⁿ_{j=0 ^ j≠i} (x-x_j) wielomian interpolacyjny ma postać w(x) =a₀ϕ₀(x)+ a₁ϕ₁(x)+…+ a_nϕ_n(x) Współczynnik wielomianu L wyznacza się z równania XA=Y, gdzie X jest macierzą diagonalną postaci:

[(x₀-x₁)….(x₀-­­­­x_n)…. 0 ]

X=[ …. …. ]

[ 0 (x_n-x₀)…(x_n-x_n-1)] zatem: A=X^-1Y stąd: a_i=y_i/((x_i-x₀)…. (x_i-x_i-1)(x_i-x_i+1)… (x_i-x_n)) Inter. wielom. Nie jest zbyt efektywna ponieważ macierz X jest macierza pełną (błędy przy odwracaniu) Macierz X nie jest zawsze dobrze uwarunkowana(może być osobliwa) Wzór interpolacyjny ma postać: w(x)=sumaⁿ_i=0 y_iL_i(x) gdzie y_i=f(x_i) jest znaną rzędną i-tego węzła interpolacji, L_i(x) jest i-tym wielomianem Lagrange’a zdefiniowanym jako L_i(x)= Πⁿ_{j=0 i j≠i} (x-x_j)/(x_i-x_j) więc możemy zauważuc, że L_i(x_k)=0 dla k≠i oraz k=1,2,3..n a dla pozostałych przypadków L_i(x_i)=1.

Jak wyznacza się współczynniki

APROKSYMACJA FUNKCJI

Aproksymacja funkcji polega na wyznaczeniu dla danej funkcji f(x) takich funkcji F(x), które w określonym sensie najlepiej przybliżają funkcje f(x). Stosujemy je najczęściej , gdy nie znamy funkcji f(x), lub gdy jest podana bardzo złożonym wzorem, który chcemy jakoś uprościć. Rodzaje aproksymacji: 1. Interpolacyjna 2.jednostajna 3. Średniokwadratowa.

APROKSYMACJA INTERPOLACYJNA podobnie jak w przypadku interpolacji żądamy spełnienia warunku aby funkcja f(x) i szukana funkcja F(X) przyjmowały te same wartość dla tych samych argumentów, czyli muszą przechodzić przez te same punkty węzłowe. Czasem chcemy też, żeby pochodne równały się w tych punktach, jeśli pochodne zostaną tam zadane.

APROKSYMACJA JEDNOSTAJNA. W przypadku takiej aproksymacji funkcję f(x) przybliżamy taką funkcją F(x) , która daje najmniejsze maksimum różnicy między F(x) a f(x) w całym przedziale <a,b>: max_xϵ<a,b>|F(x)-f(x)|=min Wielomian o dowolnie małym odchyleniu od funkcji f(x) na przedziale <a,b> możemy wyznaczyć zgodnie z twierdzeniem weierstrassa, jednak nie jest powszechnie znana metoda która pozwala jednoznacznie wyznaczyć wielomian. R=<całka od a do b> [F(x)-f(x)]^2 dx, gdzie R jest sumaryczną odchyłką często nazywaną residuum. Zatem w przypadku tej aproksymacji minimalizowane jest pole odchylenia funkcji F(x) od f(x). Gdy f(x) zadana jest tylko w określonych punktach (x_i, f(x­_i)) używamy kryterium : R=<suma> [F(x_i)-f(x_i)]^2 =min

METODA NAJMNIEJSZYCH KWADRATÓW

Zakładamy, że mamy wyznaczyć współczynnik p_i gdzie i=0(1)k wielomianu aproksymującego zadaną funkcję z=z(u) wg. Kryterium najmniejszych kwadratów z=p₀ϕ₀(u)+ p₁ϕ₁(u)+…+ p_kϕ_k(u) Aproksymowana funkcja zadana jest punktami (u_i,z_i) gdzie i=1(1)n, na podstawie zmienności funkcji przyjmujemy: φ(u)=[ ϕ₀(u), ϕ₁(u),…, ϕ_n(u)], kryterium optymalnego doboru parametrów ma postać: R<suma i=0 do n> [p₀ϕ₀(u_i)+ p₁ϕ₁(u_i)+…+ p_kϕ_k(u_i)-z_i]²=min

{∂R/∂p₀=2 <suma i=0 do n> [p₀ϕ₀(u_i)+ p₁ϕ₁(u_i)+…+ p_kϕ_k(u_i)-z_i] ϕ₀(u_i)

{…

{∂R/∂p_k=2 <suma i=0 do n> [p₀ϕ₀(u_i)+ p₁ϕ₁(u_i)+…+ p_kϕ_k(u_i)-z_i] ϕ_k(u_i)

Po przyrównaniu pochodnych do 0 otrzymujemy układ k+1 równań liniowych:

{p₀<suma>ϕ₀²(u_i)+p₁<suma>ϕ₀²(u_i) ϕ₁(u_i)+…+p_k <suma>ϕ₀(u_i) ϕ_k(u_i)=<suma>ϕ₀(u_i)z_i

{..

{p₀<suma>ϕ₀(u_i) ϕ_k(u_i)+p₁<suma>ϕ₁(u_i) ϕ_k(u_i)+…+p_k <suma> ϕ_k²(u_i)=<suma>ϕ_k(u_i)z_i

Układ ten można zapisać w postaci U^TUP=U^TZ prowadzając uogólnioną macierz danych wejściowych U i wektor prawych stron:

[ϕ₀(u₁) ϕ₁(u₁) … ϕ_k(u₁)] [z₁]

U=[ ϕ₀(u₂) ϕ₁(u₂) … ϕ_k(u₂)] Z= [z₂]

[.. ] […]

[ϕ₀(u_n) ϕ₁(u_n) … ϕ_k(u_n)] [z_n]

Przy założeniu det[U^TU]≠0 równanie można przekształcić P=(U^TU)^-1U^TZ Prawa strona jest znana i przy spełnieniu w.w. warunku można obliczyć współczynniki p_i wielomianu aproksym.

ROZWIĄZYWANIE RÓWNAŃ NIELINIOWYCH

METODA PRZESZUKIWANIA

przedział <a,b> dzielimy na n równych części i obliczmy wartości y₁ = f(x_i) Szukamy podprzedziału, dla którego f(a), f(b)<0 jeżeli funkcja ani razu nie zmienia znaku uważamy że nie ma pierwiastka w przedziale <a, b>. Jeżeli znajdziemy przedział izolacji pierwiastka <a_n, b_n>

to zawężamy go dalej aż osiągniemy zadowalającą dokładność tj: |f(x₀)\ <= <epsilon> dla a_n<=x₀<=b_n Zaletą tej metody jest jej prostota. Wady to: -krok może okazać się za duży i można pominąć jakiś pierwiastek -nie znajduje się pierwiastków w których funkcja osiąga ekstrema lokalne -zwiększenie dokładności osiąga się przez zmniejszenie kroku, co zwieksza czas obliczeń.

METODA BISEKCJI

jeżeli znajdzimy przedział izolacji pierwiastka <a_n, b_n> to zwężanie przedziału izolacji możemy wykonać metoda bisekcji: -jako pierwsze wyznaczamy ciągi przybliżeń przyjmujemy: x1=a, x2=b ,

-kolejne przybliżenie to x₁=0.5(x_i-1+x_k) gdzie k ϵ <1, i-2> , i>2

przy czym tak dobieramy k żeby |x_i-x_i-1| =~ |x_i-x_k| oraz F(x_i-1)F(x_k)<0

f’(b)=tgα tgα=f(b)/Δ Δ=f(b)²/tgα=f(b)/f’(b)

METODA SIECZNYCH

Przybliżenie przyjęte w tej metodzie to punkt przecięcia siecznej z osią odciętych. Jeżeli nie uzyska się żadanej dokładności do dalszych obliczeń jako granicę przedziału przejmuje się przybliżoną wartość pierwiastka.

Mamy: m=(f(b)-f(a))/(b-a) i y=mx+u i m=(0-f(a))/(c-a) stąd uzyskujemy (f(b)-f(a))/(b-a)=-f(a)/(c-a) zatem : c= a+(b-a) f(a)/(f(a)-f(b))

Zazwyczaj zaczynamy od metody przeszukiwania aż znajdziemy przedział izolacji perwiastka f(a) f(b)<0 następnie wykonujemy np. metodę siecznych

Metoda stycznych (Newtona)
Polega na zastąpieniu łuku krzywej y=fx w danym przedziale, w ktorym znajduje się dokładnie 1 pierwiastek rownania fx=0 odcinkiem stycznej poprowadzonej do tej krzywej w pewnym punkcie tego przediału.......
Rozwiązanie równania F(x) = 0 w przedziale [a,b] jest

przybliżone ciągiem miejsc zerowych stycznych do

funkcji F(x).
Równanie stycznej w punkcie xi−1:

y-F(x_i−1) =F' (x_i−1)(x-x_i−1)
Wybór pierwszego przybliżenia x1 = a lub x1 = b:
F'(x) * F′′(x) > 0 x­₁=b
F'(x) * F′′(x) < 0 x­₁=a
Jeżeli druga pochodna w przedziale izolacji nie ma stałego

znaku, to proces iteracyjny może być rozbieżny
1. f'(x)>0 , f''(x)<0
2. f'(x)<0 , f''(x)>0
3. f'(x)>0 , f''(x)>0
4 .f'(x)<0 , f''(x)<0

Kwadratury interpolacyjne

Rozpatrujemy funkcję f(x) ciągłą i ograniczoną w przedziale

domkniętym [a, b].

Przedział [a, b] dzielimy na skończoną liczbę podprzedziałów,

wyróżniając na osi x zbiór punktów:
a= x­­₀<x₁<x₂<...<x_i+1<...<x_n­=b
Punkty xi, i = 0, 1, ..., n tworzą siatkę o stałym kroku (z reguły):
x_i+1-x_i=h=const
Z własności całki oznaczonej wynika, że:

Istota metody kwadratur interpolacyjnych:
---Przybliżenie funkcji podcałkowej f(x) w przedziale [xi, xi+1]

lub przedziale odpowiednio poszerzonym wzorem

interpolacyjnym

---W(x) – wielomian interpolacyjny
---Wyprowadzenie wzorów przybliżających wartość

całki w przedziale [a, b]
---Wielomian interpolacyjny
I. wzór Newtona

Całkowanie numeryczne
polega na obliczaniu wartości całki oznaczonej w przedziale (a,b) opierając sie na interpretacji geometrycznej calki oznaczonej mozna wyprowadzic wzror na obliczenie dowolnj calki z zadana z gory dokladnoscia.
Metoda prostokatow teoria
Załóżmy, że chcemy obliczyć całkę z funkcji f(x) w przedziale <x_p; x_k>. Definicja całki oznaczonej Riemana, mówi nam, że wartość całki równa jest sumie pól obszarów pod wykresem krzywej w zadanym przedziale całkowania. Sumę taką możemy obliczyć w przybliżeniu dzieląc obszar całkowania na n równych części. Dla każdej takiej części możemy wyznaczyć prostokąt, który w przybliżeniu będzie odpowiadał polu obszaru pod wykresem krzywej. Jak widać na schemacie poniżej, dla funkcji rosnącej wartości tych przybliżeń będą większe niż w rzeczywistości - nadmiar powoduje część prostokąta znajdująca się ponad wykresem krzywej - dwa pierwsze prostokąty na schemacie. Natomiast dla funkcji malejącej wartości przybliżeń będą mniejsze niż rzeczywiste pole pod wkresem - niedomiar powoduje część pola znajdująca się nad wyznaczonym prostokątem - ostatni prostokąt na schemacie.
Jak już wpomnieliśmy przedział całkowania <x_p; x_k> podzielimy na n równych części. Szerokość każdej z nich wynosić będzie zatem:
dx = ( x_k - x_p ) / n.
Taka też będzie szerokość każdego prostokąta przybliżającego nam wartość całki w zadanym przedziale. Wysokość każdego z prostokątów wynosić będzie:
f( x_i ) dla i = 1, 2, ..., n, gdzie x_i = x_p + i*dx.
Całkę w zadanym przedziale uzyskamy dodając do siebie pola wszystkich tych prostokątów, wynosić będzie ona zatem:
dx * f( x₁ ) + dx * f( x₂ ) + ... + dx * f( x_n ) = dx * (f( x₁ ) + f( x₂ ) + ... + f( x_n )).
Warto zauważyć, iż im większa liczba przedziałów n z tym większą dokładnością wyznaczymy interesującą nas całkę.

Metoda trapezów
Załóżmy, że chcemy obliczyć całkę z funkcji f(x) w przedziale <x_p; x_k>. Definicja całki oznaczonej Riemana, mówi nam, że wartość całki równa jest sumie pól obszarów pod wykresem krzywej w zadanym przedziale całkowania. Sumę taką możemy obliczyć w przybliżeniu dzieląc obszar całkowania na n równych części. Dla każdej takiej części możemy wyznaczyć trapez, który w przybliżeniu będzie odpowiadał polu obszaru pod wykresem krzywej.

Jak już wspomnieliśmy przedział całkowania <x_p; x_k> podzielimy na n równych części. Szerokość każdej z nich wynosić będzie zatem:
dx = ( x_k - x_p ) / n.
Taka też będzie wysokość każdego z trapezów. Podstawy i-tego trapezu będą wynosić odpowiednio:
f( x_i-1 ) oraz f( x_i ) dla i = 1, 2, ..., n, gdzie x_i = x_p + i*dx.
Pole i-tego trapezu zgodnie ze wzorem wynosić będzie:
P_i = ((f( x_i-1 ) + f( x_i )) / 2) * dx
Całkę w zadanym przedziale uzyskamy dodając do siebie pola wszystkich wyznaczonych trapezów, wynosić będzie ona zatem:
((f( x₀ ) + f( x₁ )) / 2) * dx + ((f( x₁ ) + f( x₂ )) / 2) * dx + ... + ((f( x_n-1 ) + f( x_n )) / 2) * dx =
= (dx/2) * (f( x₀ ) + 2f( x₁ ) + 2f( x₂ ) + ... + 2f( x_n-1 ) + f( x_n )) =
= dx * (f( x₀ ) / 2 + f( x₁ ) + f( x₂ ) + ... + f( x_n-1 ) + f( x_n ) / 2)
W praktyce w pętli dodajemy do siebie wszystkie wartości funkcji od 1 do n-1, a potem dwie wartości brzegowe podzielone przez dwa. Całość mnożymy przez dx i otrzymujemy w ten sposób wynik.
Takie postępowanie daje wyniki lepsze niż całkowanie metodą prostokątów, ale i tutaj otrzymany wynik nie będzie zawsze idealny - zakładamy przecież, że funkcja w obrębie przedziałów jest liniowa, co w ogólności nie musi być prawdą. Na schemacie powyżej widać, że metoda dość dobrze (ale nie idealnie) odwzorowywuje naszą przykładową funkcję w dwóch pierwszych przedziałach, natomiast w ostatnim przedziale widać wyraźnie różnicę pomiędzy polem pod wykresem a wyznaczonym trapezem. Warto zauważyć, iż im większa liczba przedziałów n z tym większą dokładnością wyznaczymy interesującą nas całkę.

Wzór Simpsona
metod przybliżania wartości całki oznaczonej funkcji rzeczywistej
Metoda ma zastosowanie do funkcji stablicowanych w nieparzystej liczbie równo odległych punktów (wliczając końce przedziału całkowania). Metoda opiera się na przybliżaniu funkcji całkowanej przez interpolację wielomianem drugiego stopnia. Geometrycznie metoda ta odpowiada zastąpieniu w każdym z kolejnych k przedziałów zmiennej x łuku wykresu funkcji y=f(x) łukiem paraboli przeprowadzonej przez trzy kolejne węzły interpolacji (punkty wykresu o znanych współrzędnych) odpowiadające początkowi, środkowi i końcowi kolejnego przedziału.

Oznacza to:
f(x) w przedziale [xi, xi+1] jest przybliżona wielomianem

interpolacyjnym ograniczonym do trzech pierwszych

składników
Przedział [a, b] dzieli się na parzystą ilość podprzedziałów
Pole pojedynczego trapezu krzywoliniowego:

Wzór Simpsona:

TEORIA SIMPSON
Załóżmy, że chcemy obliczyć całkę z funkcji f(x) w przedziale <x_p; x_k>. Definicja całki oznaczonej Riemana, mówi nam, że wartość całki równa jest sumie pól obszarów pod wykresem krzywej w zadanym przedziale całkowania. Sumę taką możemy obliczyć w przybliżeniu dzieląc obszar całkowania na n równych części. W metodach prostokątów i trapezów zakładaliśmy, że przybliżenie funkcji w przedziale jest funkcją liniową (przybliżenie odwzorowywało funkcję na odcinek w obrębie przedziału). W metodzie Simpsona w każdym takim przedziale będziemy przybliżać funkcję dla, której obliczamy całkę przy pomocy paraboli.
Kwadratury Gaussa
Rozpatrujemy funkcję f(x) ciągłą i ograniczoną w przedziale domkniętym [a, b].
Pierwszy krok:
Sprowadzenie całki

do postaci znormalizowanej:
Normalizacja
Podstawienia:

Czyli:

Metoda łamanych Eulera

Zakładamy że ∫_xi^xi+1F[x,y(x)]dx=hiF(xi,yi) gdzie h_i=x_i+1-x_i , czyli że f-cja F[x,y(x)]=const na całym przedziale <x_i,y_i> zatem (rys1) y_xi+1=y_i+h_iF(xi,yi) ˄ y(x₀)=y₀ zazwyczaj h_i=h=const

Ulepszenia metody Eulera: Modyfikacja 1

Zakładamy że styczna do łuku rozpiętego na punktach A(xi,yi), B(x_i+1,y_i+1) w punkcie x*=0,5(x_i,x_i+1) jest w przybliżeniu równoległa do cięciwy AB (rys 2)poszczególne kroki prowadzące do rozw mają postać x*=0,5(x_i+x_i+1)˄y*yi+0,5hF(xi,yi), m*=F(x*,y*)˄y_i+1=yi+hm*

Modyfikacja2 – metoda Eulera-Cauchego

Zakładamy że współczynnik kier. siecznej AB jest średnią arytmetyczną współczynników kier. stycznych do poszukiwanej funkcji w punktach A i B (rys3 )poszczególne kroki prowadzące do rozw mają postać x*=x_i+1=x_i+h˄y*=yi=hF(xi,yi), m*=F(x*,y*)˄y_i+1=yi+0,5[F(xi,yi)+m*]h

Błędy metody i błędy zaokrągleń

Dokładność metody łamanych wyraźnie rośnie przy mniejszym kroku h (rys 4)Jednakże wtedy wymagana jest realizacja znacznie większej liczby kroków co zwiększa tzw. „błąd zaokrąglenia” istnieje wiec pewien optymalny krok h* ale niestety nie ma ogólnej metody pozwalającej na jednoznaczne obliczenie tego kroku

Wzory Rungego-Kutty

W ogólnym przypadku metoda R-K polega na doborze współczynników a₁,a₂...,b₁,b₂... oraz A₁,A₂..., tak aby : y_i+1=y_i+hϵ_jm_jA_j gdzie:

m₁=F(x_i,y_i), m₂=F(x_i+a₁h,y_i+b₁hm₁), m₃=F(x_i+a₂h,y_i+b₂hm₂+b₃hm₃), m₄=….

Ulepszona metoda łamanych –modyfikacja 1

to: a₁=0,5 A₁=0 a₂=1 b₁=0,5

Ulepszona metoda łamanych – modyfikacja 2

A₁=1 A₁=0,5 A₂=0,5 b₁=1

Graficzna ilustracja do metody RK4 (rys5) poszczególne wsp. m₁,m₂,m₃,m₄, odpowiadają stycznym 1(a),1(b),1(c),1(d).

Wzory RK 4 rzędumają postać:

m₁=F(xi,yi), m₂=F(xi+0,5h, yi+0,5hm₁), m₃=F(xi+h,yi+hm₃), m₄=F(xi+h, yi+hm₃), y_i+1=yi+h(m₁+2m₂+2m₃+m₄)/6

Wzory RK3 rzędu mają postać:

m₁=F(xi,yi), m₂=F(xi+0,5h, yi+0,5hm₁), m₃=F(xi+h, yi+hm₁+2hm₂), y_i+1=yi+h(m₁+hm₂+m₃)/6

Przeanalizujemy metodę RK z innego pkt widzenia: Y_i+1=y_i+∫_xi^xi+1F[x,y(x)]dx˄F(x,y)=y’˄y(x₀)=y₀ zgodnie z tymi warunkami dla pierwszych dwóch pkt całkowania równanie różniczkowe możemy zapisać: y1=y0+∫_x0^x1F[x,y(x)]dx