Я новичок в сообществе и прошу прощения, если я не предоставил информацию, как предполагалось.
Я пытаюсь узнать python, пришедший из Matlab. У меня есть очень простой код для начальных целей:
from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np
def plot_hc(radius, pd):
q = np.linspace(0.008, 1.0, num=500)
y = hc(q, radius, pd)
plt.loglog(q, y)
plt.show()
def hc_formfactor(q, radius):
y = (1.0 / q) * (radius * j1(q * radius))
y = y ** 2
return y
def g_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd)) * exp(
-((z - radius) / (sqrt(
2) * pd)) ** 2)
def ln_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
-(log(z / radius) / (sqrt(2) * pd)) ** 2)
# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1
@vectorize
def hc(x, radius, pd):
global d
if Dist == 1:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: g_distribution(z, radius, pd), va, vb)
elif Dist == 2:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
else:
d = 1
def fun(z, x, radius, pd):
if Dist == 1:
return hc_formfactor(x, z) * g_distribution(z, radius, pd)
elif Dist == 2:
return hc_formfactor(x, z) * ln_distribution(z, radius, pd)
else:
return hc_formfactor(x, z)
y = integrate.quad(lambda z: fun(z, x, radius, pd), va, vb)[0]
return y/d[0]
if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)
Как и некоторые, я должен использовать для l oop, но это еще больше снизило скорость. Код, используемый для l oop, выглядит следующим образом:
from numpy import vectorize
from scipy import integrate
from scipy.special import j1
from math import sqrt, exp, pi, log
import matplotlib.pyplot as plt
import numpy as np
def plot_hc(radius, pd):
q = np.linspace(0.008, 1.0, num=500)
y = hc(q, radius, pd)
plt.loglog(q, y)
plt.show()
def hc_formfactor(q, radius):
y = (1.0 / q) * (radius * j1(q * radius))
y = y ** 2
return y
def g_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd)) * exp(
-((z - radius) / (sqrt(
2) * pd)) ** 2)
def ln_distribution(z, radius, pd):
return (1 / (sqrt(2 * pi) * pd * z / radius)) * exp(
-(log(z / radius) / (sqrt(2) * pd)) ** 2)
# Dist=1(for G_Distribution)
# Dist=2(for LN Distribution)
Dist = 1
def hc(q, radius, pd):
if Dist == 1:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: g_distribution(z, radius, pd), va,vb)
elif Dist == 2:
nmpts = 4
va = radius - nmpts * pd
vb = radius + nmpts * pd
if va < 0:
va = 0
d = integrate.quad(lambda z: ln_distribution(z, radius, pd), va, vb)
else:
d = 1
def fun(z, q, radius, pd):
if Dist == 1:
return hc_formfactor(q, z) * g_distribution(z, radius, pd)
elif Dist == 2:
return hc_formfactor(q, z) * ln_distribution(z, radius, pd)
else:
return hc_formfactor(q, z)
y = empty([len(q)])
for n in range(len(q)):
y[n] = integrate.quad(lambda z: fun(z, q[n], radius, pd), va, vb)[0]
return y / d[0]
if __name__ == '__main__':
plot_hc(radius=40, pd=0.5)
Если я запускаю ту же программу для тех же значений в Matlab, это очень быстро. Я не знаю, какую ошибку я сделал здесь. Пожалуйста, помогите:)