科学的python数组语法

Question

所以我用 python 编写了这段代码。我不打算解释它，因为它是一个简单的语法修复，我似乎看不到，所以解释它的用途是没有用的。

我遇到的问题是，对于给定的 d，例如 15，我得到正确的“cuentas”值和正确的“e”值。

我想要做的是遍历一组 d 并获取每个 cuentas 和每个 e 的值，以便绘制 d 与 e。

我的问题是我似乎不知道如何在 python 中创建矩阵。

在matlab中我曾经写过两个不同的循环，比如theese：

for i=1:1:N
    for j=1:9

      a[i,j]= and so on

a[i,j] 将是一个包含 N 行和 9 列的矩阵，我可以访问和操作。

在我下面的代码中，我会故意将 # cmets 放在我想要遍历距离的地方

import numpy as np
import matplotlib.pyplot as plt
N=100000
cos=np.zeros(N)
phi=np.zeros(N)
teta=np.zeros(N)
a=np.zeros(N)


xrn=np.zeros(N)
yrn=np.zeros(N)
zrn=np.zeros(N)

x=np.zeros(N)
y=np.zeros(N)
z=np.zeros(N)

lim1=14.7
lim2=3.35
lim3=-lim1
lim4=-lim2

#d=np.array([15,20,25,30,35,40,45,50,55])
d=15

    #for j in range(9):
for i in range(N):
    cos[i]=np.random.uniform(-1,1)
    teta[i]=np.random.uniform(-np.pi,np.pi)
    phi[i]=np.random.uniform(0,2*np.pi)

# a[i]=d[j]/cos[i]*np.cos(phi[i])
a[i]=d/cos[i]*np.cos(phi[i])


xrn[i]=np.random.uniform(-1,1)*lim1
yrn[i]=np.random.uniform(-1,1)*lim2

x[i]=a[i]*np.sin(teta[i])*np.cos(phi[i])+xrn[i]
y[i]=a[i]*np.sin(teta[i])*np.sin(phi[i])+yrn[i]

#cuentas[j]=0
cuentas=0

#for j in range(9):
for i in range(N):
    if a[i]>0 and x[i] < lim1 and x[i]>lim3 and y[i] < lim2 and y[i]>lim4:
        cuentas=cuentas+1
#e[j]=cuenta[j]/N
e=cuentas/N

非常感谢那些阅读！

Answer 1

短版：

这与 Python 中的 MATLAB 代码完全相同

a = np.zeros([N, 9])
for i in range(N):
    for j in range(9):
        a[i,j]= and so on

唯一的主要区别是您需要事先定义数组，如果您希望您的代码具有合理的性能，您也应该在 MATLAB 中这样做。

但是，如果事先不知道大小，可以在 Python 中使用列表，然后在最后转换为 numpy 数组。由于处理列表和矩阵/数组的内部机制，这将比大型数组的 MATLAB 示例快得多：

