Algorytmy numeryczne

1. Algorytmy numeryczne • Podstawowe algorytmy numeryczne: • poszukiwanie miejsc zerowych funkcji; • iteracyjne obliczanie wartości funkcji; • interpolacja funkcji metodą Lagrange’a; • różniczkowanie funkcji; • całkowanie funkcji metodą Simpsona; • rozwiązywanie równań liniowych metodą Gaussa. Kompletne algorytmy, bez uzasadnienia matematycznego Literatura: Knuth D.E.: Sztuka programowania –tom 1,2,3 Wydawnictwa Naukowo-Techniczne, 2001; Praca zbiorowa pod redakcją Jerzego Klamki: Laboratorium metod numerycznych, Skrypt nr 1305 Politechniki Śląskiej w Gliwicach, Gliwice 1987; Wróblewski P.: Algorytmy, struktury danych i techniki programowania; Stożek E.: Metody numeryczne w zadaniach; Inne książki z metod numerycznych…..

2. Algorytmy numeryczne – miejsca zerowe Idea algorytmu: Systematyczne przybliżanie się do miejsca zerowego funkcji za pomocą stycznych do krzywej. Iteracyjnie powtarzamy stop, jeśli Symbol  oznacza stałą, gwarantującą zatrzymanie algorytmu, z0 jest wartością początkową, f i f’ trzeba znać jawnie.

3. Algorytmy numeryczne – miejsca zerowe #include <iostream.h> #include <math.h> const double epsilon=0.0001; double f(double x) //funkcja f(x)=3x2-2 {return 3*x*x-2;} doublefp(double x) //pochodna f’(x)=(3x2-2)’=6x {return 6*x;} double zero(double x0,double(*f)(double),double(*fp)(double)) { if(f(x0)<epsilon)return x0; else return zero(x0-f(x0)/fp(x0),f,fp); } intmain() { cout << "Zero funkcji 3x*x-2 wynosi "<<zero(1,f,fp)<<endl; //wynik 0,816497 } Zostały użyte wskaźniki do funkcji, ale funkcji można użyć też bezpośrednio.

4. Algorytmy numeryczne – wartość funkcji Idea algorytmu: Załóżmy, że dysponujemy jawnym wzorem funkcji w postaci uwikłanej: F(x,y)=0. Za pomocą metody Newtona można obliczyć jej wartość dla pewnego x w sposób iteracyjny: Wartość początkowa y0 powinna być jak najbliższa wartości poszukiwanej y i spełniać warunek: F(x,y)F’y(x,y)>0. Zalety: Czasami iloraz znacznie się upraszcza! Przykład:

5. Algorytmy numeryczne – wartość funkcji #include <iostream.h> #include <math.h> const double epsilon=0.00000001; double wart(doublex,doubleyn) { double yn1=2*yn-x*yn*yn; //fabs(x)=|x|, wartość bezwzględna dla danych double if( fabs(yn-yn1)<epsilon) return yn1; else return wart(x,yn1); } intmain() { cout << "Wartość funkcji y=1/x dla x=7 wynosi "<<wart(7,0.1); //funkcja f(x)=3x2-2 }

6. Algorytmy numeryczne – wartość funkcji Efektywne obliczanie wartości wielomianów – schemat Hornera. Ma on zastosowanie np. w algorytmach kryptograficznych, bo trzeba tam wykonywać obliczenia na bardzo dużych liczbach całkowitych. Można takie liczby traktować jako współczynniki wielomianów. Podstawa x=10 12 9876 0002 6000 0000 0054 Podstawa x=10000 Wtedy operacje na dużych liczbach zastępujemy algorytmami działającymi na wielomianach.

7. Algorytmy numeryczne – wartość funkcji W(b) – „metoda siłowa” #include <iostream.h> constn=5; // stopień wielomianu int oblicz(int b, int w[], introzm) { int res=0, pot=1; for(int j=rozm-1;j>=0;j--) { res+=pot*w[j]; //sumy cząstkowe pot*=b; //następna potęga b } returnres; } intmain() { intw[n]={1,4,-2,0,7}; // współczynniki wielomianu cout << oblicz(2,w,n) << endl; }

8. Algorytmy numeryczne – wartość funkcji W(b) – schemat Hornera #include <iostream.h> const n=5; // stopień wielomianu intoblicz_wielomian2(inta,int w[],introzm) // schemat Hornera { intres=w[0]; for(int j=1;j<rozm;res=res*a+w[j++]); return res; } intmain() { intw[n]={1,4,-2,0,7}; // współczynniki wielomianu cout << oblicz_wielomian2(2,w,n) << endl; }

