Q&A
具体的になにどううまくいかないのかを記載ください。
質問内容
下記作成中のコードにて，極座標(r-θ平面)上にベクトル図を作図しようと考えています．
r,θを変数とし，各座標の値を代入して，その場の速度を算出してベクトル図を描ければと考えましたが上手くいきませんでした．
下記エラーコードが出たため，その対策も調べ確認してみましたが，上手くいかずご質問させていただきました．
アドバイスなど頂けるとありがたいです．
作成コード
import cmath import numpy as np import scipy from scipy.integrate import quad import matplotlib.pyplot as plt #積分1 def M(func, a, b, **kwargs): def real_func(r): return scipy.real(func(r)) def imag_func(r): return scipy.imag(func(r)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) #複素積分計算1 def F_rc(func, a, b, **kwargs): def real_func(r): return scipy.real(func(r)) def imag_func(r): return scipy.imag(func(r)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) #複素積分計算2 def G_rc(func, a, b, **kwargs): def real_func(r): return scipy.real(func(r)) def imag_func(r): return scipy.imag(func(r)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) #複素積分計算3 def F_r2(func, a, b, **kwargs): def real_func(r): return scipy.real(func(r)) def imag_func(r): return scipy.imag(func(r)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) #複素積分計算4 def G_r2(func, a, b, **kwargs): def real_func(r): return scipy.real(func(r)) def imag_func(r): return scipy.imag(func(r)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) #Vr,Vtheta変換 def pol2cart(th,v_th,v_r): return v_r*np.cos(th) - v_th*np.sin(th), v_r*np.sin(th) + v_th*np.cos(th) RadiusColumn = 0.12 NumberRadii = 10 NumberThetas = 20 theta = np.linspace(0,2*np.pi,NumberThetas,endpoint=False) radius = np.linspace(RadiusColumn, 10 * RadiusColumn, NumberRadii)[:,None] f = plt.figure() ax = f.add_subplot(111, polar=True) r0 = 0.14 r1 = 0.11 r2 = 0.096 r3 = 0.0384 vr2 = 1.65 a_a = 0.14 +0.06*1j beta2 = cmath.pi / 9 time = 360 n = 1 Q = 1 ganma = - Q / cmath.tan(beta2) omega = 1.1e-06 +1.6e-05 * 1j A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1) A_F_rc = F_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / radius)**(n + 1), radius, r3) A_G_rc = G_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (radius / r)**(n - 1), r3, radius) A_F_r2 = F_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r2)**(n + 1), r3, r2) A_G_r2 = G_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r2 / r)**(n - 1), r2, r2) vr_tilde = (-1j *((0.5 * A_F_rc[0]) + (0.5 * A_G_rc[0]))) - (r2 / radius)**(n + 1) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))) - a_a *((radius / r2)**(n-1) - (r2 / radius)**(n+1) * cmath.exp(-2 * 1j * beta2)) vtheta_tilde = ((0.5 * A_F_rc[0]) - (0.5 * A_G_rc[0])) - 1j * (r2/radius)**(n + 1) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))) - 1j * a_a * ((radius / r2)**(n - 1) - (r2 / radius)**(n + 1) * cmath.exp(-2 * 1j * beta2)) VelocityRadius = Q/(2* cmath.pi * radius) + vr_tilde * cmath.exp(1j * (omega * time - n * theta)) VelocityTheta = ganma/(2* cmath.pi * radius) + vtheta_tilde * cmath.exp(1j * (omega * time - n * theta)) TotalVelocity = np.linalg.norm([VelocityRadius.real, VelocityTheta.real],axis=0) VelocityX,VelocityY = pol2cart(theta,VelocityTheta.real ,VelocityRadius.real) ax.quiver(theta,radius, VelocityX/TotalVelocity, VelocityY/TotalVelocity) plt.show()
エラーコード
ValueError Traceback (most recent call last) <ipython-input-1-ee1d8f65e98c> in <module> 86 87 A_M = M(lambda r: r2 / (r * cmath.sin(beta2)), r2, r1) ---> 88 A_F_rc = F_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / radius)**(n + 1), radius, r3) 89 A_G_rc = G_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (radius / r)**(n - 1), r3, radius) 90 A_F_r2 = F_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r2)**(n + 1), r3, r2) <ipython-input-1-ee1d8f65e98c> in F_rc(func, a, b, **kwargs) 22 def imag_func(r): 23 return scipy.imag(func(r)) ---> 24 real_integral = quad(real_func, a, b, **kwargs) 25 imag_integral = quad(imag_func, a, b, **kwargs) 26 return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) ~/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst) 336 337 # check the limits of integration: \int_a^b, expect a < b --> 338 flip, a, b = b < a, min(a, b), max(a, b) 339 340 if weight is None: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
出力結果
回答1件
0
自己解決
for文に書き換えたら上手く動きました．
for ra in range(len(radius)): for th in range(len(theta)): A_M = M(lambda r: r2 / (r * (cmath.sin(beta2))**2), r2, r1) A_F_rc = F_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / radius[ra])**(n + 1), r3, radius[ra]) A_G_rc = G_rc(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (radius[ra] / r)**(n - 1), radius[ra], r2) A_F_r2 = F_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r / r2)**(n + 1), r3, r2) A_G_r2 = G_r2(lambda r: cmath.exp((-1j* (omega * cmath.pi * (r**2 - r2**2) / Q - n * ganma * cmath.log(r / r2) / Q ))) * (r2 / r)**(n - 1), r2, r2) vr_tilde = (-1j *((0.5 * A_F_rc[0]) + (0.5 * A_G_rc[0]))) - (r2 / radius[ra])**(n + 1) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))) - a31_a32 *((radius[ra] / r2)**(n-1) - (r2 / radius[ra])**(n+1) * cmath.exp(-2 * 1j * beta2)) vtheta_tilde = ((0.5 * A_F_rc[0]) - (0.5 * A_G_rc[0])) - 1j * (r2/radius[ra])**(n + 1) * (-1j *((0.5 * A_F_r2[0]) + (0.5 * A_G_r2[0]))) - 1j * a31_a32 * ((radius[ra] / r2)**(n - 1) - (r2 / radius[ra])**(n + 1) * cmath.exp(-2 * 1j * beta2)) VelocityRadius =Q/(2* cmath.pi * radius[ra]) + vr_tilde * cmath.exp(1j * (omega * time - n * theta[th])) VelocityTheta = ganma/(2* cmath.pi * radius[ra]) + vtheta_tilde * cmath.exp(1j * (omega * time - n * theta[th])) TotalVelocity = np.linalg.norm([VelocityRadius.real, VelocityTheta.real],axis=0) ax.quiver(theta[th], radius[ra], VelocityRadius.real * cmath.cos(theta[th]) - VelocityTheta.real * cmath.sin(theta[th]), VelocityRadius.real * cmath.sin(theta[th]) + VelocityTheta.real * cmath.cos(theta[th])) print(math.degrees(math.atan2(VelocityRadius.real, VelocityTheta.real))) plt.show()
投稿2020/08/29 14:54
総合スコア2
あなたの回答
tips
プレビュー
