2
我一直在使用來自「從scipy.optimize導入根」的其他問題的解決方案中需要兩個方程，f（x ，y）和g（x，y），到目前爲止我還沒有發現任何障礙，直到現在，整個話題都是潛在的流動，而這個特殊的問題是關於一個表面上的渦旋+穩定速度， 接下來的代碼是關於找到點P的座標（Xp，YP），其中速度爲零，並且在表面上具有渦旋（渦旋的強度= -550），並且該渦旋位於壁的左側。 U：穩定速度 CV：渦流強度 H：只是速度的部件中的一個是零與程序發現的結果的渦流和表面使用根優化後的結果不是他們應該是什麼
import numpy as np
from scipy.optimize import root
from math import pi
cv = -550.0
U = 10.0
h = 18.0
'''
denom1 = (X + h) ** 2 + Y ** 2
denom2 = (X - h) ** 2 + Y ** 2
###########################################
# f(x,y)
###########################################
f_1_a1 = - cv * Y/denom1
f_1_a2 = cv * Y/denom2
# f(x, y)
f_1 = f_1_a1 + f_1_a2
dfx_1 = (- cv * Y) * ((-1) * (2) * (X + h))/(denom1 ** 2)
dfx_2 = (cv * Y) * ((-1) * (2) * (X - h))/(denom2 ** 2)
# df_x
d_f_1_x = dfx_1 + dfx_2
dfy_1 = (- cv)/denom1
dfy_2 = (- cv * Y) * (- 1) * (2) * (Y) /(denom1 ** 2)
dfy_3 = (cv)/denom2
dfy_4 = (cv * Y) * (-1) * (2 * Y) /(denom2 ** 2)
# df_y
d_f_1_y = dfy_1 + dfy_2 + dfy_3 + dfy_4
###########################################
# g(x,y)
###########################################
g_a1 = - U
g_a2 = cv * (X + h)/denom1
g_a3 = - cv * (X - h)/denom2
# g(x, y)
f_2 = g_a1 + g_a2 + g_a3
dgx_1 = cv/denom1
dgx_2 = cv * (X + h) * (-1) * (2) * (X + h)/(denom1 ** 2)
dgx_3 = (- cv)/denom2
dgx_4 = (- cv) * (X - h) * (-1) * (2) * (X - h)/denom2
dgx = dgx_1 + dgx_2 + dgx_3 + dgx_3 + dgx_4
# dg_x
d_f_2_x = dgx
dgy_1 = cv * (X + h) * (-1) * (2) * (Y)/(denom1 ** 2)
dgy_2 = (- cv) * (X - h) * (-1) * (2 * Y)/(denom2 ** 2)
dgy = dgy_1 + dgy_2
# dg_y
d_f_2_y = dgy
'''
def Proof(X, Y):
denom1 = (X + h) ** 2 + Y ** 2
denom2 = (X - h) ** 2 + Y ** 2
###########################################
# f(x,y)
###########################################
f_1_a1 = - cv * Y/denom1
f_1_a2 = cv * Y/denom2
# f(x, y)
f_1 = f_1_a1 + f_1_a2
dfx_1 = (- cv * Y) * ((-1) * (2) * (X + h))/(denom1 ** 2)
dfx_2 = (cv * Y) * ((-1) * (2) * (X - h))/(denom2 ** 2)
# df_x
d_f_1_x = dfx_1 + dfx_2
dfy_1 = (- cv)/denom1
dfy_2 = (- cv * Y) * (- 1) * (2) * (Y) /(denom1 ** 2)
dfy_3 = (cv)/denom2
dfy_4 = (cv * Y) * (-1) * (2 * Y) /(denom2 ** 2)
# df_y
d_f_1_y = dfy_1 + dfy_2 + dfy_3 + dfy_4
###########################################
# g(x,y)
###########################################
g_a1 = - U
g_a2 = cv * (X + h)/denom1
g_a3 = - cv * (X - h)/denom2
# g(x, y)
f_2 = g_a1 + g_a2 + g_a3
dgx_1 = cv/denom1
dgx_2 = cv * (X + h) * (-1) * (2) * (X + h)/(denom1 ** 2)
dgx_3 = (- cv)/denom2
dgx_4 = (- cv) * (X - h) * (-1) * (2) * (X - h)/denom2
dgx = dgx_1 + dgx_2 + dgx_3 + dgx_3 + dgx_4
# dg_x
d_f_2_x = dgx
dgy_1 = cv * (X + h) * (-1) * (2) * (Y)/(denom1 ** 2)
dgy_2 = (- cv) * (X - h) * (-1) * (2 * Y)/(denom2 ** 2)
dgy = dgy_1 + dgy_2
# dg_y
d_f_2_y = dgy
print "The values of u and v are:"
print f_1
print f_2
print "The derivates are:"
print dgx, dgy
print d_f_1_x, d_f_1_y
def fun_imp1(x):
X = x[0]
Y = x[1]
denom1 = (X + h) ** 2 + Y ** 2
denom2 = (X - h) ** 2 + Y ** 2
###########################################
# f(x,y)
###########################################
f_1_a1 = - cv * Y/denom1
f_1_a2 = cv * Y/denom2
# f(x, y)
f_1 = f_1_a1 + f_1_a2
dfx_1 = (- cv * Y) * ((-1) * (2) * (X + h))/(denom1 ** 2)
dfx_2 = (cv * Y) * ((-1) * (2) * (X - h))/(denom2 ** 2)
# df_x
d_f_1_x = dfx_1 + dfx_2
dfy_1 = (- cv)/denom1
dfy_2 = (- cv * Y) * (- 1) * (2) * (Y) /(denom1 ** 2)
dfy_3 = (cv)/denom2
dfy_4 = (cv * Y) * (-1) * (2 * Y) /(denom2 ** 2)
# df_y
d_f_1_y = dfy_1 + dfy_2 + dfy_3 + dfy_4
###########################################
# g(x,y)
###########################################
g_a1 = - U
g_a2 = cv * (X + h)/denom1
g_a3 = - cv * (X - h)/denom2
# g(x, y)
f_2 = g_a1 + g_a2 + g_a3
dgx_1 = cv/denom1
dgx_2 = cv * (X + h) * (-1) * (2) * (X + h)/(denom1 ** 2)
dgx_3 = (- cv)/denom2
dgx_4 = (- cv) * (X - h) * (-1) * (2) * (X - h)/denom2
dgx = dgx_1 + dgx_2 + dgx_3 + dgx_3 + dgx_4
# dg_x
d_f_2_x = dgx
dgy_1 = cv * (X + h) * (-1) * (2) * (Y)/(denom1 ** 2)
dgy_2 = (- cv) * (X - h) * (-1) * (2 * Y)/(denom2 ** 2)
dgy = dgy_1 + dgy_2
# dg_y
d_f_2_y = dgy
a_1 = f_1
a_2 = f_2
b_1 = d_f_1_x
b_2 = d_f_1_y
c_1 = d_f_2_x
c_2 = d_f_2_y
f = [ a_1,
a_2]
df = np.array([[b_1, b_2],
[c_1, c_2]])
return f, df
sol = root(fun_imp1, [ 1, 1], jac = True, method = 'lm')
print "x = ", sol.x
print "x0 =", sol.x[1]
print "y0 =", sol.x[0]
x_1 = sol.x[0]
x_2 = sol.x[1]
Proof(x_1, x_2)
之間的距離。 起初我以爲這是衍生品的問題，但我沒有發現任何問題。我的一位朋友曾經說過，當渦旋強度過高時（如超過150），渦旋強度有時會有不同的表現。
新增信息：
這是流線的情節：
使用此代碼後：
import numpy as np
import matplotlib.pyplot as plt
vortex_height = 18.0
h = vortex_height
vortex_intensity = -550.0
cv = vortex_intensity
permanent_speed = 10
U1 = permanent_speed
Y, X = np.mgrid[-21:21:100j, -21:21:100j]
U = (- cv * Y)/((X + h)**2 + (Y ** 2)) + (cv * Y)/((X - h)**2 + (Y ** 2))
V = - U1 + (cv * (X + h))/((X + h)**2 + (Y ** 2)) - (cv * (X - h))/((X - h)**2 + (Y ** 2))
speed = np.sqrt(U*U + V*V)
plt.streamplot(X, Y, U, V, color=U, linewidth=2, cmap=plt.cm.autumn)
plt.colorbar()
plt.savefig("stream_plot.png")
plt.show()
而且我用得到的結果該計劃是：
>>>
x = [ 1.32580109e-01 3.98170636e+02]
x0 = 398.170635755
y0 = 0.132580109151
The values of u and v are:
-8.2830922107e-05
-10.1246349802
The derivates are:
-2.20709329055 0.000624761030349
-0.000624761030349 6.22388943399e-07
>>>
其中u和v應該是：
U = 0.0
V = 0.0
代替：
U = -8.2830922107e-05（這一個是可以接受的） v = -10.1246349802（這是絕對錯誤的）
而當我將其更改爲'hybr'時
sol = root(fun_imp1, [ 1, 1], jac = True, method = 'hybr')
我得到這個：
>>>
C:\Python27\lib\site-packages\scipy\optimize\minpack.py:221: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)
x = [ -4.81817071e+02 1.96057929e+06]
x0 = 1960579.2949
y0 = -481.817070593
The values of u and v are:
2.53176901102e-12
-10.0000000052
The derivates are:
-7.14899730857e-05 5.25462578799e-15
-5.25462578799e-15 -3.87401132188e-18
>>>
我有類似的東西一次，但我不記得很清楚，我認爲在其他情況下是由於受到手功能不好的推導，以及在目前的問題中，我沒有跟蹤這方面的任何錯誤。
'sol.success'的值是什麼？你是否得到了同樣的解決方案而不通過雅可比行列式？什麼是結果而不是它應該是什麼？你是否繪製了流場圖，看看結果是否有意義？請用這個信息編輯問題 – gg349
'hybr'求解器應該有'sol.success == False';你應該檢查結果。 'lm'求解器是不同的，因爲它解決了最小二乘的問題，即使它不是真正的根也能成功返回。 –