[C++] Algorytm Newtona

[trzynastykrol]

Panowie, dzień dobry.

Zastanawiam się jak napisać algorytm Newtona w C++. Chodzi mi o wyznaczanie dla danych węzłów pojedynczych i wartości w węzłach wielomianu interpolacyjnego. Miałbym przy tym korzystać ze wzoru Newtona oraz wartości wielomianu w punkcie \(\displaystyle{ x=c}\).

Jakieś pomysły?

[Mariusz M]

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);
     }

Chodzi o coś takiego ?