lagrange approximation -c++

Question

I updated the code. What i am trying to do is to hold every lagrange's coefficient values in pointer d.(for example for L1(x) d[0] would be "x-x2/x1-x2" ,d1 would be (x-x2/x1-x2)*(x-x3/x1-x3) etc.

My problem is

1) how to initialize d ( i did d[0]=(z-x[i])/(x[k]-x[i]) but i think it's not right the "d[0]"

2) how to initialize L_coeff. ( i am using L_coeff=new double[0] but am not sure if it's right.

The exercise is: Find Lagrange's polynomial approximation for y(x)=cos(π x), x ∈−1,1 using 5 points (x = -1, -0.5, 0, 0.5, and 1).

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>

using namespace std;

const double pi=3.14159265358979323846264338327950288;

// my function
double f(double x){

return (cos(pi*x));

}


//function to compute lagrange polynomial
double lagrange_polynomial(int N,double *x){
//N = degree of polynomial

double z,y;
double *L_coeff=new double [0];//L_coefficients of every Lagrange L_coefficient
double *d;//hold the polynomials values for every Lagrange coefficient
int k,i;
//computations for finding lagrange polynomial
//double sum=0;

for (k=0;k<N+1;k++){
        for ( i=0;i<N+1;i++){
            if (i==0) continue;
            d[0]=(z-x[i])/(x[k]-x[i]);//initialization
            if (i==k) L_coeff[k]=1.0;
            else if (i!=k){
            L_coeff[k]*=d[i];

                      }

           }

cout <<"\nL("<<k<<") = "<<d[i]<<"\t\t\tf(x)= "<<f(x[k])<<endl;
}

}

int main()
{

    double deg,result;
    double *x;


    cout <<"Give the degree of the polynomial :"<<endl;
    cin >>deg;

    for (int i=0;i<deg+1;i++){
    cout <<"\nGive the points of interpolation : "<<endl;
    cin >> x[i];
    }

    cout <<"\nThe Lagrange L_coefficients are: "<<endl;
    result=lagrange_polynomial(deg,x);



    return 0;
}

Here is an example of lagrange polynomial

Answer 1

As this seems to be homework, I am not going to give you an exhaustive answer, but rather try to send you on the right track.

• How do you represent polynomials in a computer software? The intuitive version you want to archive as a symbolic expression like 3x^3+5x^2-4 is very unpractical for further computations.
• The polynomial is defined fully by saving (and outputting) it's coefficients.

What you are doing above is hoping that C++ does some algebraic manipulations for you and simplify your product with a symbolic variable. This is nothing C++ can do without quite a lot of effort.

You have two options:

• Either use a proper computer algebra system that can do symbolic manipulations (Maple or Mathematica are some examples)
• If you are bound to C++ you have to think a bit more how the single coefficients of the polynomial can be computed. You programs output can only be a list of numbers (which you could, of course, format as a nice looking string according to a symbolic expression).

Hope this gives you some ideas how to start.

Edit 1

You still have an undefined expression in your code, as you never set any value to y. This leaves prod*=(y-x[i])/(x[k]-x[i]) as an expression that will not return meaningful data. C++ can only work with numbers, and y is no number for you right now, but you think of it as symbol.

You could evaluate the lagrange approximation at, say the value 1, if you would set y=1 in your code. This would give you the (as far as I can see right now) correct function value, but no description of the function itself.

Maybe you should take a pen and a piece of paper first and try to write down the expression as precise Math. Try to get a real grip on what you want to compute. If you did that, maybe you come back here and tell us your thoughts. This should help you to understand what is going on in there.

And always remember: C++ needs numbers, not symbols. Whenever you have a symbol in an expression on your piece of paper that you do not know the value of you can either find a way how to compute the value out of the known values or you have to eliminate the need to compute using this symbol.

P.S.: It is not considered good style to post identical questions in multiple discussion boards at once...

Edit 2

Now you evaluate the function at point y=0.3. This is the way to go if you want to evaluate the polynomial. However, as you stated, you want all coefficients of the polynomial.

Again, I still feel you did not understand the math behind the problem. Maybe I will give you a small example. I am going to use the notation as it is used in the wikipedia article.

Suppose we had k=2 and x=-1, 1. Furthermore, let my just name your cos-Function f, for simplicity. (The notation will get rather ugly without latex...) Then the lagrangian polynomial is defined as

f(x_0) * l_0(x) + f(x_1)*l_1(x)

where (by doing the simplifications again symbolically)

l_0(x)= (x - x_1)/(x_0 - x_1) = -1/2 * (x-1) = -1/2 *x  + 1/2
l_1(x)= (x - x_0)/(x_1 - x_0) = 1/2 * (x+1)  = 1/2 * x  + 1/2

So, you lagrangian polynomial is

  f(x_0) * (-1/2 *x  + 1/2) + f(x_1) * 1/2 * x  + 1/2
= 1/2 * (f(x_1) - f(x_0)) * x     +     1/2 * (f(x_0) + f(x_1))

So, the coefficients you want to compute would be 1/2 * (f(x_1) - f(x_0)) and 1/2 * (f(x_0) + f(x_1)).

Your task is now to find an algorithm that does the simplification I did, but without using symbols. If you know how to compute the coefficients of the l_j, you are basically done, as you then just can add up those multiplied with the corresponding value of f.

So, even further broken down, you have to find a way to multiply the quotients in the l_j with each other on a component-by-component basis. Figure out how this is done and you are a nearly done.

Edit 3

Okay, lets get a little bit less vague.

We first want to compute the L_i(x). Those are just products of linear functions. As said before, we have to represent each polynomial as an array of coefficients. For good style, I will use std::vector instead of this array. Then, we could define the data structure holding the coefficients of L_1(x) like this:

std::vector L1 = std::vector(5); 
// Lets assume our polynomial would then have the form 
// L1[0] + L2[1]*x^1 + L2[2]*x^2 + L2[3]*x^3 + L2[4]*x^4

Now we want to fill this polynomial with values.

// First we have start with the polynomial 1 (which is of degree 0)
// Therefore set L1 accordingly:
L1[0] = 1;
L1[1] = 0; L1[2] = 0; L1[3] = 0; L1[4] = 0;
// Of course you could do this more elegant (using std::vectors constructor, for example)
for (int i = 0; i < N+1; ++i) {
    if (i==0) continue;  /// For i=0, there will be no polynomial multiplication
    // Otherwise, we have to multiply L1 with the polynomial
    // (x - x[i]) / (x[0] - x[i])
    // First, note that (x[0] - x[i]) ist just a scalar; we will save it:
    double c = (x[0] - x[i]);
    // Now we multiply L_1 first with (x-x[1]). How does this multiplication change our 
    // coefficients? Easy enough: The coefficient of x^1 for example is just
    // L1[0] - L1[1] * x[1]. Other coefficients are done similary. Futhermore, we have
    // to divide by c, which leaves our coefficient as
    // (L1[0] - L1[1] * x[1])/c. Let's apply this to the vector:
    L1[4] = (L1[3] - L1[4] * x[1])/c;
    L1[3] = (L1[2] - L1[3] * x[1])/c;
    L1[2] = (L1[1] - L1[2] * x[1])/c;
    L1[1] = (L1[0] - L1[1] * x[1])/c;
    L1[0] = (      - L1[0] * x[1])/c; 
    // There we are, polynomial updated.
}

This, of course, has to be done for all L_i Afterwards, the L_i have to be added and multiplied with the function. That is for you to figure out. (Note that I made quite a lot of inefficient stuff up there, but I hope this helps you understanding the details better.)

Hopefully this gives you some idea how you could proceed.