Metody numeryczne

Jeśli chciałbyś wykorzystać zamieszczone w tym dziale dane (tekst, rysunki, skrypty JS, procedury i funkcje) w celu zamieszczenia na stronie www lub w publikacjach (prace, artykuły) musisz się uprzednio ze mną skontaktować podając w jakim celu chcesz je wykorzystać.
[e-mail]

W dziale Programy znajduje się zbiór 5 programów z metod numerycznych do ściągnięcia.

Jeżeli zastanawiasz sie co to są metody numeryczne odpowiedź znajduje się poniżej. Jeżeli już wiesz co to są metody numeryczne możesz pominąć ten wstęp i przejść do przeglądania zawartości stron.
Czasami zdarza się tak że nie wszystkie zadania matematyczne da się rozwiązać w sposób analityczny lub ich rozwiązanie zajęło by bardzo dużo czasu, to stosuje się w tedy metody numeryczne. Dzięki nim można uzyskać rozwiązanie z zadaną dokładnością już po kilku sekundach. W większości przypadków są to metody przybliżone ponieważ opis matematyczny, który zastępują jest bardzo złożony. Zadaniem tych metod numerycznych jest uproszczenie tego zapisu lub zastąpienie go innym (lub kilkoma) prostszym, który można zapisać w programie i który daje rozwiązanie z zadanym przybliżeniem. Stosuje się je także w przypadku gdy proces rozwiązywania zadania (o prostym zapisie matematycznym) trwałby zbyt długo (nawet z zastosowaniem komputera).
Zdarza się jednak że opis matematyczny do danego zadania nie jest skomplikowany w tedy stosuje się metody dokładne, które pozwalają uzyskać rozwiązanie z dokładnością na jaką pozwala nam komputer.
W tym dziale będę zamieszczał procedury i funkcje wraz w ich opisem wykorzystywane do obliczeń w matematyce, fizyce, mechanice i innych zbliżonych dziedzinach.
Serwis przeznaczony jest dla osób znających chociaż podstawy języka Pascal ponieważ nie będę wyjaśniał zagadnień i poleceń tego języka.

Całkowanie numeryczne - Metoda trapezów

Bardzo prosta i szybka funkcja do obliczania całek oznaczonych (w tym przypadku pola pod daną funkcją). Metodę tę rozbiłem na dwie funkcje - tę właściwą i funkcję obliczającą wartość funkcji podcałkowej. Umożliwia to umieszczenie w niej kilku różnych funkcji i ich wybór np. za pomocą case. Metoda opiera się na skończonej liczbie pod przedziałów, a ilustruje ją poniższy obrazek. Im większa liczba pod przedziałów tym dokładniejsze będą obliczenia - należy pamiętać aby nie przesadzić ponieważ obliczenia będą trwać zbyt długo, a uzyskany wynik nie będzie wcale dokładniejszy, wręcz przeciwnie można w ten sposób pogorszyć dokładność. Poza tym należy dobrać tak n aby długość pod przedziału była liczą wymierną - najlepiej n=10,100,1000 itd. W ten sposób wyeliminuje się błędy powstałe w czasie obliczeń wartości funkcji (błąd zmiennej).

function calka(n: integer;a,b: real): real;//n - liczba pod przedziałów, a,b - granice całkowania b>a
var i: integer; //i - zmienna iteracyjna
  x0,fx1,fx2,d,wyn: real; //x0,fx,fx2 - pomocnicze zmienne
begin
  d:=(b-a)/n; //wyznaczanie wielkości pod przedziału
  wyn:=0; //zerowanie zmiennej wyniku
  x0:=a; //pierwszy punkt obliczeń
  fx2:=Abs(f(x0)); //obliczanie wartości w "drugim" punkcie i jeżeli funkcja ma wartość ujemną następuje zmiana na dodatnią
  for i:=1 to n do
  begin
    fx1:=fx2; //wartość funkcji w pierwszym punkcie pod przedziału
    x0:=x0+d; //zwiększenie wartości x0 o wielkość pod przedziału
    fx2:=Abs(f(x0)); //wartość funkcji w drugim punkcie pod przedziału
    wyn:=wyn+((fx1+fx2)*d*0.5);
  end;
  calka:=wyn;
