python - Numpy/Scipy求解带积分的联立方程

Question

我正在尝试使用numpy和scipy来解决以下两个方程式：

P(z) = sgn(-cos(np.pi*D1) + cos(5*z)) * sgn(-cos(np.pi*D2) + cos(6*z))

1. 0 = 2/2pi ∫ P(z,D1,D2) * cos(5z) dz + z/L

2. 0 = 2/2pi ∫ P(z,D1,D2) * cos(6z) dz - z/L

D1和 D2（积分限制为0-> 2pi）。

我的代码是：

def equations(p, z):
    D1, D2 = p
    period = 2*np.pi

    P1 = lambda zz, D1, D2: \
            np.sign(-np.cos(np.pi*D1) + np.cos(6.*zz)) * \
            np.sign(-np.cos(np.pi*D2) + np.cos(5.*zz)) * \
            np.cos(6.*zz)
    P2 = lambda zz, D1, D2: \
            np.sign(-np.cos(np.pi*D1) + np.cos(6.*zz)) * \
            np.sign(-np.cos(np.pi*D2) + np.cos(5.*zz)) * \
            np.cos(5.*zz)

    eq1 = 2./period * integrate.quad(P1, 0., period, args=(D1,D2), epsabs=0.01)[0] + z
    eq2 = 2./period * integrate.quad(P2, 0., period, args=(D1,D2), epsabs=0.01)[0] - z
    return (eq1, eq2)


z = np.arange(0., 1000., 0.01)

N = int(len(z))

D1 = np.empty([N])
D2 = np.empty([N])
for i in range(N):
    D1[i], D2[i] = fsolve(equations, x0=(0.5, 0.5), args=z[i])

print D1, D2

不幸的是，它似乎并没有收敛。我对数值方法不太了解，希望有人可以帮我一下。

谢谢。

附言我也在尝试以下应等效的方法：

import numpy as np
from scipy.optimize import fsolve
from scipy import integrate
from scipy import signal

def equations(p, z):
    D1, D2 = p
    period = 2.*np.pi

    K12 = 1./L * z
    K32 = -1./L * z + 1.

    P1 = lambda zz, D1, D2: \
            signal.square(6.*zz, duty=D1) * \
            signal.square(5.*zz, duty=D2) * \
            np.cos(6.*zz)
    P2 = lambda zz, D1, D2: \
            signal.square(6.*zz, duty=D1) * \
            signal.square(5.*zz, duty=D2) * \
            np.cos(5.*zz)

    eq1 = 2./period * integrate.quad(P1, 0., period, args=(D1,D2))[0] + K12
    eq2 = 2./period * integrate.quad(P2, 0., period, args=(D1,D2))[0] - K32

    return (eq1, eq2)


h = 0.01
L = 10.
z = np.arange(0., L, h)

N = int(len(z))

D1 = np.empty([N])
D2 = np.empty([N])
for i in range(N):
    D1[i], D2[i] = fsolve(equations, x0=(0.5, 0.5), args=z[i])
    print
    print z[i]
    print ("%0.8f,%0.8f" % (D1[i], D2[i]))
    print

PSS：
我很好地完成了您所写的内容（我想我理解了！）。谢谢。不幸的是，我在这一领域确实没有太多技能，也不知道如何做出适当的猜测，所以我只猜测为0.5（我还在最初的猜测中添加了少量噪音，以尝试和改进它） 。我得到的结果似乎有数字错误，而且我不确定为什么，我希望您能指出正确的方向。因此，从本质上讲，我进行了FFT扫描（针对每个占空比变化进行了FFT，并查看了下图所示的5处的频率分量），发现线性部分（z / L）呈锯齿状。



PSSS：
谢谢您，我注意到了您建议的一些技术。我尝试复制您的第二张图，因为它似乎非常有用。为此，我保持D1（D2）固定并扫过D2（D1），并针对各种z值执行了此操作。 fmin并不总是找到正确的最小值（这取决于初始猜测），因此我扫除了 fmin的初始猜测，直到找到正确的答案。我也收到类似的答案。 （我认为是对的吗？）

