Generating noise in particular way

Hi there,

I would like to generate noise in a certain way. Assume I have my FEM-space V_h.
Then my noise is to be generated via

\sum_{k=1}^{N_h} \xi_k e_k

where e_k is a basis for the FEM space and \xi_i are normally distributed variables. The issue is to obtain a basis for the FEM space. This method has to be generalizable, so I wonder if there is some “generate basis”-funtion for the FEM-space? There seems to be something like it,
https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/cpp/fem/BasisFunction.html
but I don’t understand the documenation, and more importantly, how do I get this to work with fenics? I need some sort of “integrable” expression, right? Furthermore, I want the covariance matrix R of the noise to be given by

(R^T R)^{-1}=M^{-1}

where M is the mass matrix of the FEM space Vh.

Perhaps this works for you:

from fenics import *
import numpy as np

# mesh, FunctionSpace
mesh = UnitIntervalMesh(5)
Vh = FunctionSpace(mesh, 'P', 1)

# allocate function u in Vh
u = Function(Vh)
Nh = u.vector()[:].size
xi = np.random.rand(Nh, ) # or any other random array function

# set the entries on the nodes of u to (the random array) xi
u.vector()[:] = xi

# generating mass matrix
z = TestFunction(Vh)
y = TrialFunction(Vh)
m = y*z*dx
Mass = assemble(m)
# Mass.array() gives you the desired matrix
M = Mass.array()
L = np.linalg.cholesky(M)
R = np.transpose(L)

I do not know if the last two lines gives you the covariance matrix R. I supposed that R should be upper triangular? The Cholesky-decomposition gives a lower triangular matrix L such that L L^T = M. Define R := L^T, which is upper triangular, then we have R^T R = L^{T^T} L^T = L L^T = M.

Thank you for your reply, I appreciate it a lot!

I am sorry, I am new to fenics (not helped by the complete lack of documentation), but I don’t see how that random generation would generate according to my formula above… it is quite important that it takes a sum of normal random coefficients times an eigenbasis function, that is to say

\sum_{k=1}^{N_h} \xi_k e_k

where then e_k is the eigenbasis.
Furhtermore, I want to use this as a RHS in an elliptic equation. I can’t get this to a form that fenics accepts or understand.

Ok. So let u = \sum_{k = 0}^{N_h} u_k e_k, where e_k is your basis of V_h. With u = Function(V) I get exactly such function in FEniCS. Here u.vector()[0] corresponds to u_0, u.vector()[5] to u_5 and so on. With u.vector()[:] I get an array U = [u_0, u_1, \ldots, u_{N_h}]. What I want to do now is to set u_k = \xi_k for k = 0, \ldots, N_h. Thats why I am generating a random-array xi. Now I have to assign the nodal values of u with the of \xi. Translating this is u.vector()[k] = xi[k] for all k. Or, in short notation u.vector()[:] = xi. Here, I assign an array with an array.
Hope that helps. Otherwise I have not understood your question.

Okey! Great! I think I get it now! Just one question, is u=Function(V) any basis element or is it an eigenbasis?

None of both. V is your functionspace of choice (P1, P2, DG, …), u = Function(V) is a Function from the functionspace V, which is 0 everywhere.

So what is u.vector()[0] then?

The value of u corresponding to the 0-th basisfunction of the functionspace V.