Creating a Lagrange Polynomial with Unknown 'x.'

I'm trying to create a Lagrange polynomial 'P(x)' for a set of points, but without the x-coordinate of the point I'm trying to find. For those who don't know, a Lagrange polynomial is a method of interpolating a set of data, such that it creates a polynomial that goes through every point in the set of data. The polynomial is formed from each of the x and y values of the data points, as well as the variable 'x,' which represents the x-coordinate of a point you are trying to find the y-value for.

However, in my situation, I have the y-value (y = 43) of a point I am trying to find an x-value for (which will later be done with a root-finding method). I initially thought to use symbolic calculations, as can be seen below:

y0 = 42.997791034585575 ;
y1 = 42.999274955274394 ;
y2 = 43.000758746700800 ;
y3 = 43.002242408852400 ;
x0 = 3.3335 ;
x1 = 3.3336 ;
x2 = 3.3337 ;
x3 = 3.3338 ; % I have other calculations to get these variables but I left that out for simplicity

syms x 

L0 = ((x - x1)/(x0 - x1))*((x - x2)/(x0 - x2))*((x - x3)/(x0 - x3)) ;
L1 = ((x - x0)/(x1 - x0))*((x - x2)/(x1 - x2))*((x - x3)/(x1 - x3)) ;
L2 = ((x - x0)/(x2 - x0))*((x - x1)/(x2 - x1))*((x - x3)/(x2 - x3)) ;
L3 = ((x - x0)/(x3 - x0))*((x - x1)/(x3 - x1))*((x - x2)/(x3 - x2)) ;

P1 = y0 * L0 + y1 * L1 + y2 * L2 + y3 * L3 % Symbolic expression of P(x)

P1 = 

P2 = sym2poly(P1) % Converting symbolic expression of P(x) to vector expression

P2 = 1×4

-2.0651 14.1836 -10.8771 -1.8562

P3 = P2 ;
P3(4) = P3(4) - 43 ;