9. Algorytmy numeryczne – interpolacja Idea algorytmu: • Interpolacja funkcji metodą Lagrange’a służy do przybliżania funkcji (np. za pomocą wielomianu określonego stopnia), gdy: • nie znamy funkcji tylko dysponujemy fragmentem wykresu, tzn. znamy jej wartości dla skończonego zbioru argumentów • albo też wyliczanie na podstawie wzoru jest zbyt czasochłonne Interpolacja funkcji f(x) za pomocą F(x).

10. Algorytmy numeryczne – interpolacja Idea algorytmu: Na rysunku na podstawie 7 par punktów udało się obliczyć wielomian F(x). Wielomian interpolacyjny konstruuje się za pomocą tzw. wyznacznika Vandermonde’a, który pozwala na wyliczenie współczynników. Jeśli jednak szukamy tylko wartości funkcji w określonym punkcie z, to prostsza jest metoda Lagrange’a: Powyższy wzór łatwo się tłumaczy na kod C++ za pomocą dwóch zagnieżdżonych pętli „for”.

11. Algorytmy numeryczne – interpolacja #include <iostream.h> #include <math.h> constint n=3; // stopień wielomianu interpolującego // wartośći funkcji (y[i]=f(x[i])) double x[n+1]={3.0, 5.0, 6.0, 7.0}; double y[n+1]={1.732, 2.236, 2.449, 2.646}; doubleinterpol(double z, double x[n], double y[n]) // zwraca wartość funkcji w punkcie 'z' { double wnz=0,om=1,w; for(int i=0;i<=n;i++) { om=om*(z-x[i]); w=1.0; for(int j=0;j<=n;j++) if(i!=j) w=w*(x[i]-x[j]); wnz=wnz+y[i]/(w*(z-x[i])); } returnwnz=wnz*om; } intmain() { double z=4.5; cout << "Wartość funkcji sqrt(x) w punkcie " << z; cout << " wynosi " << interpol(z,x,y) <<endl; }

12. Algorytmy numeryczne – różniczkowanie Idea algorytmu: Do numerycznego różniczkowania służy tzw. wzór Stirlinga (nie wyprowadzamy!). Pozwala on wyznaczyć pochodne f’ i f’’ w punkcie x0 dla pewnej funkcji f(x), której wartości znamy w postaci tabelarycznej: Tablica różnic centralnych

13. Algorytmy numeryczne – różniczkowanie Idea algorytmu: Różnice  są obliczane w identyczny sposób w całej tabeli, np.: Przyjmując upraszczające założenie, że zawsze będziemy obliczali pochodne dla punktu centralnego x=x0 wzór Stirlinga ma postać: Punktów kontrolnych może być więcej niż 5. W algorytmie poniżej jest 5 punktów, co prowadzi do tablicy różnic centralnych niskiego rzędu. • Uwaga! • Im mniejsze h, tym większy błąd! • Za duże h jest niezgodne z ideą metody. • Metoda nie jest dobra do punktów na krańcach przedziałów zmienności argumentu funkcji.

14. Algorytmy numeryczne – różniczkowanie #include <iostream.h> #include <math.h> constintn=5; // rząd obliczanych różnic centralnych wynosi n-1 double t[n][n+1]= { {0.8, 4.80}, // pary (x[i], y[i]) dla y=5x*x+2*x {0.9, 5.85}, //inicjacja tablicy: wpisane są dwie pierwsze kolumny, {1, 7.00}, // a nie wiersze! {1.1, 8.25}, {1.2, 9.60} }; structPOCHODNE{double f1,f2;}; POCHODNE stirling(double t[n][n+1]) // funkcja zwraca wartości f'(z) i f''(z) gdzie z jest elementem // centralnym: tutaj t[2][0], tablica 't' musi być uprzednio centralnie // zainicjowana, jej poprawność nie jest sprawdzana { POCHODNE res; double h=(t[4][0]-t[0][0])/(double)(n-1); // krok argumentu 'x' for(int j=2;j<=n;j++) for(int i=0;i<=n-j;i++) t[i][j]=t[i+1][j-1]-t[i][j-1]; res.f1=((t[1][2]+t[2][2])/2.0-(t[0][4]+t[1][4])/12.0)/h; res.f2=(t[1][3]-t[0][5]/12.0)/(h*h); returnres; } intmain() { POCHODNE res=stirling(t); cout << "f'=" << res.f1 << ", f''=" << res.f2 <<endl; }

