1

I am trying to solve numerically the following implicit equation, using Python:

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

where c is a constant. How would you suggest me to do that? I implemented the classical Newton-Raphson's method, but it doesn't seem to converge.

EDIT:

Here is the code:

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

3 回答 3

2

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

也许您的意思是非零解决方案？或者是其他东西？

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

转换为更好的形式

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

变换坐标使用z=y*cwhich 给出方程：

 A*z = sin(z)

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

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

如果你画出这个图表，你会发现|A| >= 1只有z=0解决方案。

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

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

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

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

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

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 并将这些值用作割线或二等分法的起点。

于 2012-04-11T08:19:43.030 回答
2

If Newton raphson is not working, perhaps you should try something simple like Bisection . Time complexity is still in the order of O(n) where n is the number of bits of precision.

于 2012-04-10T21:16:53.553 回答
1

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

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

于 2012-04-10T21:02:45.813 回答