【发布时间】:2025-12-29 05:45:07
【问题描述】:
我正在尝试从我正在运行的数值模拟中获取一些数据。
我有一个函数可以解决我的 ODE 并将我需要的值保存在一个向量中,rv。我想运行我的函数 1000 次,并为每个时间步获得 rv 的平均值,但我不知道该怎么做。
我的第一个想法是调用函数1000 次并将数据保存在一个数组中。
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)
但是当我这样做时,我收到以下错误消息:
NameError: name 'rv' is not defined
我将附上我的代码,因为我不太擅长解释。抱歉弄得这么乱
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)
【问题讨论】:
-
你在循环之前定义了 rv 吗?
-
第三行的
:没有问题? -
您的变量
rv未在给定代码的某处声明。仅存在r、v和R。您可能希望使用 this 作为datavector = np.append(datavector, mean(r, v)) -
是的,来自 circle(N,dt,x0,y0,phi0,r,T,nu,v,Omega,R):我得到了 rv(我试过打印出来)。但我无法将其中的 1000 个收集到一个数组中。
-
我已经在函数圈中定义了 rv 但没有在 for i in range(1000):-loop 中定义。我希望我的函数给我 rv
标签: python