python - 为什么我的 python 代码的运行时间这么长,我该怎么做才能让它运行得更快?
问题描述
我刚刚为我的论文写了一些 python 代码。由于某种原因,计算最终值需要花费大量时间。经过仔细检查,它与代码中的 Sk 数组有关,但对我来说,这似乎不是一个特别昂贵的计算。有人可以解释一下我在编码问题方面做错了什么,我可以做些什么来获得相同的数字结果,同时减少我的运行时间?我附上了下面的代码。
#Finding Optimal Values For Helix
import numpy as np
import math
import matplotlib.pyplot as plt
import sympy as sym
from sympy.abc import L
#Initalization of variables
omega_list = []
Rd = 1 #Inner Radius
td = 1 #Thickness Radius
h = 1 #height
l = 1 # Length
lamb = 2*sym.pi*h #Wavelength
mu = 1 #Don't touch Viscosity
q = 0.09*lamb #Don't touch
tau = sym.Matrix([[0],[0],[1]])
x = np.linspace(0.001, 10, num = 5 )
#For Loop For Optimization
for h in x:
# print(i)
#Distance to center of Mass
r1 = Rd*sym.cos(sym.abc.s)-(6*h**2*Rd)/(3*Rd**2+4*sym.pi**2*h**2)
r2 = Rd*sym.sin(sym.abc.s)+(6*h**2*Rd)/(3*Rd**2+4*sym.pi**2*h**2)
r3 = h*sym.abc.s-(3*sym.pi*h**2*Rd**2+6*sym.pi*h**2)/(3*Rd**2+4*sym.pi**2*h**2)
print('again')
#Frenet Frame Components
t = np.array([[-Rd*sym.sin(sym.abc.s)],
[Rd*sym.cos(sym.abc.s)],
[h]])/sym.sqrt(Rd**2+h**2)
n = np.array([[-Rd*sym.cos(sym.abc.s)],
[-Rd*sym.sin(sym.abc.s)],
[0]])
b = np.array([[sym.sin(sym.abc.s)],
[-sym.cos(sym.abc.s)],
[Rd/h]])/sym.sqrt(Rd**2+h**2)*h
#Skew Matrix
Sk = np.array([[0, -r3, r2],
[r3, 0, -r1],
[-r2, r1, 0]])
#Frenet Frame Assembly
R = np.hstack((t,n,b))
#Drag Matrix
C = np.array([[l*(2*sym.pi*mu)/sym.log(2*q/td), 0, 0],
[0,l*(4*sym.pi*mu)/sym.log(2*q/(td+0.5)) , 0],
[0, 0, l*(4*sym.pi*mu)/sym.log(2*q/td)]])
#Ressistance Matrix Assembly
R11 = -1*R@C@R.transpose()
R12 = R@C@R.transpose()@Sk
R21 = -1*Sk@R@C@R.transpose()
R22 = Sk@R@C@R.transpose()@Sk
R11 = sym.Matrix([R11[0,:], R11[1,:], R11[2,:]])
R12 = sym.Matrix([R12[0,:], R12[1,:], R12[2,:]])
R21 = sym.Matrix([R21[0,:], R21[1,:], R21[2,:]])
R22 = sym.Matrix([R22[0,:], R22[1,:], R22[2,:]])
#Integration of Ressistance
R11 = R11.integrate((sym.abc.s, 0, 2*sym.pi))
R12 = R12.integrate((sym.abc.s, 0, 2*sym.pi))
R21 = R21.integrate((sym.abc.s, 0, 2*sym.pi))
R22 = R22.integrate((sym.abc.s, 0, 2*sym.pi))
print('wtf^2')
#Move Around to Isolate the Angular Velocity Omega
ROmega = -1*R21@Sk+R22
print('wtf')
Omega = ROmega.inv()@tau
print('good')
Omega_final = (Omega[0,0]**2+Omega[1,0]**2+Omega[2,0]**2)
print('all good')
omega_list.append(Omega_final)
print('ik')
#Find the Velocity
Sk@Omega
#Plot the Graphs
print("Omega:", omega_list)
print(x)
#plt.plot([x], [omega_final])
for i in omega_list:
plt.plot(x, omega_list)
plt.show()
解决方案
我不确定这是否是故意的,但您的代码基本上计算了一些矩阵乘法作为 s 的函数。R21 和 R22 不是 s 的函数。而 omega 是 s 与 R21 和 R22 的函数的乘积。我假设您正在尝试将 omega 绘制为 s 的函数。
正如@hpaulj 所说,如果您对速度计算不感兴趣,则可能没有必要使用 SymPy。这是一个纯数字版本。我的回答很糟糕,但很快。我相信有人可以轻松地出现并使其更具可读性。
# Finding Optimal Values For Helix
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
# Initalization of variables
omega_list = []
Rd = 1 # Inner Radius
td = 1 # Thickness Radius
h = 1 # height
l = 1 # Length
lamb = 2 * np.pi * h # Wavelength
mu = 1 # Don't touch Viscosity
q = 0.09 * lamb # Don't touch
tau = np.array([[0], [0], [1]])
x = np.linspace(0, 10, num=5)
x[0] += 1/1000
# For Loop For Optimization
print(x)
for h in x:
print(f"Beginning loop where h={h}")
# Distance to center of Mass
r1 = lambda s: Rd * np.cos(s) - (6 * h ** 2 * Rd) / (3 * Rd ** 2 + 4 * np.pi ** 2 * h ** 2)
r2 = lambda s: Rd * np.sin(s) + (6 * h ** 2 * Rd) / (3 * Rd ** 2 + 4 * np.pi ** 2 * h ** 2)
r3 = lambda s: h * s - (3 * np.pi * h ** 2 * Rd ** 2 + 6 * np.pi * h ** 2) / (
3 * Rd ** 2 + 4 * np.pi ** 2 * h ** 2)
# Frenet Frame Components
t = lambda s: np.array([[-Rd * np.sin(s)],
[Rd * np.cos(s)],
[h]]) / np.sqrt(Rd ** 2 + h ** 2)
n = lambda s: np.array([[-Rd * np.cos(s)],
[-Rd * np.sin(s)],
[0]])
b = lambda s: np.array([[np.sin(s)],
[-np.cos(s)],
[Rd / h]]) / np.sqrt(Rd ** 2 + h ** 2) * h
# Skew Matrix
Sk = lambda s: np.array([[0, -r3(s), r2(s)],
[r3(s), 0, -r1(s)],
[-r2(s), r1(s), 0]])
# Frenet Frame Assembly
R = lambda s: np.hstack((t(s), n(s), b(s)))
# Drag Matrix
C = lambda s: np.array([[l * (2 * np.pi * mu) / np.log(2 * q / td), 0, 0],
[0, l * (4 * np.pi * mu) / np.log(2 * q / (td + 1/2)), 0],
[0, 0, l * (4 * np.pi * mu) / np.log(2 * q / td)]])
# Ressistance Matrix Assembly
R21 = lambda s: -1 * Sk(s) @ R(s) @ C(s) @ R(s).transpose()
R22 = lambda s: Sk(s) @ R(s) @ C(s) @ R(s).transpose() @ Sk(s)
# Integration of Resistance
R21 = np.array([[quad(lambda s: R21(s)[i, j], 0, 2*np.pi)[0] for j in range(3)] for i in range(3)])
R22 = np.array([[quad(lambda s: R22(s)[i, j], 0, 2*np.pi)[0] for j in range(3)] for i in range(3)])
print("Finished integration")
# Move Around to Isolate the Angular Velocity Omega
ROmega = lambda s: -1 * R21 @ Sk(s) + R22
Omega = lambda s: np.linalg.inv(ROmega(s).astype(float)) @ tau
Omega_final = lambda s: (Omega(s)[0, 0] ** 2 + Omega(s)[1, 0] ** 2 + Omega(s)[2, 0] ** 2)
# Cannot find the Velocity as a function of s
# Adding lambda functions to the array does not work
lam = Omega_final
x_vals = np.linspace(0, 2 * np.pi, 100)
omega_list.append(np.vectorize(lam)(x_vals))
# Plot the Graphs
for (i, omega) in enumerate(omega_list[0:1]):
lam = omega
x_vals = np.linspace(0, 2*np.pi, 100)
plt.plot(x_vals, omega, label=f"h={x[i]}")
plt.legend()
plt.show()
这将在约 1 秒内运行。
如果你想要一个符号版本,它是可能的,但它更长。这些表达式非常大,但我们可以稍微加快速度:如果我们扩展线性组合,矩阵乘法后的积分可以更快。这就是为什么expand(...).evalf()
使用了这么多次。simplify()
需要太长时间。
# Finding Optimal Values For Helix
import matplotlib.pyplot as plt
import sympy as sym
import numpy as np
# Initalization of variables
omega_list = []
Rd = 1 # Inner Radius
td = 1 # Thickness Radius
h = 1 # height
l = 1 # Length
lamb = 2 * sym.pi * h # Wavelength
mu = 1 # Don't touch Viscosity
q = sym.Rational("0.09") * lamb # Don't touch
tau = sym.Matrix([[0], [0], [1]])
x = list(range(0, 10, 2))
x[0] += sym.S(1)/1000
s = sym.symbols("s", real=True)
# For Loop For Optimization
# h took to long to invert the matrix if it was a symbol
for h in x:
print(f"Beginning loop where h={h}")
# Distance to center of Mass
r1 = Rd * sym.cos(s) - (6 * h ** 2 * Rd) / (3 * Rd ** 2 + 4 * sym.pi ** 2 * h ** 2)
r2 = Rd * sym.sin(s) + (6 * h ** 2 * Rd) / (3 * Rd ** 2 + 4 * sym.pi ** 2 * h ** 2)
r3 = h * s - (3 * sym.pi * h ** 2 * Rd ** 2 + 6 * sym.pi * h ** 2) / (
3 * Rd ** 2 + 4 * sym.pi ** 2 * h ** 2)
# Frenet Frame Components
t = sym.Matrix([[-Rd * sym.sin(s)],
[Rd * sym.cos(s)],
[h]]) / sym.sqrt(Rd ** 2 + h ** 2)
n = sym.Matrix([[-Rd * sym.cos(s)],
[-Rd * sym.sin(s)],
[0]])
b = sym.Matrix([[sym.sin(s)],
[-sym.cos(s)],
[Rd / h]]) / sym.sqrt(Rd ** 2 + h ** 2) * h
# Skew Matrix
Sk = sym.Matrix([[0, -r3, r2],
[r3, 0, -r1],
[-r2, r1, 0]])
Sk = sym.simplify(Sk)
# Frenet Frame Assembly
# R = np.hstack((t, n, b))
R = t.row_join(n).row_join(b)
R = sym.simplify(R)
# Drag Matrix
C = sym.Matrix([[l * (2 * sym.pi * mu) / sym.log(2 * q / td), 0, 0],
[0, l * (4 * sym.pi * mu) / sym.log(2 * q / (td + sym.S(1)/2)), 0],
[0, 0, l * (4 * sym.pi * mu) / sym.log(2 * q / td)]])
C = sym.simplify(C)
# Ressistance Matrix Assembly
R21 = -1 * sym.expand(sym.expand(Sk @ R) @ C) @ R.transpose()
R21 = sym.expand(R21).evalf()
R22 = sym.expand(sym.expand(sym.expand(Sk @ R) @ C) @ R.transpose()) @ Sk
R22 = sym.expand(R22).evalf()
# Integration of Resistance
R21 = sym.integrate(R21, (s, 0, 2 * sym.pi))
R22 = sym.integrate(R22, (s, 0, 2 * sym.pi))
print("Finished integration")
# Move Around to Isolate the Angular Velocity Omega
ROmega = -1 * sym.expand(R21 @ Sk) + R22
ROmega = sym.expand(ROmega).evalf()
print(f"ROmega = {ROmega}")
Omega = ROmega.inv() @ tau
print(f"Omega = {Omega}")
Omega_final = (Omega[0, 0] ** 2 + Omega[1, 0] ** 2 + Omega[2, 0] ** 2).evalf()
print(f"Omega_final = {Omega_final}")
omega_list.append(Omega_final)
print('ik')
# Find the Velocity
print(f"Sk @ Omega = {sym.expand(Sk @ Omega)}")
# Omega list is massive
print(f"omega_list = {omega_list}")
# Plot the Graphs
for (i, omega) in enumerate(omega_list):
lam = sym.lambdify((s,), omega)
s_vals = np.linspace(0, 2*np.pi, 100)
plt.plot(s_vals, lam(s_vals), label=f"h={x[i]}")
plt.legend()
plt.show()
这会产生相同的情节,但在我的机器上可能需要一分钟。
推荐阅读
- javascript - 已修复且始终位于 Chrome 扩展程序中的顶部弹出窗口
- javascript - 只是因为将 HTML 嵌套 div 放在新行上而不断出错
- sql-server - 是否必须在 SQL 后端和 VB.Net 前端建立表关系?
- apache-spark - 将字典传递给 pyspark udf
- c# - 在 postgreSQL 中使用 odbcdatareading 加载数据表
- javascript - 无法从函数内部删除 eventListener
- ios - 从方法向现有 NSDictionary 添加值
- import - 在原型中从父文件夹导入到子文件夹
- python - 按列号熊猫选择不相邻的列
- php - PHP 从 PHP 脚本生成的 url 下载图像