2018年12月23日 星期日

Section 3.4 Hermite Polynomial

Section 3.4 Hermite Polynomial; 

The Hermite C code needs the input data. 

2       0.5            -0.25
2.75    0.36363636     -0.13223140496
4       0.25           -0.0625

C code of Hermite Interpolation 

/*
*      HERMITE INTERPOLATION ALGORITHM 3.3
*
*      TO OBTAIN THE COEFFICIENTS OF THE HERMITE INTERPOLATING
*      POLYNOMIAL H ON THE (N+1) DISTINCT NUMBERS X(0), ..., X(N)
*      FOR THE FUNCTION F:
*
*      INPUT:   NUMBERS X(0), X(1), ..., X(N); VALUES F(X(0)), F(X(1)),
*               ..., F(X(N)) AND F'(X(0)), F'(X(1)), ..., F'(X(N)).
*
*      OUTPUT:  NUMBERS Q(0,0), Q(1,1), ..., Q(2N + 1,2N + 1) WHERE
*
*               H(X) = Q(0,0) + Q(1,1) * ( X - X(0) ) + Q(2,2) *
*                      ( X - X(0) )**2 + Q(3,3) * ( X - X(0) )**2 *
*                      ( X - X(1) ) + Q(4,4) * ( X - X(0) )**2 *
*                      ( X - X(1) )**2 + ... + Q(2N + 1,2N + 1) *
*                      ( X - X(0) )**2 * ( X - X(1) )**2 * ... *
*                      ( X - X(N - 1) )**2 * (X - X(N) ).
*/

#include<stdio.h>
#include<math.h>
#define true 1
#define false 0

double F(double);
void INPUT(int *, double *, double [][26], int *);
void OUTPUT(FILE **);
void OUTPUT2(FILE **, int, int, double *, double [][26]);

main()
{
   double Q[26][26],X[13],Z[26];
   double XX,S;
   int I,J,K,N,OK;
   char A; 
   FILE *OUP[1];

   INPUT(&OK, X, Q, &N);
   if (OK) {
      /* STEP 1 */
      for (I=0; I<=N; I++) {
         /* STEP 2 */
         Z[2*I] = X[I];
         Z[2*I+1] = X[I];
         Q[2*I+1][0] = Q[2*I][0];
         /* STEP 3 */
         if (I != 0) 
            Q[2*I][1] = (Q[2*I][0] - Q[2*I-1][0]) / (Z[2*I] -Z[2*I-1]);
      }    
      /* STEP 4 */
      K = 2 * N + 1;
      for (I=2; I<=K; I++)
         for (J=2; J<=I; J++)
            Q[I][J] = ( Q[I][J - 1] - Q[I - 1][J - 1] ) / ( Z[I] - Z[I - J] );
      /* STEP 5 */
      OUTPUT(OUP);
      OUTPUT2(OUP, K, N, X, Q);
      printf("Do you wish to evaluate this polynomial?\n");
      printf("Enter Y or N\n");
      scanf("\n%c", &A);
      if ((A == 'Y') || (A == 'y')) {
         printf("Enter a point at which to evaluate\n");
         scanf("%lf", &XX);
         S = Q[K][K] * (XX - Z[K-1]);
         for (I=2; I<=K; I++) {
            J = K - I + 1;
            S = (S + Q[J][J]) * (XX - Z[J-1]);
         }
         S = S + Q[0][0];
         fprintf(*OUP, "x-value and interpolated-value\n");
         fprintf(*OUP, "  %12.10e %12.10e\n", XX, S);
      }      
   }
   fclose(*OUP);
   return 0;
}
/* Change functions F and FP if program is to calculate the
   first column and part of the second column of Q  */
double F(double X)
{
   double f; 

   f =  1.0/X;
   return f;
}

/*  Note:  Function FP is the derivative of F  */
double FP(double X)
{
   double f; 

   f =  -1.0/(X*X);
   return f;
}

void OUTPUT(FILE **OUP)
{
   int FLAG;
   char NAME[30];

   printf("Choice of output method:\n");
   printf("1. Output to screen\n");
   printf("2. Output to text file\n");
   printf("Please enter 1 or 2\n");
   scanf("%d", &FLAG);
   if (FLAG == 2) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("for example:   A:OUTPUT.DTA\n");
      scanf("%s", NAME);
      *OUP = fopen(NAME, "w");
   }
   else *OUP = stdout;
   fprintf(*OUP, "HERMITE INTERPOLATING POLYNOMIAL\n\n");
}

void INPUT(int *OK, double *X, double Q[][26], int *N)
{
   int I, FLAG;
   char A;
   char NAME[30];
   FILE *INP; 

         printf("This is Hermite interpolation.\n");
         *OK = true;

         printf("Has a text file been created with the data in three ");
         printf("columns?\n");
         printf("Enter Y or N\n");
         scanf("\n%c", &A);
         if ((A == 'Y') || (A == 'y')) {
            printf("Input the file name in the form - ");
            printf("drive:name.ext\n");
            printf("for example:   A:DATA.DTA\n");
            scanf("%s", NAME);
            INP = fopen(NAME, "r");
            *OK = false;
            while (!(*OK)) {
               printf("Input the number of data points minus 1\n");
               scanf("%d", N);
               if (*N > 0) {
                  for (I=0; I<=*N; I++) 
                     fscanf(INP, "%lf %lf %lf", &X[I], &Q[2*I][0],
                            &Q[2*I+1][1]);
                  fclose(INP);
                  *OK = true;
               }
               else printf("Number must be a positive integer\n");
            }
         }
         else {
            printf("Please create the input file in three column ");
            printf("form with the X values, F(X), and\n");
            printf("F'(X) values in the corresponding columns.\n");
            printf("The program will end so the input file can ");
            printf("be created.\n");
            *OK = false;
         }
}

void OUTPUT2(FILE **OUP, int K, int N, double *X, double Q[][26])
{
   int I; 

   fprintf(*OUP, "The input data follows:\n");
   fprintf(*OUP, "  X, F(X), F''(X)\n");
   for (I=0; I<=N; I++) 
      fprintf(*OUP, "  %12.10e %12.10e %12.10e\n", X[I], Q[2*I][0], Q[2*I+1][1]);
   fprintf(*OUP, "\nThe Coefficients of the Hermite Interpolation Polynomial\n");
   fprintf(*OUP, "in order of increasing exponent follow:\n\n");
   for (I=0; I<=K; I++) fprintf(*OUP, " %12.10e\n", Q[I][I]);
} 

沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...