在使用根优化之后，结果不是应该达到的结果。

Question

我一直在使用“来自scipy.optimize导入根”的根函数来解决其他需要两个方程的问题，f(x，y)和g(x，y)，到目前为止我还没有发现任何障碍，整个主题是潜在的流动，这个特别的问题是一个涡旋+一个稳定的表面速度，下一个代码是关于求一个点P，(Xp，YP)的坐标，其中速度是零，在一个表面上有一个涡旋(涡旋强度= -550)，这个漩涡在墙的左边。稳定速度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
>>>

我曾经得到过类似的东西，但我记不太清楚，我认为在另一种情况下，是因为手工推导函数的错误，而在目前的问题中，我没有在这方面跟踪任何错误。

Answer 1

您使用的是method='lm'，根据文档，它只在最小二乘意义上解决了方程。使用method="hybr"，您可以得到sol.success == False。

很可能，您的jacobian是不正确的，因为根是在jac=False中找到的。

编辑：您的Jacobian似乎是错误的，至少：

x = np.array([3, 3.])
dx = np.array([1.3, 0.3])
eps = 1e-5
dx = 1e-5 * dx / np.linalg.norm(dx)

df_num = (np.array(fun_imp1(x + dx/2)[0]) - np.array(fun_imp1(x - dx/2)[0])) / eps
df_cmp = fun_imp1(x)[1].dot(dx)/eps

print df_num
print df_cmp

版画

[-1.43834392 -0.69055079]
[   -1.43834392 -1024.60208799]

总是用数值微分来检查雅可比是非常有用的。