[C++] Interpolacja Hermite'a

[zdolowany123]

Kod: Zaznacz cały

//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ą conajmniej 1e-5. Normę błędu liczyć jako L_2. Dla uproszczenia, dokonywać podziału na
//  1,2,4,8,16,...,2^n podprzedziałów

//projekt 4: interpolacja liniowa Hermite'a C^2

#include <iostream>
#include <cmath>
using namespace std;

double transform(double x, double a, double b)
{return (x-a)/(b-a);}

double f(double x)
{return 1/(1+pow(x,2));}

double fd(double x)
{return -2*x/pow(1+pow(x,2),2);}

double fdd(double x)
{return (8*pow(x,2)/pow(1+pow(x,2),3))-(2/pow(1+pow(x,2),2));}

double H1(double x)
{return 1-10*pow(x,3)+15*pow(x,4)-6*pow(x,5);}

double H2(double x)
{return x-6*pow(x,3)+8*pow(x,4)-3*pow(x,5);}

double H3(double x)
{return pow(x,2)/2-3*pow(x,3)/2+3*pow(x,4)/2-pow(x,5)/2;}

double H4(double x)
{return 10*pow(x,3)-15*pow(x,4)+6*pow(x,5);}

double H5(double x)
{return -4*pow(x,3)+7*pow(x,4)-3*pow(x,5);}

double H6(double x)
{return pow(x,3)/2-pow(x,4)+pow(x,5)/2;}

double Hermit(double x, double a, double b)
{
       double H,eta;
       eta=transform(x,a,b);
       
       H=f(a)*H1(eta)+fd(a)*H2(eta)+fdd(a)*H3(eta)+f(b)*H4(eta)+fd(b)*H5(eta)+fdd(b)*H6(eta);
       return H;
}

double norma();

int main()
{    
     double norm; //norma(błąd)
     double x,h,h_norm,a,b; //wezel, krok, krok do normy,  końce podprzedziału
     int iprz=1; //ilosc podprzedzialow
     int iprz_norm=2; //ilosc podprzedzialow
     int i=0; //iterator
     
     do
     {
         norm=0.0;
         a=0.0;
         b=5.0;
         x=a; //pierwszy wezel
         h=(b-a)/iprz; //krok
         b=a+h; //koniec pierwszego podprzedzialu
         h_norm=h/iprz_norm; //krok dla liczenia bledu(normy)
         
         for (int j=0;j<iprz;j++)
         {
            if (j==iprz-1) iprz_norm++;
            for (int k=0;k<iprz_norm;k++)
            {
                norm=norm + h_norm * (Hermit(x,a,b)-f(x)) * (Hermit(x,a,b)-f(x)); //policz blad(norme)
//                printf("iprz=%f\nx=%f\n",a,x);getchar();
                if (j==iprz-1)
                {
                    h_norm=(b-a)/(iprz_norm-1);
                    x=x+h_norm;
                }
                else
                {
                    h_norm=(b-a)/iprz_norm;
                    x=x+h_norm;
                }
            }
            a=b; //poczatek podprzedzialu
            b=b+h; //koniec podprzedzialu, h=krok
            x=a; //kolejny wezel
         }
         cout<<endl<<"ilosc podprzedzialow: "<<iprz<<endl;
         cout<<"norma: "<<norm<<endl;
         iprz*=2; // l. podprzedzialow to kolejne potegi dwojki
         i++;
     }
     while(norm>=0.0001);
     
     cout<<endl<<"ilosc podprzedzialow: "<<iprz/2<<endl;
     cout<<"norma: "<<norm<<endl;
     cout<<"ilosc iteracji: "<<i<<endl;
     
     getchar();
     return 0;
}

Może mi ktoś powiedzięc skąd się bierze h1, h2, h3 itd?