end;

function f(x: real):real;
begin
  f:=exp(-sqr(x));
end;

Pod poniższym linkiem znajduje się przystosowana dla języka JavaScript procedura całkowania. Można w niej dokonać obliczenia powyższej całki.
Przykład
Jak widać dla potrzeb metod numerycznych nie jest wymagane specjalne oprogramowanie obsługujące różne języki programowania. Wystarczy przeglądarka z obsługą JavaScript np. IE 4 i nowszy lub Netscape

Całkowanie numeryczne - Złożony wzór paraboliczny

Całkowanie to odbywa się przez przybliżenie funkcji podcałkowej na danym odcinku łukiem paraboli (w powyższym przykładzie przybliżaliśmy odcinkiem prostej). Najważniejszym zadaniem w tym przypadku jest wyznaczenie równania paraboli przechodzącej przez 3 punkty naszej funkcji, a następnie jego z całkowanie na danym odcinku. Dokonać można tego w dwojaki sposób: po przez rozwiązanie układu równań i wyznaczenie współczynników a, b, c równania paraboli przechodzącej przez trzy punkty, a następnie jej z całkowanie w danym przedziale lub po przez z całkowanie w danym przedziale wielomianu interpolacyjnego Lagrangea.
Uzyskujemy w ten sposób poniższy wzór:

function calka2(n: integer;a,b: real): real;//n - liczba pod przedziałów, a,b - granice całkowania b>a
var i: integer; //i - zmienna iteracyjna
  x0,fx1,fx2,fx3,wyn,h: real; //x0,fx,fx2,fx3 - pomocnicze zmienne, wyn- wynik
begin
  h:=(b-a)/n; //wyznaczanie wielkości pod przedziału
  wyn:=0; //zerowanie zmiennej wyniku
  x0:=a; //pierwszy punkt obliczeń
  fx3:=Abs(f(x0)); //obliczanie wartości w "trzecim" punkcie i jeżeli funkcja ma wartość ujemną następuje zmiana na dodatnią
  for i:=1 to n do
  begin
    fx1:=fx3; //wartość funkcji w pierwszym punkcie pod przedziału
    fx2:=Abs(f(x0+h/2)); //wartość funkcji w drugim punkcie pod przedziału
    fx3:=Abs(f(x0+h)); //wartość funkcji w trzecim punkcie pod przedziału
    wyn:=wyn+((fx1+4*fx2+fx3)*(h/6));
    x0:=x0+h; //zwiększenie wartości x0 o wielkość pod przedziału
  end;
  calka2:=wyn;
end;

function f(x: real):real;
begin
  f:=exp(-sqr(x));
end;

Przykład
Porównując obie metody zauważyć można iż większą dokładność daje wykorzystanie złożonego wzoru parabolicznego:
dla naszej funkcji i danych: a=0, b=2, n=10 otrzymujemy odpowiednio:
metoda 1 - 0.88183881
metoda 2 - 0.88208098
przy czym przybliżone rozwiązanie wynosi - 0.88208125

Metody numeryczne

Jeśli chciałbyś wykorzystać zamieszczone w tym dziale dane (tekst, rysunki, skrypty JS, procedury i funkcje) w celu zamieszczenia na stronie www lub w publikacjach (prace, artykuły) musisz się uprzednio ze mną skontaktować podając w jakim celu chcesz je wykorzystać.
[e-mail]

W dziale Programy znajduje się zbiór 5 programów z metod numerycznych do ściągnięcia.

Obliczanie pierwiastków równań nieliniowych - metoda bisekcji

