python - 迭代三个变量
问题描述
好的,这会很长。我正在尝试计算极地网格中解析解中的水位。它取决于 r 和 theta 以及这个变量 j。我要做的基本上是根据给定的方程计算网格中特定 r、theta 点的水位。这个方程有一部分对 j 的无限值求和。方程式的一部分在这里方程式图像
部分代码如下:
"""
Plot idealized solution for water levels in a quarter circle domain with
constant bathymetry
"""
import numpy as np
import matplotlib.pyplot as plt
#from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
# establish parameters
Ho = 300 #m
g = 9.81 #m/s2
r1 = 1000 #m
r2 = 10000 #m
rr = np.arange(r1,r2,10)
#radial size of domain
phi = np.pi/2
#theta is the angle in radians at a specific location within the domain
#theta = np.pi/4
theta = np.arange(0, phi, np.pi/360)#varies
Theta = theta[1:180]
zeta = [0] * len(rr) * len(Theta)
#converting from wind speed to wind shear stress
U = 10
Cd = (1/1000) * ((3/4) + (U/15))
Roair = 1.225 #kg/m3
Rowater = 997 #kg/m3
W = (Roair/Rowater) * Cd * (U**2)
#wind shear in m^2/s^2 in the 0 direction (W to E)
Wo = np.sqrt((W**2)/2)
#wind shear in m^2/s^2 in the phi direction
Wphi = np.sqrt((W**2)/2)
zeta = np.zeros((len(rr), len(Theta)))
#determines the bathymetry
a_star = []
n = 0
kappa = (1-n)**(0.5)
for t in range(len(Theta)):
a_star.append ( ( (np.sin(phi)) / (g*Ho*kappa* np.sin(kappa*Theta[t])) ) )
#first half of equation 19 that does not depend on j
for r in range(len(rr)):
for t in range(len(Theta)):
zeta[r,t] = ( (a_star[t] * (rr[r]**(1-n)))*(Wo*np.cos(((1-n)**(0.5))*Theta[t]) + Wphi*np.cos(((1-n)**(0.5))*(Theta[t]-phi))) )
#second half of equation 19 for j=0
ajbj = []
for t in range(len(Theta)):
j = 0
Djo = np.sin(( ( (1-n)**(0.5) ) * phi ) ) / ( (1-n)**(0.5) * (phi) )
Ejo = (np.sin(phi)) / (phi)
ajbj.append ( (r2**(1-n)) * (-a_star[t] * Djo))
for r in range(len(rr)):
zeta[r,t] = zeta[r,t] + (ajbj[t])*(Wo+Wphi)
#second half of equation 19 for j=1,2,3 (summation)
sj = []
tj = []
Dj = []
Ej = []
r1EogH = []
astarD = []
tjr1r2 = []
sjr2 = []
aj = []
bj = []
jj = [1,2,3]
for j in range(len(jj)):
sj.append(- (n/2) + np.sqrt( ( (n/2)**2) + ( (jj[j]*np.pi / phi)**2) ) )
tj.append (- (n/2) - np.sqrt( ( (n/2)**2) + ( (jj[j]*np.pi / phi)**2) ) )
Dj.append ( (2* ((-1)**jj[j]) * ((1-n)**(0.5)) * phi * np.sin( ((1-n)**(0.5)) * phi )) / ( (1-n) * (phi**2) - (jj[j]**2) * (np.pi**2) ) )
Ej.append ( (2* ((-1)**jj[j]) * phi * np.sin(phi) ) / ( (phi**2) - (jj[j]**2) * (np.pi**2) ) )
r1EogH.append ( ( (r1**(1-n)) * Ej[j] ) / ( g * Ho ) )
tjr1r2.append ( tj[j] * (r1**tj[j]) * (r2**sj[j]) )
sjr2.append ( sj[j] * (r2**tj[j]) )
for t in range(len(Theta)):
#astarD.append ( a_star[t] * Dj[j] )
aj.append ( ( a_star[t]*Dj[j] * ( ( tj[j] * (r1**tj[j]) * (r2**(1-n)) ) - ( (r2**tj[j]) * (r1**(1-n)) ) ) + ( r1EogH[j] * r2**tj[j] ) ) / ( (sjr2[j] * (r1**sj[j])) - tjr1r2[j] ) )
bj.append ( ( -1* a_star[t]*Dj[j] * ( ( sj[j] * (r1**sj[j]) * (r2**(1-n)) ) - ( (r2**sj[j]) * (r1**(1-n)) ) ) - ( r1EogH[j] * r2**sj[j] ) ) / ( (sjr2[j] * (r2**sj[j])) - tjr1r2[j] ) )
for r in range(len(rr)):
zeta[r,t] = zeta[r,t] + ( ( (aj[j] * rr[r]**(sj[j])) + (bj[j] * rr[r]**(tj[j])) ) * (Wo*np.cos( (jj[j]*np.pi*Theta[t])/phi) + Wphi*np.cos( (jj[j]*np.pi*(Theta[t]-phi))/phi)) )
x,y = np.meshgrid(Theta, rr)
X = Theta
Y = rr
fig = plt.figure()
ax = fig.add_subplot(111, polar='True')
ax.pcolormesh(X, Y, zeta) #X,Y & data2D must all be same dimensions
ax.set_thetamin(0)
ax.set_thetamax(90)
plot = ax.pcolor(zeta)
fig.colorbar(plot)
plt.show()
上述解决方案 我得到的 zeta 值非常大。它们应该更类似于以下 theta PI/4 解决方案中的值。我知道方程式很复杂。我想要的是网格中的一个点说 r=1000, theta= pi/4 计算 j=1, j=2 和 j=3 的 zeta 值,然后将它们加在一起,然后做同样的事情网格中每个点的东西。我想知道我是否只需要以不同的方式构造我的循环?还是不使用 .append 功能?有没有人有什么建议?
为 THETA PI/4 完成
"""
Plot idealized solution for water levels in a quarter circle domain with
constant bathymetry
"""
import numpy as np
import matplotlib.pyplot as plt
#from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
# establish parameters
Ho = 300 #m
g = 9.81 #m/s2
r1 = 1000 #m
r2 = 10000 #m
rr = np.arange(r1,r2,10)
zeta = [0] * len(rr)
#radial size of domain
phi = np.pi/2
#theta is the angle in radians at a specific location within the domain
theta = np.pi/4
#converting from wind speed to wind shear stress
U = 10
Cd = (1/1000) * ((3/4) + (U/15))
Roair = 1.225 #kg/m3
Rowater = 997 #kg/m3
W = (Roair/Rowater) * Cd * (U**2)
#wind shear in m^2/s^2 in the 0 direction (W to E)
Wo = np.sqrt((W**2)/2)
#wind shear in m^2/s^2 in the phi direction
Wphi = np.sqrt((W**2)/2)
#determines the bathymetry
n = 0
kappa = (1-n)**(0.5)
a_star = ( (np.sin(phi)) / (g*Ho*kappa* np.sin(kappa*theta)) )
#first half of equation 19 that does not depend on j
for r in range(len(rr)):
zeta[r] = ( a_star * (rr[r]**(1-n))*(Wo*np.cos(((1-n)**(0.5))*theta) + Wphi*np.cos(((1-n)**(0.5))*(theta-phi))) )
#plt.xlabel("rr")
#plt.ylabel("zeta")
###
#plt.plot(rr,zeta, label = 'LHS eq 19')
#second half of equation 19 for j=0
for r in range(len(rr)):
j = 0
Djo = np.sin(( ( (1-n)**(0.5) ) * phi ) ) / ( (1-n)**(0.5) * (phi) )
Ejo = (np.sin(phi)) / (phi)
ajbj = (r2**(1-n)) * (-a_star * Djo)
zeta[r] = zeta[r] + (ajbj)*(Wo+Wphi)
#plt.xlabel("rr")
#plt.ylabel("zeta")
###
#plt.plot(rr,zeta, label = 'j=0')
#second half of equation 19 for j=1,2,3 (summation)
sj = []
tj = []
Dj = []
Ej = []
r1EogH = []
astarD = []
tjr1r2 = []
sjr2 = []
aj = []
bj = []
jj = [1,2,3]
for j in range(len(jj)):
sj.append(- (n/2) + np.sqrt( ( (n/2)**2) + ( (jj[j]*np.pi / phi)**2) ) )
tj.append (- (n/2) - np.sqrt( ( (n/2)**2) + ( (jj[j]*np.pi / phi)**2) ) )
Dj.append ( (2* ((-1)**jj[j]) * ((1-n)**(0.5)) * phi * np.sin( ((1-n)**(0.5)) * phi )) / ( (1-n) * (phi**2) - (jj[j]**2) * (np.pi**2) ) )
Ej.append ( (2* ((-1)**jj[j]) * phi * np.sin(phi) ) / ( (phi**2) - (jj[j]**2) * (np.pi**2) ) )
r1EogH.append ( ( (r1**(1-n)) * Ej[j] ) / ( g * Ho ) )
astarD.append ( a_star * Dj[j] )
tjr1r2.append ( tj[j] * (r1**tj[j]) * (r2**sj[j]) )
sjr2.append ( sj[j] * (r2**tj[j]) )
aj.append ( ( astarD[j] * ( ( tj[j] * (r1**tj[j]) * (r2**(1-n)) ) - ( (r2**tj[j]) * (r1**(1-n)) ) ) + ( r1EogH[j] * r2**tj[j] ) ) / ( sjr2[j] * (r1**sj[j]) - tjr1r2[j] ) )
bj.append (- ( astarD[j] * ( ( sj[j] * (r1**sj[j]) * (r2**(1-n)) ) - ( (r2**sj[j]) * (r1**(1-n)) ) ) - ( r1EogH[j] * r2**sj[j] ) ) / ( sjr2[j] * (r2**sj[j]) - tjr1r2[j] ) )
for r in range(len(rr)):
zeta[r] = zeta[r] + ( ( (aj[j] * rr[r]**(sj[j])) + (bj[j] * rr[r]**(tj[j])) ) * (Wo*np.cos( (jj[j]*np.pi*theta)/phi) + Wphi*np.cos( (jj[j]*np.pi*(theta-phi))/phi)) )
plt.xlabel("rr")
plt.ylabel("zeta")
plt.title("Wind In at 45 degrees")
#
plt.plot(rr,zeta, label = 'Ho=100m')
plt.legend(loc='upper right')
解决方案
不幸的是,我不是数学天才。在我尝试以某种方式重构它之前,我尝试运行您的代码以查看它得到了什么输出。但是,由于提供的代码段中未定义sj
from之类的内容,因此存在错误。sj.append()
我不知道如何阅读您图片中提供的方程式(数学课程是多年前的,自大学以来就没有使用过)。
至于建议,在没有完全理解你的问题的情况下,我能给你的最好的就是考虑制作函数。这使您/其他人可以更清楚地看到代码。另一个优点是当某些东西发生变化时(似乎不太可能用一个方程,但也许你想在未来的某个时候改变你的输入类型),这些变化是在更小的部分进行的。
def calculate_j(point_or_other_necessary_input):
""" do the maths for j """
return some_equation_or_value_for_j
def calculate_t(point_or_other_necessary_input):
""" do the maths for t """
return some_equation_or_value_for_t
def calculate_r(point_or_other_necessary_input):
""" do the maths for r """
return some_equation_or_value_for_r
def process_grid(zeta_grid):
""" process each point in the provided grid"""
for point in zeta_grid:
# or whatever makes sense to combine the values/equations
calculated_sum = calculate_j(point) + calculate_t(point) + calculate_r(point)
# store value in some reasonable way, probably another grid
return calculated_points
print(process_grid(zeta))
经过进一步评估以及您编辑的更多信息后,我不得不问您是否真的打算这样做:
jj = [1,2,3]
for j in range(len(jj)):
我在猜测,但并不确定,你想j
成为1
那个时候2
,然后3
,不是吗?目前,您确实得到了0
, 1
,2
因为您正在迭代列表的长度,而不是列表中的项目。这在您获取 jj 的值的第一部分并不重要,jj[j]
但它在您使用的 r 部分中,j
而不是在以下部分jj[j]
中:
zeta[r] = zeta[r] + ( ( (aj[j] * rr[r]**(sj[j])) + (bj[j] * rr[r]**(tj[j])) )
如前所述,这可能是您的预期行为,即仅使用 0,1,2 而不是 1,2,3。至于为什么值可能太大,您是否打算为所有 rr 计算 r 3 次?您将它嵌套在第一个 for 循环中,这意味着它在循环的每次迭代中都在执行for j in ...
。所以它会zeta[r] = zeta[r] + ...
大量地添加到自身 ( ) 中。(同样,不是数学家)但不是运行 len(rr)^3 次或 3^len(rr) 次吗?比可能预期的多一些倍。
推荐阅读
- pandas - 如何在 Dask DataFrame 上正确执行多索引切片?
- jmeter - Jmeter - 同一线程组中不同http请求的计数变化
- kotlin - 如何调用具有接收器的类的函数?
- javascript - 如果索引为0,Vue.js如何添加活动类
- macos - Android Studio - 安装 APK 时出错 + 管道损坏错误(在 Mac 上)
- scala - 如何使用 Spark-Sql 高效读取 Hive 表
- angularjs - 量角器运行时错误失败:脚本超时:量角器:进程退出,错误代码为 1
- swift - 我如何使用 CGAffineTransform 进行正确的拖放
- azure - 从源 xxxx 访问 XMLHttpRequest ' 被 CORS 策略阻止:对预检的响应。预检请求不允许重定向
- c - 如何按城镇值的排序顺序将指向城市数组的指针插入城镇数组?