Dobry wieczór. Jestem w trakcie pisania programu obliczającego wartości współczynników wielomianu interpolacji Lagrange'a. Napisałem kod, który dla dwóch (dość małych) zestawów danych działał prawidłowo. Niestety, w przypadku trzeciego, trochę większego zestawu, ostatnie 10 wartości wynosi -nan, a pozostałe to jedynki. Tablice z danymi wejściowymi są prawidłowe, nigdzie też nie wyłapałem dzielenia przez zero. Zmiana z double na long double skutkuje tym, że wszystkie wyniki wynoszą -nan. Skończyły mi się już pomysły, stąd też moja prośba o pomoc w znalezieniu błędu.
Oto kod programu:
#include <iostream>
#include <iomanip>
double lj(const double xn[], int n, int idx, int x) //Oblicza wartość fj
{ //(x-x0)(x-x1)...(x-x(i-1))(x-x(i+1))...(x-xn)/(xj-x0)(xj-x1)...(xj-x(i-1))(xj-x(i+1))...(xj-xn)
double wynik=1; //wartość funkcji
//złożoność O(n)
for(int t=0; t<n; ++t)
{
if(t!=idx)
{
wynik*=(x-xn[t])/(xn[idx]-xn[t]);
}
}
return wynik;
}
template<int N>
struct x_points
{
constexpr x_points() : array() //konstruktor constexpr
{
for(int i=0; i<N; ++i) //węzły obliczane są w czasie kompilacji
{
array[i]=-1+i*(1/32.);
}
}
double array[N];
};
double* valuesOfFunction(const double x[], int size) //wartości w węzłach
{
double *tab=new double[size];
for(int i=0; i<size; ++i)
{
tab[i]=1/(1+5*x[i]*x[i]);
}
return tab;
}
int main()
{
const int precision=4;
const int zero_val=0;
const int x_max=1;
const int x_min=-1;
constexpr int size=((x_max-x_min)/(1/32.))+1; //rozmiar tablicy, liczba węzłów, wartości oraz współczynników
constexpr auto xv=x_points<size>();
const double *xn=xv.array; //węzły
/*for(int i=0; i<size; ++i)
{
std::cout<<xn[i]<<std::endl;
}*/
const double *fx=valuesOfFunction(xn, size); //wartości w węzłach
/*for(int i=0; i<size; ++i)
{
std::cout<<fx[i]<<std::endl;
}*/
//constexpr int size=std::size(xn); //rozmiar tablicy, liczba węzłów, wartości oraz współczynników
double temp[size]; //tablica dla wartości w węzłach kolejnych wielomianów
double a[size]; //współczynniki wielomianu w postaci naturalnej
int xn_zero_value_index=-1; //indeks miejsca, w którym tablica xn zawiera element zerowy
for(int i=0; i<size; ++i)
{
if(xn[i]==0)
{
xn_zero_value_index=i;
}
a[0]+=fx[i]*lj(xn, size, i, zero_val); //wyraz wolny
}
if(xn_zero_value_index!=-1)
{
temp[xn_zero_value_index]=a[0];
}
for(int j=0; j<size; ++j) //oddzielnie ze względu na brak konieczności interpolacji wartości
{ // (fj-a0)/xj
if(xn[j]!=0)
{
temp[j]=(fx[j]-a[0])/xn[j]; //wartości y1(x) w kolejnych węzłach
}
a[1]+=temp[j]*lj(xn, size, j, zero_val); //suma kolejnych wyrazów wielomianu interpolacyjnego w zerze
}
//temp[xn_zero_value_index]=a[1];
//pozostałe współczynniki wg wzoru:
// yi(xj)=(y{i-1}(xj)-a{i-1})/xj
for(int k=2; k<size; ++k) //pętla po kolejnych współczynnikach
{
if(xn_zero_value_index!=-1)
{
temp[xn_zero_value_index]=a[k-1];
}
for(int l=0; l<size; ++l) //pętla po kolejnych wartościach w węzłach
{
if(xn[l]!=0)
{
temp[l]=(temp[l]-a[k-1])/xn[l];
}
a[k]+=temp[l]*lj(xn, size, l, zero_val);
}
}
for(int m=0; m<size; ++m) //wyniki
{
std::cout<<std::setprecision(precision)<<std::fixed<<"a"<<m<<" = "<<a[m]<<std::endl;
}
return 0;
}