Metoda ta pozwala w prosty sposób wyznaczyć z zadaną dokładnością pierwiastki równania nieliniowego. Przed przystąpieniem do obliczania należy poznać najpierw przebieg funkcji by w przybliżeniu określić granice przedziału, w którym znajduje się szukany pierwiastek (można je oszacować z większym lub mniejszym przybliżeniem). Metoda ta może być stosowana w każdym przypadku, w którym funkcja w granicach podanego przedziału zmienia znak (po przejściu przez miejsce zerowe). Metoda ta zawsze jest zbieżna (zbieżność liniowa tej procedury gwarantuje odnalezienia pierwiastka równania jednak kosztem ilości iteracji). Poniższy rysunek obrazuje proces poszukiwania pierwiastków. Niżej procedura poszukująca pierwiastki zadanego równania.

procedure rownanie;
var it: integer; //zmienna iteracyjna
  d,x1,x2,x0,a,b: real;
  koniec: boolean;
Label kon;
begin
//przed wywołaniem procedury lub w tym miejscu należy wczytać granice przedziału a i b oraz dokładność obliczeń d
  if (a>=b) or (f(a)*f(b)>0) then goto kon; //sprawdza poprawność podania granicy przedziału i w razie niezgodności przerywa obliczenia
  it:=0; x1:=a; x2:=b;
  repeat //obliczanie pierwiastka
    it:=it+1;
    x0:=(x1+x1)/2; //wyznaczenie środka przedziału x1 i x2
    if Abs(f(x0))<d then koniec:=true; //jeżeli wartość funkcji w punkcie x0 jest mniejsza niż zadana dokładność - został znaleziony pierwiastek równania, kończymy obliczenia
//przesunięcie jednej granicy przedziału do punktu x0 - poniższe warunki sprawdzają naszą funkcję i przesuwają jedną z granic przedziału x1,x2 do punktu x0
    if f(x0)>0 then
    begin
      if f(x1)>0 then x1:=x0;
      if f(x2)>0 then x2:=x0;
    end
    else
    begin
      if f(x1)<0 then x1:=x0;
      if f(x2)<0 then x2:=x0;
    end;
  until koniec;
//po zakończeniu obliczeń w tym miejscu można odczytać wyniki obliczeń: pierwiastek równania x0, liczba iteracji potrzebnych do wyznaczenia pierwiastka z zadaną dokładnością it oraz wartość funkcji w punkcie x0 f(x0)
  kon:
end;

function f(x: real):real;
begin
  f:=2*x*x+1-exp(x*ln(2));
end;

Przy wykorzystaniu tej procedury wyznaczyłem 2 z 3 rozwiązań równania 2*x*x+1-exp(x*ln(2))=0:
dla a=0,1 i b=1 oraz d=1e-8
x0=0,39928074777
f(x0)=-6,072e-9
it=25
dla a=6 i b=7 oraz d=1e-8
x0=6,35234489
f(x0)=6,02e-10
it=31
Trzecim rozwiązaniem tego równania jak łatwo można zauważyć jest x0=0.
Sprawdź sam: Przykład

Obliczanie pierwiastków równań nieliniowych - metoda iteracji prostej

Metoda iteracji prostej należy do tzw. metod iteracyjnych jednopunktowych. Oparta jest na twierdzeniu Banacha o punkcie stałym. Ciąg kolejnych przybliżeń dla tej metody jest następujący: x_k+1=q(x_k) (k=1,2,3,...). Przy czym q(x) jest równaniem równoważnym otrzymanym po zastąpieniu równania wyjściowego f(x)=0 => x=q(x). Pierwszy punkt przybliżenia przyjmujemy w jak najbliższym sąsiedztwie pierwiastka równania. Funkcja ta powinna być tak dobrana aby proces był jak najszybciej zbieżny. Ilustrację graficzną przedstawia poniższy rysunek prezentujący dwa warunki gdy metoda jest zbieżna i rozbieżna.

procedure rownanie;
var it: integer; //zmienna iteracyjna
  xn,x0,eps,bl: real;
