[C++] Implementacja interpolacji Newtona

[szw1710]

Cieszę się. Zobacz do Kompendium Analizy na moje posty o tych ilorazach i wzorze Newtona.

[karpiuch]

Muszę przyznać, że dawno któryś z tematów matematycznych mnie zaciekawił. Wczoraj późniejszym wieczorem właśnie przeglądałem posty w Kompendium na ten temat. Wieloman Newton'a to taki Pana konik?

[szw1710]

Z jego użyciem zrobiłem wiele swoich badań naukowych. Tak więc może go nie badam, ale używam, więc jako tako znam.

[Mariusz M]

Ja też się trochę tym bawiłem
Oto co napisałem

Kod: Zaznacz cały

#include<stdio.h>
#include<stdlib.h>

double horner(int,double* ,double ,double* );
void invhorner(int,double*,double,double*);
void newton(int,double*,double*,double*);

int main()
{
    FILE* in;
    FILE* out;
    double *x, *y, *a;
    int i,n;
    if((in=fopen("in.txt","r"))==NULL)
    {
        printf("Nie mozna otworzyc pliku do czytania 
");
        exit(1);
    }
    if((out=fopen("out.txt","w"))==NULL)
    {
        printf("Nie mozna otworzyc pliku do pisania 
");
        exit(1);
    }
    while(!feof(in))
    {
        fscanf(in,"%d",&n);
        x=(double*)malloc((n+1)*sizeof(double));
        y=(double*)malloc((n+1)*sizeof(double));
        a=(double*)malloc((n+1)*sizeof(double));
        for(i=0; i<n; i++)
            fscanf(in,"%lf",&x[i]);
        for(i=0; i<n; i++)
            fscanf(in,"%lf",&y[i]);
        newton(n,x,y,a);
        for(i=0; i<n; i++)
            fprintf(out,"a[%d]=%.15lf
",i,a[i]);
        fprintf(out,"
");
        free(x);
        free(y);
        free(a);
    }
    fclose(in);
    fclose(out);
    return 0;
}

double horner(int n,double* a,double c,double* b)
{
    int i;
    double r;
    r=a[n];
    for(i=n-1; i>=0; i--)
    {
        b[i]=r;
        r=r*c+a[i];
    }
    return r;
}

void invhorner(int n,double* b,double c,double* a)
{
    int k;
    double* u;
    u=(double*)malloc((n+1)*sizeof(double));
    for(k=0; k<=n; k++)
        u[k]=b[k];
    a[n+1]=b[n];
    for(k=0; k<=n; k++)
        if((k-1)<0) a[k]=-c*u[k];
        else a[k]=u[k-1]-c*u[k];
    free(u);
}

void newton(int n,double* x,double* y,double* a)
{
    int i,j;
    double* b;
    double** dd;
    b=(double*)malloc((n+1)*sizeof(double));
    dd=(double**)malloc((n+1)*sizeof(double));
    for(i=0; i<=n; i++)
        dd[i]=(double*)malloc((n+1)*sizeof(double));
    for(i=0; i<=n; i++)
    {
        b[i]=0.0;
        a[i]=0.0;
    }
    for(i=0; i<=n; i++)
        for(j=0; j<=n; j++)
            dd[i][j]=0.0;
    for(j=0; j<=n; j++)
        dd[j][0]=y[j];
    for(j=1; j<=n-1; j++)
        for(i=0; i<=n-j; i++)
            dd[i][j]=(double)(dd[i+1][j-1]-dd[i][j-1])/(x[i+j]-x[i]);

    b[0]=1.0;
    for(i=0; i<=n; i++)
    {
        for(j=0; j<=i; j++)
            a[j]+=dd[0][i]*b[j];
        invhorner(i,b,x[i],b);
    }
    free(b);
    for(i=0; i<=n; i++)
        free(dd[i]);
    free(dd);
}