### Treść zadania ###

Na ile podprzedziałów należy podzielić przedział 0 do 5 aby funkcja f(x)=1/(1+x^2)
była interpolowana z dokładnością co najmniej 1e-5.

Normę błędu liczyć jako L2.

Dla uproszczenia dokonywać podziału na 1,2,4,8,16,...,2^n podprzedziałów.

Interpolacja liniowa

UWAGA: Program kompilowany w Dev-C++

#include <cstdlib>

#include <iostream>
#include <string>
#include <cmath>
#include <vector>

using namespace std;

//=================================================================================

struct wspolrzedne {
double x;
double y;
};

// Nasza funkcja
double f(double x)
{
return 1/(1+x*x);
}

// Wzor funkcji liniowej przechodzącej przez dwa punkty [x0, f(x0)] i [x1, f(x1)]
double p(double x, double x0, double x1)
{
return (
f(x0) + ( (f(x1)-f(x0))/(x1-x0) ) * (x-x0)
);
}

/* ### UWAGA: Funkcja pomocnicza. Nie jest wykorzystywana w programie.

Całkowanie numeryczne - metoda Simpsona

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.

xp - x początkowe
xk - x końcowe
n - liczba punktów podziałowych dla przedziału od xp do xk
*/
double calka_simpsona_z_f(double xp, double xk, int n)
{
double x;
double s = 0;
double st = 0;
double h = (xk - xp) / n; // odleglosc miedzy kolejnym wezlem

for(int i=1; i<=n; i++)
{
x = xp+i*h;
st += f(x-h/2);
if(i<n) {
s += f(x);
}
}
return h / 6 * (f(xp) + f(xk) + 2 * s + 4 * st);
}

/*
### UWAGA: Funkcja pomocnicza. Nie jest wykorzystywana w programie.

Całkowanie metodą Simpsona funkcji liniowej p przechodzącej przez dwa punkty x0 i x1
*/
double calka_simpsona_z_p(double xp, double xk, int n, double x0, double x1)
{
double x;
double s = 0;
double st = 0;
double h = (xk - xp) / n; // odleglosc miedzy kolejnym wezlem

for(int i=1; i<=n; i++)
{
x = xp+i*h;
st += p(x-h/2, x0, x1);
if(i<n) {
s += p(x, x0, x1);
}
}
return h / 6 * (p(xp, x0, x1) + p(xk, x0, x1) + 2 * s + 4 * st);
}

/*
Funkcja interpolująca. Jest złożona z wielu funkcji liniowych.
x - argument dla ktorego liczymy wartosc
tab - tablica wspolrzednych punktow wezlowych

*/
double finterp(double x, vector<wspolrzedne> &tab)
{
double y = 0; // przechowuje obliczoną wartosc funkcji interpolujacej

for(int i=0; i<tab.size(); i++)
{
// szukamy między jakimi punktami węzłowymi znajduje się podany jako argument x
// a następnie obliczamy wartosc funkcji liniowej przechodzącej przez te punkty węzłowe dla podanego x
if(x <= tab[i+1].x)
{
y = p(x, tab[i].x, tab[i+1].x);
break; // po znalezieniu odpowiedniego przedziału przerywamy pętle
}
}
return y;
}

/*
Funkcja g(x) = (f(x) - finterp(x))^2
Jest to funkcja pod całkowa do wzoru na norme euklidesową - L2
*/
double g(double x, vector<wspolrzedne> &tab)
{
double tmp = f(x)-finterp(x, tab);
return pow(tmp, 2);
}

/*
Całka Simpsona z funkcji g
Całka ta zostanie wykorzystana do liczenia normy L2

Zalecana liczba podprzedzialow n to iloczyn liczby podprzedziałów całej funkcji
i ilości podziałów w danym podprzedziale.
*/
double calka_simpsona_z_g(double xp, double xk, int n, vector<wspolrzedne> &tab)
{
double x;
double s = 0;
double st = 0;
double h = (xk - xp) / n; // odleglosc miedzy kolejnym wezlem

for(int i=1; i<=n; i++)
{
x = xp+i*h;
st += g(x-h/2, tab); // suma wartosci dla wczesniejszych węzłów
if(i<n) { // zabezpieczamy sie przed wyjsciem poza przedzial
s += g(x, tab); // suma wartosci dla następnych węzłów
}
}
return h / 6 * (g(xp, tab) + g(xk, tab) + 2 * s + 4 * st);
}

//=================================================================================

int main(int argc, char *argv[])
{
cout<<"### (c) Webook.pl (2010/2011) - Interpolacja liniowa ###nn";
//cout.setf(ios::showpoint); // uzupelnia liczby zerami po przecinku
cout.precision(30); // wyswietla liczby zmienno przecinkowe do 30 miejsc po przecinku

double eps = pow(10.,-5); // wymagana dokładność

double a = 0.; // początek przedziału
double b = 5.; // koniec przedziału

int licznik = 2;
double n = pow(2., licznik); // liczba podprzedzialow od ktorej startujemy

double h; // odleglosc miedzy kolejnym wezlem na osi x
double L2; // zmienna przechowujaca norme euklidesowa

wspolrzedne tmp_wsp; // zmienna pomocnicza do tworzenia struktury zawierajacej wspolrzedne punktu [x,y]

vector<wspolrzedne> tab; // tworzymy vector złożony ze struktury dwóch zmiennych x i y

// powtarzamy pętlę aż zostanie wywołany break
while(1)
{
tab.clear(); // kasujemy stan vectora przy rozpoczeciu kazdej kolejnej iteracji

h = (b-a)/n; // liczymy odleglosc miedzy kolejnym wezlem na osi x

cout<<"=== Krok "<<licznik-1<<" ================n"
<<"n = "<<n<<" (liczba podprzedzialow)n"
<<"h = "<<h<<" (odleglosc miedzy kolejnym wezlem)nn";

cout<<"Wspolrzedne punktow wezlowych:n";

// mamy n+1 punktow węzłowych dlatego zaczynamy od 0
for(double i=0; i<=n; i++)
{
tmp_wsp.x = h*i;
tmp_wsp.y = f(h*i);
tab.push_back(tmp_wsp); // wstawiamy na koniec vectora nasze wspolrzedne

// podgląd wspolrzednych
cout<<"Punkt "<<i<<" # x: "<<tab[i].x<<", y: "<<tab[i].y<<"n";
}

// każdy podprzedział jest dzielony na kolejne 1000 podprzedzialow w czasie calkowania
L2 = sqrt(calka_simpsona_z_g(a, b, n*1000, tab));

cout<<"nNorma L2 = "<<L2<<"nnn";

if(L2 <= eps)
{
break; // przerywamy pętle gdy norma jest mniejsza lub rowna zadanej dokładności
}
else
{
n = pow(2., ++licznik); // zwiekszamy liczbe podprzedzialow
}

}

// po wyjściu z pętli wypisujemy odpowiedź.
cout<<"Wymagana liczba podprzedzialow to: "<<n;

/*
W ostatnim kroku zostanie wypisane:

Norma L2 = 9.4489457156375221e-006

Wymagana liczba podprzedzialow to: 512
*/

// testowe sprawdzenie jak działają funkcje:
//cout<<"g: "<<g(1, tab);
//cout<<"finterp: "<<finterp(1.25, tab);

cout<<"nn";
system("PAUSE");
return EXIT_SUCCESS;
}