Hot answers tagged

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


5

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

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


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


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


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


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


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

As I understand, your function (call it $f(x)$) has only one simple pole at $x=a$ and possibly more than one zero, and $a$ is not known, except that it is in $[0,2\pi]$. In addition, you do not know $f(x)$ in analytical form, which means that you cannot differentiate it analytically. In view of all this, I would suggest the following: Determine $a$ by ...


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

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

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


1

Numerical recipes in c provides closed form expression for real roots of quadratic and cubic which presumably have decent precision. Since the algebraic solution of the quartic involves solving a cubic and then solving two quadratics maybe a closed form quartic w good precision isn't out of the question.


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