时间戳

首页 > 解决方案 > 超越方程的牛顿法

python - 超越方程的牛顿法

问题描述

超越方程：tan(x)/x + b = 0，其中 b 是任意实数。我需要介绍 n 并给我这个方程的 n 个解。

我的代码（Python）：

    from math import tan, cos, pi, sqrt, sin,exp
    import numpy as np 
    from matplotlib.figure import Figure
    import matplotlib.pyplot as plt

    def f(x,b):
        return tan(x)/x + b

    def f1(x,b):
        return (x/(cos(x)*cos(x)) - tan(x))/(x**2)

    e = 0.00000000001

    def newtons_method(x0, f, f1, e):
        x0 = float(x0)
        while True:
            x1 = x0 - (f(x0,b) / f1(x0,b))
            if abs(x1 - x0) < e:
                return x1
            x0 = x1

    result = []
    n = int(input("Input n: "))
    b = float(input("Input b: "))
    for i in range(2,4*n,1):
        result.append(newtons_method(i, f , f1, e))
    lambda_result = sorted(list(set(result)))
    print(len(lambda_result))

我使用步骤 1 更改初始近似值。根以 ~pi 的周期重复，因此第二个参数为 4*n。我通过集合排除重复的解决方案。如果 n 是 50，那么他只能找到 18 个解。需要修复什么才能使其正常工作？请帮帮我。

标签： pythonmathnumerical-methods

解决方案

在您的交叉帖子https://math.stackexchange.com/q/2942400/115115中，Yves Daoust 强烈建议您将牛顿方法基于函数

f(x)=sin(x) + b*x*cos(x)

或者

f(x)=sin(x)/x + b*cos(x)

因为这些函数没有极点或其他奇点（如果 x 不接近 0）。

至少对于大的c，解决方案接近初始值(i+0.5)*pi for i in range(n)。结果不需要对结果进行排序或缩短。

可以在余弦的根处使用正弦项具有交替符号。这为应用regula falsi (illinois)、Dekker's fzeroin、Brent's method等包围方法提供了一个完美的情况......这些方法几乎与 Newton 方法一样快，并且可以保证在区间内找到根。

唯一的复杂因素是区间(0,pi/2)，因为 . 将有一个非零根b<-1。必须删除x=0对于b接近-1.

当初始点足够接近根时，牛顿法仅可靠地归零到根。如果该点距离较远，接近函数的极值，则切线的根可能非常远。因此，要成功应用牛顿法，需要从一开始就找到好的根近似。为此，可以使用全局收敛的定点迭代或所考虑函数的结构简单近似。

  • 使用收缩定点迭代：周围的解k*pi也是等价方程的根x+arctan(b*x)=k*pi。这给出了近似解x=g(k*pi)=k*pi-arctan(b*k*pi)。由于即使对于小 ，圆弧切线也相当平坦k，因此这提供了一个很好的近似值。

  • 如果b<-1有一个正根k=0，那就是在区间内(0,pi/2)。以前的方法在这种情况下不起作用，因为圆弧切线1在此区间内具有斜率。根是由于方程的高阶非线性项，因此需要对方程的等价形式之一进行非线性近似。近似tan(x)=x/(1-(2*x/pi)^2)给出相同的极点，+-pi/2并且两者之间足够接近。将此近似值插入给定方程并求解，得到初始根近似值x=pi/2*sqrt(1+1/b)。

在实现中，必须移动根集b<-1以包含附加的第一个解决方案。

from math import tan, cos, pi, sqrt, sin, atan

def f(x,b):
    return sin(x)/x+b*cos(x)

def f1(x,b):
    return cos(x)/x-(b+x**-2)*sin(x)

e = 1e-12

def newtons_method(x0, f, f1, e):
    x0 = float(x0)
    while True:
        x1 = x0 - (f(x0,b) / f1(x0,b))
        if abs(x1 - x0) < e:
            return x1
        x0 = x1

result = []
n = int(input("Input n: "))
b = float(input("Input b: "))
for i in range(n):
    k=i; 
    if b >= -1: k=k+1
    x0 = pi/2*sqrt(1+1/b) if k==0 else  k*pi-atan(b*k*pi)
    result.append(newtons_method(x0, f , f1, e))
lambda_result = sorted(list(set(result)))
print(len(result), len(lambda_result))

推荐阅读