Using physical parameter as a Gaussian random variable in a simple Poisson problem

Question

I want to vary the input parameter of a physical dynamic mechanics problem, as a Gaussian Random variable and view the resulting Probability Density Function (PDF). I used the Finite Element Method to create deterministic codes in C++ and Python language. However, before I go there, I want to begin by applying this technique to a simple Poisson problem deterministic FEM code written in python or C++. I have mentioned my thoughts below. If you have done this before, please, help me through your own approach.

In a Poisson problem, I need to vary $k$, i.e., the spatial variability parameter, as an RV with some statistical distribution. Let's assume a mean of 1 and a std. deviation of 0.3, or any other convenient value. Now, I want to run this code multiple times for different values of the parameter - $k$. So, let's say, I generate an array of 100 values using randn() function of Matlab/python. Now, maybe a shell script or python script is needed to use this array with a for-loop and run the FEM code 100 times for each value of $k$. How will the deterministic FEM code be mentioned in the loop? How will it accept those values of $k$? Does anyone have a generic script that can help?

I would also really appreciate some suggestions/ comments on post-processing. I believe the output .vtk or .pvd format files can be used to create PDF in Paraview software. I am also a little confused about the generation of PDF. I believe I'll have to choose a single point on the 2D mesh and for each point, I have output values in $x$ and $y$ direction, since it is a 2D problem. Therefore, I'll have to choose 1 direction too. So, the PDF of output on 1 point, in 1 direction, is what I am looking for. Is that correct? I am not sure.

I basically want to do a non-intrusive polynomial chaos expansion using quadrature or sampling. However, as stated above, I am starting with Monte Carlo, taking a step at a time, so I understand how to do pre and post-processing.

Answer 1

In the loop over your instances of 100 random values for $k$, you could write the value for $k$ into a file, and call your finite element code; the finite element code could then open the file, read the value for $k$ from it and then use it for whatever computation you may want to do with it.

As for your second question, visualizing uncertainty is difficult. Why don't you just output the solution at a bunch of points for each value of $k$, and then create histograms from this for individual points by hand? This should be easy enough within matlab.