People who code: we want your input. Take the Survey
11

This is a root-finding problem for an analytic function over a single dimension. One standard technique is to approximate $F(k)$ as a polynomial $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ using a Chebyshev approximation, and to compute the roots of the polynomial semi-analytically, e.g. by setting up the companion matrix and computing its eigenvalues. (...


10

According to Wolfram Alpha, $x^5+3(x-1)=0$ has no closed-form solution, so you can forget about a nice closed-form expression. :) I see nothing wrong with Newton's method; it should be quick and accurate, and with some analysis like the one you are sketching I think you can identify a safe starting point and prove global convergence. It might even be faster ...


8

It all boils down to building a certain matrix from the polynomial coefficients and computing its eigenvalues. John Boyd did a lot of work in this area, these are some relevant papers: Boyd, John P. "A Fourier companion matrix (multiplication matrix) with real-valued elements: Finding the roots of a trigonometric polynomial by matrix eigensolving." ...


8

I will complement @Richard Zhang 's answer (+1) with a python implementation of his suggested approach. The MATLAB package Chebfun has been partially ported in python. Actually there are two versions available: chebpy and pychebfun. Here's an implementation of the root finding procedure with pychebfun (the approach is similar with chebpy) from scipy....


5

Here is a simple (Matlab) Newton method as a first attempt to help get started. It finds 1087 roots with error below $10^{-11}$. f = @(x) ((2*x)./(x.^2-1)) - tan(x); fp = @(x)-tan(x).^2+2.0./(x.^2-1.0)-x.^2.*1.0./(x.^2-1.0).^2.*4.0-1.0; x0 = 0; for jj = 1 : 1200 %number of iterations to find some roots x0 = x0 + (jj-1)*(jj/10^4); %take previous ...


5

As the other Answer already touches on the possibility of a symbolic root-solver being applied to this particular equation (by transforming into a polynomial form, albeit of degree $\ge 5$), I'll make some remarks about numerical root-finding. When we want to find real roots of a continuous function, the Intermediate Value Theorem is a basic tool. On a ...


5

Complementing the answers, I would like to say that directly discretizing the differential equation might work really well. I used that approach in a previous question. I used Finite Differences and used an eigensolver to find the roots. That, besides a perturbation analysis followed by Newton's method. The suggestion there was to use Chebyshev ...


4

In terms of computational effort, it is useless. I mean, you are still having a nonlinear problem with a bigger Jacobian (which is the worst part to be computed quickly). For next thoughts: dense problems are always worse that sparse ones. Try to solve a sparse system and compare it with the time needed to solve a dense one. Convergence does depend on the ...


4

You can rewrite with integer exponents through a change of variable. This looks like $d=x^r$ where $r$ is the LCM of the denominator of the rationalized exponents. In your case we have $r=4$ and: $$6.72427 - 0.125568 x^4 + 2.51327 x^{13}=0$$ Which can be solved in the usual ways. In Mathematica there is a clear distinction between symbolic and numerical ...


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

There are some unknowns in what you are doing but for simplicity, suppose we want to find $u(t)$ as discrete times $t_1, t_2, \cdots, t_n$. Let $\textbf{F} = [F(t_1), F(t_2), \cdots, F(t_n)]^T$ and $\textbf{u} = [u(t_1), u(t_2), \cdots, u(t_n)]^T$ be column vectors representing $F$ and $u$ evaluated at the desired times. From your problem statement, you wish ...


3

You have demonstrated that the polynomial has only a single positive root (I assume here that you are only interested in positive real roots). Using Cauchy's upper bound for polynomial roots $$1+ \max\left\{\left|\frac{a_{n-1}}{a_{n}}\right|, \left|\frac{a_{n-2}}{a_{n}}\right|, \ldots, \left|\frac{a_{0}}{a_{n}}\right|\right\}$$ you can obtain an upper bound ...


3

