1

I am trying to get some data from a numerical simulation I am running. I have a function that solves my ODE and saves the values I need in a vector, rv. I would like to run my function 1000 times and get the mean value of rv for each timestep but I can't figure out how to do it.

My first thought was to call the function 1000 times and save the data in an array.

for i in range(1000):
    datavector = np.array(0)
    circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R)
    datavector = np.append(datavector, rv)

But when I do that I got the following error message: NameError: name 'rv' is not defined

I'll attach my code since I'm not very good at explaining. Sorry for it being so messy

import math
import numpy as np
import matplotlib.pyplot as plt

def circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R):


kB = 1.38*10**-23 # Boltzmann constant [J/K]plt #math constants
DT = kB*T/(6*math.pi*nu*r)  # Translational diffusion coefficent [m^2/s]
DR = kB*T/(8*math.pi*nu*r**3) # Rotational diffusion coefficent [rad^2/s]
n = 0 #iteration constant


x=x0 #vectors and initial values
y=y0
phi=phi0
phiv= np.array(phi) #vector containing phi-values
xv= np.array(x) #vector containing x-values
yv= np.array(y) #vector containing y-values
rv=np.array(np.sqrt(x**2+y**2))
xss=np.array(x) #vector containing start and (soon) end value
yss=np.array(y) # same as above

phiK=math.sqrt(2*DR*dt) #constants used in the iteration
xyK=math.sqrt(2*DT*dt)
Odt=Omega*dt

while n < N: #from 0 -> N-1
    phi = phi + Odt + phiK*np.random.normal(0,1) #eq (9)
    x = x + v*math.cos(phi)*dt + xyK*np.random.normal(0,1) #eq (10)
    y = y + v*math.sin(phi)*dt + xyK*np.random.normal(0,1) #eq (11)
    if (x**2+y**2) > R**2: #if the particle is outside the boundary
        if abs(x) > R: #need to make sure arccos(x/R) is meaningful
            xn = np.sign(x)*R
            theta = np.sign(y)*np.arccos((xn/R))#angle to particle
        else:
            theta = np.sign(y)*np.arccos((x/R))#angle to particle  
        rp = np.array([x,y]) #r for the particle
        rr = np.array([np.cos(theta)*R,np.sin(theta)*R]) #r for the boundary closest to the particle
        #d = rp - rr #smallest distance from particle to boundary
        q = (2*np.linalg.norm(rr)/np.linalg.norm(rp)) - 1
        x = q*np.array(rp[0]) #x- and y-value forced inside the boundary
        y = q*np.array(rp[1])

    phiv = np.append(phiv, phi) #adding all phi-values to a vector
    xv = np.append(xv, x) # adding all x-values to a vector
    yv = np.append(yv, y) # adding all y-values to a vector'
    rv = np.append(rv,(np.sqrt(x**2+y**2)))
    n=n+1 #iteration 


return(rv)
#print(rv)

for i in range(2): #run the program a number of times
    datavector = np.array(0)
    circle(1E5,1E-3,0*np.random.uniform(-1E-5,1E-5),0*np.random.uniform(-1E-5,1E-5),np.random.uniform(-2*np.pi,2*np.pi),1E-6,300,1E-3,5*1E-6,0,2E-5)
    datavector = np.append(datavector, rv)
np.savetxt('testet.txt', datavector)
xaxablyat
  • 27
  • 1
  • 7
  • Have you defined rv before the loop ? – Sumit S Chawla May 17 '18 at 08:58
  • 1
    And you don't have problems with that `:` on 3rd line? – zipa May 17 '18 at 08:58
  • 1
    Your variable `rv` is not declared somewhere in the given code. There only exists `r`, `v` and `R`. You might want to do the mean using [this](https://stackoverflow.com/questions/7716331/calculating-arithmetic-mean-average-in-python) as `datavector = np.append(datavector, mean(r, v))` – Lord Salforis May 17 '18 at 09:00
  • Yes, from circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R): i get rv (I've tried printing it out). But I am not able to collect 1000 of them into an array. – xaxablyat May 17 '18 at 09:10
  • I have defined rv in the function circle but not in the for i in range(1000):-loop. I want my function to give me rv – xaxablyat May 17 '18 at 09:11
  • The : was a misstake from me. It isn't in the code, I've edited the post now. – xaxablyat May 17 '18 at 09:12
  • My code is pretty messy so I didn't include it but my function gives me the vector rv (it has nothing to do with the argument r or v in the function, I'm just bad at naming stuff). – xaxablyat May 17 '18 at 09:22
  • Gives you, how? Does it return the vector or store it in a global variable? – Arndt Jonasson May 17 '18 at 09:27
  • @xaxablyat: *I want my function to give me rv*. Then use `rv = circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R)` and **return** it from `circle` function. – Serge Ballesta May 17 '18 at 09:43
  • Thanks for all the help – xaxablyat May 17 '18 at 10:22

1 Answers1

1

name 'rv' is not defined means your variable has not been defined before operations have been done upon it.

From the information you provided, rv is never defined. On top of that, you are re-initialising your result vector datavector at every iteration, while it should only be initialised BEFORE your main loop.

I assume rv is the return value from circle.

So a correction of your code should look like:

datavector = []
for i in range(1000):
    rv = circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R)
    datavector.append(rv)
syltruong
  • 2,563
  • 20
  • 33
  • thank you, it was that i hadn't put rv = circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R) that was the problem! I ended up going with datav=np.ndarray(shape=(1000,1E5)) for i in range(1000): #run the program a number of times rv=circle(1E5,1E-3,0*np.random.uniform(-1E-5,1E-5),0*np.random.uniform(-1E-5,1E-5),np.random.uniform(-2*np.pi,2*np.pi),1E-6,300,1E-3,5*1E-6,0,2E-5) datav[i]=(rv) #print(datav) datav=np.transpose(datav) np.savetxt('testet.txt', datav) – xaxablyat May 17 '18 at 10:21