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

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

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$$

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.

• Thanks a lot, it makes much more sense now. Could you explain with a bit more detail how you derived the formula? I don´t understand what you equated and I wasn´t really familiar with CDF , what´s the difference from the probability density function?
– Roy
Commented Aug 5, 2017 at 17:43
• The method I used above is called the inverse transform method. en.wikipedia.org/wiki/Inverse_transform_sampling. While a PDF is essentially the probability that P(X=a) where X is a random variable. The CDF is the probability that P(X<a). Hence the CDF is the area under the PDF upto point a, therefore we integrate the PDF to get the CDF. Commented Aug 5, 2017 at 17:49
• You dont really need to know how to generate your own distributions using the inverse transform method. Most common distributions are built into most of the common programming languages. You can also think of the uniform distribution as being shifted. Given R from [0,1), we shift it by A units. Then we need to define a width of the distribution which is (B-A). Hence when R=1, X=B and R=0, X=A, and all the numbers in between A and B for other values of R. Commented Aug 5, 2017 at 17:59
• Thanks again @Sada93. If I have a most likely value of $K$ instead of lower and upper limits, say $1e8$, would it be better to use the first method and also if I want to have a normal distribution ?
– Roy
Commented Aug 6, 2017 at 23:46
• If you have an idea for what the mean ($\mu$) and spread (standard deviation $\sigma$) is you can find the limits, A and B using the equations above. The choice between uniform and normal is a fundamental difference in how you want to model your system. Choose uniform if every possibility between A and B is equally likely. Choose normal if the mean value is most likely. The first method you described can be modeled as both uniform and normal by the choice of $\epsilon$, I would still recommend using the inbuilt functions. Please upvote and accept the answer if you found it usefull. Commented Aug 7, 2017 at 4:34