使用 Python 查找隐式函数的根

Question

我正在尝试使用 Python 对以下隐式方程进行数值求解：

y*sin(c)+sin(y*c)=0

其中c 是一个常数。你会建议我怎么做？我实现了经典的 Newton-Raphson 方法，但它似乎没有收敛。

编辑：

代码如下：

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

import math

alpha = input('Define the alpha angle: ')

print >> f, 'alpha =', alpha

lambda_0 = input('Define an initial value: ')

print >> f, 'y_0 = ', y_0

gamma = math.pi - math.radians(alpha)

# def f(x,y):
#   p = 2.0 * x
#   t = p * y
#   val  = t * math.cos(t) - math.sin(t)
#   val /= (p * math.cos(t) + math.sin(p))
#   return val

toll = 1e-12

itmax = 50

diff = 1.0

it = 0

while diff > toll and it <= itmax:
    p = 2.0 * gamma
    t = p * y_0
    y_1  = t * math.cos(t) - math.sin(t)
    y_1 /= (p * math.cos(t) + math.sin(p))

    print >> f, 'p = ', p

    print >> f, 't = ', t

    print >> f, 'y at iteration ', it, ' = ', y

    diff = abs(y - y_0)
    y_0 = y
    it += 1

    print >> f, 'diff = ', diff

    print >> f, 'y_0 = ', y_0

    print >> f, 'it = ', it

f.close()

Answer 1

这是您的解决方案：y=0

也许您的意思是非零解决方案？还是别的什么？

编辑：添加一些更有用的东西

转换为更好的形式

假设您需要第一个积极的解决方案，您可以执行以下操作：

使用z=y*c 转换坐标，得到等式：

 A*z = sin(z)

在哪里A = -sin(c)/c。

这是确定通过斜率 A 原点的直线与正弦曲线的交点。

如果你画这个图，你会看到|A| >= 1 只有z=0 解决方案。

让你的方法收敛到正确的根的部分问题是选择根附近的起始值。在这种情况下，我们可以近似这些根。

临界斜率情况 (|A|~1)

对于接近 1 的|A|，我们可以看到在 0 附近会有一个正解和负解。我们可以使用正弦的低阶泰勒级数来近似这些。

A*z = z - z^3/6 + ...

对此的近似解是 z=0、z=R 和 z=-R，其中 R=sqrt( 6 (1-A) )。 这表明对于接近 1 的|A|，数值估计的良好起点是：

y=(1/c) sqrt( (6/c)( 1 + sin(c) ) )  

小坡度情况（|A|~0）

对于小 A，我们希望在 pi 附近有一个解决方案。在这种情况下，我们会进一步更改变量z=p+pi。

 A (p + pi) = sin( p + pi ) = -sin(p)

再次扩展罪恶给我们：

A p  + A pi = -p + ...

这给出了一个近似解：

p = - A pi / (A+1)

简化为

z = pi/(A+1)

这意味着我们应该寻找一个起点

y = pi/(c(A+1))

使用近似值

我可能会通过使用线性插值混合这两个值来选择一个起点。

您也可以在 A = 2/pi 附近对“中间”值进行类似的扩展，并在线性插值中使用三个点。

从这些 apprimxations 开始可能足以让牛顿法收敛到所需的值。但是，如果您确实需要确保收敛到正确的根，您可能需要放弃 Newton-Raphson 并将这些值用作割线或二等分法的起点。

Answer 2

如果 Newton raphson 不起作用，也许您应该尝试一些简单的方法，例如 Bisection 。时间复杂度仍然是 O(n) 量级，其中 n 是精度位数。

Answer 3

This 是随机 c (8) 的曲线。

据我所知，它有无限的解决方案，所以找到一个不会太难。看一下写得很好的 Raphson-Newton 的 scipy （你的代码可能是错误的）。看这里：http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.newton.html