用积分来解拟合方程。

I am trying to use numpy and scipy to solve the following two equations:

我尝试使用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

for D1 and D2 (integral limits are 0 -> 2pi).

对于D1和D2(积分限是0 -> 2)。

My code is:

我的代码是:

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

Unfortunately, it does not seem to converge. I don't know much about numerical methods and was hoping someone could give me a hand.

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

Thank you.

谢谢你！

P.S. I'm also trying the following which should be equivalent:

附注:我也在尝试以下的方法:

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: I implemented what you wrote (I think I understand it!), very nicely done. Thank you. Unfortunately, I really don't have much skill in this field and don't really know how to make a suitable guess, so I just guess 0.5 (I also added a small amount of noise to the initial guess to try and improve it). The result I'm getting have numerical errors it seems, and I'm not sure why, I was hoping you could point me in the right direction. So essentially, I did an FFT sweep (did an FFT for each dutycycle variation and looked at the frequency component at 5, which is shown below in the graph) and found that the linear part (z/L) is slightly jagged.

PSS:我实现了你写的东西(我想我懂了!)，做得很好。谢谢你！不幸的是，我在这个领域没有太多的技巧，也不知道如何做一个合适的猜测，所以我猜我只是猜了0.5(我还在最初的猜测中加入了少量的噪音来尝试改进它)。结果是我得到了数字错误，我不知道为什么，我希望你能告诉我正确的方向。所以本质上，我做了一个FFT扫描(对每个dutycycle变化做了一个FFT，并观察了5的频率分量，在图中显示)，并发现线性部分(z/L)略有锯齿。

PSSS: Thank you for that, I've noted some of the techniques you've suggested. I tried replicated your second graph as it seems very useful. To do this, I kept D1 (D2) fixed and swept D2 (D1), and I did this for various z values. fmin did not always find the correct minimum (it was dependent on the initial guess) so I swept the initial guess of fmin until I found the correct answer. I get a similar answer to you. (I think it's correct?)

PSSS:谢谢你，我已经注意到你提出的一些技巧。我试过复制你的第二张图，因为它看起来很有用。为了做到这一点，我将D1 (D2)固定并扫过D2 (D1)，我这样做是为了不同的z值。fmin并不总是能找到正确的最小值(它依赖于最初的猜测)，所以我把fmin的初始猜测扫了一遍，直到找到正确的答案。我得到了一个类似的答案。(我认为这是正确的吗?)

Also, I would just like to say that you might like to give me your contact details, as this solution as a step in finding the solution to a problem I have (I'm a student doing research), and I will most certainly acknowledge you in any papers in which this code is used.

同时,我想说,你想要给我你的联系方式,这个解决方案如一步找到解决问题的方案给我(我是一个学生做研究),和我最肯定会感谢你论文中使用这段代码。

#!/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()

1 个解决方案

#1

This is a very ill-posed problem. Let's recap what you are trying to do:

这是一个非常病态的问题。让我们回顾一下你想做的事情:

You want to solve 100000 optimization problems
你想解决100000个最优化问题。
Each optimization problem is 2 dimensional, so you need O(10000) function evaluations (estimating O(100) function evaluations for a 1D optimization problem)
每一个优化问题都是二维的，所以你需要O(10000)函数评估(估计一个一维优化问题的O(100)函数评估)
Each function evaluation depends on the evaluation of two numerical integrals
每个函数的评价依赖于两个数值积分的评价。
The integrands contain jumps, i.e. they are 0-times contiguously differentiable
被积体包含跳跃，即它们是0次连续可微的。
The integrands are composed of periodic functions, so they have multiple minima and maxima
这些积分是由周期函数组成的，所以它们有多个最小值和最大值。

So you are off to a very hard time. In addition, even in the most optimistic estimate in which all factors in the integrand that are < 1 are replaced by 1, the integrals can only take values between -2*pi and 2*pi. Much less than that in reality. So you can already see that you only have a chance of a solution for

所以你的日子很不好过。此外，即使在最乐观的估计中，被积函数的所有因子< 1被1替换，积分也只能取-2*和2*pi之间的值。比现实中要少得多。所以你可以看到你只有一个解决方案。

I1 - z = 0
I2 + z = 0

for very small numbers of z. So there is no point in trying up to z = 1000.

对于非常小的z数，所以在z = 1000时是没有意义的。

I am almost certain that this is not the problem you need to solve. (I cannot imagine a context in which such a problem would appear. It seems like a weird twist on Fourier coefficient computation...) But in case you insist, your best bet is to work on the inner loop first.

我几乎可以肯定，这不是你需要解决的问题。(我无法想象这样一个问题会出现在哪里。这似乎是傅里叶系数计算的一个奇怪的转折……但如果你坚持，最好的办法是先在内部循环。

As you noted, the numerical evaluation of the integrals is subject to large errors. This is due to the jumps introduced by the sgn() function. Functions such as scipy.integrate.quad() tend to use higher order algorithms which assume that the integrands are smooth. If they are not, they perform very badly. You either need to hand-pick an algorithm that can deal with jumps or, much better in this case, do the integrals by hand:

正如你所注意到的，积分的数值计算是有很大误差的。这是由于sgn()函数引入的跳转。像scipy.integrate.quad()这样的函数倾向于使用更高阶的算法，该算法假定这些积分是平滑的。如果不是，他们的表现非常糟糕。你需要手动选择一个算法来处理跳跃或者，在这种情况下，用手来做积分:

The following algorithm calculates the jump points of the sgn() function and then evaluates the analytic integrals on all pieces:

下面的算法计算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))

This gets rid of the problem of the numerical integration. It is very accurate and fast. However, even the integrals will not be perfectly contiguous for all parameters, so in order to solve your optimization problem, you are better off using an algorithm that doesn't rely on the existence of any derivatives. The only choice in scipy is scipy.optimize.fmin(), which you can use like:

这就消除了数值积分的问题。它非常准确和快速。然而，即使是积分也不会对所有的参数都是完全连续的，所以为了解决优化问题，你最好使用一个不依赖于任何衍生工具的算法。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)