a = []
for i in range(N):
    a.append([])
    for j in range(9):
        a[-1].append( and so on
a = np.array(a)

[-1] 表示（“a 的最后一个元素”，append 将括号内的任何内容放在列表末尾。所以。a[-1].append(foo) 表示“将foo 放入任何内容位于a 的最后一个元素中。

长版：

您的 MATLAB 代码在 Python 中的工作方式大致相同，但您需要考虑一些显着差异。

首先，分配给大于现有数组/矩阵的索引在 MATLAB 中有效，但在 numpy 中无效。所以如果你有一个大小[5, 5]数组/矩阵，在MATLAB中你可以分配给元素[5, 6]，但你不能在numpy中。这意味着在 MATLAB 中您可以从一个空数组开始，而在 numpy 中您必须事先设置数组大小。请注意，MATLAB 矩阵实际上无法调整大小，实际上每次通过循环创建一个新矩阵并将所有数据复制到其中。这非常慢，这就是 MATLAB 警告您预分配数组的原因。 Numpy 只是不假装能够调整数组的大小，因此您需要更明确地进行复制、预分配或使用列表（可调整大小）。

其次，类似地，MATLAB 不需要您在使用矩阵之前定义它，而 numpy 则需要。这是因为传统上 MATLAB 具有三种数据结构（矩阵、元胞数组和结构），每种都有自己的索引风格。因此，MATLAB 可以仅根据您对它的索引方式来确定您想要创建哪种数据结构。 Python只有一种索引风格，所以它不能做出这种猜测。

第三，在 MATLAB 中使用具有某些（但不是全部）功能的单个数组大小创建一个具有该大小的每个维度的 2D 方阵，而在 numpy 中它创建一个 1D 数组。我不确定您的代码是否是您所期望的。坦率地说，我不知道为什么 MATLAB 会这样工作。

第四，numpy 数组可以有任意数量的维度，0（标量）、1（向量）、2、3、4 等。另一方面，MATLAB 矩阵必须至少有两个维度。这可能会导致一些意想不到的差异，例如转置对 numpy 向量什么都不做。

至于您的 Python 代码，如果不说出了什么问题，我无法告诉您如何解决它。但希望我已经给了你足够的信息让你自己去做。

Answer 2

您可以通过以下方式使用 numpy 在 python 中创建矩阵：

n=5
k=4
a=np.zeros([n,k])
for i in range(n):
    for j in range(k):
        a[i,j]=i+j
print(a)

结果是

[[ 0.  1.  2.  3.]
 [ 1.  2.  3.  4.]
 [ 2.  3.  4.  5.]
 [ 3.  4.  5.  6.]
 [ 4.  5.  6.  7.]]

Answer 3

所以我采纳了你的答案并且成功了！

如果有人想知道有多少粒子会通过两个长方形探测器，这是一个蒙特卡罗模拟。当我从检测器 1 中抛出粒子时，它们通过它是微不足道的，我计算通过检测器 2 的数量。

修正后的代码是

N=100000
"La dirección viene dada por v=[rsin(teta)*cos(phi),rsin(teta)sin(phi),rcos(teta)]"
"Los vectores que vamos a usar debemos inicializarlos como un vector de ceros"
cos=np.zeros([10,N])
phi=np.zeros([10,N])
teta=np.zeros([10,N])
a=np.zeros([10,N])

xrn=np.zeros(N)
yrn=np.zeros(N)
zrn=np.zeros(N)

x=np.zeros([10,N])
y=np.zeros([10,N])
z=np.zeros([10,N])

lim1=14.7
lim2=3.35
lim3=-lim1
lim4=-lim2
"d son las disversas distancias a las que colocamos la fuente con respecto al detector"

d=np.array([0.00001,15,20,25,30,35,40,45,50,55])

"e es la eficiencia geométrica simulada"
e=np.zeros(10)

"Debemos definir el coseno como números aleatorios en vez de el ángulo teta, debido a que queremos"
"que se distribuyan uniformemente por toda la esfera"
for j in range(10):
    for i in range(N):
        cos[j,i]=np.random.uniform(0,1)

        phi[j,i]=np.random.uniform(0,2*np.pi)

        a[j,i]=d[j]/cos[j,i]

        xrn[i]=np.random.uniform(-1,1)*lim1
        yrn[i]=np.random.uniform(-1,1)*lim2

        x[j,i]=a[j,i]*np.sin(math.acos(cos[j,i]))*np.cos(phi[j,i])+xrn[i]
        y[j,i]=a[j,i]*np.sin(math.acos(cos[j,i]))*np.sin(phi[j,i])+yrn[i]


cuentas=np.zeros(10)
for j in range(10):
    for i in range(N):
        if a[j,i]>0 and x[j,i] < lim1 and x[j,i]>lim3 and y[j,i] < lim2 and y[j,i]>lim4:
            cuentas[j]=cuentas[j]+1

    e[j]=cuentas[j]/N

谢谢大家！