另外，我想说的是，您可能想给我您的联系方式，因为该解决方案是寻找解决我遇到的问题（我是一名正在研究的学生）的步骤，我肯定会承认您在使用此代码的任何论文中都可以找到您。

#!/usr/bin/env python

import numpy as np
from scipy.optimize import fsolve
from scipy import integrate
from scipy import optimize
from scipy import signal


######################################################
######################################################

altsigns = np.ones(50)
altsigns[1::2] = -1

def get_breaks(x, y, a, b):
    sa = np.arange(0, 2*a, 2)
    sb = np.arange(0, 2*b, 2)

    zx  = (( x + sa) % (2*a))*np.pi/a
    zx2 = ((-x + sa) % (2*a))*np.pi/a
    zy  = (( y + sb) % (2*b))*np.pi/b
    zy2 = ((-y + sb) % (2*b))*np.pi/b

    zi = np.r_[np.sort(np.hstack((zx, zx2, zy, zy2))), 2*np.pi]
    if zi[0]:
        zi = np.r_[0, zi]
    return zi

def integrals(x, y, a, b):
    zi = get_breaks(x % 1., y % 1., a, b)
    sins = np.vstack((np.sin(b*zi), np.sin(a*zi)))
    return (altsigns[:zi.size-1]*(sins[:,1:] - sins[:,:-1])).sum(1) / np.array((b, a))

def equation1(p, z, d2):
    D2 = d2

    D1 = p
    I1, _ = integrals(D1, D2, deltaK1, deltaK2)

    eq1 = 1. / np.pi * I1 + z
    return abs(eq1)

def equation2(p, z, d1):
    D1 = d1

    D2 = p
    _, I2 = integrals(D1, D2, deltaK1, deltaK2)

    eq2 = 1. / np.pi * I2 - z + 1
    return abs(eq2)

######################################################
######################################################

z = [0.2, 0.4, 0.6, 0.8, 1.0]#np.arange(0., 1., 0.1)
step = 0.05

deltaK1 = 5.
deltaK2 = 6.

f = open('data.dat', 'w')

D = np.arange(0.0, 1.0, step)
D1eq1 = np.empty([len(D)])
D2eq2 = np.empty([len(D)])
D1eq1Err = np.empty([len(D)])
D2eq2Err = np.empty([len(D)])
for n in z:
    for i in range(len(D)):
        # Fix D2 and solve for D1.
        for guessD1 in np.arange(0.,1.,0.1):
            D2 = D
            tempD1 = optimize.fmin(equation1, guessD1, args=(n, D2[i]), disp=False, xtol=1e-8, ftol=1e-8, full_output=True)
            if tempD1[1] < 1.e-6:
                D1eq1Err[i] = tempD1[1]
                D1eq1[i] = tempD1[0][0]
                break
            else:
                D1eq1Err[i] = -1.
                D1eq1[i] = -1.

        # Fix D1 and solve for D2.
        for guessD2 in np.arange(0.,1.,0.1):
            D1 = D
            tempD2 = optimize.fmin(equation2, guessD2, args=(n, D1[i]), disp=False, xtol=1e-8, ftol=1e-8, full_output=True)
            if tempD2[1] < 1.e-6:
                D2eq2Err[i] = tempD2[1]
                D2eq2[i] = tempD2[0][0]
                break
            else:
                D2eq2Err[i] = -2.
                D2eq2[i] = -2.

    for i in range(len(D)):
        f.write('%0.8f,%0.8f,%0.8f,%0.8f,%0.8f\n' %(D[i], D1eq1[i], D2eq2[i], D1eq1Err[i], D2eq2Err[i]))
    f.write('\n\n')
f.close()

Answer 1

这是一个非常不适的问题。让我们回顾一下您想做的事情：


您要解决100000个优化问题
每个优化问题都是二维的，因此您需要O（10000）函数评估（为一维优化问题评估O（100）函数评估）
每个函数求值取决于对两个数值积分的求值
被积数包含跳跃，即它们是0次连续可微分
整数由周期函数组成，因此它们具有多个最小值和最大值


