Newton Method For Transcendental Equation
Solution 1:
In your cross-post https://math.stackexchange.com/q/2942400/115115 it was strongly recommended by Yves Daoust to base your Newton method on the function
f(x)=sin(x) + b*x*cos(x)
or
f(x)=sin(x)/x + b*cos(x)
as these functions do not have poles or other singularities (if x is not close to 0).
The solutions are, at least for large c
, close to the initial values (i+0.5)*pi for i in range(n)
. The result does then not require sorting or shortening of the result.
One could use that at the roots of the cosine the sine term has alternating sign. This makes a perfect situation to apply a bracketing method like regula falsi (illinois), Dekker's fzeroin, Brent's method,... These methods are almost as fast as the Newton method and are guaranteed to find the root inside the interval.
The only complication is the interval (0,pi/2)
as there will be a non-zero root for b<-1
. One has to remove the root-finding process from x=0
which is non-trivial for b
close to -1
.
Newton's method only zero's in to the root reliably when the initial point is close enough to the root. If the point is farther away, close to an extremum of the function, the root of the tangent may be very far away. Thus to successfully apply the Newton method, one needs to find good root approximations from the start. To that end one may use globally convergent fixed-point iterations or structurally simple approximations of the function under consideration.
Using the contracting fixed-point iteration: the solution around
k*pi
is also a root of the equivalent equationx+arctan(b*x)=k*pi
. This gives the approximate solutionx=g(k*pi)=k*pi-arctan(b*k*pi)
. As the arcus tangent is rather flat even for smallk
, this gives a good approximation.If
b<-1
there is a positive root fork=0
, that is in the interval(0,pi/2)
. The previous method does not work in this case as the arcus tangent has a slope around1
in this interval. The root is due to the higher order, non-linear terms of the equation, so one needs a non-linear approximation of one of the equivalent forms of the equation. Approximatingtan(x)=x/(1-(2*x/pi)^2)
gives the same poles at+-pi/2
and is sufficiently close in between. Inserting this approximation into the given equation and solving gives the initial root approximation atx=pi/2*sqrt(1+1/b)
.
In the implementation one has to shift the root set for b<-1
to include the additional first solution.
from math import tan, cos, pi, sqrt, sin, atan
def f(x,b):
returnsin(x)/x+b*cos(x)
def f1(x,b):
returncos(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))
ifabs(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==0else 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))
Post a Comment for "Newton Method For Transcendental Equation"