Label kon;
begin
  //W tym miejscu procedury należy wprowadzić wielkość dopuszczalnego błędu bl oraz wartość początkową x0 (lub przy wywołaniu procedury);
  it:=0; eps:=1e4; //zerowanie zmiennej iteracyjnej, eps- zmienna przyrostu pomiędzy kolejnymi przybliżeniami rozwiązania - wartość przybliżona
  repeat
    xn:=q(x0); //kolejne wyższe (dokładniejsze) przybliżenie rozwiązania
    if Abs(xn-x0)>eps then //sprawdzamy czy przyrost pomiędzy kolejnymi przybliżeniami był większy niż poprzedni
    begin
      writeln('Proces iteracji jest rozbieżny');
      goto kon;
    end
    else eps:=Abs(xn-x0); //wyznaczamy aktualny przyrost pomiędzy kolejnymi przybliżeniami rozwiązania
    x0:=xn;// przepisanie zmiennej o wyższym (dokładniejszym) przybliżeniu tak aby przy kolejnym wykonaniu pętli wyznaczyć jeszcze dokładniej rozwiązanie równania
    it:=it+1;
  until eps<=bl;// zakończenie pętli
  //W tym miejscu odczytujemy: x0- pierwiastek równania, f(x0) oraz it
  kon:
end;

function f(x: real):real;
begin
  f:=2*x*x+1-exp(x*ln(2));
end;

function q(x: real):real;
begin
  q:=sqrt(0.5*exp(x*ln(2))-0.5);
end;

Za pomocą tej metody wyznaczyłem jedno z rozwiązań naszego równania:
dla pierwszego przybliżenia x0=1 oraz d=1e-8
x0=0,39928076
f(x0)=5,322e-9
it=32
Przy wyznaczaniu drugiego z rozwiązań równania proces iteracji był rozbieżny.

Obliczanie pierwiastków równań nieliniowych - metoda Newtona (stycznych)

Metoda ta należy do jednych z najszybciej zbieżnych metod służących do wyznaczania pierwiastków równań nieliniowych. Stopień zbieżności tej metody wynosi 2. Zbieżność tej metody została osiągnięta dzięki dodatkowemu nakładowi pracy jakim jest wyznaczenie pochodnej badanego równania. W niektórych przypadkach wyznaczenie pochodnej jest sprawą bardziej złożoną i zastosowanie tej metody nie jest opłacalne.
Na poniższym rysunku pokazany jest sposób w jaki dążymy do wyznaczenia pierwiastka równania.

  

procedure rownanie;
var it: integer; //zmienna iteracyjna
  xn,x0,eps,bl: real;
Label kon;
begin
  //W tym miejscu procedury należy wprowadzić wielkość dopuszczalnego błędu bl oraz wartość początkową x0
  it:=0; eps:=1e4; //zerowanie zmiennej iteracyjnej, eps- zmienna przyrostu pomiędzy kolejnymi przybliżeniami rozwiązania - wartość przybliżona
  repeat
    xn:=x0-f(x0)/fp(x0); //kolejne wyższe (dokładniejsze) przybliżenie rozwiązania
    if Abs(xn-x0)>eps then //sprawdzamy czy przyrost pomiędzy kolejnymi przybliżeniami był większy niż poprzedni
    begin
      writeln('Proces iteracji jest rozbieżny');
      goto kon;
    end
    else eps:=Abs(xn-x0); //wyznaczamy aktualny przyrost pomiędzy kolejnymi przybliżeniami rozwiązania
    x0:=xn;// przepisanie zmiennej o wyższym (dokładniejszym) przybliżeniu tak aby przy kolejnym wykonaniu pentli wyznaczyć jeszcze dokładniej rozwiązanie równania
    it:=it+1;
  until eps<=bl;// zakończenie pętli
  //W tym miejscu odczytujemy: x0- pierwiastek równania, f(x0) oraz it
  kon:
end;

function f(x: real):real;
begin
  f:=2*x*x+1-exp(x*ln(2));
