- #1
trmcclain
- 14
- 0
I'm hoping this will be the last time I call for help, but in any case, here it goes.
I thought I had a handle on this before, but in all of my attempts, my code diverges within a few iterations.
My problem is creating a spatial distribution of particles given a probability density. I've tried to start out with an easier function, the normal distribution [itex]f(x)=e^{-x^2}[/itex]
The method I've tried to implement, which has been failing, is the Newton-Raphson method:
[itex]x_{new}=x_{guess}-\frac{U-G(x)}{f(x)}[/itex] where U is a random number between 0 and 1 and [itex]G(x)=\int_{-∞}^x f(x)\,dx[/itex]
My code is in C++ and is as follows:
dr = rand() % 100;
U=dr/100.;
do{
Gfunc=(1./R)*(sqrt(PI)/2.)*(erf(rg)+1);
func=exp(-(rg*rg));
whole=(U-Gfunc)/func;
rn=rg-whole;
if(rn > 1E-5)
rg=rn;
}while(rn>1E-5);
(Sorry if there's a way to make code pretty on the forums)
**Note, I do use srand(time(NULL)) above this loop.**
I thought I had a handle on this before, but in all of my attempts, my code diverges within a few iterations.
My problem is creating a spatial distribution of particles given a probability density. I've tried to start out with an easier function, the normal distribution [itex]f(x)=e^{-x^2}[/itex]
The method I've tried to implement, which has been failing, is the Newton-Raphson method:
[itex]x_{new}=x_{guess}-\frac{U-G(x)}{f(x)}[/itex] where U is a random number between 0 and 1 and [itex]G(x)=\int_{-∞}^x f(x)\,dx[/itex]
My code is in C++ and is as follows:
dr = rand() % 100;
U=dr/100.;
do{
Gfunc=(1./R)*(sqrt(PI)/2.)*(erf(rg)+1);
func=exp(-(rg*rg));
whole=(U-Gfunc)/func;
rn=rg-whole;
if(rn > 1E-5)
rg=rn;
}while(rn>1E-5);
(Sorry if there's a way to make code pretty on the forums)
**Note, I do use srand(time(NULL)) above this loop.**