From this, I was given a polynomial (rounded to two decimal places) of (This is slightly different to what the above code will give as I assigned values to each of the x and y-coordinates). I need to then 'lower' the polynomial so that the x-coordinate of the point I am trying to find lies on the x-axis, which I did by converting the symbolic polynomial to a vector expression (a 1x4 matrix in the form of [-2.06 14.15 -10.82 -1.88] and subtracting 43 (the y-coordinate of the point) from the fourth term.

This worked perfectly and gave me exactly what I needed, however I later realised we weren't allowed to use symbolic calculations in this process. I really don't know where to go from here. I tried assigning each of L0, L1, L2 and L3 to function handles:

L0 = @(x) (((x) - x1)/(x0 - x1))*(((x) - x2)/(x0 - x2))*(((x) - x3)/(x0 - x3)) ;
L1 = @(x) (((x) - x0)/(x1 - x0))*(((x) - x2)/(x1 - x2))*(((x) - x3)/(x1 - x3)) ;
L2 = @(x) (((x) - x0)/(x2 - x0))*(((x) - x1)/(x2 - x1))*(((x) - x3)/(x2 - x3)) ;
L3 = @(x) (((x) - x0)/(x3 - x0))*(((x) - x1)/(x3 - x1))*(((x) - x2)/(x3 - x2)) ;

Unfortunately this didn't work as when I went to form 'P1' I got the error that operations using '*' could not be performed on function handles. I suppose I COULD hypothetically just skip forming L0, L1, L2 and L3 and then P1 by just putting all together in one line, but this would be very messy and harder for a reader to interpret.

Thank you for being so patient as to read through this really long question (I just wanted to give as many details as I could to be as clear as possible), and any ideas are welcome.

Answer 1

[John D'Errico]

I can understand where you went wrong, as it is a common mistake. What worked for symbolic variables, did not work for function handles. Thus, you CANNOT add function handles. However, you CAN add the result of a set of function handles. So you could have done this::

L = @(x) L0(x) + L1(x) + L2(x) + L3(x);

you should see the I am evaluating each function handle, at the unknown value x, and then adding the outputs.

The result is perfectly legal syntax in MATLAB, and now you can evaluate the function handle L(x) at any point x.

[Isabelle Davies]

Great! Thank you so much, it gave the correct answer! I there anyway to show what the polynomial is before evaluating it for a certain value of x?

[Torsten]

To compare your results, you can use polyfit and polyval to construct the polynomial from your interpolation values x_i and y_i and evaluate it for new points x.

[Isabelle Davies]

No sorry, I didn't mean to compare results, but once I have this:

L0 = @(x) (((x) - x1)/(x0 - x1))*(((x) - x2)/(x0 - x2))*(((x) - x3)/(x0 - x3)) ;
L1 = @(x) (((x) - x0)/(x1 - x0))*(((x) - x2)/(x1 - x2))*(((x) - x3)/(x1 - x3)) ;
L2 = @(x) (((x) - x0)/(x2 - x0))*(((x) - x1)/(x2 - x1))*(((x) - x3)/(x2 - x3)) ;
L3 = @(x) (((x) - x0)/(x3 - x0))*(((x) - x1)/(x3 - x1))*(((x) - x2)/(x3 - x2)) ;

P1 = @(x) y0 * L0(x) + y1 * L1(x) + y2 * L2(x) + y3 * L3(x)

P2 = P1 ; 
P2 = @(x) P2(x) - 43 ;

Is there any way to explicitly show that ? Or if not, to plot the polynomial? In terms of plotting, I tried the below (after changing all '*' to '.*' to account for x being a vector)

x = -2:0.1:8 ;
plot(x, polyval(P2(x), x))
xlim([-2 8])
ylim([-40 60])

I'm not sure what went wrong, but I got this graph:

I feel like plotting a graph has something to do with polyval, but I'm not sure where I'm going wrong.

[Torsten]

y0 = 42.997791034585575 ;
y1 = 42.999274955274394 ;
y2 = 43.000758746700800 ;
y3 = 43.002242408852400 ;
x0 = 3.3335 ;
x1 = 3.3336 ;
x2 = 3.3337 ;
x3 = 3.3338 ; % I have other calculations to get these variables but I left that out for simplicity

L0 = @(x) (x - x1)/(x0 - x1).*(x - x2)/(x0 - x2).*(x - x3)/(x0 - x3) ;
L1 = @(x) (x - x0)/(x1 - x0).*(x - x2)/(x1 - x2).*(x - x3)/(x1 - x3) ;
L2 = @(x) (x - x0)/(x2 - x0).*(x - x1)/(x2 - x1).*(x - x3)/(x2 - x3) ;
L3 = @(x) (x - x0)/(x3 - x0).*(x - x1)/(x3 - x1).*(x - x2)/(x3 - x2) ;
P1 = @(x) y0 * L0(x) + y1 * L1(x) + y2 * L2(x) + y3 * L3(x)

P1 = function_handle with value:
    @(x)y0*L0(x)+y1*L1(x)+y2*L2(x)+y3*L3(x)

P2 = P1 ; 
P2 = @(x) P2(x) - 43 ;
x = 2:0.1:4 ;
plot(x, P2(x),[x0 x1 x2 x3],[P2(x0) P2(x1) P2(x2) P2(x3)],'*r')
xlim([3 4])
ylim([-8 8])
grid on

[Torsten]

Or:

y0 = 42.997791034585575 ;
y1 = 42.999274955274394 ;
y2 = 43.000758746700800 ;
y3 = 43.002242408852400 ;
x0 = 3.3335 ;
x1 = 3.3336 ;
x2 = 3.3337 ;
x3 = 3.3338 ; % I have other calculations to get these variables but I left that out for simplicity
syms x
L0 =  (x - x1)/(x0 - x1)*(x - x2)/(x0 - x2)*(x - x3)/(x0 - x3) ;
L1 =  (x - x0)/(x1 - x0)*(x - x2)/(x1 - x2)*(x - x3)/(x1 - x3) ;
L2 =  (x - x0)/(x2 - x0)*(x - x1)/(x2 - x1)*(x - x3)/(x2 - x3) ;
L3 =  (x - x0)/(x3 - x0)*(x - x1)/(x3 - x1)*(x - x2)/(x3 - x2) ;
P1 =  y0 * L0 + y1 * L1 + y2 * L2 + y3 * L3;
P2 = P1;
P2 = P2 - 42;
P22 = vpa(expand(P2),4)

P22 = 

xnum = 2:0.1:4 ;
plot(xnum, double(subs(P2,x,xnum)),[x0 x1 x2 x3],[double(subs(P2,x,x0)), double(subs(P2,x,x1)) double(subs(P2,x,x2)) double(subs(P2,x,x3))],'*r')
xlim([3 4])
ylim([-8 8])
grid on