end;

function fp(x: real):real;
begin
  fp:=4*x-exp(x*ln(2))*ln(2);
end;

Przy wykorzystaniu tej procedury wyznaczyłem 2 z 3 rozwiązań tego równania:
dla x0=1 oraz bl=1e-8
x0=0,399280756662
f(x0)=5,44161e-16//bardzo duża dokładność
it=6
dla x0=7 oraz bl=1e-8
x0=6,352344892434
f(x0)=5,256906e-14
it=5

Podsumowanie

Najlepiej w porównaniu tych trzech metod wypada metoda bisekcji mimo iż jest wolno zbieżna przy jej wykorzystaniu wiele się nie napracujemy. Na jej korzyść świadczy to że wyznaczymy wszystkie pierwiastki równania pod warunkiem że nie stanowić będą jednocześnie ekstremum dla tych funkcji.
Metoda stycznych jest najszybszą metodą jednak zanim uruchomimy program musimy mieć policzoną pochodną, a to nie zawsze jest łatwe. Metoda stycznych nie zawsze jest metodą zbieżną.
Najgorzej wypadła metoda iteracji prostej - dużo iteracji i mała dokładność i no fakt że nie udało się nam wyznaczyć za jej pomocą jednego z pierwiastków.

Metody numeryczne

Jeśli chciałbyś wykorzystać zamieszczone w tym dziale dane (tekst, rysunki, skrypty JS, procedury i funkcje) w celu zamieszczenia na stronie www lub w publikacjach (prace, artykuły) musisz się uprzednio ze mną skontaktować podając w jakim celu chcesz je wykorzystać.
[e-mail]

W dziale Programy znajduje się zbiór 5 programów z metod numerycznych do ściągnięcia.

Rozwiązywanie układów równań liniowych - Metoda eliminacji Gaussa

Rozwiązywanie układów równań liniowych jest najważniejszym zagadnieniem matematyki. Tylko proste i o liczbie równań nie większej niż 3 można bez problemu rozwiązać ręcznie. Większe i bardziej złożone zagadnienia powierza się komputerom (układy kilkunastu, setek a nawet tysięcy równań). Opracowano wiele metod służących do numerycznego rozwiązywania (metody dokładne jak i przybliżone). Ja przedstawię najprostszą a zarazem dokładną metodę rozwiązywania układów równań liniowych jaką jest metoda Eliminacji Gaussa.

Metoda ta jak sama nazwa wskazuje polega na eliminowaniu poszczególnych zmiennych oraz równań. Proces eliminacji prowadzony jest w znany sposób przez wyznaczanie "wartości" jednej zmiennej z jednego równania i wstawianie jej do pozostałych równań. Następnie prowadzony jest proces odwrotny, w którym wyznaczane są dopiero wartości poszczególnych niewiadomych układu równań.
Metoda w swojej normalnej postaci jest skuteczna prawie zawsze gdy elementy macierzy współczynników na jej głównej przekątnej nie są równe zeru. W przeciwnym przypadku stosowane są jej modyfikacje z częściowym i pełnym wyborem elementu macierzy. Poniższe ilustracje przedstawiają proces eliminacji oraz proces odwrotny z ich opisem.

Weźmy pod uwagę układ równań z n-niewiadomymi. Na początek wyznaczymy "wartość" pierwszej niewiadomej x₁.

 

Tak wyznaczoną wartość x₁ wstawiamy do pozostałych równań (w przykładzie podstawiłem do drugiego równania) i w ten sposób wyeliminowaliśmy jedną z niewiadomych tego układu równań, a tym samy jedno równanie mniej (nie jest ono usuwane - pozostaje w pamięci by móc wyznaczyć później z jego pomocą niewiadome układu równań). Powstaje w ten sposób nam układ z n-1 niewiadomymi o współczynnikach z indeksem 1 (indeks oznacza nam kolejne działania eliminacji) i n-1 równaniami (rysunek niżej).

