Transcript Podstawowe algorytmy numeryczne: • poszukiwanie miejsc zerowych funkcji; • iteracyjne obliczanie wartości funkcji; • interpolacja funkcji metodą Lagrange’a; Kompletne algorytmy, bez uzasadnienia matematycznego • różniczkowanie funkcji; • całkowanie funkcji.

Podstawowe algorytmy numeryczne:
• poszukiwanie miejsc zerowych funkcji;
• iteracyjne obliczanie wartości funkcji;
• interpolacja funkcji metodą Lagrange’a;
Kompletne
algorytmy, bez
uzasadnienia
matematycznego
• różniczkowanie funkcji;
• całkowanie funkcji metodą Simpsona;
• rozwiązywanie równań liniowych metodą Gaussa.
Literatura:
1. Knuth D.E.: Sztuka programowania –tom 1,2,3 Wydawnictwa NaukowoTechniczne, 2001;
2. Praca zbiorowa pod redakcją Jerzego Klamki: Laboratorium metod
numerycznych, Skrypt nr 1305 Politechniki Śląskiej w Gliwicach, Gliwice
1987;
3. Wróblewski P.: Algorytmy, struktury danych i techniki programowania;
4. Stożek E.: Metody numeryczne w zadaniach;
5. Inne książki z metod numerycznych…..
Idea algorytmu:
Systematyczne przybliżanie
się do miejsca zerowego
funkcji za pomocą
stycznych do krzywej.
zi  zi 1 
stop, jeśli
f ( zi 1 )
f '( zi 1 )
;
f ( zi )  
Iteracyjnie
powtarzamy
Symbol  oznacza
stałą, gwarantującą
zatrzymanie
algorytmu,
z0 jest wartością
początkową,
f i f’ trzeba znać
jawnie.
#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;}
double fp(double x) //pochodna f’(x)=(3x2-2)’=6x
{return 6*x;}
Zostały użyte
wskaźniki do funkcji,
ale funkcji można
użyć też bezpośrednio.
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);
}
int main()
{
cout << "Zero funkcji 3x*x-2 wynosi "<<zero(1,f,fp)<<endl;
//wynik 0,816497
}
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:
yn 1  yn  FF' y((xx, y, y)) ;
stop, jesli yn 1  yn   .
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:
y  1x ,
F ( x, y )  0  x  1y ,
Fy' ( x, y ) 
1
y2
,
yn 1  2 yn  x ( yn ) 2 .
#include <iostream.h>
#include <math.h>
const double epsilon=0.00000001;
double wart(double x,double yn)
{
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);
}
int main()
{
cout << "Wartość funkcji y=1/x dla x=7 wynosi "<<wart(7,0.1);
//funkcja f(x)=3x2-2
}
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
x21  2x20  (9x19  8x18  7 x17  6x16 )  (2x12 )  (6x11 )  (5x1  4).
12 9876 0002 6000 0000 0054
Podstawa x=10000
(12x5 )  (9876x 4 )  (2x3 )  (6000x 2 )  54.
Wtedy operacje na dużych liczbach zastępujemy algorytmami działającymi na
wielomianach.
W ( x)  an xn  an1xn1  ... a1x1  a0 .
W(b) – „metoda siłowa”
#include <iostream.h>
const n=5;
// stopień wielomianu
int oblicz(int b, int w[], int rozm)
{
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
}
return res;
}
int main()
{
int w[n]={1,4,-2,0,7};
// współczynniki wielomianu
cout << oblicz(2,w,n) << endl;
}
W (b)  ...anb  an1 b  an2 b  ... a1 b  a0 .
W(b) – schemat
Hornera
#include <iostream.h>
const n=5;
// stopień wielomianu
int oblicz_wielomian2(int a,int w[],int rozm)
// schemat Hornera
{
int res=w[0];
for(int j=1;j<rozm;res=res*a+w[j++]);
return res;
}
int main()
{
int w[n]={1,4,-2,0,7};
// współczynniki wielomianu
cout << oblicz_wielomian2(2,w,n) << endl;
}
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).
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:
n
F ( z )  ( z  x0 )(z  x1 )...(z  xn )
j 0 ( z  x j )
yj
n
 ( x j  xi )
