Hot answers tagged

9

If you can stand to use complex arithmetic, simultaneous iteration methods might be preferable for computing all the roots of your polynomial. The simplest simultaneous iteration method, the (Weierstrass-)Durand-Kerner method, is effectively equivalent to applying Newton-Raphson to the Vieta relations relating the coefficients and roots of a polynomial, ...


8

This equation is not polynomial. Assuming both $K$ and $C$ are positive (as in your linked problem), then the solution of $C x \ln_2(x) - K = 0$ can be found in terms of the Lambert $W$ function or the Wright $\omega$ function: $$x = \frac{D}{W_0(D)} = \frac{D} {\omega\left(\ln\left( D \right)\right)}$$ where $D = \ln(2)K/C$. Corless, et al. 1996 is a ...


7

There are many good open-source implementations of root-finding methods out there already. One example is boost, whose implementation of Newton-Raphson and related methods you can find here. If you read its source code, you will be able to see what issues the authors wanted to address, and how they dealt with issues like convergence, user-specified tolerance ...


7

No. Although the question is phrased in terms of floating-point arithmetic, this is at heart a question about binary search. Binary search is an optimal algorithm among algorithms that rely on a decision tree (e.g., https://stackoverflow.com/q/4571478/491171). Any data structures and algorithms textbook should discuss this. Intuitively speaking, if at each ...


7

As pointed out by David Ketcheson, one method is to use the companion matrix and find its eigenvalues (this is what Matlab does for the roots function). However, if you want to code everything by yourself, you can try to use the Sturm sequences, which you use as a first step to find an interval with only one zero. Then, you can apply one standard methods ...


6

If you're using Matlab, you may want to try the Chebfun system (disclaimer: I used to be an active developer of this project). It can find all the roots of a one-dimensional function in a closed or open interval to machine precision. The main idea behind the Chebfun root-finder is to use a combination of recursive bisection and the Colleague Matrix, an ...


5

I suspect the main problem is the magnitude of the values. If you divide through by $\cosh$ to make all the numbers smaller, then ApproxFun doesn't seem to have a problem finding all the roots. function lhs(ε, α) A, B = sqrt((α-1)/(2ε)), sqrt((α+1)/(2ε)) return sqrt(α*α-1)*(cos(A) - 1/cosh(B)) - sin(A)*tanh(B) end u = Fun(x -> lhs(1e-3, x), [1.0,...


5

See @Kirill's answer, but depending on what sort of hardware you have available, you might be able to leverage parallelism to solve the problem more quickly. For example, say you had a way to evaluate $f$ for $P$ operands simultaneously - you could then repeatedly divide your bracket into $P+1$ intervals instead of $2$. Whether this would actually be "more ...


5

As you can see in this plot of $\log|f(\lambda)|$, $$ f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1), $$ the roots $\lambda_k$ are really regular, and are approximately equal to $-k$ (starting from $k\geq0$, $\lambda_0=1.23219$ is an exception). So the way to get the $k$-th eigenvalue ($k\geq1$), is to bracket the root in $[-k-\frac12,-k+\...


4

Since you are assuming $\eta$ is normal what i would do is try to compute the expectation as fast as possible for every $\theta$. So I would compute the expectation using any numerical integration I might know. That is, if $n$ and $N$ are the Gaussian density and distribution respectively, $$ f( \theta ) = \int_{ -\infty }^\infty R(\theta, \eta ) n( \eta ) ...


4

Let us parameterize the curve you are looking for by $(x(t),y(t))$ and let us for a moment assume that at all points on this curve $\nabla f(x(t),y(t)) \neq 0$, i.e., the curve never intersects another isocontour line. Then you know that at each point the tangent to the curve $(x(t),y(t))$ is parallel to the isocontour levels of $f$, i.e., that it is ...


4

One technique is to use a library with arbitrary large integers and use fixed point arithmetic. Just scale up x by some integer factor, then compute your polynomial in fixed-point, keeping only the sign. From two values, that you know positive and negative, you can then perform a binary search. At the very end, divide back by the fixed point factor and get ...


4

This is a classical problem in numerical methods research: evaluating the zeros of special functions. Many years of research have gone into devising efficient methods. The canonical starting point for anything related to this is the classic book by Abramowicz and Stegun. Scanned versions of the book can easily be found on the internet (I'm not sure if these ...


4

I suggest you sample the function and fit a rational function to the samples. You can find $a$ among the poles of the rational function, and the root among the real zeros. There are many good methods for finding rational function fits, and they would lead to extremely accurate approximations for this type of function with just a small number of samples. Some ...


3

Because the functions are strictly increasing, this can be reduced to the one-dimensional case. In each one-dimensional problem there are more efficient standard methods available, more efficient than bisection, and it also means one doesn't have to come up with a non-standard method oneself. After all, one-dimensional root-finding is a "solved" standard ...


3

The assumptions you have are no more specific than what you need to assume to make something like Newton's method work. In fact, you don't even assume enough to make the problem unique: you only assume that $F(x)$ is non-decreasing, when of course you need it to be monotonically increasing to make finding an answer unique. As a consequence, there is really ...


3

Assuming your jacobian isn't too large, it's a good idea to solve for your newton direction by QR factorization when your jacobian has a large condition number. It is numerically stable since it breaks the problem into two separate systems with the minimum error amplification factor allowed by the original problem. See here for details. This enables you to ...


3

You could consider $d=1$ and $$ f(x) = e^x, \qquad y=0. $$ I don't think a function's derivative being positive everywhere implies that the function's range is all of $\mathbb{R}$. It does imply, though, that solutions are unique. If $(x_1,y)$, $(x_2,y)$, $x_1\neq x_2$ are two solutions, then $$ 0 = (x_2-x_1)\cdot(f(x_1) - f(x_2)) = \int_0^1 (x_2-x_1)\cdot ...


3

Let's lay out some specs for what the halfIntervalSearch function is supposed to do, and go from there to changes necessary in your code. I'd assume from the name that it is supposed to find (and return) a root of a (hardcoded) function f using the method of interval bisection. The mathematics behind this is that a continuous function which changes sign in ...


3

This is a standard problem in numerical analysis. The standard approach if you can't compute derivatives of the function is to use a "bisection search" or, if you want to be fancy, its improvement using the golden ratio or Bent's method. See https://en.wikipedia.org/wiki/Bisection_method and https://en.wikipedia.org/wiki/Brent%27s_method If you can ...


2

There exists a clever algorithm for doing this, based on ideas due to Cayley, using the Bézout matrix. I learned this from [1], which describes an algorithm for a generic version of this problem. If we want to find all the critical points of the function $\|f(s)-f(t)\|^2$ (viewing $f$ as a function $\mathbb{R}\to\mathbb{R}^2$), a sufficient condition is ...


2

The asymptotic order of convergence of Brent's method tends to be either 1.618 and 1.689, of which the Alefeld-Potra-Shi's methods lie directly in between. The main difference of the Alefeld-Potra-Shi's methods are the tightness of the bounds, which Brent's method may fail to give during intermediate iterations. Overall, I do not recommend Alefeld-Potra-Shi'...


2

If you look at the full output of your script, sol has fields ipvt and qtf. Both of these fields come from MINPACK (http://www.netlib.org/minpack/lmdif.f), which is the backend here (https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize/_root.py#L246; https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy....


2

Any reasonable convergence criterion must be invariant to scaling of the function. A decent stopping criterion is therefore if $\|f(x_k)\|≤ \epsilon\|f(x_0)\|$ where $x_0$ is the starting point of the iteration. It's slightly annoying that the stopping criterion now depends on the starting point, but because Newton's iteration converges so quickly close to ...


2

It looks like the eigenvalues depend smoothly on your $r$ parameter. Since the problem seems to be easy to solve for small values of $r$, you can compute the eigenvalues for a set of radii $r_1,r_2,r_3$, and then extrapolate each of these eigenvalues using a quadratic function to some $r^\ast>r_3$. This extrapolation should then be an excellent starting ...


2

Assuming it is not too expensive to evaluate the function, I would recommend using the chebfun toolbox (matlab), as shown here: https://www.chebfun.org/docs/guide/guide03.html It builds an approximation of the function using Chebyshev polynomials and then finds the polynomial’s roots using a highly efficient root finder. It will find all roots in one fell ...


1

This is a technique called continuation. It typically works by using Newton's method to find one root, then you take steps along the root curve by picking a nearby point as the initial guess for the next Newton step. This is usually done to solve problems of the form $$F_\lambda(x_ \lambda)=0$$ where $\lambda$ is some parameter, but your problem can be ...


1

You can use backtracking to generate your polynomials. This can be parallelized in pretty much any language you choose. You may also want to call the Sturm sequence finder in the FLINT library from your script for speed.


1

As I said in my comment, I don't think you can necessarily count on the value of x in the last call to your function being the value that is returned as the solution. But, that being said, here is a way to do what you want. The solution in the code I show below relies on two components. First you define a MATLAB class that derives from handle. This makes it ...


1

A full log of the Dekker method run as in [1] (but with function value equality relaxed based on the tolerances) reads initial: x[ 0] = -4.000000000000000, f(x[ 0]) = -25.000000000000000; a=x[ 0], b=x[ 0], c=x[ 0] initial: x[ 1] = 1.333333333333333, f(x[ 1]) = 0.481481481481481; a=x[ 0], b=x[ 0], c=x[ 0] secant: x[ 2] = 1....


Only top voted, non community-wiki answers of a minimum length are eligible