Proces eliminacji powtarzamy do momentu aż pozostanie nam jedno równanie i jedna niewiadoma - czyli n-1 razy. Efektywność tego procesu zapewnia nam poniższy wzór (dla potrzeb programu wyrazy wolne b₁,...,b_n zapisujemy w n+1 kolumnie macierzy współczynników).

 

Teraz gdy otrzymaliśmy już jedno równanie z jedną niewiadomą nic nie stoi na przeszkodzie by ją wyznaczyć. Kolejne niewiadome od przedostatniej do pierwszej (kolejność odwrotna) wyznaczamy podstawiając do kolejnych równań (od przedostatniego do pierwszego) co dopiero wyznaczone wcześniej niewiadome.

W czasie dokonywania eliminacji jak i przy wyznaczaniu zmiennych występuje nam dzielenie. Zmusza nas ono do sprawdzenia czy współczynnik będący dzielnikiem jest różny od zera. W przeciwnym przypadku (gdy taka sytuacja nastąpi) musimy przerwać proces i zakończyć obliczenia. Jest to jedyna wada tej metody - elementy macierzy przed i po procesie eliminacji na jej głównej przekątnej muszą być różne od zera.
Po procesie eliminacji warunek ten sprawdzamy obliczając wyznacznik macierzy współczynników równania. Wyznacza się go mnożąc ze sobą elementy macierzy występujące na jej głównej przekątnej.

procedure elimgaussa;
var i,j,k,n: integer;
  sum,det: real;
  a: array [1..10,1..11] of real; //Macierz współczynników plus 1 kolumna dla wyrazów wolnych
  x: array [1..10] of real; //tablica niewiadomych
Label kon;
begin
{W tym miejscu procedury należy wczytać ilość równań n oraz dane do macierzy a (ewentualnie wczytać je wcześniej)}

  for k:=1 to n-1 do
  for i:=k+1 to n do
  for j:=k+1 to n+1 do
  if a[k,k]<>0 then //ten warunek sprawdza czy współczynniki nie są równe 0
  a[i,j]:=a[i,j]-(a[k,j]*a[i,k])/a[k,k] //jeśli nie
  else goto kon; //gdy tak przerywa obliczenia
{Gdy współczynnik będzie równy 0 nie wystąpi dzielenie przez 0 i program nam się nie wysypie.}

  det:=1;
  for i:=1 to n do det:=det*a[i,i]; //wyznaczamy wyznacznik macierzy współczynników

{Obliczony wyznacznik sprawdzamy czy nie jest zbyt mały. Poniższy warunek daje nam zabezpieczenie przed wystąpieniem błędów obliczeń - złe uwarunkowanie układu równań, lub jak wiadomo gdy wyznacznik jest równy zeru układ nie posiada rozwiązań. Jeśli jest wszystko ok to wyznaczamy niewiadome układu równań.}

  if Abs(det)>1e-8 then
  begin
    x[n]:=a[n,n+1]/a[n,n];
    for i:=n-1 downto 1 do
    begin
      sum:=0;
      for j:=i+1 to n do
      sum:=sum+a[i,j]*x[j];
      x[i]:=(a[i,n+1]-sum)/a[i,i];
    end;
  end;
//W tym miejscu możesz odczytać wyniki obliczeń
  kon:
end;

Ta niewielka procedura w zupełności pozwala na prowadzenie wielu obliczeń. Nawet gdy niektóre współczynniki są równe 0 wystarczy zamienić równania miejscami i spokojnie można kontynuować obliczenia. Przy okazji wyznaczamy także wyznacznik macierzy współczynników. Metoda ta należy do metod dokładnych i szybkich, a stosując dodatkowe modyfikacje jest bardzo efektywna.

Metoda eliminacji Gaussa z częściowym wyborem elementu podstawowego

