Метод Ньютона для трансцендентного уравнения

Question

Трансцендентное уравнение: 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. Корни повторяются спериод ~ пи, поэтому второй аргумент 4 * п.Я исключаю повторные решения через множество.Если n равно 50, то он находит только 18 решений.Что нужно исправить, чтобы оно заработало?Помогите мне, пожалуйста.

Answer 1

В вашем кросс-посте https://math.stackexchange.com/q/2942400/115115 Ив Дауст настоятельно рекомендовал основывать ваш метод Ньютона на функции

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), фцероина Деккера, метода Брента , ... Эти методы почти так же быстры, как метод Ньютона, игарантированно найти корень внутри интервала.

Единственное осложнение - интервал (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))