7
$\begingroup$

Assume, that we have points $x_i$ with $i=1,...,N+1$. We construct the Lagrange basis polynomials as \begin{align} L_j(x) = \prod_{k\not = j} \frac{x-x_k}{x_j-x_k} \end{align}

Now according to my computation and the results by Yves Daoust here, the derivative of $L_i$ can be computed as \begin{align} L'_j(x) = L_j(x)\cdot \sum_{k\not = j}\frac{1}{x_k-x_j} \end{align}

I try to reproduce the numerical results of a paper, and for this results the authors use the derivative matrix $D$ with $D_{ij} =L_i'(x_j)$.

The authors use the Legendre Gauss Lobatto quadrature points, plotted below. enter image description here Now I can easily construct the basis polynomials, also plotted below, together with the quadrature points. enter image description here

Given the concrete choice of basis points, the plots suggest, that $L'_j(x_j)=0$ except for the first and last basis function. Additionally it seems, that $L'_j(x_i)\not =0$ for $i\not = j$. But using the derived formula from above, we obtain \begin{align} L'_j(x_i) = L_j(x_i)\cdot \sum_{k\not = j}\frac{1}{x_k-x_j}=0 \end{align} for $i\not = j$ since $L_j(x_i) =\delta_{ij}$, which contradicts my observation.

I used the following Matlab functions, to construct the matrix $D$.

function y=dl(i,x,z)
n = length(x);
y = 0;
for m=1:n
    if not(m==i)
        y = y + 1/(x(m)-x(i));
    end
end
size(y)
y = y*l(i,x,z);

end



function y=l(i,x,z)

n = length(x);
% computes h_i(z)

y = 1;
for m=1:n
    if not(m==i)
        y = y.*(z-x(m))./(x(i)-x(m));
    end
end
end

Where dland lare $L'$ and $L$. $D$ is then constructed by

D = zeros(M+1,M+1);
for i=1:M+1
    for j=1:M+1
        D(i,j) =dl(i,X,X(j));
    end
end

which gives the matrix

D =

   10.5000         0         0         0         0         0         0
         0   -0.0000         0         0         0         0         0
         0         0    0.0000         0         0         0         0
         0         0         0   -0.0000         0         0         0
         0         0         0         0    0.0000         0         0
         0         0         0         0         0    0.0000         0
         0         0         0         0         0         0  -10.5000

Here I agree with the zeros on the diagonal except for the first and last element. The zeros everywhere else agree with the formula, but they don't agree with my observation.

If I use finite differences (which is no solution for the real implementation) in the following way: FD(i,j) =(l(i,X,X(j)+eps)-l(i,X,X(j)-eps))/(2*eps); I obtain the output

FD =

  -10.5000   -2.4429    0.6253   -0.3125    0.2261   -0.2266    0.5000
   14.2016    0.0000   -2.2158    0.9075   -0.6164    0.6022   -1.3174
   -5.6690    3.4558    0.0000   -2.0070    1.0664   -0.9613    2.0500
    3.2000   -1.5986    2.2667         0   -2.2667    1.5986   -3.2000
   -2.0500    0.9613   -1.0664    2.0070   -0.0000   -3.4558    5.6690
    1.3174   -0.6022    0.6164   -0.9075    2.2158   -0.0000  -14.2016
   -0.5000    0.2266   -0.2261    0.3125   -0.6253    2.4429   10.5000

which is zero on the diagonal (except first and last) and nonzero everywhere else.

So here are my questions:

Did I make an error in the implementation?

Is my formula derived under the assumption that $x$ is not a basis point?

Can I modify my code to obtain the results that my finite difference code produces?

Is my assumption correct, that rather the results of the finite difference computation are "correct"?

EDIT

It seems, like the derivation of the formula goes wrong for $x=x_i$. If I compute the derivative without simplification, I get \begin{align} L_j'(x) = \sum_{l\not = j} \frac{1}{x_j-x_l}\prod_{m\not = (j,l)} \frac{x-x_m}{x_j-x_m} \end{align} Or in code

function y = alternative_dl(j,x,z)

y = 0;
n = length(x);
for l=1:n
    if not(l==j)
        k = 1/(x(j)-x(l));
        for m=1:n
            if not(m==j) && not(m==l)
                k = k*(z-x(m))/(x(j)-x(m));
            end
        end
        y = y + k;
    end
end

end

Which agrees with the finite difference computation.

So it seems to me, that simplifying the above formula includes some "hidden division by zero" if $x=x_i$.

$\endgroup$
  • $\begingroup$ I found your post and the last expression for derivative really useful. Can you look at this post and help me out? I need to find second order derivative and have no idea how. Thanks in advance! math.stackexchange.com/questions/1572576/… $\endgroup$ – Luka Bulatovic Dec 12 '15 at 20:00
4
$\begingroup$

I think your implementation is correct. You copied the wrong formula.

$$L'_i(x)=L_i(x) \sum_{m=0,\ m\neq i}^n\frac1{x-x_m}$$

This formula would not give you zeros everywhere else.

And this is the same as your new derived formula.

$\endgroup$
  • $\begingroup$ But your formula contains $L_i(x)$ on the right hand side. Now if I evaluate $L'_i(x_j)$ I get $L_i'(x_j) = L_i(x_j)\cdot \sum...=0\cdot \sum$ since $L_i(x_j)=0$ if $i\not = j$. $\endgroup$ – k1next Jan 15 '15 at 11:39
  • $\begingroup$ Okay I see, that I copied the formula wrong, but the problem I have in the comment above still exists?! $\endgroup$ – k1next Jan 15 '15 at 11:42
  • $\begingroup$ When you plug $x_j$ into it, there is also a $x_j-x_j$ in the denominator on the right hand side. I guess that's what you called "hidden division by zero"? In fact, you can expand this formula to get the other formula in the post you referred to: $L_i(x)'=\frac{\sum_{k=0}^{n} \prod_{l=0, l \neq k}^n (x-x_l)}{\prod_{k=0, k \neq j}^n (x_j-x_k)}$, which does not contain zero in the denominator. $\endgroup$ – KittyL Jan 15 '15 at 11:44
  • $\begingroup$ Here is a simple 3 point formula. $L_1'(x)=\frac{(x-x_0)+(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}(\frac{1}{x-x_0}+\frac{1}{x-x_2})$. If you evaluate it at $x_2$, you will see you cannot just say it is zero using the second formula. $\endgroup$ – KittyL Jan 15 '15 at 11:47
  • $\begingroup$ Okay, but then I can't use the formula in my implementation, since I end up computing $0/0$ since $(x-x_2)=0$ if $x=x_2$ and $1/(x-x2)=1/0$ if $x=x_2$? $\endgroup$ – k1next Jan 15 '15 at 12:16
6
$\begingroup$

Just to complement the answer, here is the formula for second derivative: $$ L^{"}_i(x) =\sum_{l\ne i}\frac{1}{x_i-x_l}\left( \sum_{m\ne(i,l)}\frac{1}{x_i-x_m}\prod_{k\ne(i,l,m)}\frac{x-x_k}{x_i-x_k}\right) $$ through recursion, one can compute further higher derivatives.

$\endgroup$
  • $\begingroup$ Look at this guy... what a nicccceeee answer. $\endgroup$ – bremen_matt Feb 22 '18 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.