Ta modyfikacja polega na wybraniu przed każdym krokiem eliminacji takiego równania, w którym pierwszy element jest największy co do wartości, w porównaniu z pozostałymi równaniami. Następnie sprawdzamy czy element ten jest równy zeru, jeśli tak to przerywamy dalsze obliczenia. Jeśli znaleźliśmy takie równanie to zamieniamy je z naszym bieżącym i po tym kroku dokonujemy kolejnej eliminacji. Zamieniając równania miejscami dokonujemy zamiany znaków w jednym z nich aby wyznacznik był równy z wyznacznikiem układu wejściowego Poniżej przedstawiam fragment kodu, który należy wstawić do poprzedniej procedury.

for k:=1 to n-1 do //po tym wierszu wklejamy poniższy kod

max:=0; w:=k;
for i:=k to n do
if max<Abs(a[i,k]) then //sprawdzamy czy dany współczynnik jest większy od max
begin //i jeśli tak to zapamiętujemy numer równania i oraz wartość tego współczynnika
  max:=Abs(a[i,k]);
  w:=i;
end;
if max=0 then goto kon;
if w>k then //jeżeli znaleziony został większy współczynnik to zamieniamy wiersze miejscami
for i:=k to n+1 do
begin
  pom:=a[k,i];
  a[k,i]:=-a[w,i];
  a[w,i]:=pom;
end;

//Dokonujemy jeszcze małej poprawki
if a[k,k]<>0 then //ten warunek z poprzedniej procedury możemy już usunąć ponieważ teraz wcześniej dokonujemy tego sprawdzenia.

Tak zmodyfikowana metoda pozwala na prowadzenie obliczeń w przypadkach gdy na głównej przekątnej macierzy współczynników występują (lub mogą wystąpić w czasie obliczeń) zerowe współczynniki.

Metody numeryczne

Jeśli chciałbyś wykorzystać zamieszczone w tym dziale dane (tekst, rysunki, skrypty JS, procedury i funkcje) w celu zamieszczenia na stronie www lub w publikacjach (prace, artykuły) musisz się uprzednio ze mną skontaktować podając w jakim celu chcesz je wykorzystać.
[e-mail]

W dziale Programy znajduje się zbiór 5 programów z metod numerycznych do ściągnięcia.

Różniczkowanie numeryczne - Różnica centralna

Jak wiadomo wartością pochodnej funkcji f(x) w danym punkcie jest współczynnik kierunkowy prostej stycznej do danej funkcji w tym właśnie punkcie.  Przedstawia  to  poniższy  rysunek. Zaznaczyłem  na  nim  także  odcinek  łączący  2 punkty  f(x-dx) oraz f(x+dx) - jest to tzw. różnica centralna. W metodzie tej wykorzystuje się przybliżenie że odcinek ten jest równoległy do stycznej w punkcie f(x_i), a więc kąt pochylenia tego odcinak i stycznej są takie same. W wyniku tej metody otrzymujemy szereg punktów odpowiadających wartością (w przybliżeniu) do wartości pochodnej w tych punktach F'(x) (do wartości dokładnej).
W przypadku gdy funkcja f(x) jest przedstawiona w postaci dyskretnej tzn. za pomocą pewnej ilości n punktów (x,f(x)) to w przypadku tej metody otrzymamy n-2 punkty wartości pochodnej tej że funkcji. Jest to związane z tym iż dla pierwszego jak i ostatniego punktu nie mamy poprzedniego jak i następnego punktu, którego moglibyśmy wykorzystać do naszych obliczeń.

function f(x: real):real;
begin
  f:=10*x-exp(-x);
end;

procedure pochodna;
var i,n: integer; // n- liczba punktów
  a,b,h,x,fi: real; // a,b- granice przedziału, h- przyrost (dx), fi- wartość pochodnej
begin
// tutaj należy wczytać zmienne a, b oraz n ewentualnie wywołać procedurę z tymi parametrami
  x:=a; // przyporządkowanie wartość x pierwszego punktu
  h:=(b-a)/n; //obliczamy wartość przyrostu
  for i:=1 to n+1 do
  begin
    fi:=(f(x+h,nr)-f(x-h,nr))/(2*h); // obliczanie pochodnej funkcji f w punkcie x