You are still left with the problem of finding suitable starting values for all z, which is still a tricky business. Good Luck!

您仍然面临着为所有z找到合适的起始值的问题，这仍然是一个棘手的问题。好运！

Edit

编辑

To check if you still have numerical errors, plug the result of the optimization back in the equations and see if they are satisfied to the required accuracy, which is what I did above. Note that I used (0.25, 0.25) as a starting value, since starting at (0.5, 0.5) didn't lead to convergence. This is normal for optimizations problems with local minima (such as yours). There is no better way to deal with this other than trying multiple starting values, rejecting non-converged results. In the case above, if equations2(res, z) returns anything higher than, say, (1e-6, 1e-6), I would reject the result and try again with a different starting value. A very useful technique for successive optimization problems is to use the result of the previous problem as the starting value for the next problem.

为了检查你是否仍然有数值错误，将优化的结果插入公式中，看看它们是否满足要求的精度，这就是我上面所做的。注意，我使用(0.25,0.25)作为起始值，从(0.5,0.5)开始，并没有导致收敛。对于局部最小值(如您的)的优化问题，这是正常的。除了尝试多个起始值，拒绝非聚合的结果之外，没有更好的方法来处理这个问题。在上面的例子中，如果equations2(res, z)返回任何高于(比如，1e-6, 1e-6)的值，那么我就会拒绝这个结果，然后再尝试一个不同的初始值。对于连续的优化问题，一个非常有用的技术是将前面问题的结果作为下一个问题的起始值。

Note however that you have no guarantee of a smooth solution for D1(z) and D2(z). Just a tiny change in D1 could push one break point off the integration interval, resulting in a big change of the value of the integral. The algorithm may well adjust by using D2, leading to jumps in D1(z) and D2(z). Note also that you can take any result modulo 1, due to the symmetries of cos(pi*D1).

但请注意，你不能保证D1(z)和D2(z)的平滑解决方案。D1的一个微小变化可以将一个断点从积分区间上推掉，从而导致积分值的大变化。该算法可以使用D2进行适当的调整，导致D1(z)和D2(z)的跳跃。还要注意，由于cos(pi*D1)的对称性，您可以使用任何结果模块1。

The bottom line: There shouldn't be any remaining numerical inaccuracies if you use the analytical formula for the integrals. If the residuals are less than the accuracy you specified, this is your solution. If they are not, you need to find better starting values. If you can't, a solution may not exist. If the solutions are not contiguous as a function of z, that is also expected, since your integrals are not contiguous. Good luck!

底线:如果你使用积分的解析公式，不应该有任何剩余的数值不准确。如果残差小于您指定的精度，这就是您的解决方案。如果不是，那么您需要找到更好的启动值。如果不能，解决方案可能不存在。如果解不是连续的，作为z的函数，这也是期望的，因为你的积分不是连续的。好运！

Edit 2

编辑2

It appears your equations have two solutions in the interval z in [0, ~0.46], and no solutions for z > 0.46, see the first figure below. To prove this, see the good old graphical solution in the second figure below. The contours represent solutions of Eq. 1 (vertical) and Eq. 2 (horizontal), for different z. You can see that the contours cross twice for z < 0.46 (two solutions) and not at all for z > 0.46 (no solution that simultaneously satisfies both equations). If this is not what you expected, you need to write down different equations (which was my suspicion in the first place...)

你的方程在区间z(0， ~0.46)中有两个解，z > 0.46没有解，见下面第一个图。为了证明这一点，请参见下面第二个图中的老图形解决方案。等高线代表的是Eq. 1(垂直)和Eq. 2(水平)的解，对于不同的z，你可以看到z < 0.46(两个解)的等高线，而不是z > 0.46(没有同时满足两个方程的解)。如果这不是你所期望的，你需要写下不同的方程(这是我首先怀疑的…)

Here is the final code I was using:

下面是我使用的最终代码:

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

1 个解决方案

#1

更多相关文章

随机推荐