超越方程的牛顿方法



超越方程: tan(x(/x + b = 0,其中 b 是任何实数。 我需要引入 n 并给我这个方程的 n 个解。

我的代码(蟒蛇(:

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 个解。需要修复什么才能使其工作?请帮帮我。

在你的交叉帖子 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(伊利诺伊州(,Dekker的fzeroin,Brent的方法提供了一个完美的情况,...这些方法几乎与牛顿方法一样快,并且保证在区间内找到根。

唯一的复杂情况是间隔(0,pi/2)因为b<-1将有一个非零根。人们必须从x=0中删除寻根过程,这对于接近-1b来说并非易事。


牛顿的方法只有零在起点足够接近根时可靠地进入根。如果点更远,接近函数的极值,则切线的根可能很远。因此,要成功应用牛顿方法,需要从一开始就找到良好的根近似。为此,可以使用全局收敛的定点迭代或所考虑函数的结构简单近似。

  • 使用收缩定点迭代:围绕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))

最新更新