15. Algorytmy numeryczne – całkowanie Idea algorytmu: Przedstawiamy skomplikowaną funkcję w postaci przybliżonej , prostszej obliczeniowo. Na danym etapie i trzy kolejne punkty funkcji podcałkowej są przybliżane parabolą. Stosujemy to do każdego obszaru złożonego z 3 kolejnych punktów. • Odstępy h muszą być jednakowe. • Całka globalna jest sumą całek cząstkowych.

16. Algorytmy numeryczne – całkowanie #include <iostream.h> #include <math.h> constint n=4; // ilość punktów= 2n+1 double f[2*n+1]={41, 29, 19, 11, 5, 1, -1, -1, 1}; doublesimpson(double f[2*n+1], double a, double b) // funkcja zwraca całkę funkcji f(x) w przedziale [a,b], // której wartości są podane tabelarycznie w 2n+1 punktach { double s=0,h=(b-a)/(2.0*n); for(int i=0;i<2*n;i+=2) // skok cp dwa punkty! s+=h*(f[i]+4*f[i+1]+f[i+2])/3.0; return s; } intmain() { cout << "Wartość całki =" << simpson(f,-5,3) << endl; // 82.667 }

17. Algorytmy numeryczne – całkowanie #include <iostream.h> #include <math.h> constint n=4; // ilość punktów= 2n+1 doublefun(double x) { return x*x-3*x+1; } // funkcja x*x-3*x+1 w przedziale [-1,3] doublesimpson_f(double(*f)(double), // wskaźnik do f(x) double a, double b, int N) // funkcja zwraca całkę znanej w postaci wzoru funkcji f(x) // w przedziale [a,b], N - ilość podziałów { double s=0,h=(b-a)/(double)N; for(int i=1;i<=N;i++) s+=h*(fun(a+(i-1)*h)+4*fun(a-h/2.0+i*h)+fun(a+i*h))/3.0; return s; } intmain() { cout << "Wartość całki =" << simpson_f(fun,-5,3,8) << endl; // 82.667 } Całkowanie funkcji znanej w postaci analitycznej też jest możliwe

18. Algorytmy numeryczne – układy równań liniowych Idea algorytmu: Aby zastosować algorytm Gaussa, musimy najpierw zapisać układ w postaci uporządkowanej. Przykład: Dalej, stosujemy metodę eliminacji Gaussa, sprowadzając macierz do macierzy trójkątnej. Eliminacja x z wierszy 2 i 3 Eliminacja y z wierszy 1 i 3

19. Algorytmy numeryczne – układy równań liniowych Idea algorytmu: Stosujemy redukcję wsteczną, by wyznaczyć zmienne: • Niebezpieczeństwa: • Eliminacja zmiennych może prowadzić do dzielenia przez zero. • Zamieniamy wtedy wiersze, chyba, że poniżej wiersza i-tego nie ma tej zmiennej różnej od zera. • Wtedy układ nie ma rozwiązania.

20. Algorytmy numeryczne – układy równań liniowych #include <iostream.h> #include <math.h> constint N=3; double x[N]; double a[N][N+1]={ {5 , 0, 1, 9}, {1 , 1,-1, 6}, {2, -1, 1, 0}}; intgauss(double a[N][N+1], double x[N]) { int max; doubletmp; for(int i=0; i<N; i++) // eliminacja { max=i; for(int j=i+1; j<N; j++) if(fabs(a[j][i])>fabs(a[max][i])) max=j; for(intk=i; k<N+1; k++) // zamiana wierszy wartościami { tmp=a[i][k]; a[i][k]=a[max][k]; a[max][k]=tmp; } if(a[i][i]==0) return 0; // Układ sprzeczny! for(j=i+1; j<N; j++) for(k=N; k>=i; k--) // mnożenie wiersza j przez współczynnik "zerujący": a[j][k]=a[j][k]-a[i][k]*a[j][i]/a[i][i]; } cd.

21. Algorytmy numeryczne – układy równań liniowych // redukcja wsteczna for(int j=N-1; j>=0; j--) { tmp=0; for(int k=j+1; k<=N; k++) tmp=tmp+a[j][k]*x[k]; x[j]=(a[j][N]-tmp)/a[j][j]; } return 1; // wszystko w porządku! } intmain() { if(!gauss(a,x)) cout << "Układ (1) jest sprzeczny!\n"; else { cout << "Rozwiązanie:\n"; for(int i=0;i<N;i++) cout << "x["<<i<<"]="<<x[i] << endl; } }

22. Algorytmy numeryczne Warto poczytać o szybkiej transformacie Fouriera.