因此，您的处境非常艰难。此外，即使在最乐观的估计中，将被乘数中所有<1的因子都替换为1，积分也只能采用-2*pi和2*pi之间的值。比实际情况要少得多。因此，您已经看到只有解决方案

I1 - z = 0
I2 + z = 0

对于非常小的z。因此，尝试使用 z = 1000没有任何意义。

我几乎可以肯定，这不是您需要解决的问题。 （我无法想象会出现这样的问题的上下文。在傅里叶系数计算上看起来像是一个怪异的曲折……）但是，如果您坚持认为，最好的选择是首先在内部循环上工作。

如您所指出的，积分的数值评估会出现较大的误差。这是由于 sgn()函数引入的跳转引起的。诸如 scipy.integrate.quad()之类的函数倾向​​于使用高阶算法，这些算法假定被积物是光滑的。如果不是，它们的性能会很差。您要么需要手动选择一个可以处理跳跃的算法，要么在这种情况下更好的是手动进行积分：

以下算法计算 sgn()函数的跳转点，然后评估所有部分的解析积分：

altsigns = np.ones(50)
altsigns[1::2] = -1

def get_breaks(x, y, a, b):
    sa = np.arange(0, 2*a, 2)
    sb = np.arange(0, 2*b, 2)

    zx  = (( x + sa) % (2*a))*np.pi/a
    zx2 = ((-x + sa) % (2*a))*np.pi/a
    zy  = (( y + sb) % (2*b))*np.pi/b
    zy2 = ((-y + sb) % (2*b))*np.pi/b

    zi = np.r_[np.sort(np.hstack((zx, zx2, zy, zy2))), 2*pi]
    if zi[0]:
        zi = np.r_[0, zi]
    return zi

def integrals(x, y, a, b):
    zi = get_breaks(x % 1., y % 1., a, b)
    sins = np.vstack((np.sin(b*zi), np.sin(a*zi)))
    return (altsigns[:zi.size-1]*(sins[:,1:] - sins[:,:-1])).sum(1) / np.array((b, a))

这消除了数值积分的问题。这是非常准确和快速的。但是，即使是积分，对于所有参数而言也不是完全连续的，因此，为了解决优化问题，最好使用不依赖于任何导数的算法。 scipy中唯一的选择是 scipy.optimize.fmin()，您可以像这样使用：

def equations2(p, z):
    x, y = p
    I1, I2 = integrals(x, y, 6., 5.)

    fact = 1. / pi
    eq1 = fact * I1 + z
    eq2 = fact * I2 - z
    return eq1, eq2

def norm2(p, z):
    eq1, eq2 = equations2(p, z)
    return eq1**2 + eq2**2  # this has the minimum when eq1 == eq2 == 0

z = 0.25

res = fmin(norm2, (0.25, 0.25), args=(z,), xtol=1e-8, ftol=1e-8)
print res
# -> [ 0.3972  0.5988]

print equations2(res, z)
# -> (-2.7285737558280232e-09, -2.4748670890417657e-09)

您仍然面临为所有z查找合适的起始值的问题，这仍然是一项棘手的事情。祝好运！

编辑

要检查是否仍然存在数值误差，请将优化结果重新插入方程式中，看看它们是否满足要求的精度，这就是我上面所做的。请注意，我使用 (0.25, 0.25)作为起始值，因为从 (0.5, 0.5)开始并不会导致收敛。对于局部最小值（例如您的最小值）的优化问题，这是正常现象。除了尝试多个起始值，拒绝未收敛的结果外，没有更好的方法来处理此问题。在上述情况下，如果 equations2(res, z)返回的值比 (1e-6, 1e-6)高，我将拒绝该结果，然后使用其他起始值重试。对于连续优化问题，一种非常有用的技术是将前一个问题的结果用作下一个问题的起始值。

但是请注意，您不能保证 D1(z)和 D2(z)的平滑解决方案。 D1中的微小变化可能会使积分间隔偏离一个折点，从而导致积分值发生较大变化。通过使用 D2可以很好地调整算法，从而导致 D1(z)和 D2(z)的跳转。另请注意，由于 cos(pi*D1)的对称性，您可以取任何取模1的结果。

