MATLAB: Cant get integral from Lagrange interpolation polynom

Question

I have writed a function to get integral from interpolation Lagrange polynom.

function fxn = lagrange()
clc; clear;
k = [0 1 2 3 4 5 6 7 8 9 10];
xpoints = 0.3 * k;
ypoints = 1 - exp(-xpoints);
L = @(x) 0;
for i = 1 : 11
    Li = @(x) 1;
    Ld = @(x) 1;
    for j = 1 : 11
        if i ~= j
            Li = @(x) Li(x) * (x - xpoints(j));
            Ld = @(x) Ld(x) * (xpoints(i) - xpoints(j));
        end
    end
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
end
fprintf("integral from L(x): %f",integral(L, 0, 3));
fxn = @(x) L(x);
end

But when i run this code I have this message appears:

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To
perform elementwise multiplication, use '.*'.
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)Li(x)*(x-xpoints(j)) (line 12)
            Li = @(x) Li(x) * (x - xpoints(j));
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in lagrange>@(x)L(x)+(Li(x)/Ld(x))*ypoints(i) (line 16)
    L = @(x) L(x) + (Li(x)/Ld(x))*ypoints(i);
Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);
Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);
Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);
Error in lagrange (line 23)
fprintf("integral from L(x): %f",integral(L, 0, 3));

What i do wrong? How can i fix this?

Answer 1

Ok. As I said, this is a misuse of recursively built function handles. It is easy enough to build a Lagrange polynomial. First, you need to learn how to use a function. I'll create xpoints and ypoints as variables in my base workspace.

k = 0:10;
xpoints = 0.3 * k;
ypoints = 1 - exp(-xpoints);

So, we have two vectors of points, as (x,y) pairs. Now, build a function that does a Lagrange interpolation. So you will pass it any value x or a vector of values x. It will return the interpolant. There is no need for this massive recursive set of function handles. On the other hand, it seems you want to do this using function handles. So let me see how I might write this as a function, first.

function y = LagrangeInterp(x,xpoints,ypoints)
  y = zeros(size(x));
  N = length(xpoints);
  for j = 1:N
    L = ones(size(x));
    for m = setdiff(1:N,j)
      L = L.*(x - xpoints(m))/(xpoints(j) - xpoints(m));
    end
    y = y + L;
  end
end

Note there is no need to put clear and clc inside the beginning of a function. As well, you should write a function to take input variables, as they will be passed in. Now, we can use this to evaluate the N-point Lagrange interpolant.

Also, note the nice trick of

for m = setdiff(1:N,j)

to create a loop that variable that takes on all values of 1:N, EXCEPT that it skips j. That is how I would write a Lagrange interpolant. Save the function as an m-file on your search path. We can now use it. At the command line, you can now create a function handle, as simply as:

Lagfun = @(x) LagrangeInterp(x,xpoints,ypoints);

We can evaluate it at any point, or list of points, as simply as

Lagfun(1.35)
ans =
         0.740759739304464

We can plot it. It actually looks reasonable. (Always dangerous with Lagrange interpolating polynomials.)

fplot(Lagfun,[0,3])
hold on 
plot(xpoints,ypoints,'ro')

And since it is a function handle, we can call integral with it.

integral(Lagfun,0,3)
ans =
          2.04978706821567
          
integral(@(x) 1-exp(-x),0,3)
ans =
          2.04978706836786

Not too far off. We can even get tricky, and extract the polynomial itself in symbolic form, just by passing in a symboic variable as x.

syms X
vpa(Lagfun(X),5)
ans =
8.1804*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 2.1) - 3.3233*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 1.8)*(X - 2.1) + 0.044345*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 0.9)*(X - 1.8)*(X - 2.1) - 9.1364*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) - 0.43532*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) + 1.9096*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) + 6.8486*X*(X - 0.3)*(X - 0.6)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) + 0.94753*X*(X - 0.3)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) - 0.12096*X*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)*(X - 2.1) - 4.9144*X*(X - 0.3)*(X - 0.6)*(X - 1.2)*(X - 2.4)*(X - 2.7)*(X - 1.5)*(X - 3.0)*(X - 0.9)*(X - 1.8)
vpa(expand(Lagfun(X)),5)
ans =
- 6.3837e-8*X^10 + 1.6007e-6*X^9 - 0.000020676*X^8 + 0.00018853*X^7 - 0.0013729*X^6 + 0.0083158*X^5 - 0.041654*X^4 + 0.16666*X^3 - 0.5*X^2 + 1.0*X

So that seems to do what you need. But could we have written things using function handles more aggressively and recursively? Hmm. Let me take a shot, and see if I can fix the code you wrote. (Shiver.)

function L = lagrange(xpoints,ypoints)
L = @(x) zeros(size(x));
N = length(xpoints);
for j = 1 : N
  Li = @(x) ones(size(x));
  Ld = 1;
  for m = setdiff(1:N,j)
    Li = @(x) Li(x) .* (x - xpoints(m));
    % see that LD is not a function of x
    Ld = Ld * (xpoints(j) - xpoints(m));
  end
  L = @(x) L(x) + (Li(x)/Ld)*ypoints(j);
end

Does it work? Again, save lagrange as an m-file on your search path.

fxn = lagrange(xpoints,ypoints)
fxn =
  function_handle with value:
    @(x)L(x)+(Li(x)/Ld)*ypoints(j)
    
fxn(1.35)
ans =
         0.740759739304464
         
integral(fxn,0,3)
ans =
          2.04978706821567
          

So it does work, and is still properly vectorized. See that I used zeros and ones in my creation of L and Li. That was important. My point is though, that this is a nest of function handles, that each call function handles with the same name that you created just a millisecond before in nested fashion. Something like that will be pure hell to debug. It is a birds nest of recursive function handles, something that still makes me shiver.