Weights and interpolation -- Replicating Python code with hand calculations

[confused_engineer]

Hello everyone. I have a Python code which calculates, given a continuos uniform random variable U(-1,1), the order of a interpolation polynomial and a set of points the evolution of a function of this random variable. i.e.
v0 = cp.Uniform(-1,1)
t = np.linspace(0, 10, 10)
order=1
.
.
.
plt.plot (t, v0*t)

and it returns the Nodes, Weights, Mean and Variances in each of the 10 points.

I want to replicate this by hand.

I have no issue with the zeros, since I just have to find the ones of (½)*x*(5*x²-3), but I don't know how to calculate the weights. I am following the book Numerical Methods for Stochastic Computations A Spectral Method Approach by Dongbin Xiu (page 40), where it suggests that I should write w_j^(N)=Integrate l_j^(N)wdx.
Does this means that I shoud integrate (½)*x*(5*x²-3) multiplied by the random variable U(-1, 1)? Am I understanding this wrong? I don't know how to integrate a random variable...

 

[Mark44]

Hello everyone. I have a Python code which calculates, given a continuos uniform random variable U(-1,1), the order of a interpolation polynomial and a set of points the evolution of a function of this random variable. i.e.
v0 = cp.Uniform(-1,1)
t = np.linspace(0, 10, 10)

The line above creates an array with 10 elements, spaced along the interval [0, 10]. If you want the elements to be 0, 1, 2, ..., 10, do this:

t = np.linspace(0, 10, 11)

A partition of an interval with 10 subintervals is defined by 11 points.

order=1
.
.
.
plt.plot (t, v0*t)

and it returns the Nodes, Weights, Mean and Variances in each of the 10 points.

I don't see how this is done, based on the code you're showing.

I want to replicate this by hand.

I have no issue with the zeros, since I just have to find the ones of (½)*x*(5*x²-3), but I don't know how to calculate the weights. I am following the book Numerical Methods for Stochastic Computations A Spectral Method Approach by Dongbin Xiu (page 40), where it suggests that I should write w_j^(N)=Integrate l_j^(N)wdx.

What is $l_j^{(N)}$?
Also, is w your weighting function -- $w(x) = \frac 1 2 x(5x^2 - 3)$?

Does this means that I shoud integrate (½)*x*(5*x²-3) multiplied by the random variable U(-1, 1)?

U(-1, 1) isn't the random variable. This notation means that the random variable X is uniformly distributed in the interval [-1, 1].

Am I understanding this wrong? I don't know how to integrate a random variable...

BTW, I used BBCode code tags for the Python code I showed above, and I used LaTeX scripting for the equations.
We have tutorials on both of these:
Code tags - https://www.physicsforums.com/threads/read-this-before-posting-code-geshi-syntax-highlighter.773407/
LaTeX primer - https://www.physicsforums.com/help/latexhelp/

 

[confused_engineer]

I apologize if I haven't expressed myself propplerly.
Here it is the full code:

import chaospy as cp
import numpy as np
import odespy as od
import matplotlib.pyplot as plt
np.set_printoptions(threshold=1000,edgeitems=1000)
def s(t, v0):
return v0*t
v0 = cp.Uniform(-1,1)
#a = cp.Beta(2,2)

#v0 = cp.Normal(0,1)
#a = cp.Normal(0,1)
dist = v0
#t = np.linspace(0, 10, 1000)
t = np.linspace(0, 10, 10)
print(t)
order=1
P, norms = cp.orth_ttr(order, dist, retall=True)
nodes, weights = cp.generate_quadrature(order+1, dist, rule="G")
samples_u = [s(t, *n) for n in nodes.T]
u_hat, coef = cp.fit_quadrature(P, nodes, weights, samples_u, retall=True)
print(nodes.shape)
print(weights.shape)
print(len(samples_u))
sensitivity_Q = cp.Sens_t(u_hat, dist)
#plt.plot(t, cp.E(u_hat, dist))
#plt.plot(t, cp.Var(u_hat, dist))
E = cp.E(u_hat, dist)
Var = cp.Var(u_hat, dist)
Analit_Var = np.array((1./240.)*t**2*(20+3*t**2))
print('Polynomials')
print(P)
print('Samples')
print(samples_u)
print('Norms')
print(norms)
print('Nodes')

print(nodes)
print('Weights')
print(weights)

print('Mean')
print(E)
print('Variance')
print(Var)
print('u_hat')
print(u_hat)
print(coef)
#print(sensitivity_Q)
#plt.plot(t,E)
#plt.plot(t,Var)
#plt.plot (t, [1, 2, 4, 8, 16, 32, 64, 128, 129, 130]) #esto funciona
plt.plot (t, Analit_E, 'bs')
plt.plot(t, Analit_Var, 'g^')
plt.plot(t, Analit_E, t, Analit_Var)
plt.show()

I wanted to replicate what this program returns. The proper base for a uniform variable are the legrendre polynomials. If I need three zeros, the legendre polynomial which has three zeros is (1/2)*x*(5*x³-3), and I find the same zeros in this polynomial as the ones the code returns.

Now, I want to calculate the weights, and for that I want to use the expression found in the book which I mentioned, where l_j^(N) is the corresponding Legendre polynomial. I will rewrite the expression: Integrate between -1 and 1: l_j^(N)*wdx.
This weight calculation is what I don't know how to do.

Edit: the Python anwser is
Weights
[ 0.27777778 0.44444444 0.27777778]

 

[eys_physics]

I don't think you have understood the text you refer to correctly. Firstly, $l_j^N$ is the Lagrange interpolation polynomial, and not the Legendre polynomial. Now, what is discussed in the text is the Gauss-quadrature formula of the form
$$\int_{-1}^{1} W(x) f(x) dx =\sum_{k=1}^n w_k f(x_k) $$
where $W(x)$ is the weight function and $x_k$ are the zeroes of the corresponding orthogonal polynomial, and $w_k$. In your case one has $W(x)=1$ meaning that the polynomial is the Legendre polynomial of order $n$. The easiest way to compute the weights is to use the formula (see http://mathworld.wolfram.com/Legendre-GaussQuadrature.html)
$$w_k=\frac{2}{nP_{n-1}(x_k)P'_n(x_k)}$$

Your last post also contains a typo. The polynomial should be
$$P_3(x)=(1/2)x(5x^2-3)$$

I hope this helps.

 

[confused_engineer]

I don't think you have understood the text you refer to correctly. Firstly, $l_j^N$ is the Lagrange
I hope this helps.

Indeed it does. Thank you.
I only have one question left. Since I am working with a random variable, what shoud I use as f(x)? if I used the expected value, I would get a zero, since I am dealing with a uniform random variable.

 

[eys_physics]

I have to admit that I'm very far from an expert on this. But, I believe you can use a probability distribution of the form
$$f(x)=\begin{cases} 1/(b-a) & a\leqslant x \leqslant b \\
0 & x<a \text{ or } x>b\end{cases}$$

Probably, it is also useful for you to read the Wikipedia page https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)

 

[confused_engineer]

I have to admit that I'm very far from an expert on this. But, I believe you can use a probability distribution of the form
$$f(x)=\begin{cases} 1/(b-a) & a\leqslant x \leqslant b \\
0 & x<a \text{ or } x>b\end{cases}$$

Probably, it is also useful for you to read the Wikipedia page https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)

Don't worry. You have already helped me a lot.
Thanks