How to implement a Monte Carlo Simulation for studying uncertainty in dynamic stiffness?

Question

I have a dynamic mechanical system of the form: $$\mathbf{M{\ddot q}+{\Omega_c}G{\dot q}+{Kq}=F(t)}$$

and I need to introduce some uncertainty for the $K$ parameter, which is a symetric matrix of $q$ degrees of freedom.

I did a bit of research of Monte Carlo simulation but I´m a bit confused of how different textbooks approach the problem.

For example , my actual stiffness value is of $1e8$, and I wish to generate values uniformly distributed between $K_l = 1e6$ and $K_u = 1e12$.

I found that I can do this by applying: $$\lambda_i = \mu_{\lambda_i} + \epsilon_i$$ where $\mu_{\lambda_i}$ is the mean value of $\lambda_i$ and $\epsilon$ is a random parameter, and also $$\sigma_{\lambda_i}^2 = E[{(\lambda_i-\mu_{\lambda_i})^2}]$$

I would be using Matlab for running the simulation , with the command rand I can generate a randomly distributed numbers between 0 and 1.

The problem I find with this is that I don't know know to input the boundaries of the $K$ stiffnes values I want.

Then I found another text book that explains how to generate random variables using: $$ X_i = A + (B-A)R_i$$ Where $R_i$ is the random number which can be generated by Matlab with rand command to generate a new random variable $X_i$, with lower and upper limits $A$ and $B$. So in matlab I would do a for loop like this:

for i=1:num
    Ki=Kl+(Ku-Kl)*rand;
end

My main question would actually be, what is the difference between these two methods and for the first one, would $\epsilon$ just be a random 0-1 number adding to the average designated value?

Thanks

Answer 1

In generating uniform pseudo random numbers the method I have always used is the second method, where you have,

$$R_i \sim U(0,1)$$

And if you want a uniformly distributed number $X_i$ between A, B. The formula is,

$$X_i = A+(B-A)R_i$$

It's straight forward to derive this formula by looking at the CDF for a uniform distribution between A,B.

First consider the PDF,

\begin{cases} \frac{1}{B-A} &, A\leq x\leq B\\ 0 &, Otherwise \end{cases}

The resulting CDF is,

$$\int_A^x\frac{1}{B-A}dx$$ \begin{cases} 0&,x<A\\ \frac{x-A}{B-A}&,A\leq x<B\\ 1 &, x\geq B \end{cases}

equating the CDF to $R$ we get,

$$R = \frac{x-A}{B-A}$$ $$x = R_i(B-A)+A $$

Because we equated the CDF to $R$, $R_i \sim U(0,1)$

The mean and variance of this distribution is,

$$\mu = \frac{1}{2}(A+B)$$ $$\sigma^2 = \frac{1}{12}(B-A)^2$$

Coming back to your question,

Yes you could apply $\lambda = \mu + \epsilon$ where $\mu$ is $1e8$ but now determining epsilon is complicated because you want a uniform distribution in both negative and positive direction to achieve $1e6<\lambda<1e12$. Therefore, $\epsilon$ will not be a random number between 0-1. The code you have written looks correct. A good test is to generate n samples and calculate the mean and variance. Compare this value with the theoretical values obtained from the equations above.

It is standard practice however to use a uniform random distribution $R\sim U(0,1)$. Using this one uniform random number you can generate any distribution you want. If you can determine the CDF, equate it and take the inverse as shown above. There are more "complicated" techniques to generate things like a normal distribution from a uniform distribution between 0,1 but it is possible.