Hot answers tagged

20

I would strongly advice against using closed form solutions since they tend to be numerically very unstable. You need to take extreme care in the way and order of your evaluations of the discriminant and other parameters. The classical example is the one for the quadratic equation $ax^2+bx+c=0$. Calculating the roots as $$x_{1,2} = \frac{-b \pm \sqrt{b^2-...


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

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

If $q:=|f'(x^*)|<1$, where $x^*$ is the solution, the fixed point iteration you talk about is locally linearly convergent with convergence rate $q$. Thus if $q$ is small or zero, the method is competitive with Newton's method. Far away from the solution, convergence is difficult to predict in the absence of global information (such as a Lipschitz ...


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 ...


6

I don't know about the theoretical complexity. But the practically fastest algorithm for finding all zeros of a polynomial given as a finite power series is the algorithm by Traub [1]. It factors a degree $n$ polynomial in typically $O(n^2)$ time. See also [2]. For factoring a polynomial accessible only through its function values, see, e.g., my paper [3]....


6

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 ...


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+\...


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

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

The Feigenbaum fractal is a good example of how strange fixpoint iteration can be: http://en.wikipedia.org/wiki/Feigenbaum_fractal http://en.wikipedia.org/wiki/File:Logistic_Bifurcation_map_High_Resolution.png The second link plots the behavior of fixpoint iteration applied to the logistic map as one of the parameters varies. For certain values it ...


4

The Banach fixed-point theorem describes the standard situation when a fixed-point iteration is globally convergent. Especially the uniqueness part of the theorem indicates that you can only expect local convergence if the solution is not unique. Most situations of local convergence can be explained by this theorem, at least in theory. This is even true for ...


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 ...


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

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

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 ) ...


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

As mentioned in [1], Reif [2] gives an $O(n \log^2 n (\log n + \log \epsilon))$ algorithm for the case when all roots are real, which is $O(n \log^3 n)$ as long as $\epsilon = n^{O(1)}$. He analyzes only the case of exact arithmetic, but after going through the algorithm I believe it is (or could be modified to be) backward stable. His method is an ...


3

See these: Solving Quartics and Cubics for Graphics, originally published in Graphics Gems V. The original code is here. See also this and this. A Universal Method Of Solving Quartic Equations.


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 ...


2

You may consider useful this reference: A Homotopy for Solving Large, Sparse and Structured Fixed Point Problems. R. Saigal. Mathematics of Operations Research, Vol. 8, No. 4 (Nov., 1983), pp. 557-578.


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

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 ...


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