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))