i 0 ,i  j
Powyższy wzór łatwo się tłumaczy na kod C++ za pomocą dwóch
zagnieżdżonych pętli „for”.
.
#include <iostream.h>
#include <math.h>
const int 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};
double interpol(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]));
}
return wnz=wnz*om;
}
int main()
{
double z=4.5;
cout << "Wartość funkcji sqrt(x) w punkcie " << z;
cout << " wynosi " << interpol(z,x,y) <<endl;
}
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:
...(x0  2h, f ( x0  2h)), ( x0  h, f ( x0  h)), ( x0 , f ( x0 )), ( x0  h, f ( x0  h)),...
Tablica różnic
centralnych
Idea algorytmu:
Różnice  są obliczane w identyczny sposób w całej tabeli, np.:
f ( x0  32 h)  f ( x0  h)  f ( x0  2h).
Przyjmując upraszczające założenie, że zawsze będziemy obliczali pochodne
dla punktu centralnego x=x0 wzór Stirlinga ma postać:
f ' ( x) 
f ' ' ( x) 

1
1
1 f ( x  2 h ) f ( x  2 h )
h
2
1
h2

2

3
3
1
1
1  f ( x  2 h )  f ( x  2 h )
6
2

5
5
1
1
1  f ( x  2 h )  f ( x  2 h )
30
2


 ...
f ( x)  121  4 f ( x)  901  6 f ( x)  ...
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.
#include <iostream.h>
#include <math.h>
const int n=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}
};
struct POCHODNE{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);
return res;
}
int main()
{ POCHODNE res=stirling(t);
cout << "f'=" << res.f1 << ", f''=" << res.f2 <<endl;
}
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.
x2

x0
f ( x)dx  h
f ( x0 )  4 f ( x1 )  f ( x2 )
3
.
#include <iostream.h>
#include <math.h>
const int n=4; // ilość punktów= 2n+1
double f[2*n+1]={41, 29, 19, 11, 5, 1, -1, -1, 1};
double simpson(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;
}
int main()
{
cout << "Wartość całki =" << simpson(f,-5,3) << endl; // 82.667
}
#include <iostream.h>
#include <math.h>
const int n=4; // ilość punktów= 2n+1
double fun(double x)
{
return x*x-3*x+1;
}
// funkcja x*x-3*x+1 w przedziale [-1,3]
double simpson_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;
}
Całkowanie funkcji
znanej w postaci
analitycznej też
jest możliwe
int main()
{
cout << "Wartość całki =" << simpson_f(fun,-5,3,8) << endl; // 82.667
}
Idea algorytmu:
Aby zastosować algorytm Gaussa, musimy najpierw zapisać układ w postaci
uporządkowanej. Przykład:
5x  z  9
xz y 6
2x  y  z  0
5x  0 y  z  9
1x  1y  1z  6
2 x  1y  1z  0
 5 0 1  x   9 

   
 1 1  1 y    6 
 2  1 1  z   0 

   
Dalej, stosujemy metodę eliminacji Gaussa, sprowadzając macierz do
macierzy trójkątnej.
5x  0 y  z  9
1x  1y  1z  6
2 x  1y  1z  0
Eliminacja x z
wierszy 2 i 3
5x  0 y  z  9
0 x  1y  1,2 z  4,2
0 x  1y  0,6 z  3,6
Eliminacja y z
wierszy 1 i 3
5x  0 y  z  9
0 x  1y  1,2 z  4,2
 0 x  0 y  0,6 z  0,6
Idea algorytmu:
Stosujemy redukcję wsteczną, by wyznaczyć zmienne:
z  0,6 / 0,6  1
y  1,2 z  4,2  3
x  (9  z ) / 5  2
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.
#include <iostream.h>
#include <math.h>
const int N=3;
double x[N];
double a[N][N+1]={ {5 , 0, 1, 9}, {1 , 1,-1, 6}, {2, -1, 1, 0}};
int gauss(double a[N][N+1], double x[N])
{ int max;
double tmp;
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(int k=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.
// 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!
}
int main()
{
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;
}
}
Warto poczytać o szybkiej transformacie Fouriera.