Box-Mueller算法出现对数函数负/零参数运行时错误的解决咨询

Hey there, let's break down why you're hitting that runtime error and how to fix it quickly.

The Root Cause #

Your error stems from a logical mistake in the do while loop condition meant to filter valid points for the Box-Mueller algorithm. Let’s look at that problematic line:

Do while (w>=1.0.and.w<0.0)

This condition is impossible to satisfy—w can’t be both ≥1 and <0 at the same time! The loop never runs, so you end up using invalid values of w (either your initial w=2 or occasionally w=0 if r1 and r2 both land on 0) when calculating log(w). When w=0, log(w) throws the "negative or zero argument" error; when w>1, -2*log(w) becomes negative, which would cause another error when taking the square root later.

The Fix #

The Box-Mueller algorithm requires generating points strictly inside the unit circle, meaning w = r1² + r2² must be 0 < w < 1. Update your loop condition to keep generating new r1 and r2 until w falls into this valid range:

Do while (w >= 1.0 .or. w <= 0.0)

This way, if w is too large (≥1) or zero (the only invalid non-positive value possible here), the loop will re-run to get valid values.

Corrected Code Snippet #

Here’s the fixed section of your program:

Do j=0,n
    ! Reset w to an invalid value to ensure we enter the loop
    w = 2.0
    Do while (w >= 1.0 .or. w <= 0.0)
        Call random_number(r1)
        Call random_number(r2)
        r1 = 2.0*r1 - 1.0
        r2 = 2.0*r2 - 1.0
        w = r1*r1 + r2*r2
    End do
    w = sqrt((-2.0*log(w))/w)
    r1 = r1*w
    r2 = r2*w
    If(mod(j,2)==0) then
        w=r1
    Else if(mod(j,2)==1) then
        w=r2
    End if
    fx(j+1)=fx(j)+w*sqrt(2.0*D*dt)
    ! Rest of your code...
End do

Bonus Optimization (Optional) #

Right now, you generate a new pair of Gaussian numbers every j iteration but only use one each time. You can reuse the second number to cut down on random_number calls:

! Track if we have an unused Gaussian number from the last generation
logical :: has_spare = .false.
real :: spare_gauss

Do j=0,n
    if (.not. has_spare) then
        w = 2.0
        Do while (w >= 1.0 .or. w <= 0.0)
            Call random_number(r1)
            Call random_number(r2)
            r1 = 2.0*r1 - 1.0
            r2 = 2.0*r2 - 1.0
            w = r1*r1 + r2*r2
        End do
        w = sqrt((-2.0*log(w))/w)
        r1 = r1*w
        r2 = r2*w
        spare_gauss = r2
        has_spare = .true.
        w = r1
    else
        w = spare_gauss
        has_spare = .false.
    end if
    ! Rest of your code for updating fx(j+1)...
End do

This is just a performance tweak—your core error is fixed with the loop condition change.

Why This Works #

By fixing the loop condition, you ensure w is always in the (0,1) range before calculating log(w). This eliminates invalid arguments to the logarithm routine and guarantees the square root in the Box-Mueller formula uses a positive value.

内容的提问来源于stack exchange，提问作者SchrodingersCat