Suppose you want to minimize $$\Phi(x)=\frac{1}{2}||Ax-b||^2$$ The gradient is $$\frac{\partial \Phi}{\partial x} = A^T(Ax-b)$$ The step size to guarantee convergence is $$\alpha=||A^TA||^{-1}$$ Why? The direct solution to the problem is: $$x_{opt}=(A^TA)^{-1}A^Tb$$ This can be achieved iteratively if we look at the update on the estimate $x_k$. Suppose we ...


3

As this is a transcendent equation, finding all roots is not an option. You cannot (in general) find an expression that gives a closed-form for the roots of the equation. So you need to do some analysis first. Plotting the function, or plotting the left-hand-side and right-hand-side of the equation (and hence, the intersections of both curves are the roots) ...


2

You can at least take a quick look at the zeros using gnuplot by plotting contours at zero of the equation $$0 = \frac{2x}{x^2-1} - \tan(x)$$ set terminal png set output "test.png" set xlabel "x" set ylabel "y" set contour set cntrparam levels discrete 0 set view map unset surface set isosamples 1000,1000 splot 2*x/(x**2-1) - tan(x) Since the equation ...


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


2

Newton's method can refer either to a method for solving $f(x)=0$ where $f: R^{n} \rightarrow R^{n}$, or to a method for minimizing/maximizing a function $g: R^{n} \rightarrow R$ by solving the system of equations $\nabla g(x)=0$. Your function $h$ maps $R^{2}$ to $R$ and you want to find a zero of the function. This is typically done by minimizing $\min h(...


2

You can modify the iteration, for example by including a projection. That's because if $x\in D$ is indeed a fixed point of the function, i.e. $x=f(x)$, then it is also a fixed point of $x=\Pi(f(x))$ where $\Pi$ is some kind of operation so that $$ \Pi(x) = \begin{cases} x & \text{if $x\in D$} \\ \text{some $y\in D$} & \text{otherwise}. \end{cases} $...


2

A few options are KINSOL (Fortran and C interfaces are available; page 18 in the KINSOL documentation says you can use the banded solvers from LAPACK) scipy.optimize.root (particularly the large scale nonlinear solvers) fsolve in MATLAB allows you to specify the Jacobian as a sparse matrix. MATLAB codes from the book "Iterative Methods for Linear and ...


1

I think the main use of it is abstraction and the advantages that come with it: Applying a function to each element of an array comes up in many situations. Wrapping this pattern into a function raises the level of abstraction: you do not need to think in loops, your intent is probably clearer. If you are familiar with Python, comprehensions do this, as ...


1

There's a good discussion of globalization strategies for Newton's method root finding in Numerical Methods for Unconstrained Optimization and Nonlinear Equations by J. E. Dennis, Jr. and Robert B. Schnabel. See section 6.5 in particular. Dennis and Schnabel advocate treating the problem of finding a root of $f(p)=0$ as a nonlinear least squares problem: $\...


1

Adding another answer to the second part of your question: "Besides a constant step size, is there an easy variable-step size I could implement and play with?" An easy way to implement some variable step size would be the following algorithm: Consider your cost function $\Phi(x)$ you would like to minimize. Choose an initial step size $\alpha$ ...


1

The sum of squares should add up to what MATLAB says it adds up to, that is for $F:\mathbb{R}^m\to\mathbb{R}^n$ $$ r(x) = \sum\limits_{i=1}^{n}(F_i(x))^2$$ It is then easy to compute the gradient of $r$ as $$\nabla r(x) = J(x)^TF(x) $$, where $J$ is the Jacobian of $F$. Let us now try to reverse engineer if this is indeed what MATLAB means. Since it is not ...


1

You can use the companion matrix to find the roots via eigenvalue calculation: Companion matrix Note that you have to consider Sturm's theorem in order to find the number of complex/real roots without explicitly calculating them (that's the only way): Sturm's theorem


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

One of the most basic techniques in numerical analysis, when solving a complicated problem, is to construct an approximately-similar easy problem, solve that to obtain an approximate solution, and restart the process using that solution as a new starting point. This is a very general technique, and it’s not even restricted to root finding methods. Newton’s ...


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