Random Number generator does not refresh for each iteration of the loop due to operation is done on it

I want the loop to run the rand function until while condition is met. The moment you have multiple operations on the rand function the code just breaks and the loop just doesnt work properly. Took me awhile to pinpoint the problem. When I change it from C*(1-cos(rand*pi)) to B*(1-cos(rand*pi)), the loop works but doesnt display the right results.

the only reason I could come up with would be, B is a number of its own while C is B with something else which somehows break the rand function to properly run in the while loop until conditions are met.

Part shown here is the problem:

while B > A

C=B/(m0*c2);

B = B / (1+(C*(1-cos(rand*pi))));

s = s + 1;

end

Complete program:

clc
clear
clearvars

%Constants%
%=========%

h=6.6261e-34;
c=2.9979e8;
c2=c^2;
m0= 9.10956e-31;


%Variables%
%=========%

A=15.000e3;
n=100000;
D=224.060e3;
Energydist= normrnd(D,10000,[1,n]);


Sc=zeros(1,n);
Ef=zeros(1,n);


%===========================================================%


for i=1:n
    
    B=Energydist(i);
    s=0;
    
    
    while B > A 
        
        C=B/(m0*c2);
        B = B / (1+(C*(1-cos(rand*pi))));
        
        s = s + 1;
        
        
    end
    
    
    Sc(i)= s;
    Ef(i)= B;
    
end

MaxScattered=(Energydist/A);

Scattered = sum(Sc)
Scattering_efficiency = Scattered / n

tiledlayout(2,2)

ax1=nexttile;
histfit(Sc,50,'poisson')
xlabel('Number of Maximum Scatters Achieved')
ylabel('Frequency')

ax2=nexttile;
histfit(MaxScattered,50,'poisson')
xlabel('Maximum Number of scatter expected')
ylabel('Frequency')
ax3=nexttile;
plot(Ef,Sc)
xlabel('Final Energy of electron (eV)')
ylabel('Number of scatter')

ax4=nexttile;
plot(Energydist,MaxScattered)
xlabel('Initial Energy of electron (eV)')
ylabel('Maximum Number of scatter expected')

[Voss]

What is this code supposed to be doing, and what are the right results?

I don't think the problem is the usage of rand. Also, whether I use C*(1-cos(rand*pi)) or B*(1-cos(rand*pi)), the loop runs, though it may just iterate once for each i.

Answer 1

[Walter Roberson]

D=224.060e3;
Energydist= normrnd(D,10000,[1,n]);

Those random values are on the order of 224000. Those become the initial values of B.

You calculate

        C=B/(m0*c2);
        B = B / (1+(C*(1-cos(rand*pi))));

and you want to compare the result to A.

For the moment, call cos(rand*pi) by the name R. cos(rand*pi) (here named R) will be a value between -1 and +1 (theoretically exclusive on both sides because rand can never exactly equal 0 or 1)

So in order to loop again for the moment we have

        C=B/(m0*c2);
        A == B / (1+(C*(1-R)))

which is the condition

        -B/((79228162514264337593543950336*B*(R - 1))/6486512639245165 - 1) == 15000

or

        -B/(12214292474343.21*B*(R - 1) - 1) == 15000

when you restrict R (the value of the cos(pi*rand)) to the range 0 to 1, this has a solution only if B is in the range of -8E-11 to -4E-14 . For all positive B, the new calculated value for B is far far less than A, so the loop stops.

