[Math] Lagrange polynomial and derivative problem

derivativeserror functionMATLABpolynomials

I need to code this function in matlab

$$L'(x) = \sum_{k = 0}^{n} f_k l_k(x) $$

Where's $l_k$ looks like this

$$l_k(x) = \sum^{n}_{j=0, j \neq k} \frac{\prod^{n}_{i=0}(x – x_i)}{(x – x_k)(x-x_j)\prod^{n}_{i=0, i \neq k} (x_k – x_i)} $$

Well I got an array of $f_k, k = 1,\dots,n$ (there are results of error function $erf(x)$ in some points). So I wrote two functions.

First – to find $ l_k$ in some point (in arguments: res_arr array with $f_k$, $x$ – some point, k it's step in $\sum_{k}$)

function lDiffout = lDiff(res_arr, x, k)

n = size(res_arr,2);
result = 0;
p1 =1;
p2 =1;

%   П(x - x_i)
    for i = 1:n 
        p1 = p1*(x - res_arr(i));
    end;
%    П(x_k - x_i)   
    for i = 1:n 
        if(k ~= i)
          p2 = p2*(res_arr(k) - res_arr(i));
        end
    end;


for j=1:n
     p_res = p1/p2;  %П(x - x_i)/П(x_k - x_i) 
    if(j ~= k)
       p_rd = (x - res_arr(k))*(x - res_arr(j));% (x - x_k)(x - x_j)
       result = result + p_res/p_rd; % += П(x - x_i)/ (x - x_k)(x - x_j)П(x_k - x_i)
    end
end;

lDiffout = result;
end

Second to find result
(here is start, wstep, final in arguments it's because of x it's vector (in my situation from 0 to 2 with setp 0.2))

 function erdifOut = erfDiff(start, wstep, final, arr_res);
    clc
     cnt = 1;
      for xloop = start:wstep:final 
          erdifOut(cnt) = arr_res(cnt)*lDiff(arr_res, xloop, cnt);
          disp(['in x=', num2str(xloop), ' result is L = ',  num2str(erdifOut(cnt))]);
          cnt = cnt + 1;
      end;
    end

But it does not works correctly for some reason. My teacher said that i need to do something with $l_k$ function, make it simpler or something. Please how I can do it? Or maybe something wrong with the code?

Best Answer

In the following code, I assume all vectors are row vectors.

Try to avoid so many loops. MATLAB has plenty of built-in statements to do what you want to do.

To compute $\prod_{i=0}^n x - x_i$, try the following code:

% xi = [xi1 xi2 ... ] is a vector of your xi's. x is a scalar value.
delta_x = x - xi; % Now, delta_x is a vector, because a scalar minus a vector returns the element-wise subtraction.
Pn = prod(delta_x);

This will give you the product in the numerator.

Assuming $k$ is given as a parameter, we can then write the sum a few different ways. Here is the most straightforward way:

S = 0;
for j=1:n
    dxj = delta_x(j); % this is x-xj
    dxk = delta_x(k); % this is x-xk
    y = [xi(k)-xi(1:k-1) xi(k)-xi(k+1:end)]; 
    Pd = prod(y); % This is the product in the numerator
    if j ~= k
        S = S + Pn/(dxj*dxk*Pd);
    end
end

This should give you what you want for $l_k$.

We can actually do better (see my comment):

dxk = delta_x(k); % this is x-xk
y = [xi(k)-xi(1:k-1) xi(k)-xi(k+1:end)]; 
Pd = prod(y); % This is the product in the numerator
delta_x1 = [delta_x(1:k-1) delta_x(k+1:end)];
lk = Pn/(dxk*Pd)*sum(delta_x1.^(-1));

These two code snippets are equivalent.

Related Question