I need help with Newton's DD (nested form)

[ayoesquire]

Sorry guys, please this is a bit off-the-topic but I wouldn't have asked if I didn't need your urgent help.

I am taking a MSc course in Applied Numerical Analysis and the programming language/software for the class is Matlab which is a fairly unfamiliar territory to me and I have limited time to master all of its syntax and semantics.

We were asked to derive a 6th order polynomial p(x) (where n =6) that is approximately equal to the function f(x) = log10(x) and subsequently solve for f(x) when the value of x = 1.43 using the Newton's Divided difference as follows:

p(x) = a0 + (x-x0)[a1 + (x-x1)[a2 + (x-x2)[a3 + (x-x3)[a4 +
% (x-x4)[a5]]]]]

Here's the best I could come up with for now:
********************************************************************************
1. % NewtonsDiffMtd_150517_.m (Assignment)
2. clear
3. clc
4. close all
5.
6. % Given Data points are:
7. X = [1.2 1.3 1.4 1.5 1.6 1.7]; % Values for X
8. Y = log10(X); % Y = f(X) = log(X)
9. x_i = 1.43; % value of interest to be
10. %interpolated
11.
12.
13. % First, we find the Newton Coefficient, a. That is,
14. %function a = newtonCoeff(X,Y)
15. n = length(X(:,1)); % no. of X data points to be read
16. a = Y; % Stores the value of Y data point
17. % Y data points in the array a
18. for k = 2:n % For loop iterates over each elmt
19. a(k:n) = (a(k:n) - a(k-1))./(X(k:n)- X(k-1));
20. end
21.
22. % Next, we program Newton Method Polynomial:
23. % Given by: p(x) = a0 + (x-x0)[a1 + (x-x1)[a2 + (x-x2)[a3 + (x-x3)[a4 +
24. % (x-x4)[a5]]]]]
25. % NOTE: This is achieved by deploying a backward recursion:
26.
27. p = a(n);
28. for k = 1:n-1;
29. p = a(n-k) + (x_i - X(n-k))*p;
30. end
31.
32. % Display the results:
33. disp(['Values of p: ' num2str(p)]);
*******************************************************************************

I used an algorithm provided in a book: "Numerical Methods for Engineering with MATLAB" by Jaan Kiusalaas.

Please, I would be glad if you could help me correct this so that the each values of x is displayed and that the value of interest (x_i) is used in the final result of p(x).

Thank you in advance.

[lastchance]

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


double fn( double x ) { return log10( x ); }                                        // function to approximate


double D( const vector<double> &xpt, int i, int j, double f( double x ) )           // successive differences
{
   if ( j == i ) return f( xpt[i] );
   else          return ( D( xpt, i + 1, j, f ) - D( xpt, i, j - 1, f ) ) / ( xpt[j] - xpt[i] );
}


vector<double> coeffs( const vector<double> &xpt, double f( double x ) )            // set coefficients
{
   int imax = xpt.size();
   vector<double> a( imax );
   for ( int i = 0; i < imax; i++ ) a[i] = D( xpt, 0, i, f ) ;
   return a;
}


double poly( const vector<double> &xpt, const vector<double> &a, double x )          // evaluate Newton polynomial
{
   int imax = xpt.size();
   double result = 0.0;
   for ( int i = imax - 1; i >= 0; i-- ) result = a[i] + ( x - xpt[i] ) * result;
   return result;
}


int main()
{
    // Set the interpolation points
    const int NPT = 6;
    const double xstart = 1.2, dx = 0.1;
    vector<double> xpt( NPT );
    for ( int i = 0; i < NPT; i++ ) xpt[i] = xstart + i * dx;

    // Find the coefficients
    vector<double> a = coeffs( xpt, fn );

    // Set some points to interpolate to (add the originals as a check)
    const int NCHECK = NPT * 2 - 1;
    vector<double> xcheck;
    for ( int i = 0; i < NCHECK; i++ ) xcheck.push_back( xstart + 0.5 * i * dx );
    #define SP << " " << setw( 15 ) <<
                              cout SP "x" SP "Polynomial"       SP "Actual" << '\n';
    for ( double x : xcheck ) cout SP  x  SP  poly( xpt, a, x ) SP  fn( x ) << '\n';
}

               x      Polynomial          Actual
             1.2       0.0791812       0.0791812
            1.25       0.0969099         0.09691
             1.3        0.113943        0.113943
            1.35        0.130334        0.130334
             1.4        0.146128        0.146128
            1.45        0.161368        0.161368
             1.5        0.176091        0.176091
            1.55        0.190332        0.190332
             1.6         0.20412         0.20412
            1.65        0.217484        0.217484
             1.7        0.230449        0.230449

[ayoesquire]

Hmm, brother

lastchance

, thank you for this. I'll certainly try it out. I am so grateful.

Also, I gave a friend of mine the initial Matlab code and he was able to spot the error in it. It was in line 7:

I had written:

X = [1.2 1.3 1.4 1.5 1.6 1.7];

when it should have been:

X = [1.2; 1.3; 1.4 1.5 1.6 1.7];

What was missing was the semi-colon i.e.: ";". And as you might guess, I tried it out and it worked perfectly fine.

Thank you.