Lagrange basis polynomials subroutine

Help
1

Hi Fortran expert here,

I am trying to learn the basic of FEM with John Burkardt’s 1D FEM1D_LAGRANGE_PRB.
I have some difficulty to make sense of a few line of codes.
below is the subroutine for evaluation of lagrange basis polynomials.

regarding below code:

 li(i,j) = product ( ( xi(i) - xd(1:j-1)  ) / ( xd(j) - xd(1:j-1)  ) ) &
              * product ( ( xi(i) - xd(j+1:nd) ) / ( xd(j) - xd(j+1:nd) ) )

why is there two product here? should we just make:
li(i,j) = product ( ( xi(i) - xd(1:j-1) ) / ( xd(j) - xd(1:j-1) ) ) , ( ( xi(i) - xd(j+1:nd) ) / ( xd(j) - xd(j+1:nd) ) )
basically multiply two matrix?

sorry for bring up this very basic question, there are no people around that know enough of FEM.

Thanks

Martin

 !*****************************************************************************80
!
!! LAGRANGE_VALUE evaluates the Lagrange basis polynomials.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    09 October 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) ND, the number of data points.
!
!    Input, real ( kind = 8 ) XD(ND), the interpolation nodes.
!
!    Input, integer ( kind = 4 ) NI, the number of evaluation points.
!
!    Input, real ( kind = 8 ) XI(NI), the evaluation points.
!
!    Output, real ( kind = 8 ) LI(NI,ND), the value, at the I-th point XI, 
!    of the Jth basis function.
!
  implicit none

  integer ( kind = 4 ) nd
  integer ( kind = 4 ) ni

  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  real ( kind = 8 ) li(ni,nd)
  real ( kind = 8 ) xd(nd)
  real ( kind = 8 ) xi(ni)
  
  do i = 1, ni
    do j = 1, nd
      li(i,j) = product ( ( xi(i) - xd(1:j-1)  ) / ( xd(j) - xd(1:j-1)  ) ) &
              * product ( ( xi(i) - xd(j+1:nd) ) / ( xd(j) - xd(j+1:nd) ) )
    end do
  end do

  return
end subroutine lagrange_value
2

The function is a direct implementation of the definition of a Lagrange polynomial.

The Fortran PRODUCT() function returns the product of all the elements in an entire array. It doesn’t multiply two matrices.

The wikipedia definition is for a single data point x and indexes the interpolation nodes from 0 to k.

The Fortran function

  • has NI evaluation points x=XI(i), i=1, …, NI
  • indexes the interpolation points x=XD(j), j=1,…, ND

li(i,j) is the j-th polynomial for data point xi(i). The two products are for the terms 1,…,j-1 and j+1,…ND. The arguments of the PRODUCT() functions are vectors. The Fortran standard cleverly defines a product of a zero length vector = 1 so it all just magically works.

I would suggest the best way to “get it” is to do a simple example by hand. Start with NI=1 and ND=2 or 3.

3

If you look at the formula that defines l_j(x) on the Wikipedia page, you will observe that the product is over m = 0…k, skipping m = j. In other words, it has two parts that are multiplied together. The first part covers m = 0…j-1, and the second part covers m = j+1…k. The term for m = j has to be skipped.

4

Hi DavidB,
Thanks for your help. I understand it now. showhow I am confused with another function: dot_product.

Best regards

Martin

5

Hi Mecej4,
It is a quite clever trick to split the two series, I see it after your explanation, many thanks.

Martin