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


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

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

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

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

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

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

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

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

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

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

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


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


1

It seems your main concern is bracketing the root in as few iterations as possible, since each iteration is costly. In some cases you have found the regula falsi method to be unreliable, suggesting the bounds on the root may be large or the function cannot be linearly approximated easily. To handle such cases without being significantly slower than bisection,...


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


1

The derivatives $\partial_u f(u,v)$ and $\partial_v f(u,v)$ are two vectors tangent to yor surface that in general span the tangent space to the surface at the point $f(u,v)$. Then the normal vector to the surface at that point is $n(u,v)=\partial_u f(u, v)\times\partial_v f(u,v)$ and you want it parallel to say the coordinate direction $e_z=(0,0,1)^T$. Then ...


1

It probably depends on how much accuracy you want. Call $\pi - A = d$. If you know $d$, then you could use $\sin(\pi - x) = \sin(x)$, and evaluate $\sin(d)$. It's very close to zero, and a first-order Taylor series yields $\sin(d) \approx d$. (Python's math.sin(2e-33) yields 2 \times 10^{-33} to at least 10 significant digits.) Additional terms in the ...


1

The question changed, so this a new answer. If $(z,J(x)z)\geq\epsilon>0$, then take (w.l.o.g.) $f(0)=0$, and consider a path $g:[0,1]\to\mathbb{R}^n$ such that $$ g(0) = 0, \qquad f(g(0)) = 0, \qquad f(g(t)) = y t, \quad 0\leq t\leq 1, $$ so that $g(1)$ is the solution to $f(x)=y$. Then we differentiate with respect to $t$ to get an ODE for $g(t)$: $$ g'(...


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