最重要的是：如果对积分使用解析公式，则不应有任何剩余的数值误差。如果残差小于您指定的精度，这就是您的解决方案。如果不是，则需要找到更好的起始值。如果不能，则可能不存在解决方案。如果解不是作为 z的函数连续的，那也是可以预期的，因为您的积分不是连续的。祝好运！

编辑2

看来您的方程式在间隔 z in [0, ~0.46]中有两个解，而对于 z > 0.46没有解，请参见下面的第一个图。为了证明这一点，请参见下面第二个图中的旧图形解决方案。等高线表示等式的解。 1（垂直）和等式。 2（水平），用于不同的 z。您可以看到，轮廓对于 z < 0.46交叉两次（两个解），对于 z > 0.46根本不交叉（没有一个解同时满足两个方程）。如果这不是您所期望的，则需要写下不同的等式（首先是我的怀疑……）




这是我使用的最终代码：

import numpy as np
from numpy import sin, cos, sign, pi, arange, sort, concatenate
from scipy.optimize import fmin

a = 6.0
b = 5.0

def P(z, x, y):
    return sign((cos(a*z) - cos(pi*x)) * (cos(b*z) - cos(pi*y)))

def P1(z, x, y):
    return P(z, x, y) * cos(b*z)

def P2(z, x, y):
    return P(z, x, y) * cos(a*z)

altsigns = np.ones(50)
altsigns[1::2] = -1

twopi = 2*pi
pi_a = pi/a
da = 2*pi_a
pi_b = pi/b
db = 2*pi_b
lim = np.array([0., twopi])

def get_breaks(x, y):
    zx  = arange(x*pi_a, twopi, da)
    zx2 = arange((2-x)*pi_a, twopi, da)
    zy  = arange(y*pi_b, twopi, db)
    zy2 = arange((2-y)*pi_b, twopi, db)
    zi = sort(concatenate((lim, zx, zx2, zy, zy2)))
    return zi

ba = np.array((b, a))[:,None]
fact = np.array((1. / b, 1. / a))

def integrals(x, y):
    zi = get_breaks(x % 1., y % 1.)
    sins = sin(ba*zi)
    return fact * (altsigns[:zi.size-1]*(sins[:,1:] - sins[:,:-1])).sum(1)

def equations2(p, z):
    x, y = p
    I1, I2 = integrals(x, y)

    fact = 1. / pi
    eq1 = fact * I1 + z
    eq2 = fact * I2 - z
    return eq1, eq2

def norm2(p, z):
    eq1, eq2 = equations2(p, z)
    return eq1**2 + eq2**2

def eval_integrals(Nx=100, Ny=101):
    x = np.arange(Nx) / float(Nx)
    y = np.arange(Ny) / float(Ny)
    I = np.zeros((Nx, Ny, 2))
    for i in xrange(Nx):
        xi = x[i]
        Ii = I[i]
        for j in xrange(Ny):
            Ii[j] = integrals(xi, y[j])
    return x, y, I

def solve(z, start=(0.25, 0.25)):
    N = len(z)
    res = np.zeros((N, 2))
    res.fill(np.nan)

    for i in xrange(N):
        if i < 100:
            prev = start
        prev = fmin(norm2, prev, args=(z[i],), xtol=1e-8, ftol=1e-8)
        if norm2(prev, z[i]) < 1e-7:
            res[i] = prev
        else:
            break
    return res


#x, y, I = eval_integrals(Nx=1000, Ny=1001)

#zlvl = np.arange(0.2, 1.2, 0.2)
#contour(x, y, -I[:,:,0].T/pi, zlvl)
#contour(x, y,  I[:,:,1].T/pi, zlvl)

N = 1000
z = np.linspace(0., 1., N)
res = np.zeros((N, 2, 2))

res[:,0,:] = solve(z, (0.25, 0.25))
res[:,1,:] = solve(z, (0.05, 0.95))