Please note that your C coefficient will almost certainly be so large that (C*(1-cos(rand*pi)) is about 1e+18, at which point the difference between representable floating point numbers is in the range of 512 or at least 256 -- so the 1+ is probably not going to make any difference to the calculation.

No "mistake" or malfunction of rand() is needed: your equations just work out that after 1 round, the new B value is going to be so much smaller than A that the loop will stop.

[IMAN SYAHMI BIN JASNI]

I made the plots and histogram of which left column are simulated result while on the right column is the theoretical average.

MaxScattered=(Energydist/A);

Scattered = sum(Sc)
Scattering_efficiency = Scattered / n

tiledlayout(2,2)

ax1=nexttile;
histfit(Sc,50,'poisson')
xlabel('Number of Maximum Scatters Achieved')
ylabel('Frequency')

ax2=nexttile;
histfit(MaxScattered,50,'poisson')
xlabel('Maximum Number of scatter expected')
ylabel('Frequency')
ax3=nexttile;
plot(Ef,Sc)
xlabel('Final Energy of electron (eV)')
ylabel('Number of scatter')

ax4=nexttile;
plot(Energydist,MaxScattered)
xlabel('Initial Energy of electron (eV)')
ylabel('Maximum Number of scatter expected')

[Walter Roberson]

Okay, but I do not know how that relates to the previous discussion?

[IMAN SYAHMI BIN JASNI]

It shows the issue with rand not iterate in the while loop but only just due to for loop.

[Walter Roberson]

That proves nothing of the sort. Look below, where i have it print out the random number for the first 15 iterations: the random number is clearly not stuck in any way.

I showed you mathematically why you cannot get more than one scatter -- because with the constants you are using and the initial values you are using, in order to get more scatters, your B would have to be between -8E-11 to -4E-14

format long g

%Constants%
%=========%

h=6.6261e-34;
c=2.9979e8;
c2=c^2;
m0= 9.10956e-31;


%Variables%
%=========%

A=15.000e3;
n=100000;
D=224.060e3;
Energydist= normrnd(D,10000,[1,n]);


Sc=zeros(1,n);
Ef=zeros(1,n);


%===========================================================%


for i=1:n
    
    B=Energydist(i);
    s=0;
    
    
    while B > A 
        
        C=B/(m0*c2);
        R = rand;
        B = B / (1+(C*(1-cos(R*pi))));
        if i <= 15 || s > 0
          fprintf('i = %d, scatter #%d, rand = %.15f\n', i, s, R);
        end
        s = s + 1;
        
        
    end
    
    
    Sc(i)= s;
    Ef(i)= B;
    
end

i = 1, scatter #0, rand = 0.219466075263635
i = 2, scatter #0, rand = 0.911841110212500
i = 3, scatter #0, rand = 0.307685677560821
i = 4, scatter #0, rand = 0.647390154628003
i = 5, scatter #0, rand = 0.700127070480684
i = 6, scatter #0, rand = 0.601458379302843
i = 7, scatter #0, rand = 0.678973873394930
i = 8, scatter #0, rand = 0.099393719644295
i = 9, scatter #0, rand = 0.567386988308393
i = 10, scatter #0, rand = 0.574242047152035
i = 11, scatter #0, rand = 0.240506493462864
i = 12, scatter #0, rand = 0.390816309225101
i = 13, scatter #0, rand = 0.347826238255798
i = 14, scatter #0, rand = 0.332174225017237
i = 15, scatter #0, rand = 0.092818489982529

MaxScattered=(Energydist/A);

Scattered = sum(Sc)

Scattered = 
      100000

Scattering_efficiency = Scattered / n

Scattering_efficiency = 
     1

tiledlayout(2,2)

ax1=nexttile;
histfit(Sc,50,'poisson')
xlabel('Number of Maximum Scatters Achieved')
ylabel('Frequency')

ax2=nexttile;
histfit(MaxScattered,50,'poisson')
xlabel('Maximum Number of scatter expected')
ylabel('Frequency')
ax3=nexttile;
plot(Ef,Sc)
xlabel('Final Energy of electron (eV)')
ylabel('Number of scatter')

ax4=nexttile;
plot(Energydist,MaxScattered)
xlabel('Initial Energy of electron (eV)')
ylabel('Maximum Number of scatter expected')

[IMAN SYAHMI BIN JASNI]

This is the result you should expect from the simulation vs theoretical average.

[Walter Roberson]

Okay, but that does not mean that there is a problem with rand()

• you might have a problem with your constants
• your formulas for the calculation might be wrong