// tu należy odczytywać wartość pochodnej w danym punkcie lub zapamiętać ją do tablicy np.
    x:=a+i*h; // następny punkt obliczeń
  end;
end;

Tak obliczone wartości pochodnej w punktach można poddać dalszej obróbce, a mianowicie aproksymacji lub interpolacji uzyskując w ten sposób równanie pochodnej. Równanie to w zależności od wybranej metody może być przedstawione w postaci wielomianu lub funkcji trygonometrycznych ewentualnie ich kombinacji.

Shemat Hornera - Obliczanie wartości wielomianu

Shemat Hornera jest najczęściej wykorzystywanym algorytmem do obliczania wartości wielomianu, szczególnie w przypadku wielomianów wysokiego stopnia.
W metodzie tej wykorzystuje się przekształcenie wzoru tradycyjnego w inną postać umożliwiającą szybsze wyznaczenie szukanej wartości. A mianowicie:

Jest to nic innego jak proces wyciągania x jako części wspólnej przed nawias, a następnie wyznaczanie wartości tych nawiasów.

function P(x: real;n: integer):real;// x- argument, n- stopień wielomianu
var i: integer;
  b: real;
  a: array [0..10] of real;// zakres można dowolnie zmieniać
begin
// wczytywanie współczynników wielomianu
//można wczytywać poza funkcją - a jest wtedy zmienną globalną
//wczytywane współczynniki zapisujemy w kolejności odwrotnej(a[0]=a_n, a[1]=a_n-1)
  for i:=n downto 0 do
  readln(a[n-i]);
// Obliczanie wartości wielomianu o stopniu n dla żądanego x
  b:=a[0];
  for i:=1 to n do
  b:=b*x+a[i];
  P:=b;
// obliczona wartość wielomianu
end;

Ułamki łańcuchowe - obliczanie wartości pierwiastków kwadratowych

Ułamki łańcuchowe są to ułamki o piętrowym mianowniku rozwijane w nieskończoność o charakterze cyklicznym. Poniższy rysunek pokazuje przykład takiego ułamka rozwiniętego dla pierwiastka z liczby 3, oraz zapis skrócony, dla kilku innych liczb. W dziale Cyferki >> znajduje się liczba pierwiastek z 2 rozwinięta w postaci takiego ułamka.

Poniższa funkcja pozwala na obliczanie ułamków łańcuchowych skończonych - o skończonej liczbie ogniw.

function ulamek(a: real;n: integer;p: integer):real;// a- część całkowita, n- wielkość ułamka, p- liczba cykli
var i,j: integer;
  l,lp: array [0..6] of real;// ilość mianowników można dowolnie zmieniać
begin
// wczytywanie mianowników
//można wczytywać poza funkcją - l jest wtedy zmienną globalną
  for i:=1 to n do
  begin
    readln(l[i]);
    lp[i]:=l[i];
  end;
// Obliczanie wartości ułamka o stopniu n
  l[n]:=l[n]-1;
  l[n+1]:=1;
  for j:=1 to p do
  begin
    for i:=n downto 1 do
    l[i]:=l[i]+1/l[i+1];
    l[n+1]:=l[1];
    for i:=1 to n do
    l[i]:=lp[i];
  end;
  ulamek:=a+1/l[n+1];
// obliczona wartość ułamka
end;

Przykładowe wywołanie funkcji dla pierwiastka z 3:
wynik:=ulamek(1,2,5); oraz 1 2.
Dla pierwiastka z 7: wynik:=ulamek(2,4,5); oraz 1 1 1 4.
Od liczby powtórek zależy dokładność obliczenia pierwiastka, a powinna ona wahać się w granicach 4 - 8. Mniejsze i większe wartości nie mają praktycznego sensu. W przykładzie wynosi ona 5.

---