Lagrange Polynomials construction for the stochastic collocation methods

[ma6yu]

Hello,

I am constructing the Lagrange Polynomials using lagrange_polynomial function in chaospy and have a problem. I have 5 sampling points, the problem is 2D. I am writing as:

Lagrange = cp.lagrange_polynomial(nodes)

and get: File "/usr/local/lib/python3.5/dist-packages/chaospy/orthogonal/lagrange.py", line 47, in lagrange_polynomial matrix = numpy.prod(abscissas.T[idx]**indices[idy], -1) IndexError: index 3 is out of bounds for axis 0 with size 3

The "nodes" is [[ 0.17 0.15267949 0.17 0.18732051 0.17 ] [ 14.94015203 16.68992088 16.68992088 16.68992088 19.83651501]]

Are there any restriction of the relationship between the number of samples and the dimension?

best regards, Mayu

[jonathf]

There is two things going on here:

1. there is a bug in the code. I have fixed and pushed an new release version 3.0.12.
2. Lagrange is known for being a bit unstable for some combination of nodes. Looks like that is the case for you.

Hope this helps.

[ma6yu]

hi,

I updated the code to the new version(3.0.13), it works as you mentioned above, thanks. But I still have a problem using the lagrange_polynomial function.

Now I leave my unstable nodes for now and try another case which should be stable. The case is mentioned in the following literature, Example 2 on page 4. (http://evoq-eval.siam.org/Portals/0/Publications/SIURO/Vol1_Issue1/A_Simple_Expression_for_Multivariate.pdf?ver=2018-03-30-130233-050)

my code is:

nodes = np.array([[0.,2.,1.,-2.,-3., -1.], [1.,1.,3.,-1.,2.,2,]]) Lagrange = cp.lagrange_polynomial(nodes) print(Lagrange(*nodes))

I was expecting: [[1. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. ] [0. 0. 1. 0. 0. 0. ] [0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 1. 0. ] [0. 0. 0. 0. 0. 1. ] ]

but I got: [[-10.6 -3.6 -20.8 23.6 -1.2 -17.2 ] [ -4.22 -2.22 -7.96 8.62 0.76 -6.04] [ -3.44 -1.44 -6.92 6.24 0.52 -5.08] [ -1.98 -0.98 -3.64 2.58 -2.16 -3.36] [ -3.9 -1.9 -8.2 6.9 -1.3 -6.3 ] [-12.74 -4.74 -25.32 24.54 -1.58 -20.18]]

I am checking some parameters in the lagrange_polynomial, the "matrix" in the lagrange_polynomial is as same as the "sample matrix" in the literature. I think, the alignment sequence of sort type "G" in the chaospy.bertran.bindex and chaospy.poly.basis are different. (if I print out the "indices" and "vec" in the lagrange_polynomial, I get:

indices: [[0 0] [0 1] [1 0] [0 2] [1 1] [2 0]] vec: [1, q1, q0, q0q1, q0^2, q1^2]

but even when I change the sequence of the indices as same as the vec, I still don't get the expected result of the Lagrange polynomial. Am I using the code wrong or do you have any idea what is the problem?

best, Mayu

[jonathf]

That isn't right. I'll take a look at it.

[ma6yu]

I changed

• the sequence of the indices as same as the "vec"

• determinant calculation of "matrix" and got the expected result. I think the determinant calculation is missing change of sign. So the code line 69-76 should be:

  else:
      for i in range(size):
          ## begin test ##
          if i%2 != 0:
              k = 1
          else:
              k=0
          ## end test ##
          for j in range(size):
              #coeffs[i, j] += numpy.linalg.det(matrix[1:, 1:]) #org
  
              ## begin test ##
              if k%2 == 0:
                  coeffs[i, j] += numpy.linalg.det(matrix[1:, 1:])
              else:
                  coeffs[i, j] += -numpy.linalg.det(matrix[1:, 1:])
              ## end test ##
  
              matrix = numpy.roll(matrix, -1, axis=0) 
              k += 1 #test
          matrix = numpy.roll(matrix, -1, axis=1) 

The additional part is indicated by #test/###begin test###~###end test### and one original line is commented out. Though, there might be smarter way to do it but it gives at least correct coefficients. I don't want to touch the original part of chaospy.bertran.bindex and chaospy.poly.basis about inconsistency of the sequence of the indices,(for now I am changing the sequence manually inside the lagrange_polynomial function), I wait for you to take a look.

best, Mayu

[jonathf]

Lagrange polynomials are hard, and not my strongest suite, and I have likely taken too lightly on the matter first time around. It is also a few years back since I touched the topic, so I am likely a bit rusty.

So the example provided in the doctest works, but clearly isn't good enough to capture the examples you are now providing.

Just to note, I am not getting your code to work either. Do you mind proposing it as a pull request?

I am willing to do my due diligence and write the code if needed. But as noted above, this isn't my strongest suite. If we can get your code to work, it would be preferable. Otherwise, if you have a references that is worth looking at it would be very appreciated. If not your approach, I'll likely dig into the implementation of Burkardt to see how they did it and use that as a baseline. Not sure how complicated it will be yet.

https://people.sc.fsu.edu/~jburkardt/m_src/lagrange_nd/lagrange_nd.html